Geodynamic modelling of lithosphere subduction and the topographic response to slab tearing. Application to the Gibraltar Arc.

Geodynamic modelling of lithosphere subduction and the topographic response to slab tearing. Application to the Gibraltar Arc. Tesis doctoral presentada por KITTIPHON BOONMA al Departament de Dinàmica de la Terra i de l’Oceà de la Universitat de Barcelona para optar al grado de Doctor en Ciencias de la Tierra Group of Dynamics of the Lithosphere (GDL) Geociències Barcelona (GEO3BCN) Consejo Superior de Investigaciones Científicas (CSIC) Directores: Dr. Daniel Garcia-Castellanos Dr. Ivone Jiménez-Munt Tutor: Prof. Dr. Juan José Ledo Fernández Barcelona, 2021 Geodynamic modelling of lithosphere subduction and the topographic response to slab tearing. Application to the Gibraltar Arc. Ph.D. thesis presented at the Faculty of Geology of the University of Barcelona to obtain the Degree of Doctor in Earth Sciences Ph.D student: Kittiphon Boonma Approved: Dr. Daniel Garcia-Castellanos GEO3BCN-CSIC, Supervisor Dr. Ivone Jiménez-Munt GEO3BCN-CSIC, Supervisor Prof. Dr. Juan José Ledo Fernández Universitat de Barcelona, Tutor To my beloved grandpa ACKNOWLEDGEMENTS Scientific work cannot be done without team support, both professionally and personally, all of whom I would like to express my gratitude. I would like to express my sincere gratitude to Daniel and Ivone, my main supervisors at GEO3BCN, for giving me the opportunity to carry out this Ph.D. project. Over the last 4 years, they both have been providing me with helpful guidance and support. I would also like to thank Manel, who I worked with for the second project in the thesis. Throughout the project Manel provided me with insightful advice and enlightening discussions, all of which were invaluable. I also enjoyed all of his beer T-shirts. I am very grateful to Taras who had been the cornerstone of my numerical modelling work for the past 5-6 years. His insightful discussions were invaluable to me, also his calming cups of tea. I also thank Juan José Ledo Fernández for accepting the tutoring of this thesis and the doctorate commission from the University of Barcelona. Many thanks to my colleagues at GEO3BCN. They had always been friendly and welcoming, and their mental support was something I could not do without. Many thanks to all the support staffs at GEO3BCN, all of whom greeted me with a smile everyday. I would like to sincerely thank Ajay, who had been my Ph.D. brother from the very beginning. I appreciate all of your meals, memes, and encouraging words. My life in Barcelona would not be so great without my little Thai family: Mew, Au, and Fah. I miss all of our dinner parties, food adventures, and laugh attacks. I also thank my friends in Zürich for dinner parties, badminton games, and life chats. A special thank to Nes who had been like a mother to me during my years in Switzerland, I am forever indebted. Massive thanks to my best friend, Dr. Duff, you make one hell of a cheerleader. Your encouraging words and cheerful chats always helped me take my mind off things. Finally, I am forever grateful for my partner and family, for their love and support. Words cannot describe my gratefulness. Funding This work has been supported by EU Marie Curie Initial Training Network ‘SUBITOP’ (674899-SUBITOP-H2020-MSCA-ITN-2015) and MITE (Spanish Government national research program; CGL201459516). This thesis has been developed using the facilities of the Laboratory of Geodynamic Modelling from Geo3Bcn-CSIC. ABSTRACT Lithospheric slab tearing, the process by which a subducted lithospheric plate is torn apart and sinks into the Earth’s mantle, has been proposed as a cause of significant surface vertical motions. Although this has been linked to the change implied in the isostatic balance in subduction zones, little is known about the mechanisms and rock properties determining the tear propagation and the uplift-subsidence rates involved. This thesis aims to explore the link between the tearing of subducted lithospheric slabs and the associated vertical motions. To this purpose, I first explore the mechanisms controlling the buoyancy of a subducted lithosphere and then, with this understanding, numerically simulate the process of lithospheric tearing upon continental collision, using the Betic Cordillera as a reference scenario where such tearing-uplift interaction has been proposed for this region. With a mineral-physics approach, where a lithospheric mantle can be less dense than the underlying asthenosphere, I explore the controls on lithospheric buoyancy using a 2D thermaldiffusive model of plate convergence. Five chemical compositions and tectonothermal ages were considered, namely Archon ( >2.5 Ga), Proton ( 2.5−1.0 Ga), Tecton ( <1.0 Ga), and two oceanic lithospheric plates of 30 Ma and 120 Ma. While the advection of colder rock in oceanic-like plates always results in negative buoyancy, Protons and Tectons exhibit an ability to slowly flip from negative to positive buoyancy at low convergence rates: they first favour the sinking due to advection and then become more buoyant because they are thinner and heat up faster during subduction. In contrast, the lighter density of cratons (Archons) overprints this effect and hinders delamination or subduction, regardless of the convergence rate. This may explain why Archons are more stable during the Wilson Cycle. Having gained these insights into the role of lithospheric buoyancy in subduction settings, I then set to explore the characteristics surrounding lithospheric slab tearing and the associated surface uplift. I used 3D thermo-mechanical numerical modelling to investigate the geodynamic parameters affecting the slab-tearing initiation and its lateral propagation, and to quantify the corresponding surface vertical motions. The Betics-inspired model suggests CHAPTER 1 INTRODUCTION Plate subduction is one of the key mechanisms intervening in plate tectonics and shaping the surface and internal geodynamics of the Earth. Slab pull has been identified as a main driver of plate motions and subduction resulting from the higher potential density of the lithosphere relative to the underlying asthenospheric mantle. Slab breakoff, the process by which a subducted lithospheric slab detaches from the rest of the plate, has long been suspected to have a strong effect on the vertical motions of the surface. Lithospheric slab breakoff occurs in various styles including slab tearing, where an initial weakness develops into tearing and the tear then propagates laterally along the slab (Wortel and Spakman, 1992; Davies and von Blanckenburg, 1995). Through seismic tomography, a subducted slab has been interpreted to be laying under the Gibraltar Arc, the westernmost Mediterranean region (Wortel, 2000; Spakman and Wortel, 2004; Garcia-Castellanos and Villaseñor, 2011), with a presence of a partial lateral tear which is thought to be a result of slab tearing process (Fig. 1.1). This process, deep within the mantle, has been invoked to explain the changes in plate kinematics in the Western Mediterranean and the tectonic uplift that led to the Messinian Salinity Crisis event (MSC) ( 5.96−5.33 Ma) (Garcia-Castellanos and Villaseñor, 2011; Coulson et al., 2019). However, several aspects of slab tearing process are still poorly constrained, such as the condition of its initiation, the timing/duration of the tearing, its characteristics, and the associated surface response. This thesis aims to investigate the causes and links between the slab tearing process and the surface vertical motions, and apply this to the Gibraltar Arc/Western Mediterranean region. The significance of this work is that it provide new mechanistic understanding of slab tearing and the spatio-temporal development of vertical motions associated with this process in quantitative connection with the characterization of the lithosphere and with the 2Introduction constraints imposed by the geodynamic setting. The achievements attained here may be applicable in he future to geological scenarios other than the Gibraltar Arc. 1.1 Motivation The perception that large regions of continental crust have risen at rates that cannot be explained by crustal thickening or fault activity alone, has led to the necessity to identify the mechanism responsible for such long-wavelength surface uplift (England and Molnar, 1990). Slab breakoff is among the deep-seated mechanisms invoked to justify the long-wavelength, high rates of surface uplift (Davies and von Blanckenburg, 1995). It is driven by the same force that drives slab pull and subduction, i.e. the negative buoyancy of the lithospheric slab relative to the mantle (Garcia-Castellanos et al., 2000; Boonma et al., 2019; Jiménez-Munt et al., 2019). The concept of slab breakoff, as inferred from seismic tomography, was first proposed to be involved in the geodynamical evolution of the Mediterranean by Wortel and Spakman (1992). Slab breakoff was then used to explain post-collisional magmatism and exhumation of high-pressure rocks in the European Alps by Davies and von Blanckenburg (1995). Garzanti et al. (2018), and references therein, gave a comprehensive global overview of where slab breakoff has been invoked to explain changes in plate kinematics and tectonic deformation, e.g. the Alps (Davies and von Blanckenburg, 1995; Sinclair, 1997; Fox et al., 2015), the Mediterranean region (Carminati et al., 1998; Wortel, 2000; Rosenbaum et al., 2008; van Hinsbergen et al., 2010; Chertova et al., 2014; Spakman et al., 2018), the Anatolia-Zagros orogen (¸Sengör et al., 2003; Faccenna et al., 2006), and Himalayas and Tibet (van Hinsbergen et al., 2012; Wu et al., 2014; Liang et al., 2016). These studies often ascribe short-lived, long-wavelength uplift or exhumation events or sudden pulses in sediment supply to slab breakoff. However, they often neglect the influence of the 3D geometrical configuration of each tectonic regions. How likely were the tectonic configurations in those domains to have caused the slab tearing in the first place? How much does slab breakoff contribute to the buoyancy-driven isostatic surface uplift? The Gibraltar Arc region is an ideal tectonic scenario to study this link between the slab tearing and its consequent surface vertical motions because it supposedly triggered a conspicuous uplift of the Betic and Rifean Cordilleras that had a dramatic impact on the entire Mediterranean Basin. Based on seismic tomographic imaging, Wortel (2000) suggested that slab tearing might have occurred in the Gibraltar Arc region as a consequence of the continental collision and the subsequent slab rollback declined during early Miocene. The majority 1.2 Objectives 3 of topographic growth in the Betics appear to have initiated after late Tortonian ( ∼ 7 Ma), sometimes under little amounts of tectonic fault deformation (Comas et al., 1999; Iribarren et al., 2009). The surface vertical motions observed in the Internal Betics zone after late Tortonian are best constrained from the present elevation of tectonically undeformed Miocene marine sediment in that region, often above 600 m elevation (Garcés et al., 1998; Iribarren et al., 2009). This has been interpreted as the result of a westward migration of a lateral tear within the steeply hanging Ligurian-Tethys slab seen in tomography (Garcia-Castellanos and Villaseñor, 2011). The uplift of the intramountain basins within the Betics and Rif has been linked to the closure of the Gibraltar marine gateways during the Late Miocene which led to the partial desiccation of the Mediterranean Sea, known as the Messinian Salinity Crisis event (Garcia-Castellanos and Villaseñor, 2011; Coulson et al., 2019). The Messinian Salinity Crisis (MSC) ( 5.96−5.33 Ma) marks a period of dramatic sea-level change, possibly the most abrupt environmental change on Earth since the beginning of the Tertiary. Despite the slab tearing process being linked to these changes in plate kinematics and surface uplift, the timing of its mechanism has been poorly constrained and barely tested by thermo-mechanical models, prompting questions regarding the timing of tearing initiation and the duration of tearing process. The key motivation behind this work, therefore, is to understand the characteristics of the geodynamic setting that lead to slab tearing and control the speed of the associated surface vertical motions and slab tearing process, using numerical models inspired by the Western Mediterranean. Additionally, this could shed some light into how the tectonic uplift and slab tearing influenced on the uplift of the Betic Cordillera during the onset of the Messinian Salinity Crisis event in the westernmost Mediterranean sea. 1.2 Objectives In continental collision settings, slab breakoff is one scenario which can occur. It is driven by the density contrast between the lithospheric mantle and the asthenosphere, and between the subducted oceanic lithosphere and a more buoyant continental lithosphere that follows into the subduction zone. The overarching objectives of this thesis are to shed light on the lithospheric buoyancy relative to the asthenosphere, and to understand the role of this buoyancy force in promoting the lithospheric slab tearing during subduction. The scenario we use to address these objectives is the Gibraltar Arc, with particular attention to its Late Messinian evolution, during which the Messinian Salinity Crisis (MSC) occurred. The tectonic uplift that closed off the Mediterranean Sea from the Atlantic ocean, during the first 4Introduction Fig 1.1 (A) Map of lithosphere thickness for the Betic–Rif system (Fullea et al. (2010); their Fig. 6B) with location of three seismic tomographic transects. These image the 3D geometry of the cold lithospheric slab hanging beneath the orogen, which is continuous along its westernmost part (section 1, (Spakman and Wortel, 2004); their Fig. 2.6A), and along section 2, (Garcia-Castellanos and Villaseñor (2011), their Fig. SI-4), but is interrupted beneath the eastern and central Betic (section 3). (B) Cartoon showing the geometry of the subducted lithospheric slab beneath the Betic-Rif system (from Vergés and Fernàndez (2012)). 1.3 Thesis structure 5 stage of the MSC event, is thought to be a consequence of slab tearing process deep within the mantle. Specifically, the first objective is to understand the buoyancy forces involved in a subduction process. To this end, I first look at the role of lithospheric buoyancy in subduction setting and its dependence on mantle composition and convergence rate. For this study, I developed a 2D kinematic numerical code to model the tectonic shortening and subduction process. I set out to study different mantle compositions of various tectonothermal ages and how their buoyancy forces evolve through a subduction scenario, and what implication do they have on the stability of the lithosphere during tectonic plates convergence. Understanding the role of lithospheric buoyancy will help better the understanding of the dynamics and the characteristics of lithospheric slab tearing. In order to study the characteristics of slab tearing, I first have to investigate the mechanisms involved in the slab tearing initiation, such as the roles of tectonic configurations and mantle rheologies. To achieve this objective, I use 3D thermomechanical numerical modelling to model lithospheric subduction, slab rollback, and eventually slab tearing. The 3D numerical code, I3ELVIS (Gerya, 2013) is used here because of its ability to model complex non-linear problems and its calculation routines for various geodynamical processes. The key physical parameters for the slab tearing which I will test out are the mechanical properties of the mantle (viscosity, brittle failure strength, etc.). The two main datasets of the mantle’s properties will be from Karato and Wu (1993) and Ranalli (1995). The final objective aims to study the surface vertical motions resulting from slab tearing process. This objective hopes to provide an insight into the coupling between deep (slab tearing) and shallow (surface uplift) processes. I will look at dynamic topography and how it reflects the mantle dynamics during the tearing process; the topographic evolution throughout the tearing process; and the quantification of uplift rate resulting from the tearing process. The resulting vertical motions should also help constraining the timing and uplift rate during the first stage of the MSC event. 1.3 Thesis structure This thesis is split into 7 chapters, completing with the essential geological background, the mathematics and physics behind the numerical modelling process, and the output from the investigations. Chapter 2 introduces the reader to the general concepts of plate tectonics, subduction systems, and the geodynamic processes necessary to understand the context of this work. It will 6Introduction also give and overview of the geological settings and history of the Western Mediterranean, specifically the Gibraltar Arc region, which inspires the work in this thesis. Chapter 3 describes the fundamentals of numerical methods utilised in this thesis, including the governing equations, rheologies, and boundary conditions. The work in this volume consists of both 2D and 3D numerical modelling work. The 2D kinematic numerical modelling was carried out using a thermally advective-diffusive code ‘LithBuoy’ 1 (Boonma et al., 2019), which models the buoyancy of the lithospheric mantle in subduction setting, with an implication for lithospheric delamination. The 3D thermomechanical coupled numerical modelling code ‘I3ELVIS’ (Gerya, 2013) was used to study the dynamics, the initiation, and characteristics of lithospheric slab breakoff (specifically the one-sided slab tearing) and the corresponding topographic response. Chapter 4 investigates the effect of tectonic convergence and mantle composition on the buoyancy of lithospheric mantle, using the 2D kinematic numerical code that I developed. The study looked at mantle density from a mineral physics viewpoint where the continental lithospheric mantle is lighter than the underlying asthenosphere, posing a problem of the initiation of mantle delamination. The study considered five different types of lithospheric mantle of different tectonothermal ages The material related to this chapter is published: Boonma, K. , Kumar, A., Garcia-Castellanos, D., Jiménez-Munt, I., and Fernández, M. (2019). Lithospheric mantle buoyancy: the role of tectonic convergence and mantle composition. Scientific Reports, 9(1), 17953. DOI: 10.1038/s41598-019-54374-w The numerical modelling code LithBuoy, which came from this study, contributes to a manuscript in preparation: Irene DeFelip, Juan Alcalde, Juan Carlos Afonso, Eldar Baykiev, Isabel Bernal, Kittiphon Boonma , Ramon Carbonell, Stephanie Flude, Arnau Folch, Javier Fullea, Daniel García-Castellanos, Adelina Geyer, Santiago Giralt, Armand Hernández, Ivone Jiménez-Munt, Ajay Kumar, Maria-Gema Llorens, Joan Martí, Cecilia Molina, Andrés Olivar-Castaño, Andrew Parnell, Martin Schimmel, Montserrat Torné, Sergi Ventosa (in preparation). Towards a digital twin of the Earth system: Geo-Soft-CoRe, a Geoscientific Software & Code Repository. Chapter 5 focuses on studying the characteristic of slab tearing, what could initiate the tearing, and how does the surface response to this process happening in the deep mantle. This chapter presents the results from 3D numerical modelling of the lithospheric slab tearing. 1available at https://github.com/kboonma/LithBuoy 1.3 Thesis structure 7 Several model were constructed with various parameterisation to study and how each of them affect the subduction process and slab tearing process. The manuscript related to this chapter is submitted. Boonma, K. , Garcia-Castellanos, D., Jiménez-Munt, I., and Gerya, T. (submitted, under revision). Thermomechanical modelling of lithospheric slab tearing and its topographic response in the Gibraltar Arc (westernmost Mediterranean Sea). Chapter 6 discusses various aspects from our studies such as the effect of convergence rate and mantle compositions on slab buoyancy and its tectonic relevance, the initiation of slab tearing, the dynamics of slab tearing process and the associated surface vertical motions, and finally the overall implications for the Gibraltar Arc region and the MSC event. Lastly, Chapter 7 draws final conclusions from each studies in this thesis volume and proposes outlook for future work regarding the geodynamic processes surrounding continental collisions. 8Introduction Geodynamic modelling of lithosphere subduction and the topographic response to slab tearing. Application to the Gibraltar Arc. •Understanding the stability of the lithosphere during the convergence of tectonic plates. •Implication for processes surrounding tectonic convergence setting such as mantle delamination, slab breakoff, and lithospheric underthrusting. •Interaction between advection and diffusion in a slab density dependence on temperature, pressure and chemical composition 1. Role of slab buoyancy in subduction setting •Linking deep and shallow processes. •The changes in crustal surface in response to the slab tearing process. •Implication of uplift rate for the Messinian Salinity Crisis. 3. Topographic response to slab-tearing •The role of model configuration in assisting with the initiation of slab tearing. •Understanding how mantle rheology controls the characteristics of slab tearing. •The geodynamic surrounding the slab-tearing event. 2. Geodynamics of slab-tearing LithBuoy (Boonma et al.,2019) •2D kinematic modelling of plate convergence. •Study the effect of mantle compositions and convergence rate. I3ELVIS (Gerya, 2013) •3D thermo-mechanical geodynamic codes. •Conservation laws, viscoplastic rheology, and grainsize reduction. Constraints on geodynamic setting of Gibraltar Arc •Seismic tomography data •Surface uplift data (geology & sedimentology) •Palaeogeographic reconstruction Motivation: What drives surface vertical motions in subduction settings? Fig 1.2 Thesis workflow, outlining the key studies and their main scopes CHAPTER 2 FUNDAMENTALS & GEOLOGICAL SETTING This chapter presents the fundamental geological concepts of plate tectonics, specifically convergent plate boundary and continental collision settings, and the key geodynamic processes involved in such settings. The last section of the chapter sets the scene with geological settings of the Western Mediterranean, the Gibraltar Arc region, and the Messinian Salinity Crisis (MSC) event. Plate tectonics is a theory which describes the relative motions of the Earth’s outer shell − the lithosphere, which is qualitatively defined as the Earth’s cold upper layer. The lithosphere is divided into a patchwork of thin and rigid plates (‘tectonic plates’) which move slowly relative to one another, with various velocities up to 100 mm/year. There are 15 major plates and many minor plates (Fig 2.1). The boundaries between these plates are where the majority of the world’s earthquakes, volcanic eruptions, mountain-building, and oceanic trench formation occur. 2.1 The Lithosphere The Earth has a cool and mechanically strong outermost shell called the lithosphere (Greek lithos means ‘rock’). The lithosphere is a chemical, mechanical, and thermal boundary layer which floats on top of the mechanically weak asthenosphere. The asthenosphere (Greek asthenia means ‘weak’ or ‘sick’) is a part of the mantle immediately beneath the lithosphere. The lithosphere includes, both, the crust and the uppermost part of the mantle (‘lithospheric mantle’), being the thinnest in oceanic regions (20–100 km) and the thickest in continental regions (50–250 km). The lithosphere-asthenosphere boundary (LAB) lies at the base of the lithosphere, and therefore determines the thickness of the lithosphere. The lithosphere has 16 Fundamentals & Geological setting Fig 2.4 A schematic diagram showing two continental lithospheric plates approaching each other with an intermediate oceanic lithosphere being consumed at the subduction zone (from Encyclopedia Britannica). Fig 2.5 Different possible subduction scenarios during continental collision (from Magni et al. (2013)). 2.3 Continental Collision 17 2.3.2 Lithospheric delamination The lithosphere delamination process was introduced by Bird (1979). According to Bird’s model, delamination occurs when the asthenosphere comes into direct contact with the crust and, then, the lithospheric mantle peels back from the upper-middle crust along the weaker crust, migrating the delamination point. It is a geodynamic mechanism often invoked for the evolution of collision zones. Mantle delamination process recycles the lithospheric mantle back into the asthenosphere, and it is driven by the gravitationally instability of the lithosphere itself. Mantle delamination could arise from four possible causes, convective instability, rifting, plume erosion, or continental collision. Mantle delamination could lead to regional uplift, increased heat flow, reduced seismic velocities, and/or mafic volcanism (Bird, 1979; Houseman et al., 1981; Kay and Mahlburg Kay, 1993). Several studies looked at different component that could influence the delamination process such as the viscosity contrast between lower crust and the underlying lithospheric mantle (Bird, 1979; Channell and Mareschal, 1989; Schott and Schmeling, 1998; Gö˘ gü¸s and Pysklywec, 2008; Faccenda et al., 2009; Valera et al., 2011; Gray and Pysklywec, 2012; Bajolet et al., 2012; Magni et al., 2013; Gö˘ gü¸s et al., 2016; Gö˘ gü¸s and Ueda, 2018), but only few studies shed light on the density contrast between the lithospheric mantle and the underlying asthenosphere (Channell and Mareschal, 1989; Gö˘ gü¸s and Pysklywec, 2008; Valera et al., 2008). There are also numerous geophysical (seismicity, seismic tomography, or gravity) and geological (subsidence or uplift) observations that some authors accounted for mantle instability and potential lithospheric removal to explain them, such as in the Alboran domain (Seber et al., 1996; Turner et al., 1999; Platt et al., 2006; Valera et al., 2011; Thurner et al., 2014), the Apennines (Piana Agostinetti et al., 2002; Bartolini, 2003; Carminati et al., 2003, 2004; Chiarabba and Chiodini, 2013), Eastern Anatolia (Al-Lazki et al., 2003; Aydin et al., 2005; ¸Sengör et al., 2008; Berk Biryol et al., 2011; Sen et al., 2011; Komut, 2015), Tibet (Molnar et al., 1993; Chung et al., 2005; Jiménez-Munt and Platt, 2006; Ren and Shen, 2008), or Sierra Nevada (Ruppert et al., 1998; Wernicke and Snow, 1998; Saleeby and Foster, 2004; Jones et al., 2014; Valera et al., 2014). Despite several questions regarding delamination have already been addressed, from different perspectives, there are still open questions about the conditions under which the mechanism operates. Two major factors which control the delamination processes, as inferred by geodynamic models, have been reviewed and classified (Gö˘ gü¸s and Ueda, 2018): (i) the presence and combination of the weak lower crust, the vertical decoupling channel, and mantle uprising; and (ii) the negative buoyancy force of the foundering lithosphere. The general condition that 18 Fundamentals & Geological setting is thought to lead to mantle delamination is that the lithospheric mantle must be denser than the underlying asthenosphere so as it sinks downwards, the buoyant asthenosphere makes contact with the crust, replacing the denser lithospheric mantle. 2.3.3 Lithospheric slab breakoff The concept of lithospheric mantle delamination has been conceived as a mechanism to cause tectonic uplift and magmatism in orogenic belts (Bird, 1979; Houseman et al., 1981; Kay and Mahlburg Kay, 1993). However, in the past few decades, the concept of slab breakoff has been widely accepted as the more favourable mechanism to cause collision zone magmatism and subduction-zone rocks exhumation (Wortel, 1982; Davies and von Blanckenburg, 1995; von Blanckenburg and Davies, 1995; Niu, 2017). Slab breakoff is one of the plausible scenario in continental collision setting where the subducted oceanic lithosphere detaches from the buoyant continental lithosphere that follows it into the subduction zone. Once the continental lithosphere is introduced into the subduction zone, its positively buoyant nature opposes the negative buoyancy force exerted by the hanging oceanic lithosphere, this opposing forces cause tensile stress to build up and the hanging oceanic slab eventually detach from the buoyant continental part (Wortel and Spakman, 1992; Davies and von Blanckenburg, 1995). Multiple geological and geophysical investigation have acknowledged slab detachment as a plausible explanation for observations in several regions (Andrews and Billen, 2009), also several laboratory and numerical modelling experiments have been under taken to characterise the breakoff process (Davies and von Blanckenburg, 1995; Carminati et al., 1998; Andrews and Billen, 2009; Duretz et al., 2014; Gerya et al., 2004; Faccenda et al., 2009; van Hunen and Allen, 2011; Duretz and Gerya, 2013; Boutelier and Cruden, 2017; Magni et al., 2017; Schellart, 2017; Fernández-García et al., 2019). Since the appearance of the concept, several efforts have been made in order to better understand the dynamics of slab breakoff. Numerical modelling studies, in particular, have shown that the depth of the slab detachment is influenced by the age and the rheology of the subducting plate (Fernández-García et al. (2019) and references therein). In this thesis, we will take a closer look at the slab breakoff process, specifically the ‘slab tearing’ process (Fig 2.6), e.g. what tectonic configuration could initiate the tearing, how does the tearing propagates, and what is the surface response. 2.4 Western Mediterranean 19 Fig 2.6 Schematic diagram displaying the ‘slab tearing’ process within the context of continental collision (from Wortel (2000)). 2.4 Western Mediterranean 2.4.1 Tectonic reconstructions Several studies have attempted to reconstruct the relative position of Iberia, relative to NW African/Moroccan plate through time (Fig 2.7), Specifically in the Western Mediterranean, many have proposed detailed tectonic reconstructions of the last ∼35 Myr, all of which based on similar geological and plate kinematic constraints (e.g. Rosenbaum et al. (2002); Jolivet et al. (2009); Carminati et al. (2012); Vergés and Fernàndez (2012); Van Hinsbergen et al. (2014)). There appears to be a consensus that the main geodynamic process that drives the tectonic change in the Western Mediterranean for the past ∼35 Myr could be the subduction rollback, however, there are still disagreements about the spatial-temporal tectonic evolution. There are three schools of thought in debate about the slab rollback evolution in the Western Mediterranean (Fig 2.8). The first scenario involves an initially short NW-dipping slab confined to the Baleares margin (Fig 2.8a). After the rollback started it split into two sections, one section continues towards Africa while another section rotates and rollback towards Gibraltar (Rosenbaum et al., 2002; Spakman and Wortel, 2004; de Lamotte et al., 2009; Van Hinsbergen et al., 2014). The second scenario consists of a rollback originating from a long N-NW dipping subduction zone, from Gibraltar in the west to the Baleares in the north-east (Fig 2.8b). The rollback occurs as one long section with the eastern side 20 Fundamentals & Geological setting rolled back more than the western side (Faccenna et al., 2004; Jolivet et al., 2009). The third scenario has the rollback starting from a SE dipping subduction under the North-African margin then migrates and rotates into southern Iberia and Gibraltar Arc (Fig 2.8c) (Vergés and Fernàndez, 2012). 2.4.2 Tectonic evolution of the Gibraltar Arc The convergence between Iberian peninsula and north-western African plate marks the westernmost termination of the Mediterranean orogenic belt (Carminati et al., 2012; Faccenna et al., 2014a). This convergence gave rise to a lithospheric collision that formed the Betics and Rif Cordilleras, as well as the Alboran Sea (Fig 2.9). This NW Africa-Iberia convergence started during Late Cretaceous, with convergence rate of 4.5 mm/yr (McClusky et al., 2003; Stich et al., 2006; Macchiavelli et al., 2017), and persisted until the present with decreasing convergence rates (Argus et al., 1989; Macchiavelli et al., 2017). The Western Mediterranean subduction zone can be classified, based on the tectonic evolution, into three distinct structural domains: Internal Zone, External Zone, and the Suture Zone (e.g. Jabaloy Sánchez et al. (2019) and references therein). The Internal Zone (Alboran Domain) is composed of rocks from Paleozoic to TriassicJurassic age, which had been deformed by a complex multiphase deformation during Eocene to early Miocene. The rocks in the internal zone underwent Alpine ductile deformation and metamorphism (Andrieux et al., 1971; Casciello et al., 2015; Jabaloy Sánchez et al., 2019). The External Zone (South-Iberain paleomargin) rocks were recognised as Cretaceous and Tertiary with Miocene sediments, and underwent deformation after early Miocene and before the late Miocene (Andrieux et al., 1971; Casciello et al., 2015). This zone is composed of a thin-skinned fold-and-thrust belt, which is resulted from the tectonic inversion of the South-Iberian rifted paleomargin. The Suture Zone is the tectonic boundary between the Internal and External Zones. This zone is considered an orogenic suture attesting the closure of an oceanic domain, despite the absence of ophiolite (uplifted and resurfaced Earth’s oceanic crust and the underlying upper mantle). The suture zone consists of two sediment groups, the Flysch Trough Units (Campo de Gibraltar Complex) and the Frontal Units. The Flysch Trough Units overthrusts the External Zone in the westernmost Betics and the Rif. These units were presumably deposited over an oceanic crust, as inferred from a relationship between deep-water sedimentary setting and oceanic-type crust (Durand-Delga et al., 2000; Michard et al., 2005). The Flysch Trough Units are, therefore, interpreted as an accretionary wedge formed at an oceanic-continental subduction zone, which became imbricated during the subsequent continental collision (Luján 2.4 Western Mediterranean 21 Fig 2.7 Plate reconstructions for the Iberian and NW African/Moroccan plates. Panel a-d show the reconstructions of the relative motion between the European, Iberian, and NW African/Moroccan plates. We are particularly interested in the evolution of the plate margin above in African/Moroccan plate which you can see migrate towards the southern Iberia in d. Panel e displays a predicted N-S convergence rate between the three plates. Again, here, we pay interest to the green line which represents the N-S convergence rate between NW African and Iberian plates. (from Macchiavelli et al. (2017)) 22 Fundamentals & Geological setting Fig 2.8 The three different scenarios of the rollback reconstruction of the of the past 35 Myr in the Western Mediterranean. Scenario 1 starts from an initial short subduction zone near the Baleares (Rosenbaum et al., 2002; Spakman and Wortel, 2004; Van Hinsbergen et al., 2014); Scenario 2 involves a long initial trench along the entire Gibraltar-Baleares margin (Faccenna et al., 2004; Jolivet et al., 2009); and Scenario 3 starts from a S-SE dipping initial subduction zone under African margin (Vergés and Fernàndez, 2012). Dashed lines represent proposed transform fault regions. The triangular zone between the transform zones of Scenario 1 depicts the lithosphere that rolls back toward the east Kabylides (Spakman and Wortel, 2004) (from Chertova et al. (2014)) et al., 2006). The Frontal Units appear as a narrow and discontinuous band, composing of non-metamorphic Mesozoic to Tertiary sequences, overriding the Flysch Trough Units (Wildi, 1983). The formation of Gibraltar Arc System began in the Late Oligocene. The subduction system migrated westwards, during the Miocene, until the Alboran terrane collided with the Tethyan South-Iberian and Africa margins (Platt et al., 2013; Van Hinsbergen et al., 2014). The terrane stacked over Iberian and African passive margins which led to the formation of the Inner Zone of the Betics and the Rif mountain ranges. During the same time, tectonic extensional processes also took place within the Alboran Basin (Watts et al., 1993; Comas and Soto, 1999). Several geological processes have been proposed in order to explain this extensional process such as subduction and roll-back models (Royden, 1993; Blanco and Spakman, 1993; Wortel, 2000; Gutscher et al., 2002; Faccenna et al., 2004; Vergés and Fernàndez, 2012; Platt et al., 2013; Chertova et al., 2014; Van Hinsbergen et al., 2014; Casciello et al., 2015; Spakman et al., 2018), convective removal (Platt and Vissers, 1989; Dewey et al., 1989), and lithospheric mantle delamination (Platt and Vissers, 1989; García-Dueñas et al., 1992; Seber et al., 1996; Calvert et al., 2000; Valera et al., 2011; Williams and Platt, 2018). Although, the Gibraltar Arc has been related with a subduction system, as inferred from the tomographic imaging of a hanging lithospheric slab underneath (Wortel, 2000; Garcia-Castellanos and Villaseñor, 2011; Bezada et al., 2013; Bonnin et al., 2014; Villaseñor et al., 2015; Civiero et al., 2018). In this context, Chertova et al. (2014) 2.4 Western Mediterranean 23 Fig 2.9 Schematic geological map of the westernmost Mediterranean highlighting the position of the Alboran domain enclosed between Iberia and north-west Africa. Its metamorphic rocks compose the floor of the western Alboran Sea and form the inner zones of the Betics and Rif mountain chains (from Casciello et al. (2015)). used a set of geodynamic numerical models to illustrate that slab roll-back and lithospheric slab tearing at the edges of the subduction zone could be the most plausible cause for the tectonic extension. Along with the subduction system, volcanic activity also occurred in the basin and three episodes have been proposed: (1) Late Oligocene-Early Miocene tholeiitic dykes observed exclusively onshore (Torres-Roldan et al., 1986; Marchesi et al., 2012), (2) calk-alkaline and tholeiitic volcanism in the central and eastern part of the Alboran Basin (Duggen et al., 2008), and (3) Messinian-Quaternary alkaline volcanism at the basin margins which supports the asthenospheric mantle influx (Duggen et al., 2004, 2008). The subduction activity terminated or greatly decreased during the Messinian ( ∼6− 7 Ma), giving rise to a new geodynamic scenario (Iribarren et al., 2007). The tectonic extensional processes ended offshore and the area is now experiencing tectonic shortening 24 Fundamentals & Geological setting due to the European and African plates convergence, which in turn causes tectonic inversion of several fault systems within the area (Martínez-García et al., 2013; Giaconia et al., 2015). The surface deformation in the Gibraltar Arc System may have resulted from plate subduction, as corroborated by detailed geodetic data (Fadil et al., 2006; Pérouse et al., 2010; Palano et al., 2013) and numerical modelling (Jiménez-Munt and Negredo, 2003; Spakman et al., 2018). Despite previous efforts, what drives the present-day surface deformation still remains poorly constrained. Variation in surface topography and geodetic displacements may have been results of processes at different depth: (1) shallow processes, e.g. tectonic deformation due to plates’ relative motion, sedimentation and erosion, all of which could affect the crustal or lithospheric thickness. The consequent isostatic compensation may lead to an elevation change which translates to surface uplift or subsidence in the geodetic signal; (2) deep processes, e.g. mantle flow induces stress at the lithosphere’s base causing both vertical and horizontal surface deformation (Hager et al., 1985; Faccenna et al., 2014b). The observation of a positive seismic velocity anomaly in the upper mantle of the western Mediterranean was first interpreted by Spakman (1986) to show a deeply subducted slab. Since then, the advance in seismic tomography has produced several tomographic models of the western Mediterranean, all of which exhibit a positive P-wave velocity anomalies extending from under the Rif-Gibraltar towards north-east and east direction under the eastern Betics (e.g. Spakman and Wortel (2004); Bezada et al. (2013)). Gutscher et al. (2012) interpreted a decade worth of tomographic data and provided support for a narrow east-dipping subduction zone under the young extensional marine basin surrounded by an arcuate fold-and-thrust belt (Africa–Eurasia plate boundary). The tomography of the upper mantle in the region reveals a steep and east-dipping high Pwave velocity body beneath the Gibraltar. The 3D morphology of this subducted slab was analysed by Spakman and Wortel (2004), which was named Rif-Gibraltar-Betic (RGB) slab to emphasise its lateral extent (Fig 2.10). The key features of the RGB slab are the easterly and steeply dipping slab under the Gibraltar Strait and Alboran basin, and its curved slab morphology which extends from under the north-western Rif towards the north and then curves eastwards beneath the Betics. The tomographic interpretation also show that the hanging slab is partially torn under the eastern Betics (Fig 2.10). The slab appears to be continuous with the Moroccan Rif lithosphere and also slab partly attached to Iberian lithosphere under the central Betics. The Mesozoic extensional episodes in southern Iberia were followed by Cenozoic subduction. The restoration works of seismic transects by (Pedrera et al., 2020) reveal an asymmetric rift with shallow sub-lithospheric mantle beneath the sector with maximum lithospheric thin- 2.4 Western Mediterranean 25 ning. The subduction, plausibly, nucleated in these thin lithospheric sectors, i.e. the exhumed mantle domain. The ongoing lithospheric shortening caused the deformation to migrate towards the main thrusts of the External Zone, and also caused the burial of the leading edge of the subducted Iberian margin. This continental subduction induced the crustal thickening, as portrayed in the gravity models beneath Sierra Nevada and neighbouring area, which combined with erosion and crustal extensional deformation finally exposed high-pressure rocks at the surface (Pedrera et al., 2020). The horizontal forces from the Iberia-Africa convergence, which coincide with shear forces at the base of the crust, are the key driving force behind the continental subduction in the Central Betics. In the Western Betics, the continental lithosphere of the Iberian plate has been imaged at 100 km-depth (Morales et al., 1999), while the attached oceanic lithospheric slab has been imaged down to 600 km with seismic topographies (Spakman and Wortel, 2004). Therefore, the continental subduction occurred after an oceanic subduction (Morales et al., 1999; Ruiz-Constn et al., 2012). Towards the east, the negatively buoyant oceanic lithospheric slab also induced the south-east subduction of the South-Iberian continental paleomargin beneath the future Internal Zone, which in turn causing the subduction of the adjoining part of the Iberian margin in the Central Betics (Pedrera et al., 2020). Geophysical-petrological modelling also infer a deep sub-lithospheric anomaly related with the slab beneath the Gibraltar arc system. Jiménez-Munt et al. (2019) explain the observed gravity and geoid in Gibraltar by a denser deep anomaly in agreement with the modelled mantle seismic velocities and the regional elevation. Kumar et al. (2021) found a detached cold sublithospheric anomaly beneath the eastern Betics, which agrees with the slab tearing proposed by seismic tomographic models. 2.4.3 Messinian Salinity Crisis The Messinian Salinity Crisis (MSC) is an event occurred between 5.96−5.33 Ma during which the Mediterranean Sea was repeatedly cut off from the Atlantic Ocean (Hsü et al., 1973; Krijgsman et al., 1999). The dessication by evaporation during the MSC resulted in an extensive salt precipitation across the Mediterranean basin, reaching 2−3 km thickness in some regions. Since the birth of the concept (Hsü et al., 1973), MSC is considered as the largest and the most geologically abruptive environmental changes during the Neogene. The sedimentary record of the MSC displays complex feedback between climate, geomorphology, geodynamics, plate tectonics, and biota. The closure of marine gateways across the Gibraltar Strait (i.e. the cause of MSC) has been suggested to have a tectonic origin. A global sea level changes, during MSC, has been 32 Methods where T is temperature [ ◦ C], t is time [s], α is thermal diffusivity [m 2 /s], and v is velocity (plate convergence rate)[m/s]. The first term on the right side of Equation 3.1 is related to thermal diffusion and the second term is related to thermal advection. Equation 3.1 is solved with an explicit finite-difference scheme on a rectangular grid with a resolution of 5×5 km. 3.1.1.2 Density distribution The initial density distribution is calculated by considering a depth-dependent density according to: dρ=∂ρ ∂PdP +∂ρ ∂TdT (3.2) where the pressure derivative ∂ρ ∂P and the temperature derivative ∂ρ ∂T for each lithosphere type are calculated by computing the stable mineral assemblages from major oxide compositions using Perple_X (Connolly, 2005; Afonso et al., 2008) and averaging the obtained values of the derivatives over the lithospheric mantle thickness. The considered lithospheric mantle types correspond to Arc_3 from the Slave craton, Pr_1 garnet-average, and Tc_1 garnet-average (Griffin et al., 2009), with their respective thicknesses to calculate representative density contrast at the LAB for Archean, Proterozoic, and Phanerozoic lithospheres. Density variations are recalculated for each moving particle according to Equation 3.2. 3.1.2 Numerical scheme The methodology presented here has been implemented in a 2D finite-difference code LithBuoy. The Python code divides the numerical domain into evenly distributed rectangular grids, with different nodes/area assigned as either the lithospheric mantle or the asthenospheric mantle. The simplest numerical approach is to simulate the process using a model consisting only the lithospheric mantle, to which we impose a pressure-temperature dependent density, overlaying the asthenosphere of depth-dependent density, which is almost constant since the asthenosphere is not moving. The thermal diffusive and advective components will be calculated separately during each time-step. 3.1.2.1 Thermal diffusion The thermal diffusive term from Equation 3.1 is: 3.1 2D kinematic numerical modelling 33 ∂T ∂t=∇·(α∇T)(3.3a) α=k ρCp (3.3b) where α is thermal diffusivity [m 2 /s], k is thermal conductivity [W/(m · K)], ρ is density [kg/m 3 ], and Cp is specific heat capacity [J/kg · K)]. This lead to the 2D thermal diffusion equation: ρCp ∂T ∂t=∇·(k∇T)(3.4) or in a 2D discretised form: ∇·(k∇T)≈kTi+1,j−2Ti,j+Ti−1,j (∆x)2+kTi,j+1−2Ti,j+Ti,j−1 (∆y)2(3.5a) ρCp ∂T ∂t≈ρCp Tn+1 i,j−Tn i,j ∆t!(3.5b) The fully discretised form of 2D thermal diffusion equation (Equation 3.4) becomes: ρCp Tn+1 i,j−Tn i,j ∆t!=kTi+1,j−2Ti,j+Ti−1,j (∆x)2+kTi,j+1−2Ti,j+Ti,j−1 (∆y)2(3.6) In the new time-step update, we can obtain the new temperature for each node through: Tn+1 i,j=Tn i,j+α∆tTn i+1,j−2Tn i,j+Tn i−1,j (∆x)2+α∆tTn i,j+1−2Tn i,j+Tn i,j−1 (∆y)2(3.7) where (i,j)=(x,y) are spatial components and n is the temporal component (i.e. n+1= new and n=old). 3.1.2.2 Thermal boundary conditions The numerical modelling code is implemented with ‘insulating boundary conditions’, meaning that the heat does not flux through the boundary (i.e the temperature gradient does not exist across this boundary): qx=−k∂T ∂x=∂T ∂x=0 (3.8) 34 Methods Utilising central-difference approximation together with the insulating boundary condition, we obtain: dT dx =Tn i+1,j−Tn i−1,j (∆x)2=0 (3.9a) Tn i+1,j=Tn i−1,j(3.9b) The boundary condition is imposed laterally (x-direction) only. Temperature within the lithospheric mantle is calculated with boundary conditions at the Moho ( Tmoho ) and the LAB ( TLAB ), and temperature from the adiabatic gradient in the subducting slab boundary ( Tadiabatic ) (Fig 3.1). Temperature within the asthenosphere is fixed with the adiabatic gradient dT/dy =0.5◦C/km. 3.1.2.3 Advection The terms advection and convection are often used as a synonym of each other in literature. Convection applies to the movement of a fluid (due to density gradients created by thermal gradients), whereas advection refers to the movement of some material/bulk quantity by the fluid velocity field. In this numerical code we use advection equation which is the partial differential equation that governs the motion of a conserved scalar field (temperature, pressure, and density in this case) as it is advected by a known velocity vector field. Assuming that the medium with incompressible flow ( ∇·v=0 ), the advection equation is formulated as: ∂ψ ∂t+v·∇ψ=0 (3.10) where ψ is the scalar field being advected (temperature ( T ), pressure ( P ), and density ( ρ )), andvis the velocity vector. Since our model is 2D, the equation above becomes: ∂ψ ∂t+vx ∂ψ ∂x=0 (3.11a) ∂ψ ∂t+vy ∂ψ ∂y=0 (3.11b) The discretisation of the advection equation is implemented using the 2D first-order upwind scheme in both xand y-direction. The advection equation, of temperature as an 3.1 2D kinematic numerical modelling 35 example, in the discretised form (Figure 3.1) are: If vx>0 : Tn+1 i,j=Tn i,j+∆vx(i,j) Tn i−1,j−Tn i,j ∆x(3.12a) If vx<0 : Tn+1 i,j=Tn i,j+∆vx(i,j) Tn i,j−Tn i+1,j ∆x(3.12b) If vy>0 : Tn+1 i,j=Tn i,j+∆vy(i,j) Tn i,j−1−Tn i,j ∆x(3.12c) If vy<0 : Tn+1 i,j=Tn i,j+∆vy(i,j) Tn i,j−Tn i,j+1 ∆x(3.12d) The numerical stability of the upwind scheme follows the Courant-Friedrichs-Lewy condition (Hirsch, 2007), which states: c= v∆t ∆x≤1 (3.13) 3.1.2.4 Slab Buoyancy The key output of LithBuoy code is the buoyancy force ( Fb ) of the lithospheric mantle. Specifically, the buoyancy force of the down-going part of the lithospheric mantle. Since this is a 2D exercise, areal density is formulated as: ρ=m A(3.14) where ρ is the average density [kg/m 3 ], m is the total mass of the object [kg], and A is total area of the object. The buoyancy force at each time-step is calculated as the integral of the density changes relative to the initial stage over the entire subducting slab, Fb=Zg∆ρdxdy (3.15a) where g is gravity and ∆ρ is the density difference between the slab and the asthenosphere. The buoyancy force, measured in N/m, results from two components: diffusive ( Fd ) and advective ( Fa ), such that the total Fb=Fd+Fa is directed downwards when Fb<0 (corresponding with slab pull). 36 Methods Fig 3.2 Schematic drawing of the implementation of Natural Boundary Conditions. Within the domain, internal diffusion is calculated using 5-points stencil, highlighted in red. At the boundaries, Central Difference approximation together with Insulating Boundary Condition is applied laterally, highlighted in green. 3.1.3 Workflow Currently, the code contain database for the 5 types of lithospheric mantle we will study (Archon, Proton, Tecton, 30-Myr Oceanic, and 120-Myr Oceanic). The parameters that can be input into the code are listed in Table 3.2. The schematic model box is shown in Fig 3.2. The simplified workflow (Fig 3.3) of our 2D kinematic numerical modelling code is outlined below: 1. The 2D initial temperature distribution is created according to each lithospheric mantle type, lithospheric thickness, and Moho temperature. Subsequently, the initial density and pressure distributions are then calculated. 2. The temperature nodes are advected down along the prescribed subduction plane with a prescribed constant velocity. 3. A new temperature distribution is obtained after thermal diffusion and advection routines. 4. New density and pressure distributions were calculated from the diffused and advected temperature distribution. 5. The new density distribution is then used to calculate the buoyancy force ( Fb ) of the lithospheric mantle. 6. The new temperature, density, and pressure distribution are passed onto the next time-step. 3.1 2D kinematic numerical modelling 37 LithBuoy •Initial temperature distribution T(mantle type & thickness) •Compute initial density (ρ), pressure (P) distributions Advection & Diffusion of T-nodes New T distribution Compute new ρand Pfrom the new Tdistribution Use new ρto compute buoyancy force (Fb) New T, ρ, P pass to next time-step •Mantle’s physical properties database •Thermodynamics database from Perple_X 𝜕𝜌 𝜕𝑇 ,𝜕𝜌 𝜕𝑃 , 𝑘, 𝜅 Fig 3.3 A flowchart briefly describing the workflow of the LithBuoy code (Boonma et al., 2019). 38 Methods Table 3.2 The input parameters for LithBuoy numerical code. Input parameters Shorthand Unit Model Geometry Box height hm Box width wm Discretisation in x-plane dx m Discretisation in y-plane dy m Subduction angle angle (φin Fig 3.2) degree Position (in x-direction) of subduction hinge hinge_ax m Convergence velocity vel (v in Fig 3.2) mm/year Thermal parameters Thermal conductivity (lithosphere) k_lith W/(m·K) Specific heat capacity Cp J/(kg·K) Thermal diffusivity (lithosphere) kappa_lith m2/s Thermal diffusivity (asthenosphere) kappa_asth m2/s Mantle database Moho depth d_moho m LAB depth d_lab m Moho temperature TMoho ◦C LAB temperature TLAB ◦C Temperature at the bottom of the model box Tbottom ◦C 3.2 3D thermomechanical numerical modelling 39 3.2 3D thermomechanical numerical modelling We decided to use I3ELVIS code for the 3D thermomechanical numerical modelling of slab tearing because of its ability to model complex non-linear multi-physics problems such those we face in geodynamics. It has a large rheological database and is completed with calculation routines for various geodynamical processes. The code is developed by Gerya et al. (2013) at ETH Zürich. 3.2.1 Governing equations In I3ELVIS, the physical principles implemented in the code are in Lagrangian frame of reference (denoted by D/Dt ), meaning that we, an observer, follow an individual parcel of quantity as it moves through space and time. Interlinking with Lagrangian view point is the Eulerian frame of reference (denoted by ∂/∂t ). In an Eulerian viewpoint, an observer stays or focuses on a specific location in space through which the quantity flows as time passes. 3.2.1.1 Conservation of mass (Continuity equation) The conservation of mass can be mathematically expressed in a differential form as the continuity equation as follow: ∂ρ ∂t+▽·(vρ) = 0 (3.16) where ρ is density [kg/m 3 ], t is time [s], and v is the flow velocity vector field [m/s]. In fluid mechanics, the continuity equation states that, in a steady-state system, mass leaves a system at the same rate as at which mass enters the system. Therefore, the continuity equation can also be expressed in a Lagrangian time derivative form as: Dρ Dt +ρ(▽·v) = 0 (3.17) Assuming that the continuous medium is incompressible (i.e. the density of mobile points does not change with time, ▽·v=0 ), a common assumption in the studies of lithospheric dynamics (since material volume changes due to pressure and temperature variation is insignificant), the continuity equation can be simplified to: ∂vi ∂xi =0 (3.18) 40 Methods where i (as well as j shown hereafter) is index, which denotes spatial directions i= (x,y,z) in 3D using Einstein notation. 3.2.1.2 Conservation of momentum (Stokes equation) An equation of motion is used to relate the various internal and external forces that act on continuous media to the consequent deformation. The momentum equation describes the conservation of momentum for a continuous medium in the gravity field (Gerya, 2009). In high-viscosity flow, the inertial forces become much smaller than the body force (gravity) and therefore can be neglected, which leads to the general formulation of momentum equation: ∂σ′ i j ∂xj −∂Pi ∂xi =−ρ(T,P,c)gi(3.19) where σ′ i j is deviatoric stress tensor (i.e. deviations of stresses from the hydrostatic stress state), Pis the total pressure (i.e. mean normal stress), gis gravitational acceleration. In the numerical modelling of fluid dynamics, the fluid temperature is allowed to vary from one position to another, driving the fluid flow. This also allows the Boussinesq approximation to be imposed onto our buoyancy-driven flow problem. The approximation only takes the density difference that is multiplied by gravity into consideration, and the rest is assumed to be constant. This is because in the buoyancy force term, the temperature and volatile content play vital parts. Therefore, the density in the buoyancy force term, ρg , is allowed to be varied locally as a function of temperature (T), pressure (P), and composition (c). 3.2.1.3 Heat Transport equation Heat transport equation arose from the law of energy conservation under the consideration of slowly moving system or external work with the ‘energy’ being a product between specific heat ( Cp ) and temperature ( T ). A few assumptions are also in place: (i) Cp is constant; and (ii) mass conservation equals zero. Thus, the obtained Lagrangian heat transport equation reads the following, ρCp DT Dt =∂ ∂xik∂T ∂xi+Hr+Ha+HL+Hs(3.20) where Cp is the heat capacity [ J/K ], T is temperature [K], k denotes thermal conductivity [W/(m · K], which is a function of temperature, pressure, and composition, i.e. k(T, P, c) , Hr is radioactive heat production, Ha is adiabatic heat production/consumption, HL is latent heat 3.2 3D thermomechanical numerical modelling 41 production/consumption (due to phase transformation), and Hs is the shear heat production (a product of deviatoric stress and strain rate). Hr=constant (3.21a) Ha=Tαv∇P(3.21b) HL=constant (3.21c) Hs=σi j ˙ εi j (3.21d) In adiabatic heat source ( Ha ), T denotes temperature [K] and α denotes thermal expansion [1/K]. The shear heat ( Hs ) is composed of components of the deviatoric stress tensor ( σi j ) and strain rates ( ˙ εi j ). This shows that the conservation laws is coupled with rheology laws through the strain rate ( ˙ εi j ). The combination of the set of Equation 3.21 and Equation 3.20 , together, form the Extended Boussinesq Approximations. Lagrangian temperature equation is solved on the Eulerian nodes, and temperature increments are interpolated from nodes to markers by using the subgrid diffusion operation (Gerya and Yuen, 2003, 2007; Gerya, 2009), which can ensure physical consistence between nodal and marker thermal fields. Advection of temperature is implemented through marker advection. 3.2.1.4 Density model The rocks’ densities vary with temperature T [K] and pressure P [Pa] according to the equation of state: ρP,T=ρ0[1−α(T−T0)+β(P−P0)] (3.22) where ρ0 is the reference density at P0=1 MPa and T0=298.15 K, the coefficient of thermal expansion α=2×10−5 1/K, and the coefficient of thermal compressibility β=6×10−12 1/Pa. Our models take phase transition of olivine in the mantle into account. As the dry olivine is subjected to greater pressure at depths, it first undergoes exothermic phase transition ( ∼ 410 km) and transforms into wadsleyite (Katsura and Ito, 1989). At a greater depth and pressure, the wadsleyite exothermically transforms into ringwoodite ( ∼ 520 km), which decompose (endothermically) into bridgmanite (silicate perovskite) at an even greater depth ( ∼ 660 km) (Ito et al., 1990). The eclogitization of the subducted oceanic crust (basaltic and gabbroic) is taken into account by linearly increasing the density of the crust with pressure 48 Methods Fig 3.6 A 2D sketch to show an example of the velocity and thermal boundary conditions used in our 3D model. More specific model modification details, for our purpose, are described in Chapter 5. surface and fill them with sticky air, i.e. low density (1 kg/m 3 ) and high viscosity ( 1018 Pa · s) material (but sufficiently lower viscosity compared to the lithospheric rocks). Viscosity of sticky air is prescribed as the same value of lower limit of viscosity cut-off in the model, and it is several orders of magnitude smaller than viscosity of upper crust. Large viscosity contrast along crustal surface ensures that stress in sticky air layer is sufficiently small. (3) Permeable lower boundary – is prescribed as an infinite-like external free-slip condition at a distance away from the bottom of the model. The external free-slip permits the global mass conservation in the computational domain and is implemented as ∂vx/∂y=0 ; ∂vy/∂y= −vy/∆yexternal , where ∆yexternal is the vertical distance between the bottom of the model box and the bottom of the external boundary where the free-slip condition ∂vx/∂y=0 , vy=0 is satisfied (Fig 3.6). (4) Fast erosion condition – corresponds to a fast erosion/sedimentation at the upper free surface. This condition ensures the conservation of mass in the model. An instantaneous sedimentation layer limits a trench depth to 8 km below the water level and the an instantaneous erosion layer is prescribed at 8 km above the initial continental crustal surface where rock markers change into sticky-air markers. 3.2 3D thermomechanical numerical modelling 49 3.2.3.2 Thermal boundary conditions Thermal boundary conditions have to be specified in order to numerically solve the temperature equation (Equation 3.20). In our numerical modelling, we prescribed a constant temperature condition on the top and bottom boundaries and (1) Constant temperature condition – the temperature at a boundary is assigned a given value (which may change spatially both along the boundary and in time), T=const(x,y,z,t)(3.42) In our models, the constant temperature is prescribed as following: (i) the upper boundary (i.e. crustal surface) is prescribed as 0◦ C (273 K); (ii) the lower boundary is prescribed as the infinite-like external constant temperature (Fig 3.6), which is implemented as ∂T/∂y= (Texternal −T)/∆yexternal , where Texternal is the temperature at the bottom external boundary, and ∆yexternal is the vertical distance between the bottom of the model and the bottom external boundary where the condition T=Texternal is satisfied. (2) Insulating boundary condition – the heat does not flux through a boundary which implies that no temperature gradient exists across the lateral boundary, i.e. 0=qx=−k∂T ∂x⇒∂T ∂x=0(3.43) 3.2.4 Computational strategy The I3ELVIS (3D visco-plastic rheology, no elasticity) (Gerya, 2013) thermomechanical numerical modelling code uses conservative finite-difference numerical method and a markerin-cell technique to simulate solid phase flow. The governing equations and rheologies are solved using a regularly-spaced staggered Eulerian grid (Fig 3.7). Using markers to solve advection equations can often cause numerical problems if solved by various Eulerian advection schemes (Gerya, 2009). Marker-in-cell technique means that a certain amount of markers are prescribed in Eulerian grid cells. Since Eulerian grids and nodes are fixed, moving Lagrangian markers are used to advect temperature, as well as other physical properties such as viscosity, density, heat conductivity, and etc. Computed physical variables (by solving the governing equations) on Eulerian nodes are then interpolated from nodes to markers by using weighted-distance average method (Gerya and Yuen, 2003), and using Runge-Kutta (fourth-order accuracy in space) method to move the markers. Once the markers have been displaced, physical variables will be interpolated from markers to nodes. 50 Methods The 3D numerical modelling code utilises two main types of solver: (i) direct (e.g. Gaussian elimination) and (ii) iterative method (e.g. Gauss-Seidel), both of which are commonly used to solve large systems of linear equations. The direct method produces more accurate solutions compared to iterative method, but it requires large amount of consumed memory (typically proportional to the square of the number of unknowns) and large amount of operations (typically proportional to the square or even to the cube of number of unknowns). Therefore, the direct method is more often used in 1D or 2D numerical models as they are less computationally costly. The iterative methods are commonly used in 3D numerical models, since model has one more dimension and therefore number of equations becomes very large. An initial guess is given for each unknown in iterative method, and computed residues need to fulfil a certain criterion. Despite giving less accurate results, the iterative methods have two advantages: requiring small amount of consumed memory (typically proportional to the number of unknowns) and small amount of operations (typically proportional to the number of unknowns per solution cycle). Multi-grid method is used to speed up the convergence of the Gauss-Seidel iterations for coupled solving of mass and momentum conservation equations. 3.3 Computational demands 2D numerical modelling The 2D kinematic modelling code, LithBuoy, was developed during the course of this thesis to study the dependence of lithospheric buoyancy on tectonic convergence rate and mantle composition. Both the computational code and the visualisations (2D) scripts are written in Python language. The code is open-access and available in a GitHub repository 1 . The key contents in the repository are: •params.py - main executable script/ change basic parameters •lithbuoy.py - computational sub-routines and time-stepping •mantle.py - properties of each lithospheric mantle type •sim_figs.py - plot simulation output •plotcsv.py - plot graphs/data analysis The code does not need parallelisation and can be run on a single CPU/local computer. It can also be run on a cluster. The computational time for each run is in the range of 10−120 minutes, depending on the dimensions of the model domain and its discretisation. 1https://github.com/kboonma/LithBuoy 3.3 Computational demands 51 3D numerical modelling The I3ELVIS code is highly optimised for geodynamical simulations. It is a proprietary code and can only be obtained with the authors’ authorisation. The numerical code is written in C language, and utilises OpenMP (multi-platform shared-memory multiprocessing programming) parallelisation in the computation of 3D models. The code was developed at ETH Zürich, so the parallelisation was tailored to the home cluster i.e. the 24-cores Euler cluster (Technical specification as of 2019 are: 2 × 12-core Intel Xeon E5-2680v3 processors (2.5−3.3 GHz), between 64 and 512 GB of DDR4 memory clocked at 2133 MHz). The high resolution 3D simulations are computationally expensive. The 3D numerical models in this thesis has a vertical (depth) resolution of ∼ 2 km, which is on a higher end for a large-scale geodynamic simulation, and with that comes the long simulation time. One model can require 104−106 core hours per model, which can means model readjustments can be time-consuming if under a time constraint. For example, each 3D model presented in this thesis took 1-2 months to run. Such large computational time leads to this thesis’ weakness being the lack of the full and detailed model parameterisation. Data processing (3D visualisations and analysis) can require high-performance computing (HPC) resources. The main software used for 3D visualisation is ParaView (open-source, Ayachit (2015)) and MATLAB is used for 2D visualisation. 52 Methods Fig 3.7 Elementary volume (cell) of a 3D staggered grid used for discretization of momentum, continuity, Poisson and temperature equations in the case of incompressible viscous flow with variable viscosity and thermal conductivity (from Gerya (2009)). CHAPTER 4 LITHOSPHERIC MANTLE BUOYANCY:THE ROLE OF TECTONIC CONVERGENCE AND MANTLE COMPOSITION The work in this chapter has been published in: Boonma, K., Kumar, A., Garcia-Castellanos, D., Jiménez-Munt, I., and Fernández, M. (2019). Lithospheric mantle buoyancy: the role of tectonic convergence and mantle composition. Scientific Reports, 9(1), 17953. DOI: 10.1038/s41598-019-54374-w 4.1 Introduction Plate tectonics is thought to be mainly driven by the negative buoyancy of the lithospheric mantle relative to the asthenosphere, the driving force for both oceanic plate subduction and mantle delamination (the peeling off of the subcontinental lithospheric mantle from the crust and its sinking). The lithosphere interacts differently with the underlying asthenosphere in oceanic and continental domains. Oceanic lithosphere is formed at mid-ocean ridges and is soon (generally in less than 200 Ma) reworked back into the deeper mantle through subduction. In contrast, continental lithosphere is an order of magnitude older and has grown through accretion over longer time-scales, modifying its chemical composition. The continental lithosphere is believed to be recycled back into the asthenosphere at collision zones through the removal of parts of the lithospheric mantle by delaminating from the overlying crust. Geophysical and geological observations support mantle delamination in regions such as the Tibet (Chung et al., 2005; Ren and Shen, 2008); Alboran domain 54 Lithospheric buoyancy: the role of tectonic convergence and mantle composition (Seber et al., 1996; Valera et al., 2011; Thurner et al., 2014); the Apennines (Bartolini, 2003; Carminati et al., 2004; Chiarabba and Chiodini, 2013); Eastern Anatolia (Al-Lazki et al., 2003; Gö˘ gü¸s and Pysklywec, 2008; ¸Sengör et al., 2008); and Sierra Nevada (Saleeby and Foster, 2004; Valera et al., 2014), some of which are shown in Fig 4.1. Numerical models suggest that this process is initiated by a conduit connecting the sublithospheric mantle with the weak lower crust and the negative buoyancy of the lithospheric mantle (Bird, 1979). Fig 4.1 Geographic overview of regions with continental collision. Three regions within continental collisions that can be related to our study according to their convergence-rate and lithospheric mantle thickness. Whereas most geodynamic studies have focused on the role of the viscosity contrast between the lower crust and the lithospheric mantle in the development of delamination (Bird, 1979; Schott and Schmeling, 1998; Valera et al., 2011; Bajolet et al., 2012; Gray and Pysklywec, 2012; Magni et al., 2013), or on the role of composition and thickness of the crust on the lithospheric subductability (Cloos, 1993), only few of them have scrutinized the influence of the density contrast between the lithospheric mantle and the asthenosphere (Channell and Mareschal, 1989; Afonso et al., 2007; Gö˘ gü¸s and Pysklywec, 2008; Valera et al., 2008; Afonso and Zlotnik, 2011; Valera et al., 2011, 2014). These geodynamic models generally adopt densities that are either constant over time and space or temperaturedependent only, predefining a density for the subcontinental lithospheric mantle (SCLM) about 50 kg/m 3 higher than of the asthenosphere. More recently, the role of pressure, temperature, and composition has been incorporated in calculating the buoyancy of the 4.1 Introduction 55 descending lithospheric mantle forced by plate convergence (Afonso et al., 2008; Li et al., 2016). Global geochemical analyses on mantle xenoliths, xenocrysts, and outcrops show that the mean composition of the SCLM is mainly related to the tectonothermal age of the overlying crust (Griffin et al., 1999, 2009). As the continental lithosphere undergoes cycles of melting, it gradually depletes in incompatible elements, such as Ca , Al , and Fe , relative to the primary source. The lithospheric mantle in Archean cratons (Archons), corresponding to tectonothermal ages > 2.5 Ga, is generally highly depleted, while most lithospheric mantle beneath Neo-proterozoic/Phanerozoic mobile belts (Tectons, < 1 Ga) is closer to the composition of the Primitive Upper Mantle (PUM). The SCLM beneath Mesoand Paleo-proterozoic shields and mobile belts (Protons), with tectonothermal ages of 2.5−1.0 Ga, is typically intermediate in composition. These compositional variations affect the bulk density of the SCLM and the greater the degree of depletion, the lower the density (Griffin et al., 2009). In contrast, the composition of the oceanic lithospheric mantle corresponds to that of PUM after melt extraction at mid-ocean ridges (MOR), being relatively homogeneous except beneath large oceanic plateaus. Petrological and geochemical studies show that at identical P-T conditions as at the LAB, the density of the SCLM is lower than that of the PUM (Griffin et al., 2009) (Fig 4.2), which is at odds with the aforementioned density distribution adopted in most geodynamic and static models (Cloos, 1993). In fact, the density of the sinking SCLM increases with pressure and decreases with temperature, thus P and T having competing effects on the depth-dependence of density. Whether mantle delamination or subduction are promoted by a negative buoyancy forced by plate convergence depends on the gradients of density relative to both parameters (Table 4.1). Here we explore the idea that, whenever the lithosphere is incipiently forced to sink into the asthenosphere, it can become positively or negatively buoyant depending on its composition and on how the pressure-temperature evolution, imposed by plate convergence, affects its density evolution. This may explain why the older regions of the continental lithosphere (Archons) become tectonically more stable, self-prolonging their life span and favouring the Wilson cycle. To this purpose, we calculate the buoyancy of the sinking mantle in a kinematic model that accounts for a thermodynamically-consistent density dependence on temperature, pressure, and chemical composition. 56 Lithospheric buoyancy: the role of tectonic convergence and mantle composition Table 4.1 Physical parameters of all of the lithosphere types used in this study, together with those of the Primitive Upper Mantle (PUM). Lithospheric Mantle type (LM) Continental Lithosphere TMOHO =650◦C Oceanic Lithosphere TMOHO =300◦CPUM Archon (Arc_1) Proton (Pr_1) Tecton (Tc_1) OC120 Ma OC30 Ma Lithospheric mantle thickness (km) 160 110 80 100 60 - LAB depth (km) 200 150 120 110 70 - ∆ρLAB (kg/m3)+68 +39 +19 +17 +17 - Mg# (mantle fertility) 92.0 90.5 89.9 89.7 89.7 89.3 dρ/dP(kg/m3MPa) 0.042 0.043 0.044 0.046 0.046 0.044 dρ/dT(kg/m3K) -0.1302 -0.1285 -0.1255 -0.1236 -0.1236 -0.1165 Density at 3GP/800◦C (kg/m3)3319 3343 3365 3362 3362 3379 Density at 6GP/1300◦C (kg/m3)3341 3367 3391 3389 3389 3407 4.2 Model setup To study all these different parameters that affect the density of the slab and, in turn, its buoyancy, I developed the 2D kinematic numerical modelling code, LithBuoy (Boonma et al., 2019), in which the basic principles are described in Chapter 3.1. The approach is to model plate convergence and subduction of the lithospheric mantle. The model domain has a width ( w ) of 1500 km and a height ( h ) of 600 km, with nodes arranged on a rectangular grid with a resolution of 5×5 km. The down-going lithospheric mantle is submitted to a prescribed velocity field ( v=1,4,10,20,40, and 80 mm/yr), which is directed downward with an angle of φ=30◦(Fig 3.1 and 4.3). The model’s uppermost boundary is the Moho depth (40 km for continental lithosphere and 10 km for oceanic lithosphere), which has a fixed temperature of TMoho =650◦ C for continental lithosphere and TMoho =300◦ C for oceanic lithosphere. The next thermal boundary is the lithosphere-asthenosphere boundary (LAB) where the temperature is TLAB =1300◦ C. The initial temperature distribution within the lithosphere is interpolated between the TMoho and the TLAB across the thickness of each lithosphere type (Table 4.1), while the adiabatic temperature gradient in the asthenosphere is fixed to dT/dy =0.5◦ C/km (Fig 4.3). The details on the calculations and boundary conditions are described in Section 3.1. 4.3 Results 57 Fig 4.2 Temperature (a) and density (b) profiles for the five lithosphere types considered. The temperature is calculated by solving for the steady-state heat flow equation with conductive regime in the lithosphere and adiabatic gradient in the asthenosphere. Thermal conductivity and density in the mantle are calculated by solving for stable mineral assemblages from major oxides compositions using Perple_X . It should be noted that these profiles are not the initial profiles used in the models, but they are the reference profiles. 4.3 Results As lithospheric convergence initiates, the downward advection of cold lithospheric material generates a negative temperature anomaly that extends along the lithospheric slab, from the axial plane (Fig 4.3) to a growing depth below the undisturbed LAB, depending on the convergence velocity. Simultaneously, thermal diffusion is responsible for the reduction of that negative temperature anomaly by transferring heat from the asthenosphere to the subducting slab. Therefore, heat advection and thermal diffusion have opposite effects on temperature anomaly and buoyancy. The prevalence of each depends on the convergence rate, and the composition/age of the lithospheric mantle and its thickness. CHAPTER 5 THERMOMECHANICAL MODELLING OF LITHOSPHERIC SLAB TEARING AND ITS TOPOGRAPHIC RESPONSE IN THE GIBRALTAR ARC (WESTERN MEDITERRANEAN SEA) 5.1 Introduction The perception that large regions of continental crust have risen at rates that cannot be explained by crustal thickening alone, has led to the necessity to identify the mechanism responsible for such rapid surface uplift (England and Molnar, 1990). Slab breakoff is among the deep-seated mechanisms invoked to justify the high rate of surface uplift (Davies and von Blanckenburg, 1995). It is driven by the same force that drives slab pull and subduction, i.e., a positive contrast of the potential density of the lithospheric slab relative to the mantle (e.g. Garcia-Castellanos et al. (2000); Boonma et al. (2019); Jiménez-Munt et al. (2019)). Slab breakoff is a process happening at depth within the mantle consisting of the detachment of a subducted oceanic lithospheric slab from the more buoyant continental lithosphere during a continental collision. The concept of slab breakoff, as inferred from seismic tomography, was first proposed to be involved in the geodynamical evolution of the Mediterranean by Wortel and Spakman (1992). Slab breakoff was then used to explain post-collisional magmatism and exhumation of high-pressure rocks in the European Alps by 66 Thermomechanical modelling of slab tearing and its topographic response Davies and von Blanckenburg (1995). Garzanti et al. (2018), and references therein, gave a comprehensive global overview of where slab breakoff has been invoked to explain changes in plate kinematics and deformations e.g. the Alps (von Blanckenburg and Davies, 1995), the Mediterranean region (Carminati et al., 1998; Wortel, 2000; Rosenbaum et al., 2008; van Hinsbergen et al., 2010), the Anatolia-Zagros orogen (¸Sengör et al., 2003; Faccenna et al., 2006), or Himalaya and Tibet (van Hinsbergen et al., 2012; Wu et al., 2014; Liang et al., 2016). These studies often ascribe short-lived and long-wavelength uplift or exhumation events, or sudden pulses in sediment supply to slab breakoff. However, they often neglect the influence of the 3D geometrical configuration of each tectonic regions. How likely was the tectonic configurations in those domains to have caused the slab tearing in the first place? How much does slab breakoff contribute to the buoyancy-driven isostatic surface uplift? Our study addresses these phenomena, taking the westernmost Mediterranean as a reference geological region. The Western Mediterranean underwent subduction, slab fragmentation, and rollback, under intense crustal deformation including simultaneous N-S compression and E-W extension. Seismic tomography studies (Wortel, 1982; Spakman and Wortel, 2004; Garcia-Castellanos and Villaseñor, 2011; Bezada et al., 2013; Bonnin et al., 2014; Villaseñor et al., 2015; Civiero et al., 2019) indicate that the entire western Mediterranean region overlies structurally complex remnants of subducted lithosphere including fragments of oceanic Tethyan lithosphere inherited from the Mesozoic extension between Eurasia and Africa. The subduction between the Iberia and Africa plates is possibly associated with the uplift of the intramountain basins within the Rif and the Betics (Fig 5.1). This uplift, in turn, has been linked to the closure of the marine gateways across the Gibraltar Arc during the Late Miocene and the potential partial desiccation of the Mediterranean (Messinian Salinity Crisis; MSC) (Garcia-Castellanos and Villaseñor, 2011). According to this hypothesis, once the marine gateways closed off due to the regional epeirogenic uplift, it triggered the massive salt accumulation of the MSC (5.96–5.33 Ma) (Garcia-Castellanos and Villaseñor, 2011), possibly the most abrupt environmental change on Earth since the beginning of the Tertiary. This study aims to better understand the lithospheric slab tearing process and the consequent surface uplift, through 3D thermomechanical modelling. We investigate how different scenarios and/or parameters act on the initiation of the slab tearing and its propagation along the trench, and quantify the resulting surface elevation changes robustly coupled with the deep geodynamic processes (viscous flow, temperature evolution, and dynamic topography). 5.2 Model setup 67 Fig 5.1 Gibraltar Arc region. A present-day schematic map of the Western Mediterranean region, specifically the Gibraltar Arc System (after Rosenbaum et al. (2002)). The inset is a topographybathymetry map, with our region of interest in the red box. The map displays key units of the Gibraltar Arc System: External Zone, Internal Zone, and the Flysch Zone (Suture Zone). The blue lines are the reconstruction of the Ligurian-Tethys domain between Iberia and Africa (Late Cretaceous) before the onset of the African convergence, as proposed by Vergés and Fernàndez (2012). The grey frame outlines the area and the basic features therein, which the model setup is based on. The model frame is at such angle to encapsulate the stages from the reconstruction that involve slab detachment i.e. from 18−0 Ma. The numbers in white rectangles are the ages in Ma of the transition from marine to continental conditions of intramountain basins within the Betics (Iribarren et al., 2009; Krijgsman et al., 2018; van der Schee et al., 2018). 5.2 Model setup The work in this chapter utilised the 3D numerical code I3ELVIS (Gerya, 2013). The code is based on finite-differences and marker-in-cell numerical schemes (Gerya and Yuen, 2003; Harlow and Welch, 1965). The governing physical laws such as conservation of mass, conservation of momentum, and heat equation are discretised on a staggered Eulerian grid, assuming an incompressible medium (i.e. ∇·v=0 ). For each time step, the markers 68 Thermomechanical modelling of slab tearing and its topographic response are spatially advected with Runge-Kutta scheme. The convergence of the Gauss-Seidel iterative solver is sped up with multi-grid method. In our models, a visco-plastic rheology is implemented. For further details about the code please refer to Chapter 3.2. The continental collision is modelled with an incoming continental block (Africa), overriding a subducting oceanic plate which is connected to a stationary continental block (Iberia) through a passive margin. The 3D model domain (Fig 5.2) measures 1500×780×1200 km, with grid resolution of 4.6×3.0×4.6 km, in the x, y (vertical), and z directions, respectively. The 40-km thick continental crust splits into the upper (20 km) and lower (20 km) continental crust, and thinning toward the ocean. The 8-km thick oceanic crust also splits into the upper (basaltic, 3 km) and lower (gabbroic, 5 km) oceanic crust. Partial melting and melt extraction processes are neglected in our simplified models. The initial adiabatic temperature gradient (0.5 ◦ C/km) is prescribed in the asthenospheric mantle. The continental geotherm is prescribed as a linear variation from the model surface (0 ◦ C, 273 K) to the lithosphere-asthenosphere boundary (1344 ◦ C, 1617 K) at 110 km depth. The initial thermal structure of the oceanic lithosphere is calculated using the half-space cooling model (e.g. Turcotte and Schubert (2014)) based on a slab age of 110 Ma and a thermal diffusivity of 10−6 m 2 /s. The thermal boundary in the lower boundary of the model is prescribed as the infinite-like external constant temperature, which is implemented by ∂T/∂y= (Texternal −T)/∆yexternal , where Texternal is 1707 ◦ C (1980 K) at 1080 km depth at the bottom external boundary (outside the model box), and ∆yexternal is the vertical distance between the bottom of the model box and the bottom external boundary where T=Texternal , in this case ∆yexternal =300 km. The velocity boundary conditions are free slip on all sides except the bottom boundary, which is permeable in both upward and downward directions. This permeable bottom boundary is prescribed as an infinite-like external free slip conditions at 1080 km depth. The external free slip permits the global mass conservation in the computational domain and is implemented as ∂vx/∂y=0 ; ∂vy/∂y=−vy/∆yexternal ; ∂vz/∂y=0 , where ∆yexternal is the vertical distance between the bottom of the model box and the bottom external boundary where the free slip condition ( ∂vx/∂y=0 , vy=0 , ∂vz/∂y=0 ) is satisfied, in this case ∆yexternal = 300 km. The elevation of the lithosphere is calculated dynamically as an internal free surface through a 22 km thick layer of ‘sticky air’ ( ηair =1018 Pa · s, ρair =1 kg/m 3 ) on top of the continental plate and 25 km on top of the oceanic plate (Gerya, 2009). We implemented a simplified erosion condition in our model, where instantaneous sedimentation limits a trench 5.2 Model setup 69 depth to 8 km below the water level and the instantaneous erosion is prescribed at 8 km above the initial continental crustal surface where rock markers change into sticky-air markers. All of the experiments are two-stage experiments. The first stage is a period of forced convergence (rate of 47 mm/yr) until the subducted slab reaches 200 km depth. We prescribed the initial convergence at x=1386 km within the two transform fault weak-zones, labelled ‘ridge’ in Fig 2c. The initial convergence rate is purposely fast in order to create a sufficient hanging slab with minimised thermal diffusion. After the first stage, the obtained thermomechanical state is used as an initial setup for continental collision. In the second stage, the prescribed convergence rate is either removed, so that the slab sinks due to its own weight (Mod1-reference, Mod2, Mod3, and Mod4), or reduced to a lower value of 4 mm/yr, resembling the convergence rate in the Western Mediterranean (Mod5). 70 Thermomechanical modelling of slab tearing and its topographic response Fig 5.2 Setup of Mod1-reference model. a) 3D model domain ( 1500×780×1200 km) with colours representing different compositions. The flow law abbreviations are Wt Qtz. - wet quartzite; Plag. - Plagioclase; and Dry Olv. - Dry Olivine. Convergence is imposed by applying a uniform velocity ( v ) of 47 mm/yr until the slab reaches 200 km depth, after which v is either controlled by the sinking slab or reduced to another constant velocity. b) A cross-section profile of the viscosity. c) A map view of the model. 5.2 Model setup 71 Table 5.1 Material properties used in the numerical experiments. The flow law include: A is the pre-exponential factor; E denotes activation energy; V is activation volume; n is the stress exponent; m is grain size exponent; σcr is critical stress or the the assumed diffusion-dislocation transition stress; C is the rock compressive strength at P=0 MPa; µ0 and µ0 are the initial and final internal coefficients, respectively parameters (Karato and Wu, 1993; Ranalli, 1995; Hirth and Kohlstedt, 2003; Turcotte and Schubert, 2014). The subscripts ‘diff’ and ‘disl’ indicate that those parameters are associated with diffusion and dislocation creep processes, respectively. Mod3 has higher values for mantle activation volume than other models. Mod4 has higher final internal friction coefficient for lower oceanic crust and the mantle. Other properties for all rock types include: heat capacity Cp=1000 J/(kg·K), thermal expansion α=2×10−51/K; and compressibility β=6×10−12 1/Pa. Material Density ρ0 (kg/m3) Thermal conductivity (W/m·K) at TK, PMPa Flow Law Flow law parameters Upper continental crust (Felsic) 2750 0.64+807/(T+77)Wet quartzite (Ranalli, 1995) A=1.97 ×1017 Pa n s, n=2.3 , E=1.54 ×105 J/mol, V=0 J/(mol·MPa), σcr =3×104Pa, C=3 MPa, µ0=0.3, µ1=0 Lower continental crust (Gabbro) 3000 1.18+474/(T+77)Plagioclase An75 (Ranalli, 1995) A=4.80 ×1022 Pa n s, n=3.2 , E=2.38 ×105 J/mol, V=0 J/(mol·MPa), σcr =3×104Pa, C=3 MPa, µ0=0.3, µ1=0 Upper oceanic crust (Basalt) 3000 1.18+474/(T+77)Wet quartzite (Ranalli, 1995) A=4.80 ×1022 Pa n s, n=3.2 , E=2.38 ×105 J/mol, V=0 J/(mol·MPa), σcr =3×104Pa, C=3 MPa, µ0=0, µ1=0 Lower oceanic crust (Gabbro) 3000 1.18+474/(T+77)Plagioclase An75 (Ranalli, 1995) A=4.80 ×1022 Pa n s, n=3.2 , E=2.38 ×105 J/mol, V=0 J/(mol · MPa), σcr =3×104 Pa, C=3 MPa, µ0=0.6 , µ1=0 Mod4:µ1=0.3 Lithospheric and asthenospheric mantle 3300 0.73+1293/(T+77)× exp(0.000004P) Dry olivine (Hirth and Kohlstedt, 2003) m=3 , Adi f f =1.50×1015 Pa s, Edi f f =3.75×105 J/mol, Vdi f f = 0.7 J/(mol · MPa), Adisl =1.10×1016 Pa n s, n=3.5 , Edisl =5.30× 105 J/mol, Vdisl =2.6 J/(mol · MPa), C=3 MPa, µ0=0.6 , µ1=0 Mod3:Vdi f f =0.8 J/(mol·MPa), Vdisl =3.0 J/(mol·MPa) Mod4:µ1=0.3 Weak zone 3300 0.73+1293/(T+77)× exp(0.000004P) Dry olivine (Hirth and Kohlstedt, 2003) m=3 , Adi f f =1.50 ×1015 Pa · s, Edi f f =3.75 ×105 J/mol, Vdi f f =0.7 J/(mol × MPa), Adisl =1.10 ×1016 Pa n s, n=3.5 , Edisl =5.30 ×105 J/mol, Vdisl =2.6 J/(mol · MPa), C=3 MPa, µ0=0, µ1=0 72 Thermomechanical modelling of slab tearing and its topographic response 5.3 Results All the numerical experiments were performed using 24 cores on the ETH-Zürich EULER cluster. Note that, for clarity, all times ‘t’ are expressed in Myr after the initiation of stage 2 for each model. 5.3.1 Evolution of the reference model (Mod1-reference) After the slab has reached the depth of 200 km, the prescribed convergence rate stopped. As the slab continues to sink, due to its density being higher than the mantle (Table 1), the Iberian continental margin started to bend downward and the incoming African block overrides the passive margin (Fig 5.3). At t= 3.84−4.10 Myr, the lithospheric thinning/necking started on the slab’s easternmost side (z=800 km) at 120 km depth (Fig 5.3c). Immediately after the detachment, at 4.24 Myr, the incoming continental block (Africa) came to a complete stop, which, in turn, causes a change in the slab’s downward velocity. The slab’s portion in the vicinity of the tearing appear to have lowered downward velocity, which means the attached portion of the slab continues to sink with a faster downward velocity (Fig 5.3c, d). The tearing point propagates westward reflecting in the tilted angle of the slab’s top edge as shown in the Fig 5.3e and Fig 5.4. While the slab is fully attached, the down-dip motion of the slab induced corner flows, and the large slab body induced a large flow around the slab’s edges (Fig 5.5) The slab’s easternmost part is the only region in which the continental-continental collision occurs which leads to slab tearing (Fig 5.4a1). Any previously present oceanic crust in the forearc wedge appear to have been removed by the thermal erosion. Once the tearing caused the incoming Africa block to stop completely, the tear propagates westwards. As you look westward, the subsequent tearing now is a result of the tearing process that has been set in motion from the east, and not a tearing due to continental-continental collision. The different amount of exhumed oceanic crust in the forearc wedge (Fig 5.4b2,c3) appear to be depending on how large the oceanic domain are in between the incoming African plate and the Iberian plate. We observe a larger amount of exhumed oceanic crust in the westernmost side (Fig 5.4c3). This initiation of slab tearing is observable as a sharp surface uplift along the collisional belt. The rise in elevation (Fig 5.4a, b, c) also evolve westward, reflecting the tear propagation occurring deeper in the mantle. The tearing occurs due to great stress in the bending zone created by both the buoyancy of the Iberia block (upward force) and the weight of the hanging 5.3 Results 73 slab (downward force). The slab is completely detached after t=5.75 Myr (tearing duration of ∼1.65 Myr). Fig 5.3 The evolution of the slab’s downward velocity (Mod1-reference). The slab structure shown here comes from the temperature isosurface,T= 1300◦ C. The cross-section (z=800 km) shows the lithology/composition (for rock composition colour legends please refer to Fig 5.2. The red ‘T’ illustrates the position of the slab tear. Prior necking or slab tearing, the slab subducts with little lateral velocity variation across the slab (a and b). Once the necking and the tearing has started, the higher downward velocity now shifted to side of the slab that is still attached (c and d). After the slab is completely detached (f), the slab’s downward velocity regained the lateral uniformity of downward velocity. Model parameterisation 5.3.2 Effect of a continental-oceanic collision (Mod2) The reference model (Mod1-reference) had an incoming buoyant continental block implemented to create a continental-continental collision, which then led to a one-sided slab tearing. We now move on to look at how the lack of this incoming buoyant continental block would affect the subduction zone dynamic. Rheologically, Mod2 mimics Mod1-reference but the absent of an incoming continental block creates a continental-oceanic collision. At t=3.48 Myr, the oceanic-oceanic subduction zone made contact with the continental passive margin, after which the trench continued forwards. After 0.5 Myr, high topography developed over the trench (Fig 5.6a). Here on the eastern side of the slab, the accumulation of crustal materials above the trench prevented the trench from proceeding any further and lead to 80 Thermomechanical modelling of slab tearing and its topographic response Table 5.2 Model list. CON: continental lithosphere, OC: oceanic lithosphere Model Description Incoming continent Slab detachment Slab tear propagation (cm/yr) Tearing duration (Myr) Uplift rate from tearing (mm/yr) Mod1 Reference model Yes Yes 42.6 1.65 0.23-2.16 Mod2 CON-OC collision No Yes 67.6 1.04 0.71-1.35 Mod3 Increased mantle ductile viscosity Yes Yes 37.6 1.87 0.75-1.68 Mod4 Increased mantle brittle strength Yes No - - - Mod5 Limited convergence velocity Yes No - - - Iberian and African plates (Macchiavelli et al., 2017). This velocity is much slower than the velocity resulting from the hanging slab, it means that this model slows down the slab retreat. Such slow velocity exposes the hanging slab to great amount of thermal diffusion by the surrounding asthenosphere. No slab tearing occurs in this model, but lithospheric dripping occurred instead (Fig 5.10) Lithospheric delamination needs a velocity that is fast enough to prevail the thermal diffusion and maintain the internal temperature and, therefore, density of the slab (Boonma et al., 2019). However, the convergent velocity is too low in this case which means the slab is affected by thermal diffusion more than thermal advection and therefore could not maintain its internal temperature and density, all of which led to thermal erosion and lithospheric dripping (Fig 5.10). The great amount of thermal diffusion that the slab experienced and the amount of time that the slab is hanging in the asthenosphere allow an arcuate (in plan-view) slab to develop. In the models with spontaneous subduction (Mod1-reference, Mod2, Mod3, and Mod4), the subduction and trench migration comes to a stop once the slab reached the passive margin and the tear has started. In Mod5, however, the continuous pushing of the incoming continental block creates a band of high elevation over the arcuate trench (Fig 5.10d). 5.3 Results 81 Fig 5.10 Evolution of model Mod5. The surface topography of when the trench first arrive at the continental passive margin is shown in (a) with the corresponding viscosity cross-sections ( z=300 km and z=800 km) shown in (b) and (c). The right panels (d, e, and f) are from the stage where the tear has propagated and reached the westernmost point of the slab. CHAPTER 6 DISCUSSION This thesis aims to better understand the causes and surface implications of the tearing lithospheric plates, with an application to the late Neogene events affecting the connectivity between the Mediterranean Sea and the Atlantic Ocean. In this chapter, I will synthesize the results obtained using both methodologies and discuss them in the frame of precedent studies. 6.1 Geodynamic relevance of the lithospheric buoyancy Previous studies have established that subduction is a process mostly driven by the negative buoyancy of the oceanic lithosphere with respect to the asthenosphere, and the tectonic forces such as ridge push and slab pull (Bird, 1979; Faccenna et al., 1999; Gurnis et al., 2004; Mart et al., 2005; Afonso et al., 2007; Goren et al., 2008; Valera et al., 2008; Nikolaeva et al., 2010; Zhong and Li, 2019; Candioti et al., 2020). In this context, our results highlight the role of convergence rate, and the lithospheric mantle composition, which determines the density contrast relative to the underlying asthenosphere. Our study shows that the oceanic lithosphere acquires negative buoyancy even at very low convergence rates and predicts ever-increasing slab pull during oceanic plate subduction, proportionally to the accumulated shortening and to the age of the plate. This is consistent with the observation that flat oceanic subduction is rare, requiring additional kinematic and rheological conditions to occur (Huangfu et al., 2016). It is also consistent with the strong load ( 4−10 TN/m) inferred from the flexural bending of the Pacific plate near the Tonga and Kermadec trenches, aged at 105 Ma (Garcia-Castellanos et al., 2000). 84 Discussion Fig 6.1 Combined effect of density contrast and convergence rate. Time needed to reach a downwards lithospheric buoyancy of Fb,re f =−3 TN/m as a function of the density contrast across the LAB for a velocity of plate convergence of 4, 20, and 80 mm/yr. Curves correspond to a continent of lithospheric thickness of 80 (blue) and 160 km (red). The mantle type ranges at the bottom correspond to the various types of Tecton ( 8−29 kg/m 3 ), Proton ( 12−61 kg/m 3 ), and Archon ( 43−93 kg/m 3 ) (O’Reilly et al., 2001). The density contrast for oceanic lithospheric mantle ranges from 0 to +17 kg/m 3 . The arrows indicate the density contrast values used for our model (note that OC120Ma is off the chart and, therefore, not marked). The shaded areas indicate that the −3 TN/m value is not attained. Observations from three geological settings (Alboran (Macchiavelli et al., 2017), Zagros (Hatzfeld and Molnar, 2010), and Tibet (Hatzfeld and Molnar, 2010)) are shown. The horizontal error bars ( ±5 kg/m 3 ) represent uncertainties in the average mantle composition; the vertical error bars correspond to the uncertainty in convergence rate. In contrast, the average Archon-type lithosphere always retains its positive buoyancy, whereas intermediate mantle compositions, corresponding to average Proton-type, can attain negative buoyancy depending on the convergence rate and the elapsed time. For convergence rates v<70 mm/yr, the maximum amplitude of negative buoyancy is attained after 80−250 km shortening (Fig 4.5 ). This is similar to results from models of delamination of a lithospheric root during orogeny (Schott and Schmeling, 1998), where it was concluded that lithospheric roots (equivalent to our ‘subducted slab’ portion) of at least 100−170 km 6.1 Geodynamic relevance of the lithospheric buoyancy 85 length are needed to generate sufficient negative buoyancy for delamination and detachment to proceed. Mantle compositions and densities actually vary over a wide range even within a given age and P-T conditions, such as the sub-types of mantle composition outlined by Griffin et al. (2009). In this thesis, the presented results are based on average chemical compositions representative of five lithospheric mantle types, and based on the assumption that convergence rates and compositions are constant through time. Such limitations help isolate the key parameters we want to study while allowing our model to provide a first-order approximation of how each subcontinental lithospheric mantle type responds to a certain convergence rate. The positive buoyancy of Archon’s is not warranted for lherzolite mantle compositions, as in Kaapvaal (South Africa) and Almklovdalen (western Norway) cratons, showing densities 30 kg/m 3 larger than average Archon-type and close to average garnet Proton-type (Griffin et al., 2009). Similarly, Proton’s with harzburgite compositions can have densities 32 kg/m 3 lower than the average Proton-type and will never get negative buoyancy. Moreover, metasomatic processes can change the major-element mantle composition converting harzburgite/dunite (lighter end-member) to lherzolite (heavier end-member) and back, modifying its buoyancy (O’Reilly and Griffin, 2013). The other key parameter in our results is convergence rate, which may also vary through time, especially for slow speeds ( v≤20 mm/yr) since they involve long lasting processes ( 20−150 Myr) to attain significant shortening. These buoyancy variations, caused either by lateral compositional heterogeneities, metasomatic and differentiation processes, far-field tectonic stresses, or a combination of these, can contribute to isostatic destabilization of cratons by dripping and detachment of heavier mantle portions, as can be the case of western Gondwana craton (Hu et al., 2018). An outstanding outcome of the model is the predicted change from negative to positive buoyancy through time for continental subduction. The initial density contrast across the LAB determines how soon diffusion overcomes advection such that intermediate density contrasts 20 kg/m 3≤∆ρLAB ≤50 kg/m 3 and convergence velocities ≤50 kg/m 3 eventually lead to a change from negative to positive buoyancy (Fig 4.5). This may promote a rising of the subducting continental lithosphere and a subsequent flattening below the overriding plate producing lithosphere underthrusting as proposed in the India-Eurasia collision region (Fig 4.1 and Fig 6.2). According to seismic tomography and potential field modelling (Zhao et al., 2010; Tunini et al., 2016), underthrusting beneath the western Himalaya-Tibet extends 200 km beyond the suture, whereas to the east, the convergence is accommodated by a steep subducting slab. This two-fold behaviour is compatible with a 160-km-thick lithospheric 86 Discussion mantle of Proton-type and subducting at a high velocity of 45 kg/m 3 (Fig 6.1). Clearly, our model disregards other relevant mechanisms such as the weight of oceanic lithospheric slabs subducted before collision or the role of subduction of the continental crust (Cloos, 1993; Faccenda et al., 2009; Afonso and Zlotnik, 2011; Magni et al., 2013; Gö˘ gü¸s et al., 2016; Li et al., 2016), and hence quantitative comparisons with real scenarios must be taken with caution. I have addressed only the influence of mantle composition and convergence rate on lithospheric buoyancy in a kinematic manner. It worked as a basic proof-of-concept; however, my findings do not imply that the same conclusion could be drawn from the modelling of dynamic systems (e.g. a model with active mantle convection), which should be studied further. Keeping these limitations in mind, Fig 6.1 also shows three geological scenarios (Alboran (Macchiavelli et al., 2017), Zagros (Hatzfeld and Molnar, 2010), and Tibet (Hatzfeld and Molnar, 2010), located in Fig 4.1) that can be quantitatively compared to our study. The Alboran Basin is a back-arc setting characterised by a thinner lithospheric mantle (Carballo et al., 2015) (80 km) undergoing Cenozoic mantle subduction, and possibly late Miocene delamination (Duggen et al., 2004), whereas Zagros and Tibet are in the thicker lithospheric mantle group (Tunini et al., 2016) (red lines from Fig 6.1, 160 km). Our model suggests that in order to develop negative buoyancy under higher ∆ρLAB as in cratons, the convergence rate should be higher than 80 mm/yr or remain forced for tens of million years, making it a rare phenomenon (Fig 6.2). The results show that incorporating realistic mineralogy-based densities to geodynamic models rises up an unforeseen control on the development of negative buoyancy. Whereas assigning a constant or temperature-dependent higher density to the lithosphere always results in slab pull, accounting for the effects of composition and pressure reveals that the plate convergence velocity is key to determine the development of negative buoyancy, delamination, and subduction. The model designed here provides a methodological framework for understanding the stability of the lithosphere during the convergence of tectonic plates, and suggests a simple thermodynamic explanation for the long-term preservation of older continental regions (cratons) in the Wilson cycle. As continents aggregated during the early Earth evolution, their average buoyancy relative to the asthenosphere increased (i.e. an accumulation of buoyant material), making them less prone to subduct or delaminate and hence more stable. However, even Archean and Proterozoic lithospheric plates retain a chance to become recycled into the mantle if they are forced to sink by a fast-enough plate convergence, depending on the geographical configuration of tectonic plates. 6.2 Geometry of the passive margin and slab-tearing dynamic 87 Fig 6.2 Cartoon summarizing the effect of the convergence rate. (a) shows the case with a slow convergence rate which allows for great thermal diffusion, which in turn, reducing the slab’s density and the slab becomes more buoyant. (b) shows the case with a faster convergence rate, which allows the slab to keep cooler and denser, leading to a higher density and an increase in negative buoyancy. 6.2 Geometry of the passive margin and slab-tearing dynamic Having learned about the role of lithospheric buoyancy in subduction setting, we now shift our attention to the mechanisms involved in the initiation of slab tearing, specifically, the roles of tectonic configurations in initiating slab tearing. The lateral slab tearing has been invoke to explain some observations, such as the seismicity pattern in the Vrancea slab (Mitrofan et al., 2016) or the seismic tomography in the South Iberian margin, Apennines and Hellenic arc (Wortel, 2000). Previous 3D numerical models of subduction that exhibited slab tearing, such as those by van Hunen and Allen (2011); Duretz et al. (2014); Magni et al. (2017), all utilised a straight passive margin which ended up producing slab tearing in the slab interior rather than on the slab edge. However, in this work, I performed several models that included an oblique continental passive margin which added asymmetry to the collision process. The results from our Betics-inspired models suggest that the obliquity of the continental passive margin has great influence over the initiation of slab tearing because it promotes a laterally diachronous continental collision which leads to slab tearing. In our models, the oblique passive margin gave rise to an east-to-west slab tearing. Averaging over approximately 500 km distance, the tear velocities are 42.6 cm/yr in Mod1-reference, 67.6 cm/yr in Mod2, and 37.6 cm/yr in Mod3 (Table 5.2). Our tear-propagation rates fall well within the range of previous estimations: 7−45 cm/yr from the Carpathians’ depocenter migration by Meulenkamp et al. (1996) from the evolution of the Carpathian-Pannonian system whose 88 Discussion geological model was constructed using regional chronostratigraphic sequences; and 10−80 cm/yr from 3D numerical modelling of continental collision by van Hunen and Allen (2011). A 3D stress model, as a function of rheology, tear length, and force distribution, by Yoshioka and Wortel (1995) showed a tear propagation rate in the range as wide as 2−94 cm/yr. Overall, the slab takes less than 2 Myr to completely detach (over the width of 600− 700 km), which is fast in a geological timescale (compared to the timescales needed for subduction). The factor that seems to have control over the timing of the tearing is the mantle rheology. The viscosity of the sublithospheric mantle in Mod3 is ≤1022 Pa · s whereas it is ≤1021 Pa · s in Mod1-reference. The more viscous mantle in Mod3 slowed down the sinking slab, hence the slowest tear-propagating velocity. Chertova et al. (2014) and Spakman et al. (2018) used numerical models to optimize the fit between a slab morphology as predicted by a geodynamic model, and the seismic tomographic models. The lower mantle in their numerical models has a viscosity of 2×1022 Pa · s which is higher than the value used in my models, hence their slow rollback speed and slow east-to-west slab tearing under the Betics. The sublithospheric mantle viscosity in my models allows the slab to sink faster and, therefore, leads to faster tearing velocity. Chertova et al. (2014) illustrated that in order to create an arcuate slab, such as that imaged under the Western Mediterranean, a variable rollback speed along the strike of the trench is needed. However, the slab rollback in our models was terminated early (where the slab is still attached) due to being exceeded by the fast lateral tearing velocity, hence, the lack of an arcuate slab in our models. This early termination of the slab rollback is also why the gap between the trench and the passive margin exist (Fig 5.7a-c). A slower slab tearing would allow the slab rollback to continue in the region where the slab is still attached, and the western trench could also reach the passive margin (Spakman et al., 2018). Using a geodynamic model, Peral (2020) explains the curvature of the slab by pinning the western edge to the west-Africa continental block, while the NE slab retreats free along the transform fault. The ranges of breakoff depth from our models ( ∼80 −150 km on the eastern side and ∼170−200 km on the western side) fall within similar ranges as previous numerical modelling studies: 80−240 km from Freeburn et al. (2017), 95−140 km from Schellart (2017), 100−400 km from Gerya et al. (2004), and 120−145 km from Duretz et al. (2014). The difference in breakoff depths could plausibly comes from the different rheological setting in each numerical model, as well as different tectonic setup. In the models presented here, the slab’s eastern side has a shallow tearing depth which could be caused by the weakness in the transition zone between the continental and the oceanic lithosphere. While on the western side the tearing depth is deeper as the tear, here, is not only caused by the tectonic variation 6.3 Dynamic topography 89 of the transition zone but also: (i) the negative buoyancy of the hanging and detached portion of the slab; and (ii) the high velocity mantle influx in the slab tear window (Fig 5.5). A similar pattern is reflected in the breakoff location. On the easternmost side, the breakoff tends to occur within the subducted continental lithosphere portion, such that the detached slab pinched out some continental crust. Since the slab tear depth is deeper westward, therefore, the breakoff tends to be within the subducted oceanic lithosphere portion. 6.3 Dynamic topography There are two components which are thought to be shaping the surface topography we observed today, the crustal isostatic compensation effect and dynamic topography (Forte et al., 1993). Dynamic topography is caused by the buoyancy-driven mantle convection exerting vertical stress onto the lithosphere. Dynamic subsidence is caused by downward mantle flow (downwelling), while dynamic uplift is caused by upward mantle flow (upwelling). A key objective in this thesis is to look into the coupling between slab tearing (deep process) and surface uplift (shallow process), and dynamic topography represents such interaction. Fig 6.3 shows model Mod3’s modelled density distribution and the evolution of the modelled elevation, and the two contributions, the isostatic and the dynamic components. We calculated the isostatic effect with a compensation depth of 150 km (128 km below crustal surface). This isostatic elevation is due to the density changes at crustal and lithosphere scales, without accounting the dynamics of the slab subduction. The dynamic topography then came from taking the isostatic effect away from the modelled elevation. From the dynamic topography results, we can identify the slab pull effect with a subsidence and the corner flow and mantle upwelling with an uplift. The dynamic uplift is at its peak at t=10.44 Myr (Fig 6.3a) when the incoming continental crust terminates. While there is mantle downwelling in the mantle wedge (corner flow), the sublithospheric mantle flow upward and, in the process, exhuming the subducted oceanic crustal material up toward the trench. This upward flow and exhumation gave rise to the dynamic uplift which spans over the forearc region (x= 350−440 km). Once the westward lateral tear has reached this cross-section (z=300 km), the upward flow and exhumation in the subduction channel stops (Fig 6.3c). The mantle flow now focuses on the sinking and detaching slab. The reduction of mantle convection in the sublithosphere reduces dynamic uplift greatly. When the detached slab is at a depth of 450−660 km (Fig 6.3d), the mantle convection cells re-established themselves. 96 Discussion evolution regime in the central Betics is sensitive to the initial rheological and thermal conditions. The models in this thesis differ from previous works of the region because I am not trying to model the evolution of the region. Instead, I used the Western Mediterranean as an inspiration to model the relationship between slab tearing and surface uplift because such interaction has previously been proposed in this region (Garcia-Castellanos and Villaseñor, 2011; Spakman et al., 2018). Our model brings the surface-uplift component into play, which connects the mantle dynamic with the surface uplift. The uplift rates from this work (0.23 mm/yr to 2.16 mm/yr) are consistent with the situation during the first stage of the MSC event, in which the uplift of the seaway is compensated for the erosion of seaway, allowing continuous but limited water inflow from the Atlantic into the Mediterranean Sea. The tectonic and erosion model by Garcia-Castellanos and Villaseñor (2011) proposed that a critical uplift rate of 5 mm/yr is needed to close the seaways across the Gibraltar Arc. Coulson et al. (2019) built upon Garcia-Castellanos and Villaseñor (2011)’s model by incorporating an ice-age sea level theory, which predicts a critical uplift rate of <1.5 mm/yr. Previous biostratigraphic studies in the Betics suggested that the age of marine uplift in the region was in the range of ∼10−8 Ma, younging westward (Garcés et al., 1998; Iribarren et al., 2009). The last of these basins’ emersion was previously dated at 5.3 Ma which indicates that the duration of tear propagation in the Betics (across ∼400 km distance) is ∼3−5 Myr. This assumption would give a tearing rate of 80−133 km/Myr. The newer biostratigraphic studies by Krijgsman et al. (2018) and van der Schee et al. (2018) revised the age of the uplift and provided a new age constraint on the western Betic intramontane basins to be older than 7.51 Ma (late Tortonian). The new age constraint narrows down the plausible tear propagation duration to ∼1−3 Myr, which implies a faster tearing ( ∼133 −400 km/Myr) rate than previously thought . Our tearing rate of 370−670 km/Myr, albeit very fast, still falls within the newly constrained tear propagation rate. A limitation of our models lies in the lack of an arcuate slab so a full comparison cannot be made with the interpreted 3D slab structure beneath the Gibraltar Arc System (Spakman and Wortel, 2004; Chertova et al., 2014; Spakman et al., 2018). The uplift of intramountain basins within the Betics in southern Iberia is higher on the eastern side (Garcés et al., 1998; Iribarren et al., 2009; Garcia-Castellanos and Villaseñor, 2011), where the slab is interpreted to be detached based on seismic tomography (Spakman and Wortel, 2004; Mancilla et al., 2015, 2018) (Fig 6.6b). Such uplift is not detected in the western Betics where the tear point is and the part of the same slab still remains attached (Mancilla et al., 2015) (Fig 6.6a). Our models predict a similar trend, with earlier and higher uplift on the eastern parts of the oblique margin (due to both continental-continental collision and slab tearing), and later and lower 6.5 Implications for the Western Mediterranean 97 uplift in the west, where the slab still remain attached. The study of magnetostratigraphic sequences shows that the transition from marine to continental conditions of intramountain basins within the Betics is younging westward (Fig 5.1) (Garcés et al., 1998; Iribarren et al., 2009). This trend corresponds with our models’ westward tear propagation, where the oldest uplifted region would be toward the east and the younger uplifted region toward the west. My study has been primarily concerned with the initiation of slab tearing and the consequent surface response. I found a specific model setup that produced a laterally propagating slab-tear, which I then carried on modifying. The findings of my slab-tearing study are, therefore, restricted to specific mantle rheologies and database, and my analysis has concentrated on the geodynamic processes arising from those specific setups. In this study, I used 3D high-resolution geodynamic models, which can take up to 1−2 months to run per model, which put a heavy time-constraint on the workflow. Such high computational demands led to one of the key limitations in this study, which is the lack of a more systematic exploration of the parameter field. 98 Discussion Fig 6.6 Comparison of slab structure from Mod1-reference with the seismic tomography of the Western Mediterranean. (a) and (b) are the seismic tomography from Garcia-Castellanos and Villaseñor (2011). (c) and (d) are viscosity cross-sections from model Mod1-reference with (c) sliced from z=300 km, where the subducting slab is still attached and (d) sliced from z=800 km, where the slab has just started tearing. The subsets in (c) and (d) shows the plane-view (x-z) of the surface elevation, and the red lines indicate the position of the corresponding cross-sections. The cross-sections from Mod1-reference resembles, to an extent, the seismic tomography from the Western Mediterranean, with the attached portion of the slab on the NW side (a) and the detached slab towards the NE side (b). CHAPTER 7 CONCLUSION AND OUTLOOK The presented thesis aimed to explore the link between the tearing of subducted lithospheric slabs and the associated vertical motions. To this purpose, I first studied the mechanisms controlling the buoyancy of a subducted lithosphere and then, with this understanding, numerically simulate the process of lithospheric tearing upon continental collision, using the Betic Cordillera as a reference geological region. My results support the idea that the negative buoyancy of a subducting lithosphere is primarily controlled by its composition and the subduction velocity (Chapter 4). Further conclusions are: • With a mineral-physics approach, where a lithospheric mantle can be less dense than the underlying asthenosphere, I found that Proton ( 2.5−1.0 Ga) and Tecton( <1.0 Ga) continental lithospheric mantle can attain negative buoyancy and increase their subductibility depending on the subduction velocity. • The less dense cratons (Archon continental lithospheric mantle ( >2.5 Ga)) overprints this effect and hinders subduction or delamination, regardless of the convergence rate. This may explain why Archons are more stable during the Wilson Cycle. • The results also support the idea that the subducted lithospheric mantle, having attained negative buoyancy, can also return to being positively buoyant overtime, provided that the initial density contrast across the LAB was in the range of 20−50 kg/m 3 and the convergence velocity of ≤50 mm/yr. This may promote a rising of the subducting continental lithosphere and a subsequent flattening below the overriding plate producing lithosphere underthrusting. 100 Conclusion and outlook The Betics-inspired numerical modelling study of lithospheric slab tearing aimed to study the characteristics surrounding slab tearing and the associated surface uplift. Based on my results (Chapter 5), and within the limitations of my methodologies, I conclude that: • The obliquity of the continental passive margin (relative to the trench axis) is a major influence on the initiation of slab tearing because it promotes a laterally diachronous continental collision which leads to slab tearing. • The slab tearing depth increases as it propagates along the slab, from 80−150 km to 170−200 km. • The viscosity of the sublithospheric mantle and the amount of shortening (oceanic subduction) prior tearing (which affects how large and heavy the slab would be) are the key controls on the slab-tearing dynamics. • This obliquity of the continental collision produced an east-to-west slab tearing (tearing velocity ∼37.6−67.6 cm/yr with the lower-mantle viscosity of up to 1022 Pa · s), which led to surface uplift signature of 0.5−1.5 km across the forearc region throughout the tearing process. • The uplift rate produced by a slab tearing process, 0.23−2.16 mm/yr, is compatible with the uplift rate needed to achieve an equilibrium between seaway-uplift and seaway-erosion, which could have led to the closure of marine gateways that reduced the water-flow from the Atlantic Ocean into the Mediterranean Sea during the first stage of the Messinian Salinity Crisis. Outlook To better understand the geodynamic of slab tearing, future studies could consider a proper parameterization between mantle viscosity and velocity of the slab-tearing propagation, as well as addressing the effect of horizontal dimensions of the slab on the tearing propagation. Such an systematic exploration of the full parameter field will require significantly higher computing capacity than presently available. More specific to the techniques developed in this thesis, the code LithBuoy can be improved further by implementing active mantle convection, which will add a dynamic component to the work. It will also be useful to investigate how likely the ‘triggering’ buoyancy force value ( Fb=3 TN/m) could lead to lithospheric delamination or flat-subduction in a dynamic convective model. The composition of the subcontinental lithosphere (SCLM) has 101 a wide range and the ones I used in this thesis are the averaged properties of each mantle type. I would suggest a further study which include those sub-types of mantle composition as outlined by Griffin et al. (2009). Regarding the link between tectonic uplift from slab tearing (geodynamics) and surface processes (geomorphology) that led to the Messinian Salinity Crisis event, while the results support the idea that the rates of uplift are quantitatively compatible with those needed to close the Mediterranean-Atlantic gateways, a more detailed coupling between the two sets of mechanisms is encouraged, since the feedback phenomena between the erosion of the uplifted seaways and the salinization of the Mediterranean have been previously described. The different scales of the two study fields will be problematic when trying to relate one to another. Some geodynamic codes have been coupled with high-resolution surface process codes, which link the deep process directly with the shallow (surface) process, e.g. ASPECT-FastScape coupling and I3ELVIS-DAC (Ueda et al., 2015). Although this coupling of computing codes of very different nature will be very computationally expensive, going forward, the computing facilities can only improve. This coupling should provide more constraints on the feedback loop between the uplift from deep processes and the surface processes. 102 Conclusion and outlook Fig 7.1 An illustration summarising the lithospheric tearing process and the role of lithospheric buoyancy. An oblique continental margin promotes a laterally diachronous continental collision. Continental collision hinders or slows down subduction, which creates an upward force that counteracts the slab pull force. This built-up tension eventually leads to the initiation of slab tearing, which creates a tear-gap through which mantle rushes, giving viscous support to the overlying lithosphere. The tear then propagates laterally towards another side of the slab. The extra viscous support under the lithosphere leads to vertical uplift on the surface, which is a signature of slab tearing. The surface uplift follows the general direction of slab tear propagation. The uplifted region subsides after the mantle dynamic below returns to the state prior to tearing. The torn portion of the slab now has smaller slab pull force (negative buoyancy) compared to the still-attached portion, in which the slab pull force is still being transferred up to the subduction zone. REFERENCES Afonso, J. C., Fernàndez, M., Ranalli, G., Griffin, W. L., and Connolly, J. A. D. (2008). Integrated geophysical-petrological modeling of the lithosphere and sublithospheric upper mantle: Methodology and applications. Geochemistry, Geophys. Geosystems, 9(5). Afonso, J. C., Ranalli, G., and Fernàndez, M. (2007). Density structure and buoyancy of the oceanic lithosphere revisited. Geophys. Res. Lett., 34(10):L10302. Afonso, J. C. and Zlotnik, S. (2011). The subductability of continental lithosphere: The before and after story. Front. Earth Sci., 4:53–86. Al-Lazki, A. I., Seber, D., Sandvol, E., Turkelli, N., Mohamad, R., and Barazangi, M. (2003). Tomographic Pn velocity and anisotropy structure beneath the Anatolian plateau (eastern Turkey) and the surrounding regions. Geophys. Res. Lett., 30(24). Andrews, E. R. and Billen, M. I. (2009). Rheologic controls on the dynamics of slab detachment. Tectonophysics, 464(1-4):60–69. Andrieux, J., Fontbote, J. M., and Mattauer, M. (1971). Sur un modele explicatif de l’arc de Gibraltar. Earth Planet. Sci. Lett., 12(2):191–198. Argus, D. F., Gordon, R. G., Demets, C., and Stein, S. (1989). Closure of the Africa-EurasiaNorth America Plate motion circuit and tectonics of the Gloria Fault. J. Geophys. Res., 94(B5):5585–5602. Austermann, J. and Forte, A. M. (2019). The importance of dynamic topography for understanding past sea-level changes. Past Glob. Chang. Mag., 27(1). Ávila, P. and Dávila, F. M. (2020). Lithospheric thinning and dynamic uplift effects during slab window formation, southern Patagonia. J. Geodyn., 133:101689. Ayachit, U. (2015). The ParaView Guide: A Parallel Visualization Application. Kitware, Inc., Clifton Park, NY, USA. Aydin, I., Karat, H., and Koçak, A. (2005). Curie-point depth map of Turkey. Geophys. J. Int., 162(2):633–640. Bajolet, F., Galeano, J., Funiciello, F., Moroni, M., Negredo, A. M., and Faccenna, C. (2012). Continental delamination: Insights from laboratory models. Geochemistry, Geophys. Geosystems, 13(1):1–22. 104 References Bartolini, C. (2003). When did the Northern Apennine become a mountain chain? Quat. Int., 101-102:75–80. Bercovici, D. and Ricard, Y. (2012). Mechanisms for the generation of plate tectonics by two-phase grain-damage and pinning. Phys. Earth Planet. Inter., 202-203:27–55. Bercovici, D. and Ricard, Y. (2013). Generation of plate tectonics with two-phase graindamage and pinning: Source-sink model and toroidal flow. Earth Planet. Sci. Lett., 365:275–288. Bercovici, D., Schubert, G., and Ricard, Y. (2015). Abrupt tectonics and rapid slab detachment with grain damage. Proc. Natl. Acad. Sci., 112(5):1287–1291. Berk Biryol, C., Beck, S. L., Zandt, G., and Özacar, A. A. (2011). Segmented African lithosphere beneath the Anatolian region inferred from teleseismic P-wave tomography. Geophys. J. Int., 184(3):1037–1057. Bezada, M., Humphreys, E., Toomey, D., Harnafi, M., Dávila, J., and Gallart, J. (2013). Evidence for slab rollback in westernmost Mediterranean from improved upper mantle imaging. Earth Planet. Sci. Lett., 368:51–60. Bird, P. (1979). Continental delamination and the Colorado Plateau. J. Geophys. Res. Solid Earth, 84(B13):7561–7571. Blanco, M. J. and Spakman, W. (1993). The P-wave velocity structure of the mantle below the Iberian Peninsula: evidence for subducted lithosphere below southern Spain. Tectonophysics, 221(1):13–34. Bonnin, M., Nolet, G., Villaseñor, A., Gallart, J., and Thomas, C. (2014). Multiple-frequency tomography of the upper mantle beneath the African/Iberian collision zone. Geophys. J. Int., 198(3):1458–1473. Boonma, K., Kumar, A., Garcia-Castellanos, D., Jiménez-Munt, I., and Fernández, M. (2019). Lithospheric mantle buoyancy: the role of tectonic convergence and mantle composition. Sci. Rep., 9(1):17953. Boschi, L., Faccenna, C., and Becker, T. W. (2010). Mantle structure and dynamic topography in the Mediterranean Basin. Geophys. Res. Lett., 37(20):20303. Boutelier, D. and Cruden, A. R. (2017). Slab breakoff: Insights from 3D thermo-mechanical analogue modelling experiments. Tectonophysics, 694:197–213. Bürgmann, R., Dresen, G., Urgmann, R. B. ., and Dresen, G. (2008). Rheology of the lower crust and upper mantle: Evidence from rock mechanics, geodesy, and field observations. Annu. Rev. Earth Planet. Sci., 36:531–567. Calvert, A., Sandvol, E., Seber, D., Barazangi, M., Roecker, S., Mourabit, T., Vidal, F., Alguacil, G., and Jabour, N. (2000). Geodynamic evolution of the lithosphere and upper mantle beneath the Alboran region of the western Mediterranean: Constraints from travel time tomography. J. Geophys. Res. Solid Earth, 105(B5):10871–10898. References 105 Campbell, I. H. and Taylor, S. R. (1983). No water, no granites - No oceans, no continents. Geophys. Res. Lett., 10(11):1061–1064. Candioti, L. G., Schmalholz, S. M., and Duretz, T. (2020). Impact of upper mantle convection on lithosphere hyperextension and subsequent horizontally forced subduction initiation. Solid Earth, 11(6):2327–2357. Capella, W., Spakman, W., Hinsbergen, D. J. J., Chertova, M. V., and Krijgsman, W. (2020). Mantle resistance against Gibraltar slab dragging as a key cause of the Messinian Salinity Crisis. Terra Nov., 32(2):141–150. Carballo, A., Fernandez, M., Jiménez-Munt, I., Torne, M., Vergés, J., Melchiorre, M., Pedreira, D., Afonso, J. C., Garcia-Castellanos, D., Díaz, J., Villaseñor, A., Pulgar, J. A., and Quintana, L. (2015). From the North-Iberian Margin to the Alboran Basin: A lithosphere geo-transect across the Iberian Plate. Tectonophysics, 663:399–418. Carminati, E., Doglioni, C., and Barba, S. (2004). Reverse migration of seismicity on thrusts and normal faults. Earth-Science Rev., 65(3-4):195–222. Carminati, E., Doglioni, C., and Scrocca, D. (2003). Apennines subduction-related subsidence of Venice (Italy). Geophys. Res. Lett., 30(13). Carminati, E., Lustrino, M., and Doglioni, C. (2012). Geodynamic evolution of the central and western Mediterranean: Tectonics vs. igneous petrology constraints. Tectonophysics, 579:173–192. Carminati, E., Wortel, R., Spakman, W., and Sabadini, R. (1998). The role of slab detachment processes in the opening of the western–central Mediterranean basins: some geological and geophysical evidence. Earth Planet. Sci. Lett., 160(3-4):651–665. Carnevale, G., Dela Pierre, F., Natalicchio, M., and Landini, W. (2018). Fossil marine fishes and the ’Lago Mare’ event: Has the Mediterranean ever transformed into a brackish lake? Newsletters Stratigr., 51(1):57–72. Casciello, E., Fernández, M., Vergés, J., Cesarano, M., and Torné, M. (2015). The Alboran Domain in the Western Mediterranean evolution: the birth of a concept. Bull. Soc. Geol. Fr., 186(4-5):23–34. Channell, J. E. T. and Mareschal, J. C. (1989). Delamination and asymmetric lithospheric thickening in the development of the Tyrrhenian Rift. Geol. Soc. London, Spec. Publ., 45(1):285–302. Chertova, M. V., Spakman, W., Geenen, T., van den Berg, A. P., and van Hinsbergen, D. J. J. (2014). Underpinning tectonic reconstructions of the western Mediterranean region with dynamic slab evolution from 3-D numerical modeling. J. Geophys. Res. Solid Earth, 119(7):5876–5902. Chiarabba, C. and Chiodini, G. (2013). Continental delamination and mantle dynamics drive topography, extension and fluid discharge in the Apennines. Geology, 41(6):715–718. Chopin, C. (1984). Coesite and pure pyrope in high-grade blueschists of the Western Alps: a first record and some consequences. Contrib. to Mineral. Petrol., 86(2):107–118. 112 References Jiménez-Munt, I., Torne, M., Fernàndez, M., Vergés, J., Kumar, A., Carballo, A., and GarcíaCastellanos, D. (2019). Deep Seated Density Anomalies Across the Iberia-Africa Plate Boundary and Its Topographic Response. J. Geophys. Res. Solid Earth, 124(12):13310– 13332. Jolivet, L., Faccenna, C., and Piromallo, C. (2009). From mantle to crust: Stretching the Mediterranean. Earth Planet. Sci. Lett., 285(1-2):198–209. Jones, C. H., Reeg, H., Zandt, G., Gilbert, H., Owens, T. J., and Stachnik, J. (2014). P-wave tomography of potential convective downwellings and their source regions, Sierra Nevada, California. Geosphere, 10(3):505–533. Karato, S. I. (2008). Deformation of earth materials: An introduction to the rheology of solid earth, volume 9780521844. Cambridge University Press. Karato, S.-i. and Wu, P. (1993). Rheology of the upper mantle: a synthesis. Science, 260(5109):771–778. Katsura, T. and Ito, E. (1989). The system Mg 2 SiO 4 -Fe 2 SiO 4 at high pressures and temperatures: Precise determination of stabilities of olivine, modified spinel, and spinel. J. Geophys. Res. Solid Earth, 94(B11):15663–15670. Kay, R. and Mahlburg Kay, S. (1993). Delamination and delamination magmatism. Tectonophysics, 219(1-3):177–189. Komut, T. (2015). High surface topography related to upper mantle flow beneath Eastern Anatolia. Geophys. J. Int., 203(2):1263–1273. Krijgsman, W., Capella, W., Simon, D., Hilgen, F. J., Kouwenhoven, T. J., Meijer, P. T., Sierro, F. J., Tulbure, M. A., van den Berg, B. C., van der Schee, M., and Flecker, R. (2018). The Gibraltar Corridor: Watergate of the Messinian Salinity Crisis. Mar. Geol., 403(June):238–246. Krijgsman, W., Hilgen, F. J., Raffi, I., Sierro, F. J., and Wilson, D. S. (1999). Chronology, causes and progression of the Messinian salinity crisis. Nature, 400(6745):652–655. Kumar, A., Fernàndez, M., Vergés, J., Torne, M., and Jiménez-Munt, I. (2021). Opposite symmetry in the lithospheric structure of the Alboran and Algerian basins and their margins (Western Mediterranean): Geodynamic implications. J. Geophys. Res. Solid Earth. Li, Z.-H. H., Liu, M., and Gerya, T. (2016). Lithosphere delamination in continental collisional orogens: A systematic numerical study. J. Geophys. Res. Solid Earth, 121(7):5186– 5211. Liang, X., Chen, Y., Tian, X., Chen, Y. J., Ni, J., Gallegos, A., Klemperer, S. L., Wang, M., Xu, T., Sun, C., Si, S., Lan, H., and Teng, J. (2016). 3D imaging of subducting and fragmenting Indian continental lithosphere beneath southern and central Tibet using body-wave finite-frequency tomography. Earth Planet. Sci. Lett., 443:162–175. Lithgow-Bertelloni, C. and Silver, P. G. (1998). Dynamic topography, plate driving forces and the African superswell. Nature, 395(6699):269–272. References 113 Lugli, S., Manzi, V., Roveri, M., and Schreiber, B. C. (2010). The Primary Lower Gypsum in the Mediterranean: A new facies interpretation for the first stage of the Messinian salinity crisis. Palaeogeogr. Palaeoclimatol. Palaeoecol., 297(1):83–99. Luján, M., Crespo-Blanc, A., and Balanyá, J. C. (2006). The Flysch Trough thrust imbricate (Betic Cordillera): A key element of the Gibraltar Arc orogenic wedge. Tectonics, 25(6). Macchiavelli, C., Vergés, J., Schettino, A., Fernàndez, M., Turco, E., Casciello, E., Torne, M., Pierantoni, P. P., and Tunini, L. (2017). A New Southern North Atlantic Isochron Map: Insights Into the Drift of the Iberian Plate Since the Late Cretaceous. J. Geophys. Res. Solid Earth, 122(12):9603–9626. Magni, V., Allen, M. B., van Hunen, J., Bouilhol, P., Hunen, J. V., Bouilhol, P., van Hunen, J., and Bouilhol, P. (2017). Continental underplating after slab break-off. Earth Planet. Sci. Lett., 474:59–67. Magni, V., Faccenna, C., Van Hunen, J., and Funiciello, F. (2013). Delamination vs. break-off: The fate of continental collision. Geophys. Res. Lett., 40(2):285–289. Mancilla, F. d. L., Booth-Rea, G., Stich, D., Pérez-Peña, J. V., Morales, J., Azañón, J. M., Martin, R., and Giaconia, F. (2015). Slab rupture and delamination under the Betics and Rif constrained from receiver functions. Tectonophysics, 663:225–237. Mancilla, F. d. L., Heit, B., Morales, J., Yuan, X., Stich, D., Molina-Aguilera, A., Azañon, J. M., and Martín, R. (2018). A STEP fault in Central Betics, associated with lateral lithospheric tearing at the northern edge of the Gibraltar arc subduction system. Earth Planet. Sci. Lett., 486:32–40. Marchesi, C., Garrido, C. J., Bosch, D., Bodinier, J.-L. L., Hidas, K., Padrón-Navarta, J. A., and Gervilla, F. (2012). A late oligocene suprasubduction setting in the westernmost mediterranean revealed by intrusive pyroxenite dikes in the Ronda peridotite (Southern Spain). J. Geol., 120(2):237–247. Mart, Y., Aharonov, E., Mulugeta, G., Ryan, W., Tentler, T., and Goren, L. (2005). Analogue modelling of the initiation of subduction. Geophys. J. Int., 160(3):1081–1091. Martínez-García, P., Comas, M., Soto, J. I., Lonergan, L., and Watts, A. B. (2013). Strike-slip tectonics and basin inversion in the Western Mediterranean: The Post-Messinian evolution of the Alboran Sea. Basin Res., 25(4):361–387. McClusky, S., Reilinger, R., Mahmoud, S., Ben Sari, D., and Tealeb, A. (2003). GPS constraints on Africa (Nubia) and Arabia plate motions. Geophys. J. Int., 155(1):126–138. Meissner, R. and Mooney, W. (1998). Weakness of the lower continental crust: a condition for delamination, uplift, and escape. Tectonophysics, 296(1-2):47–60. Meulenkamp, J. E., Kovᡠc, M., and Cicha, I. (1996). On Late Oligocene to Pliocene depocentre migrations and the evolution of the Carpathian-Pannonian system. Tectonophysics, 266(1-4):301–317. 114 References Michard, A., Frizon de Lamotte, D., and Chalouan, A. (2005). Comment on “The ultimate arc: Differential displacements, oroclinal bending, and vertical axis rotation in the External Betic-Rif arc” by J. P. Platt et al. Tectonics, 24(1). Mitrofan, H., Anghelache, M. A., Chitea, F., Damian, A., Cadicheanu, N., and Vi¸san, M. (2016). Lateral detachment in progress within the Vrancea slab (Romania): Inferences from intermediate-depth seismicity patterns. Geophys. J. Int., 205(2):864–875. Molnar, P., England, P., and Martinod, J. (1993). Mantle dynamics, uplift of the Tibetan Plateau, and the Indian Monsoon. Rev. Geophys., 31(4):357. Mooney, W. D. (2007). Crust and Lithospheric Structure - Global Crustal Structure. In Treatise Geophys., volume 1, pages 361–417. Elsevier. Morales, J., Serrano, I., Jabaloy, A., Galindo-Zaldívar, J., Zhao, D., Torcal, F., Vidal, F., and Lodeiro, F. G. (1999). Active continental subduction beneath the Betic Cordillera and the Alborán Sea. Geology, 27(8):735. Mulyukova, E. and Bercovici, D. (2017). Formation of lithospheric shear zones: Effect of temperature on two-phase grain damage. Phys. Earth Planet. Inter., 270:195–212. Negredo, A. M., Mancilla, F. d. L., Clemente, C., Morales, J., and Fullea, J. (2020). Geodynamic Modeling of Edge-Delamination Driven by Subduction-Transform Edge Propagator Faults: The Westernmost Mediterranean Margin (Central Betic Orogen) Case Study. Front. Earth Sci., 8:6. Nikolaeva, K., Gerya, T. V., and Marques, F. O. (2010). Subduction initiation at passive margins: Numerical modeling. J. Geophys. Res. Solid Earth, 115(3). Niu, Y. (2017). Slab breakoff: a causal mechanism or pure convenience? Sci. Bull., 62(7):456–461. Ohneiser, C., Florindo, F., Stocchi, P., Roberts, A. P., Deconto, R. M., and Pollard, D. (2015). Antarctic glacio-eustatic contributions to late Miocene Mediterranean desiccation and reflooding. Nat. Commun., 6. O’Reilly, S. Y. and Griffin, W. L. (2013). Mantle Metasomatism. pages 471–533. O’Reilly, S. Y., Griffin, W. L., Djomani, Y. H., and Morgan, P. (2001). Are lithospheres forever? Tracking changes in subcontinental lithospheric mantle through time. GSA Today. Palano, M., González, P. J., and Fernández, J. (2013). Strain and stress fields along the Gibraltar Orogenic Arc: Constraints on active geodynamics. Gondwana Res., 23(3):1071– 1088. Pedrera, A., Ruiz-Constán, A., García-Senz, J., Azor, A., Marín-Lechado, C., Ayala, C., Díaz de Neira, J. A., and Rodríguez-Fernández, L. R. (2020). Evolution of the South-Iberian paleomargin: From hyperextension to continental subduction. J. Struct. Geol., 138:104122. Peral, M. (2020). Dynamics of subduction Systems with opposite polarity in adjacent segments: application to the Westernmost Mediterranean. Phd thesis, Universitat de Barcelona. References 115 Pérouse, E., Vernant, P., Chéry, J., Reilinger, R., and McClusky, S. (2010). Active surface deformation and sub-lithospheric processes in the western Mediterranean constrained by numerical models. Geology, 38(9):823–826. Piana Agostinetti, N., Lucente, F. P., Selvaggi, G., and Di Bona, M. (2002). Crustal structure and moho geometry beneath the Northern Apennines (Italy). Geophys. Res. Lett., 29(20):60–1–60–4. Platt, J. P., Anczkiewicz, R., Soto, J.-I., Kelley, S. P., and Thirlwall, M. (2006). Early Miocene continental subduction and rapid exhumation in the western Mediterranean. Geology, 34(11):981. Platt, J. P., Behr, W. M., Johanesen, K., and Williams, J. R. (2013). The Betic-Rif Arc and Its Orogenic Hinterland: A Review. Annu. Rev. Earth Planet. Sci., 41(1):313–357. Platt, J. P., Kaus, B., and Becker, T. W. (2008). The mechanics of continental transforms: An alternative approach with applications to the San Andreas system and the tectonics of California. Earth Planet. Sci. Lett., 274(3-4):380–391. Platt, J. P. and Vissers, R. L. (1989). Extensional collapse of thickened continental lithosphere: a working hypothesis for the Alboran Sea and Gibraltar arc. Geology, 17(6):540–543. Ranalli, G. (1995). Rheology of the Earth. 2nd ed. edition. Ren, Y. and Shen, Y. (2008). Finite frequency tomography in southeastern Tibet: Evidence for the causal relationship between mantle lithosphere delamination and the north–south trending rifts. J. Geophys. Res., 113(B10):B10316. Rohling, E. J., Schiebel, R., and Siddall, M. (2008). Controls on Messinian Lower Evaporite cycles in the Mediterranean. Earth Planet. Sci. Lett., 275(1-2):165–171. Rosenbaum, G., Gasparon, M., Lucente, F. P., Peccerillo, A., and Miller, M. S. M. S. (2008). Kinematics of slab tear faults during subduction segmentation and implications for Italian magmatism. Tectonics, 27(2). Rosenbaum, G., Lister, G. S., and Duboz, C. (2002). Relative motions of Africa, Iberia and Europe during Alpine orogeny. Tectonophysics, 359(1-2):117–129. Royden, L. H. (1993). The tectonic expression slab pull at continental convergent boundaries. Tectonics, 12(2):303–325. Rozel, A., Ricard, Y., and Bercovici, D. (2011). A thermodynamically self-consistent damage equation for grain size evolution during dynamic recrystallization. Geophys. J. Int., 184(2):719–728. Ruiz-Constn, A., Pedrera, A., Galindo-Zaldívar, J., Pous, J., Arzate, J., Roldn-García, F. J., Marin-Lechado, C., and Anahnah, F. (2012). Constraints on the frontal crustal structure of a continental collision from an integrated geophysical research: The central-western Betic Cordillera (SW Spain). Geochemistry, Geophys. Geosystems, 13(8). Ruppert, S., Fliedner, M. M., and Zandt, G. (1998). Thin crust and active upper mantle beneath the Southern Sierra Nevada in the western United States. Tectonophysics, 286(14):237–252. 116 References Saleeby, J. and Foster, Z. (2004). Topographic response to mantle lithosphere removal in the southern Sierra Nevada region, California. Geology, 32(3):245–248. Schellart, W. P. (2017). A geodynamic model of subduction evolution and slab detachment to explain Australian plate acceleration and deceleration during the latest Cretaceous-early Cenozoic. Lithosphere, 9(6):976–986. Schmeling, H., Babeyko, A. Y., Enns, A., Faccenna, C., Funiciello, F., Gerya, T., Golabek, G. J., Grigull, S., Kaus, B. J., Morra, G., Schmalholz, S. M., and van Hunen, J. (2008). A benchmark comparison of spontaneous subduction models-Towards a free surface. Phys. Earth Planet. Inter., 171(1-4):198–223. Schott, B. and Schmeling, H. (1998). Delamination and detachment of a lithospheric root. Tectonophysics, 296(3-4):225–247. Seber, D., Barazangi, M., Ibenbrahim, A., and Demnati, A. (1996). Geophysical evidence for lithospheric delamination beneath the Alboran Sea and Rif–Betic mountains. Nature, 379(6568):785–790. Sen, S., Antoine, P.-O., Varol, B., Ayyildiz, T., and Sözeri, K. (2011). Giant rhinoceros Paraceratherium and other vertebrates from Oligocene and middle Miocene deposits of the Ka˘ gızman-Tuzluca Basin, Eastern Turkey. Naturwissenschaften, 98(5):407–423. Sinclair, H. D. (1997). Flysch to molasse transition in peripheral foreland basins: The role of the passive margin versus slab breakoff. Geology, 25(12):1123–1126. Sippl, C., Schurr, B., Yuan, X., Mechie, J., Schneider, F. M., Gadoev, M., Orunbaev, S., Oimahmadov, I., Haberland, C., Abdybachaev, U., Minaev, V., Negmatullaev, S., and Radjabov, N. (2013). Geometry of the Pamir-Hindu Kush intermediate-depth earthquake zone from local seismic data. J. Geophys. Res. Solid Earth, 118(4):1438–1457. Spakman, W. (1986). Subduction beneath Eurasia in connection with the Mesozoic Tethys. Geol. en Mijnb., 65(2):145–153. Spakman, W., Chertova, M. V., van den Berg, A., and van Hinsbergen, D. J. J. (2018). Puzzling features of western Mediterranean tectonics explained by slab dragging. Nat. Geosci., 11(3):211–216. Spakman, W. and Wortel, R. (2004). A Tomographic View on Western Mediterranean Geodynamics. In TRANSMED Atlas. Mediterr. Reg. from Crust to Mantle, pages 31–52. Springer Berlin Heidelberg, Berlin, Heidelberg. Stich, D., Serpelloni, E., de Lis Mancilla, F., and Morales, J. (2006). Kinematics of the Iberia-Maghreb plate contact from seismic moment tensors and GPS observations. Tectonophysics, 426(3-4):295–317. Thurner, S., Palomeras, I., Levander, A., Carbonell, R., and Lee, C.-T. (2014). Ongoing lithospheric removal in the western Mediterranean: Evidence from Ps receiver functions and thermobarometry of Neogene basalts (PICASSO project). Geochemistry, Geophys. Geosystems, 15(4):1113–1127. References 117 Torres-Roldan, R. L., Poli, G., and Peccerillo, A. (1986). An early Miocene arc-tholeiitic magmatic dike event from the Alboran Sea — Evidence for precollisional subduction and back-arc crustal extension in the westernmost Mediterranean. Geol. Rundschau, 75(1):219–234. Tunini, L., Jiménez-Munt, I., Fernandez, M., Vergés, J., Villaseñor, A., Melchiorre, M., and Afonso, J. C. (2016). Geophysical-petrological model of the crust and upper mantle in the India-Eurasia collision zone. Tectonics, 35(7):1642–1669. Turcotte, D. and Schubert, G. (2014). Geodynamics. Cambridge University Press. Turner, S. P., Platt, J. P., George, R. M. M., Kelley, S. P., Pearson, D. G., and Nowell, G. M. (1999). Magmatism Associated with Orogenic Collapse of the Betic-Alboran Domain, SE Spain. J. Petrol., 40(6):1011–1036. Ueda, K., Willett, S. D., Gerya, T., and Ruh, J. (2015). Geomorphological-thermo-mechanical modeling: Application to orogenic wedge dynamics. Tectonophysics, 659:12–30. Valera, J. L., Negredo, A. M., Billen, M. I., and Jiménez-Munt, I. (2014). Lateral migration of a foundering high-density root: Insights from numerical modeling applied to the southern Sierra Nevada. Lithos, 189:77–88. Valera, J. L., Negredo, A. M., and Jiménez-Munt, I. (2011). Deep and near-surface consequences of root removal by asymmetric continental delamination. Tectonophysics, 502(1-2):257–265. Valera, J.-L. L., Negredo, A.-M. M., and Villaseñor, A. (2008). Asymmetric delamination and convective removal numerical modeling: Comparison with evolutionary models for the Alboran Sea region. Pure Appl. Geophys., 165(8):1683–1706. van der Schee, M., van den Berg, B. C. J., Capella, W., Simon, D., Sierro, F. J., and Krijgsman, W. (2018). New age constraints on the western Betic intramontane basins: A late Tortonian closure of the Guadalhorce Corridor? Terra Nov., 30(5):325–332. van Gerven, L., Deschamps, F., and van der Hilst, R. D. (2004). Geophysical evidence for chemical variations in the Australian Continental Mantle. Geophys. Res. Lett., 31(17). van Hinsbergen, D. J., Kaymakci, N., Spakman, W., and Torsvik, T. H. (2010). Reconciling the geological history of western Turkey with plate circuits and mantle tomography. Earth Planet. Sci. Lett., 297(3-4):674–686. van Hinsbergen, D. J. J., Lippert, P. C., Dupont-Nivet, G., McQuarrie, N., Doubrovine, P. V., Spakman, W., and Torsvik, T. H. (2012). Greater India Basin hypothesis and a two-stage Cenozoic collision between India and Asia. Proc. Natl. Acad. Sci. U. S. A., 109(20):7659–64. Van Hinsbergen, D. J. J., Vissers, R. L. M., and Spakman, W. (2014). Origin and consequences of western Mediterranean subduction, rollback, and slab segmentation. Tectonics, 33(4):393–419. van Hunen, J. and Allen, M. B. (2011). Continental collision and slab break-off: A comparison of 3-D numerical models with observations. Earth Planet. Sci. Lett., 302(1-2):27–37. 118 References Vergés, J. and Fernàndez, M. (2012). Tethys–Atlantic interaction along the Iberia–Africa plate boundary: The Betic–Rif orogenic system. Tectonophysics, 579:144–172. Villaseñor, A., Chevrot, S., Harnafi, M., Gallart, J., Pazos, A., Serrano, I., Córdoba, D., Pulgar, J. A., and Ibarra, P. (2015). Subduction and volcanism in the Iberia-North Africa collision zone from tomographic images of the upper mantle. Tectonophysics, 663:238–249. von Blanckenburg, F. and Davies, J. H. (1995). Slab breakoff: A model for syncollisional magmatism and tectonics in the Alps. Tectonics, 14(1):120–131. Watts, A. B., Platt, J. P., and Buhl, P. (1993). Tectonic evolution of the Alboran Sea basin. Basin Res., 5(3):153–177. Wernicke, B. and Snow, J. K. (1998). Cenozoic Tectonism in the Central Basin and Range: Motion of the Sierran-Great Valley Block. Int. Geol. Rev., 40(5):403–410. Wildi, W. (1983). La chaine tellorifaine (Algerie, Maroc, Tunisie): stratigraphie et evolution du Trias au Miocene. Rev. Geol. Dyn. Geogr. Phys., 24(3):201–297. Williams, J. R. and Platt, J. P. (2018). A new structural and kinematic framework for the Alborán domain (Betic–Rif arc, western Mediterranean orogenic system). J. Geol. Soc. London., 175(3):465–496. Wortel, M. J. R. (2000). Subduction and Slab Detachment in the Mediterranean-Carpathian Region. Science (80-. )., 290(5498):1910–1917. Wortel, R. (1982). Seismicity and rheology of subducted slabs. Nature, 296(5857):553–556. Wortel, R. and Spakman, W. (1992). Structure and dynamics of subducted lithosphere in the Mediterranean region. In Proc. Kon. Ned. Akad. v. Wetensch, volume 95, pages 325–347. Proc Kon Ned Akad v Wetensch. Wu, F.-Y., Ji, W.-Q., Wang, J.-G., Liu, C.-Z., Chung, S.-L., and Clift, P. D. (2014). Zircon U-Pb and Hf isotopic constraints on the onset time of India-Asia collision. Am. J. Sci., 314(2):548–579. Yoshioka, S. and Wortel, R. (1995). Three-dimensional numerical modeling of detachment of subducted lithosphere. J. Geophys. Res., 100(B10). Zhao, J., Yuan, X., Liu, H., Kumar, P., Pei, S., Kind, R., Zhang, Z., Teng, J., Ding, L., Gao, X., Xu, Q., and Wang, W. (2010). The boundary between the Indian and Asian tectonic plates below Tibet. Proc. Natl. Acad. Sci. U. S. A., 107(25):11229–11233. Zhong, X. and Li, Z. H. (2019). Forced Subduction Initiation at Passive Continental Margins: Velocity-Driven Versus Stress-Driven. Geophys. Res. Lett., 46(20):11054–11064.