!
! ! ∀!!
##∃∃∃
%!&%&∋( ##∃ ∃#! )
∗! ∃ + +
∃
Acknowledgements
Firstly I want to give my thanks to my main supervisor, Dr. Abigail L. Bull, for always be- ing motivating and excited about this project, giving good guidance and for sharing helpful scripts. Thanks to my co-supervisor, Dr. Grace E. Shephard, for always giving constructive feedback, shearing GMT plotting scripts and providing me with the agegrid files used in the data assimilation. Also, thanks to my co-supervisor, Dr. Carmen Gaina, for suggestions on the thesis along the way.
To everyone at CEED, thank you for being such a welcoming bunch of people, may your beers be forever cold and your glasses always full. As for my fellow students, thanks for good laughs and fun throughout the study.
Lastly, I want to thank my family for all their support, and especially Daniel Angler Valrygg for always being there and cheering me on.
This work was performed on the Abel Cluster, owned by the University of Oslo and the Norwegian metacenter for High Performance Computing (NOTUR), and operated by the De- partment for Research Computing at USIT, the University of Oslo IT-department.
Abstract
To perform any numerical model of convection in Earth’s mantle, initial internal conditions (e.g., temperature and composition), boundary conditions (e.g., the velocity fields at the sur- face and the core-mantle boundary), and parameters such as the mantle viscosity profile must be imposed. Such conditions and parameters are crucial for simulating mantle convec- tion, however, none are fully constrained, and assumptions on their values must be made.
Using tectonic plate motion velocities from global plate reconstructions as time-dependent surface boundary conditions in numerical convection models, it is possible to simulate con- vective flow within Earth’s mantle for periods encompassing several hundreds of millions of years. Recent work has shown that such models can, to first-order, successfully reproduce the global present-day mantle velocity signature inferred from seismic tomography.
In this project, the finite-element mantle convection code CitcomS is used in parallel with the interactive plate tectonic reconstruction software, GPlates, to study the global and regional effects of the assumed initial condition, radial viscosity profile and history of subduction on the resulting 3D flow field in simulations of mantle convection. Both isochemical (i.e., purely thermal) and thermochemical (thermal and compositional) modes of convection in Earth’s mantle are investigated. In addition, we will also explore the time-dependent influence of sub- duction on regional mantle structure, including varying slab dip, subduction zone geometry and subduction polarity. This analysis is applied to investigate ambiguities regarding the mid-Mesozoic closure of the Mongol-Okhotsk Ocean and the resulting present-day observed mantle velocity structure.
This study demonstrates the influence and importance of the radial viscosity profile in nu- merical models on slab behaviour and creation of LLSVP structures in the lowermost mantle, as well as the importance of an earth-like initial temperature condition. Lastly, the location of the slab remnant from the closure of the Mongol-Okhotsk Ocean is discussed through a quantitative and qualitative comparison of the predicted present-day mantle structure to inferences from alternative seismic tomography models.
Contents
1 Introduction 1
1.1 Outline of the thesis . . . 1
1.2 The Earth’s interior and plates . . . 2
1.2.1 The mantle . . . 2
1.2.2 The core . . . 4
1.2.3 The crust and lithosphere . . . 6
1.3 Studying the Earth’s interior . . . 8
1.3.1 Seismic tomography . . . 8
1.3.2 Numerical modelling . . . 10
1.4 The Mongol-Okhotsk subduction zone . . . 11
1.4.1 Tectonic history . . . 12
1.4.2 Plate reconstructions . . . 14
1.4.3 Seismic tomography of the Mongol-Okhotsk slab . . . 15
2 Method 18 2.1 Numerical modelling of mantle convection . . . 18
2.1.1 Governing equations for mantle convection . . . 18
2.1.2 Boundary and initial conditions . . . 21
2.2 Constructing a tectonic plate reconstruction model . . . 21
2.2.1 Absolute and relative reference frames . . . 22
2.2.2 GPlates . . . 23
2.3 Numerical modelling of subduction using data assimilation . . . 26
2.3.1 A priori data files . . . 26
2.3.2 Lithosphere assimilation . . . 28
2.3.3 Slab assimilation . . . 31
2.4 Model setup . . . 31
2.4.1 Viscosity profiles . . . 31
2.4.2 Average initial temperature profiles . . . 32
2.5 Conversion to seismic wave speeds . . . 34
2.6 Numerical comparison of models . . . 36
3 Results 38 3.1 General model observations on a global scale . . . 38
3.1.1 Thermochemical versus isochemical models . . . 38
3.1.2 Models with data assimilation versus no data assimilation . . . 39
3.2 Radial viscosity profiles . . . 39
3.2.1 Main viscosity profiles V1 and V2 . . . 40
3.2.2 Viscosity profile V3 . . . 41
3.3 Average initial temperature profiles . . . 41
3.3.1 Non-Earth-like initial temperature profiles . . . 41
3.3.2 Earth-like initial temperature profiles . . . 42
3.4 The Mongol-Okhotsk slab . . . 43
3.4.1 Subduction polarity and location . . . 43
3.4.2 Slab dip . . . 47
3.5 Synthetic tomography models . . . 49
4 Discussion 53 4.1 The effect of the initial temperature structure . . . 53
4.2 The effect of the radial viscosity profile . . . 54
4.3 The Mongol-Okhotsk slab . . . 55
4.3.1 Behaviour of the Mongol-Okhotsk slab remnant . . . 56
4.3.2 Location of the slab remnant in the synthetic tomography models . . . 57
5 Conclusions and future directions 60 References 63 A Appendix 69 A.1 Model parameters . . . 69
A.2 The South American subduction zone . . . 72
A.3 Model plots . . . 73
List of Figures
1.1 The Preliminary Reference Earth Model (PREM) . . . 3
1.2 The structure of the Earth . . . 5
1.3 Global distribution of earthquakes from 1914 to 1963 . . . 7
1.4 Global seismic tomography of LLSVPs at 2800 km . . . 9
1.5 Topographic map of present-day NE Asia, with location of the Mongol-Okhotsk suture . . . 11
1.6 The tectonic plate history of the Mongol-Okhotsk Ocean . . . 13
1.7 Map of continental regions and plate velocities at 200 Ma for the Mongol-Okhotsk Ocean . . . 14
1.8 Seismic tomography of the Mongol-Okhotsk slab . . . 16
2.1 Example of a plate circuit . . . 24
2.2 Velocities and plate boundary files loaded in GPlates . . . 25
2.3 Flow chart: from reconstruction to analysis of models . . . 29
2.4 The thermal slab stencil used for the data assimilation . . . 30
2.5 Radial viscosity profiles . . . 33
2.6 Averaged horizontal initial temperature profiles . . . 33
2.7 Initial mantle temperature structure for TP1 models . . . 35
2.8 Spherical harmonic power spectrum comparison . . . 37
3.1 Vertical cross-section through the equator, for V1 and V2 models . . . 40
3.2 The Mongol-Okhotsk subduction through time, comparing models with R1 and R2 and a slab dip of 45° . . . 44
3.3 The Mongol-Okhotsk subduction through time, comparing models with R1 and R2 and a slab dip of 80° . . . 45
3.4 Vertical cross-sections of V1 synthetic tomography models . . . 50
3.5 Horizontal cross-sections of V1 synthetic tomography models . . . 52
4.1 The present-day mantle flow field near the Mongol-Okhotsk slab . . . 57
A.1 Seismic tomography of the South American slab . . . 72
A.2 Locations of the vertical cross-sections . . . 73
A.3 Locations of the horizontal cross-sections . . . 73
List of Tables
3.1 Synthetic tomography assessment table . . . 48
A.1 Fixed model parameters . . . 69
A.2 Fixed data assimilation parameters . . . 69
A.3 List of main models . . . 70
A.4 List of alternative models . . . 71
A.5 Plots of main models with viscosity profile V1 . . . 74
A.6 Plots of main models with viscosity profile V2 . . . 75 A.7 Plots of alternative models with viscosity profile V3 and R1 plate reconstruction 76 A.8 Plots of alternative models with viscosity profile V3 and R2 plate reconstruction 77
1 Introduction
The Earth is a complex and dynamic system, and by investigating the link between processes in its interior and structures at the surface we get closer to understanding the evolution of our planet. For instance, the evolution of the density structure of the Earth affects, among others, global sea level, surface geology and Earth’s moment of inertia. A complete understanding of large-scale mantle convection through time is still to come, and is a fundamental goal in solid-Earth geophysics. Numerical models, together with observations and datasets from geophysics and geochemistry, provide us with the means to test our understanding of the Earth’s dynamics and generate new hypotheses which can guide the research forward.
Although a rapid increase in computational power and storage capacity over recent years has opened the possibility for numerical models with higher resolution and more earth-like rheology, global numerical models of mantle convection still have limitations. In order to solve geodynamical problems, initial conditions, boundary conditions, and parameters such as the mantle viscosity profile must be imposed. However, none of these are fully constrained, so assumptions and/or simplifications have to be made. A sound understanding of the physics and chemistry of the Earth’s interior is crucial, as it provides us with the necessary knowledge to choose the most realistic assumptions, and ultimately move towards self-consistent models of mantle convection with Earth-like parameters.
1.1 Outline of the thesis
In this thesis, numerical modelling of mantle convection will be conducted with global plate reconstructions imposed as kinematic boundary conditions at the surface in addition to im- posed slab parameters. The effect of different initial temperature conditions and radial vis- cosity profiles on global structures such as large, hot chemically distinct structures in the lowermost mantle, and the general behaviour of cold downwellings, will be investigated in the models. Additionally, this study will provide a regional investigation on the subduction of a Mesozoic palaeo-ocean; the Mongol-Okhotsk Ocean. The Mongol-Okhotsk Ocean was relatively isolated from the subduction of nearby oceans in an otherwise tectonically complex region, making it ideal for studying the influence of model parameters and regional tectonic reconstructions. The evolution of the Mongol-Okhotsk Ocean as well as the behaviour of the slab is poorly understood. Thus, polarity of the Mongol-Okhotsk subduction will be investig- ated, by testing the effect of northern and southern subduction in the plate reconstructions.
Lastly, the temperature field of the geodynamic models is converted to seismic wave speeds, and filtered to fit the resolution of the tomography model it will be compared to, as a final assessment of the reliability of the models.
The thesis is divided into five main sections; 1) Introduction, 2) Method, 3) Results, 4) Dis-
cussion and lastly 5) Conclusions. Section 1 provides an introductory overview of the chem- istry and physics of Earth’s interior, together with a short summary of the tectonic history of the Mongol-Okhotsk Ocean and interpretations of the Mongol-Okhotsk slab in seismic tomo- graphy from previous studies. Following this, Section 2 provides details on the fundamental equations that are solved in the numerical models, construction of plate reconstructions that are used as surface kinematic boundary conditions, details on the data assimilation method used to constrain the thermal structure of subducting slabs, details on creating synthetic tomography models and an overview of the setup of the models. Sections 3 and 4 present the details of the models describing and discussing similarities between models on a global and regional scale, and relate them to the changing parameters, model type or initial conditions.
Finally, section 5 provides a summary of the findings of this thesis. At the end of the thesis the appendices A.1 and A.3 give tables containing the model parameters and plots of all the models conducted, also a short description of the South American subduction zone is given in Appendix A.2 as this slab has been briefly investigated as well.
1.2 The Earth’s interior and plates
The radial structure of the Earth’s interior is fairly well known from seismology, as changes in its properties with depth (i.e. the density and elasticity of the medium) lead to changes in seismic wave velocities, which in turn lead to reflection and refraction of the seismic waves.
By studying the arrival times of different seismic phases an Earth model such as PREM (Pre- liminary Reference Earth Model) (Dziewonski and Anderson, 1981), which includes averaged Earth properties as a function of depth (Figure 1.1), can be inferred. There are two particularly steep and large jumps in seismic velocity in PREM which relate to the main radial divisions of the Earth; the core, mantle and crust (Figures 1.1 and 1.2), each differing in chemical composition and rheology.
1.2.1 The mantle
The mantle accounts for the largest portion of the Earth, comprising 84% of Earth’s volume (e.g., Robertson, 2011) and 67% of its mass (e.g., Stacey and Davis, 2013). The mantle domain ranges from the core mantle boundary (CMB) at∼2890 km depth to the Mohoroviˇci´c discon- tinuity (Moho; boundary between the crust and the mantle) at∼40 km. The most abundant elements in the mantle are oxygen, magnesium, and silicon, and thus it is mostly composed of the silicate minerals olivine and pyroxene (e.g., Stacey and Davis, 2013). Polymorphic phase changes at 410 km and 660 km depth (Figure 1.1) are present in the mantle, with the 660 km boundary representing the boundary between the upper and lower mantle (e.g., Anderson, 1989).
The mantle is heated internally from the radioactive decay of uranium, thorium and po-
Figure 1.1:The Preliminary Reference Earth Model (PREM) showing the radially averaged P- and S-wave velocities (Vp andVs, respectively) and density (ρ) with depth (modified after Dziewonski and Anderson, 1981). a.) The velocity and density variations from the surface to the centre of the Earth. The transitions between the inner and outer core and the core-mantle boundary (CMB) are marked as solid vertical lines.
b.) Enlarged plot from 0 km (surface) to 1000 km (uppermost lower mantle), here the polymorphic phase changes at 410 km and 660 km, and the crust-mantle transition (Moho) are marked.
tassium (e.g., Turcotte and Schubert, 2002), as well as basally by thermal conduction across the CMB (e.g., Bercovici et al., 2000). At the same time, the mantle is cooled because of thermal conduction through the crust, creating a gravitationally unstable structure of cold (dense) material on top of hot (less dense) material. Because the mantle is able to transmit S-waves (shear waves) (Figure 1.1) it is a solid, but exhibits a fluid behaviour over geological time scales due to solid-state creep (e.g., Turcotte and Schubert, 2002). This fluid behaviour, together with gravitationally unstable structures within the mantle, leads to convection. In turn, this leads to mixing of the hot and cold materials, generating an average temperature throughout most of the mantle, with two transition zones from mantle temperatures to crustal temperatures at the top (UTBL: Upper Thermal Boundary Layer), and core temperatures at the bottom (LTBL: Lower Thermal Boundary Layer). The thickness of the thermal boundary layers are affected by how vigorous the convection is, closely related to the Rayleigh number (Bercovici et al., 2000) which itself is governed, among other parameters, by viscosity and temperature (see Equation 2.17). High Rayleigh numbers will result in vigorous convection, which will result in thinner thermal boundary layers, and vice versa (Bercovici et al., 2000).
However, mantle convection is only permitted for Rayleigh numbers above a certain value, the critical Rayleigh number. For the Earth’s mantle, the Rayleigh number is calculated to be between106 and109 (Turcotte and Schubert, 2002), which is well above the critical Rayleigh number on the order of103(Bercovici et al., 2000), and indicates vigorous convection.
In addition to radial heterogeneity, the mantle is also thought to be laterally heterogeneous, which is clear from seismic tomography (see section 1.3.1). A striking feature across several global whole-mantle seismic tomography models is the presence of two regions with lower- than-average shear-wave speeds in the lowermost mantle. The reduced wave speeds in these regions can possibly be explained by the existence of deep dense reservoirs of material termed LLSVPs (Large Low Shear Velocity Provinces) (e.g., Dziewonski et al., 1977).
1.2.2 The core
Below depths of∼2890 km (CMB) is the Earth’s core, with a radius of 3486 km. It is primarily composed of iron (e.g., Stacey and Davis, 2013), and is the densest part of the Earth. It is divided by the inner core boundary (at 5155 km depth) into a solid inner core and a liquid outer core, apparent from the lack of S-wave propagation through the outer core (Figure 1.1).
Because the core is enriched in iron, and has an convecting outer core, it is a good thermal conductor. In comparison, the mantle is a poor conductor, and convects slower than the core, creating a LTBL (e.g., Anderson, 1989) described in section 1.2.1. The combination of rotation (the Coriolis effect) of the Earth, convection and an electrically conducting fluid in the outer core, makes the core act like a dynamo (Stacey and Davis, 2013) generating the Earth’s magnetic field.
660
5155 Lithosphere
Figure 1.2:The main radial divisions of the Earth (not to scale). Convergent and divergent plate boundar- ies are illustrated with associated features, such as subduction of oceanic lithosphere at subduction zones, and lithospheric thickening at spreading ridges (modified from Stern, 2002).
Throughout Earth’s history the magnetic field has reversed its polarity (i.e. the magnetic north and the magnetic south poles switch). The reversals occur at irregular intervals, and because the magnetic field can be recorded as remnant magnetization in rocks when they cool below the Curie temperature, it can be seen as magnetic anomalies in the oceanic lithosphere.
The magnetic anomalies are “frozen” into the oceanic lithosphere when magma cools at mid- oceanic ridges to form symmetrical magnetic stripes mirrored on either side of the mid-ocean ridge. Remnant magnetization in continental lithosphere due to cooling of molten rocks or reheating above the Curie temperature is important for studying the movement and history of the crust.
1.2.3 The crust and lithosphere
The crust constitutes 0.5% of Earth’s mass and has a high content of silica (SiO2). It is convenient to divide the crust into continental and oceanic domains, as the oceanic crust has a higher concentration of magnesium and iron, and a shorter lifespan than the more silica- rich continental crust. The crust is a seismically defined layer from the surface to the Moho (at∼40 km depth), but there is also a thicker layer called the lithosphere (“rocky layer”) which is defined by its rheological behaviour as it, unlike the mantle, behaves rigidly over geologic time scales. The lithosphere includes the crust and extends down to the uppermost part of the mantle, the asthenosphere (“weak layer”), where the temperature and pressure is high enough for solid-state creep to occur (e.g., Anderson, 1989; Stein and Wysession, 2003).
The lithosphere is divided into several tectonic plates (Figure 1.3), with thicknesses ranging from almost zero at spreading ridges to∼200 km at stable cratons (e.g., Stein and Wysession, 2003). There are approximately 12 plates at present-day (Anderson, 2002), consisting of large plates (e.g. the Pacific and African plates) and smaller plates (e.g. Caribbean plate), which move relative to each other with (RMS) velocities of around 1 - 8 cm/yr (Zahirovic et al., 2013). The relative motion between tectonic plates is responsible for the continuous creation and consumption of lithosphere through time, and is known as plate tectonics. There are three different types of plate boundaries, distinguished by the relative motion between the two neighbouring plates; converging plate boundaries (e.g. along the western coast of South America) with relative motion towards each other, diverging plate boundaries (e.g. the Pacific mid-oceanic ridge) with relative motion away from each other, and transform plate boundaries (e.g. the San Andreas Fault) with relative motion parallel to the boundary.
Creation of oceanic lithosphere, commonly referred to as seafloor spreading, occurs at the oceanic ridges where two adjacent plates move apart, creating a divergent boundary, where magma rises to fill the gap between the plates. As the hot magma ascends, it cools and accretes to the plates to form new oceanic lithosphere. The newly formed lithosphere will continue to cool and thicken as it moves farther away from the ridge making it denser
Figure 1.3: Global distribution of earthquakes from 1914 to 1963. Location of the epicen- ters is plotted as red dots, with present-day global plate boundaries in black and continents in grey. (Earthquake data from Villaseñor et al. (1997); downloaded on 30/03/2015 from http://earthquake.usgs.gov/data/iss_summ.php)
because of thermal contraction (Figure 1.2). To accommodate for new lithosphere created at the ocean ridges, there needs to be an equal amount of consumption of lithosphere. This occurs at convergent plate boundaries, most typically where oceanic lithosphere bends and sinks down into the mantle beneath a less dense overriding plate (Figure 1.2). However, convergent sites can also include the building of major mountain belts such as the Himalayas or Andes. Gravity exerts a pull on the subducting lithosphere, known as slab pull. Because of the elastic behaviour of the lithosphere, the stress from the slab pull is transferred to the surface of the lithosphere, pulling the subducting plate towards the subduction zone (e.g., Turcotte and Schubert, 2002). This pull acts as an important driving mechanism for plate tectonics.
Present-day plate boundaries represent regions with high seismicity (Figure 1.3), especially so in subduction zones where earthquakes occur as a result of the descending lithosphere.
The earthquakes in subduction zones occur both at shallow depths as a result of the stress release between the colliding plates and bending of the plates, and at greater depths in the subducting lithosphere. If earthquake locations are plotted in a vertical cross-section through a subduction zone, the earthquake epicenters preferentially follow the subducted lithosphere and indicate that the dip of the descending slab greatly differs between subduction zones.
The earthquakes cease at depths around 700 km (e.g., Stein and Wysession, 2003), and different methods must be sought in order to study the descending lithosphere, such as seismic tomography (see Section 1.3.1).
1.3 Studying the Earth’s interior
The deepest that humans have been able to directly observe the Earth’s interior is ∼12 km depth from drilling of the continental crust (Kozlovskij and Andrianov, 1987), in other words we have only barely scratched the surface. The extreme temperatures and pressures in our planet’s interior, from ∼800 K and ∼10 GPa in the uppermost part of the mantle to∼ 5000 K and ∼360 GPa in the innermost part of the core (e.g., Stacey and Davis, 2013), make it for most part inaccessible. Although samples of mantle xenoliths have provided constraints on the mantle chemistry, geoscientists must generally turn to remote sensing techniques in order to observe and understand the evolution of the interior of our planet. One important technique is seismic tomography, which over the last few decades has provided crucial insight into the structure and dynamics of the Earth’s interior. Also, numerical models of mantle convection, although not strictly a remote sensing technique, provide means to investigate mantle convection and give further understanding of its evolution.
1.3.1 Seismic tomography
Large magnitude earthquakes (> 6 Mw) release a lot of energy, which is transferred through the Earth as seismic waves. There are two principal types of seismic waves; body waves and surface waves. The surface waves travel in the upper part of the crust, and are the main cause of destruction associated with earthquakes. Body waves, on the other hand, spread out spherically through Earth’s interior, and consist of compressional waves (P-waves) and shear waves (S-waves).
Seismic tomography is used to infer the heterogeneous structure of the mantle by studying seismic waves that have travelled through it. This is similar to the non-invasive medical imaging method CT (computed tomography), which uses X-rays sent through the human body to construct an image of the internal structure. It is convenient to look at the seismic waves as rays emanating from the center of the earthquakes, rather than a wave front, when explaining seismic tomography. As the rays emanating from the earthquake travel through the heterogeneous earth, the rays will be reflected and refracted differently, causing them to follow different paths through Earth’s interior from the source (earthquake) to the receiver (seismometer). The arrival time at the receiver will depend on the ray path and wave speed of the medium the rays are propagating through. By comparing the arrival times of individual rays with the expected arrival times derived from a reference model (e.g. PREM), it is possible to determine if the ray arrives earlier or later than expected. A ray that arrives earlier than expected would likely have travelled through a seismically fast region resulting in an positive anomaly, while a ray arriving later than expected would have travelled through a seismically slow region resulting in a negative anomaly.
The observed wave speed anomalies can be interpreted as density anomalies, making it
Figure 1.4: Horizontal cross section of global seismic tomography at 2800 km showing the LLSVPs under Africa and central Pacific (modified from Mégnin and Romanowicz, 2000)
possible to map the heterogeneous structure of the mantle, both on a local and global scale.
Seismic tomography images actively subducting slabs well, particularly in the uppermost mantle, as the downgoing lithosphere is colder than the surrounding mantle, hence also denser, and imaged as higher-than-average seismic wave speeds (positive anomalies). Seis- mic tomography of subducting slabs led to the important observation that slabs are able to penetrate into the lower mantle (e.g., Spakman et al., 1988; Sigloch et al., 2008). In fact, higher-than-average seismic wave speeds, thought to relate to cold slabs from ancient sub- duction, are seen as deep as the CMB (e.g., van der Hilst et al., 1997; Grand et al., 1997).
These ancient slabs create a higher-than-average velocity belt beneath Asia and the Americas, encircling the Pacific (Figure 1.4).
Global seismic tomography suggest that regions with lower-than-average seismic wave speeds (e.g., Dziewonski et al., 1977; Mégnin and Romanowicz, 2000) exist near the CMB beneath Africa and the central Pacific, the LLSVPs (Figure 1.4). The higher-than-average velocity belt encircling the LLSVPs creates a degree 2 pattern (i.e. two antipodal anomalies) in the lowermost mantle (∼2800 km). The nature and evolution of these lower mantle structures is currently debated; whether they are purely thermal structures (i.e. hot upwellings at the CMB) or if they also are chemically distinct from their surroundings (e.g., Bull et al., 2009) is unclear. Furthermore, questions remain as to the behaviour of the LLSVPs throughout Earth’s history; whether they are relatively stable structures over geologic time scales or short lived structures, is also a subject of discussion.
Seismic tomography provides a snapshot of the present-day mantle structure, hence it provides us with crucial constraints about the interior of our planet. However, tomography alone can not be used to infer the evolution of the mantle in the past or in the future. In order to investigate the time-dependence of the system, tomography must be used in conjunction with other methods and datasets, such as numerical modelling, plate tectonic reconstruc-
tions and deep earth geochemistry. It must be noted that because of the uneven distribution of earthquakes and seismometers, the Earth’s interior is not sampled evenly. As a result, seismic tomography models have an uneven resolution throughout the mantle, the southern hemisphere and the highest latitudes are particularly poorly resolved.
1.3.2 Numerical modelling
Numerical modelling is an excellent tool to study phenomena that are not easily studied by laboratory experiments or by collection of data directly. In the field of computational geody- namics, computers are utilized to study the Earth’s interior and its dynamics, for example by modelling global mantle convection.
Numerical models of global mantle convection are based on a set of partial differential equations that govern geodynamic processes, which are derived from basic conservation laws (see Section 2.1.1). The solution of these equations is dependent on the boundary and initial conditions of the problem. The boundary conditions relate to the behaviour of the boundaries in the model, and in 3D spherical models of mantle convection the boundaries correspond to the surface and the CMB. Tectonic plate velocities are often imposed as a kinematic boundary condition at the Earth’s surface, while the CMB has a free-slip kinematic boundary condition (Section 2.1.2). The initial condition defines the initial state of the model, and determines the values which go into the equations at the onset of the calculation. The temperature structure is one of the initial conditions that must be prescribed. In forward models of mantle convection, such as the one used in this study, the model is started at a certain age in the past and run to present-day, hence the initial condition is not known. A rheology model (i.e.
viscosity), and a set of model parameters have to be chosen as well.
Previous studies have utilized numerical models of mantle convection in the aim of ex- plaining inferences from seismic tomography. For instance, isochemical and thermochemical studies of mantle convection have shown a correlation between cold downwellings related to subduction with the high velocity belt in the lower mantle (e.g., Kellogg et al., 1999; McNamara and Zhong, 2004, 2005; Bull et al., 2009; Bower et al., 2013). The thermochemical models seem to better predict the location and geometries of LLSVPs than the isochemical models, however, Bull et al. (2009) point out that this could be due to the models being compared directly to seismic tomography without taking into account the low resolution of tomography relative to the high resolution of the numerical models. Recent studies have also looked at distinguishing ambiguities in plate reconstructions and understanding the effects of mantle convection on dynamic topography, by using global mantle convection models where the plate velocities, thermal structure of the lithosphere and slabs (down to a specific depth in the mantle) are consistent with geophysical and geological data (e.g, Shephard et al., 2013; Fla- ment et al., 2014).
1.4 The Mongol-Okhotsk subduction zone
80 ˚ 80 ˚
90 ˚
10 0˚
11 0˚
12 0˚ 13
0˚ 14
0˚
40˚
50˚
MONGOLIA
CHINA RUSSIA
LAKE BAIKAL
SEA OF OKHOTSK
..
ULAANBAATAR
ADAATSAG
−4000 −2000 0 2000 4000 6000 8000 10000 Topography ETOPO2 (m)
Figure 1.5: Topographic map of present-day NE Asia with the location of the Mongol-Okhotsk suture in purple, stretching from middle Mongolia to the Sea of Okhotsk, and political borders in grey. The location of the Adaatsag ophiolite studied by Tomurtogoo et al. (2005) is also plotted. Suture zone location was digitized from the map of Tomurtogoo et al. (2005). Global topography data from ETOPO2v2 (U.S. Depart- ment of Commerce, National Oceanic and Atmospheric Administration, National Geophysical Data Center, 2006).
The Mongol-Okhotsk Ocean has been chosen as a regional case study in this thesis be- cause of the uncertainties regarding the history of the ocean, which will presented in this section, and behaviour of the slab remnant after the final closure of the ocean. Additionally, from a modelling point of view, it is ideal for investigating the behaviour of the subducted slabs, and influence of parameters and initial conditions in numerical models, as it is a long- lived, linear subduction relatively isolated from other subduction zones.
The Mongol-Okhotsk Ocean existed in the Palaeozoic (542-251 Ma;Gradstein et al. (2004)) and Mesozoic (251-66 Ma) eras between the continental blocks of Siberia to the North and the Amuria (or Mongolia) block to the South (Figures 1.6 and 1.7)1. It was separated from Panthalassa (pre-Pacific ocean basin) to the East by an island arc and subduction zone (Figure 1.6). The existence of this palaeo ocean is evident from the Mongol-Okhotsk suture stretching
1Unless specified, the relative locations of terranes and boundaries are discussed relative to the present-day locations.
from central Mongolia, following the boundary between the Siberian and Amuria blocks, north to the Sea of Okhotsk (Figure 1.5). The suture zone contains ophiolites (obducted oceanic lithosphere) (Zonenshain et al., 1990; Tomurtogoo et al., 2005) and sediments containing marine fossils which confirms the existence of an ocean in the Palaeozoic to Mid-Late Jurassic (Halim et al., 1998, and references therein). The continuation of the suture towards the west is not clear and has lead to alternative subduction scenarios (see Section 1.4.1).
Further evidence for the existence of the Mongol-Okhotsk Ocean can be deciphered from the present-day underlying mantle, as a slab remnant in seismic tomography. The Mongol- Okhotsk slab is inferred to be visible as a seismically fast anomaly near the CMB beneath Siberia (Figure 1.8) (van der Voo et al., 1999; van der Meer et al., 2010; though perhaps shifted longitudially, see Shephard et al., 2013). See section 1.3.1 for a detailed description of seismic tomography.
1.4.1 Tectonic history
Constraining plate reconstructions is increasingly difficult as you go further back in time due to the continual production and destruction of oceanic lithosphere by plate tectonics. Today, limited oceanic lithosphere before ∼180 Ma (Early Jurassic) remains, and plate reconstruc- tions rely primarily on continental rock data, such as palaeomagnetics and the geological record. It follows that the opening timing of the Mongol-Okhotsk Ocean is debated. Estimates range from Cambrian (∼540-490 Ma) (Harland et al., 1990), Ordovician (∼485-445 Ma) (Cocks and Torsvik, 2007) to Permian (∼300-250 Ma) (Zorin, 1999; Kravchinsky et al., 2002), but Seton et al. (2012) point out that part of the timing ambiguity is related to different definitions of the Mongol-Okhotsk Ocean, i.e. at which time the associated terranes were in the config- uration needed to form the ocean. Dating of the Adaatsag ophiolite located on the western part of the suture (Figure 1.5), reveals an age of∼325 Ma (Tomurtogoo et al., 2005), indicating active seafloor spreading in the Carboniferous. This narrows the estimate of the initial ocean basin opening to Early-Middle Palaeozoic (Domeier and Torsvik, 2014).
Evidence of subduction related magmatism is found on both sides of the Mongol-Okhotsk suture (Zorin, 1999) indicating that subduction occurred underneath both the Siberian and Amurian margins. This has been interpreted as bivergent subduction where both margins were active during an overlapping time period (Zorin, 1999). Subduction beneath the southern Amurian margin is thought to have ceased in the Late Permian and Triassic (∼260-200 Ma) as the North China block accreted to the Amuria block, resulting in a passive Amurian margin during the final stage of the Mongol-Okhotsk Ocean closure. Final consumption of the ocean is therefore thought to have been along its northern margin and was completed when the combined Amuria-North China blocks collided with Siberia (Zorin, 1999; Zonenshain et al., 1990).
Figure 1.6:The tectonic plate history of the Mongol-Okhotsk Ocean at selected ages based on the R1 plate reconstruction (Seton et al., 2012). Reconstructed continental regions of Eurasia (purple), Amuria (yellow) and North China (green); other global terranes are coloured in grey and oceanic domains are white. The main Mongol-Okhotsk subduction zone for R1 is represented in red, and subduction for R2 is represented in blue. Other global subduction zones are shown in black. EUR: Eurasia (Siberia) block, AMU: Amuria block, NCH: North China block, MOK: Mongol-Okhotsk Ocean, PAN: Panthalassa and TEH: Tethys Ocean.
−180˚
−150˚
−120˚
−90˚
−60˚
−30
˚
0˚
30˚
60
˚
90˚
120˚
150˚
−180˚
−150˚
−120˚
−90˚
−60˚
−30
˚
0˚
30˚
60
˚
90˚
120˚
150˚
−180˚
−150˚
−120˚
−90˚
−60˚
−30
˚
0˚
30˚
60
˚
90˚
120˚
150˚
−180˚
−150˚
−120˚
−90˚
−60˚
−30
˚
0˚
30˚
60
˚
90˚
120˚
150˚
Figure 1.7: Map of the Mongol-Okhotsk Ocean based on the R1 plate reconstruction at 201 Ma. The velocities are plotted relative to Eurasia (purple) as this gives more intuitive velocities with respect to the subduction. Note that with northern subduction (red) the Amurian margin is passive, while with southern subduction (blue) the Mongolian margin is passive. Colouring of terranes and subduction zones is the same as in Figure 1.6.
The final closure is suggested to have occurred sometime between the Late Jurassic (∼155 Ma) and beginning of the Early Cretaceous (∼120 Ma) (Zonenshain et al., 1990; Kravchinsky et al., 2002; van der Voo et al., In review). In a palaeomagnetic study of northern China (Zhao et al., 1990), a gradual closure of the Mongol-Okhotsk Ocean was proposed as a result of counterclockwise rotation of Amuria relative to Siberia of 117°. The palaeomagnetic data from the study also indicated that the closure started in the west and ended in the east due to the coincidence of the Late Permian, Early Triassic and Late Jurassic poles of rotation (Zhao et al., 1990). This is supported by intrusions and marine fossils found within the suture which young from west to east (Zhao et al., 1990; Zonenshain et al., 1990; Halim et al., 1998;
Tomurtogoo et al., 2005). Kravchinsky et al. (2002) also suggested a rapid increase in subduc- tion velocity in Late Jurassic, due to large differences in palaeolatitude in the palaeomagnetic data from the Trans-Baikal area, resulting in the final closure of the ocean.
1.4.2 Plate reconstructions
As the tectonic history of this region is ambiguous, several plate reconstructions have been proposed. Two alternative, end-member plate reconstructions will be considered in this study;
one by Seton et al. (2012), and a second plate reconstruction with an alternative subduction history. The Seton et al. (2012) reconstruction is a global plate model from 200-0 Ma with
fully-closing plates (Gurnis et al., 2012) in arbitrary single million year increments. Here, we use an extended digital plate reconstruction back to 230 Ma (beginning of Late Triassic), consistent with the interpretations discussed in Seton et al. (2012) and with updated regional reconstructions elsewhere (e.g. Arctic from Shephard et al. (2013); South East Asia from Zahirovic et al. (2013)). The absolute reference frame is based on moving hotspots for 0-70 Ma (Torsvik et al., 2008) and True Polar Wander (TPW)-corrected palaeomagnetics for 105- 230 Ma (Steinberger and Torsvik, 2008), see section 2.2 for a detailed explanation of TPW and absolute reference frames. As the Mongol-Okhotsk subduction zone is the main region of interest in this thesis, a short summary of the plate reconstructions in this region is included (see Seton et al., 2012, for further details).
In the plate reconstruction by Seton et al. (2012), herein referred to as R1, the initial open- ing of the Mongol-Okhotsk Ocean is implemented in the Late Carboniferous (∼310 Ma) with continued seafloor spreading until 250 Ma (Early Triassic). Subduction along the northern Siberian margin started in the Late Permian (∼260 Ma) and final closure of the ocean occurred at 150 Ma (Late Jurassic). In contrast to other studies (e.g., Zorin, 1999), this reconstruction only includes a single location of subduction, rather than the contemporaneous subduction on both sides of the ocean. In this study, the geodynamical models based on R1 and R2 are started at 230 Ma, meaning that a subduction history is implemented from 230-150 Ma for the Mongol-Okhotsk Ocean, however, the initial slab depth can also be refined, implying an earlier history of subduction (see Section 2.3). Note, that as a global geodynamic model, subduction in different regions is also implemented within the 230-0 Ma time period.
The alternative plate reconstruction, herein referred to as R2, has a different subduction polarity for the Mongol-Okhotsk Ocean. Instead of northerly subduction under Siberia, R2 has subduction located adjacent to, and dipping under, Amuria (Figure 1.6) from 230 Ma to 150 Ma. This alternative southern subduction is permitted by the evidence of subduction on both margins (discussed in section 1.4.1), and allows a comparison between purely north- ward and purely southward subduction. Other than the difference in polarity, the R1 and R2 reconstructions are identical; the timing of Mongol-Okhotsk subduction as well as the abso- lute reference frame and regional plate reconstructions elsewhere on the globe are consistent between R1 and R2.
1.4.3 Seismic tomography of the Mongol-Okhotsk slab
Global seismic tomography models image a positive seismic anomaly under Siberia located between 1500 (mid-mantle) and 2880 km (CMB). This anomaly has been interpreted as the Mongol-Okhotsk slab (e.g. van der Voo et al., 1999; van der Meer et al., 2010) (Figure 1.8).
The identification and interpretation of these deep mantle structures are important in linking deep earth processes with tectonic processes at the surface.
S40RTS
GypsumS
M
M
M M
M M
1500 km 1900 km
2300 km 2700 km
S40RTS
−8000 −4000 0 4000 8000 Topography ETOPO2 (m)
Figure 1.8: Vertical and horizontal slices through seismic tomography illustrating the Mongol-Okhostk slab “M,” as inferred by van der Voo et al. (1999). Top panels, vertical slices (location shown on inset) through S40RTS (Ritsema et al., 2011) and GypsumS (Simmons et al., 2010) models and horizontal slices, through S40RTS, at selected depths between 1500 and 2700 km depth.
For example, assuming knowledge of the timing of subduction at the trench, and that slabs, once detached, sink vertically in the mantle, slab sinking rates in the mantle can be derived. Thus positive anomalies in seismic tomography can be linked to cold remnants of subducted lithosphere of a discrete subduction age. van der Meer et al. (2010) recently applied this approach on a global scale and derived a global average mantle sinking rate of
∼1.3 cm/yr. These assumptions make it possible to track remnant slabs vertically in seismic tomography, calculate the approximate age of the slabs (time of their subduction) at different depths, and project their location to subduction zones at the surface. In seismic tomography, the top of the Mongol-Okhotsk slab (i.e. the last subducted lithosphere before closure) is located at∼1500 km depth (Figure 1.8). Applying the analysis of van der Meer et al. (2010) to the Mongol-Okhotsk, a slab located at 1500 km, with a sinking rate of 1,3 cm/yr, infers a time of subduction at the surface at∼115 Ma (Early Cretaceous). This is not consistent with the plate reconstructions used in this study, which model final closure at 150 Ma. Assuming the top of the Mongol-Okhotsk slab is located at 1500 km depth, a final closure at 150 Ma would indicate a sinking rate of ∼1.0 cm/yr for the Mongol-Okhotsk slab, which is within global estimate ranges (van der Meer et al., 2010; Butterworth et al., 2014). It should be noted that slab stalling in the transition zone may also affect the average sinking rates, thus they are only used as a rough guide.
In Shephard et al. (2013), mantle structure inferred from seismic tomography was used as a constraint to refine plate kinematic models for the Circum Arctic, specifically for con- straining the timing and opening of ocean basins. Building further on this, Shephard et al.
(2014) performed numerical models of mantle convection enabling the possibility to track the slabs through time, and aid in determining the origin of the slab remnant. The authors found that the Mongol-Okhotsk slab might be located farther west (no farther east than 35°) than suggested in previous studies (van der Voo et al., 1999; van der Meer et al., 2010), and that the previously interpreted Mongol-Okhotsk slab from seismic tomography instead is the NW Panthalassa slab (Shephard et al., 2014). Whilst noting a possible component of net litho- spheric rotation (rotation of the lithosphere relative to the mantle), uncertainty in the mantle modelling parameters and subduction history, this alternative tomography-slab interpretation should be kept in mind.
2 Method
2.1 Numerical modelling of mantle convection
In order to model mantle convection we will use the 3D spherical finite element code CitcomS (Zhong et al., 2000). In this study, we will employ two different versions of CitcomS to perform different models for the same case (i.e. parameters, initial condition and boundary conditions are the same), which will be compared to each other.
2.1.1 Governing equations for mantle convection
Numerical modelling of mantle convection is performed by solving a set of equations for the conservation of mass, momentum and heat. The only significant density variation in thermal convection is caused by thermal expansion. Using the Boussinesq approximation, which states that the density difference is sufficiently small so it can be neglected except in cases where it is multiplied by the acceleration of gravitation, the problem is simplified to that of incompressible flow. This will result in the following dimensional and non-dimensional conservation equations.
Dimensional conservation equations
For an incompressible fluid the steady-state conservation of mass is given by:
∇ ·u=0 (2.1)
whereuis the velocity vector.
The steady-state conservation of momentum is given by:
∇p= ∇ ·σ+ρgzˆ (2.2)
wherep is pressure,σis the stress tensor,g is the acceleration of gravity,ρ is the density andzˆ is the unit vector in the depth direction. In equation 2.2 the terms on the left hand side represent the pressure forces, whilst those on the right hand side represent the viscous and buoyancy forces respectively. The density is dependent on temperature, and can be defined in two different ways depending on whether the calculation is isochemical (i.e., convection is dominated by thermal forces) or thermochemical (i.e., convection has both a thermal and chemical component). For the isochemical case, we have:
ρ=ρ0−αρ0(T−T0) (2.3)
and for the thermochemical case:
ρ=ρ0−αρ0(T−T0)+∆ρC (2.4)
whereρ0is the reference density at the reference temperatureT0,αis the thermal expans- ivity, T is the temperature andC∈[0, 1]is the compositional coefficient. It is also possible to separate the pressure term into hydrostatic (p) and dynamic (P) components:
p=P+ρ0g z ⇒ ∇p= ∇P+ρ0gzˆ (2.5)
the expression forρ in the isochemical case then be written as:
− ∇P+ ∇ ·(ηε)˙ =αρ0g(T−T0) ˆz (2.6)
and for the thermochemical case:
− ∇P+ ∇ ·(ηε˙)=αρ0g(T−T0) ˆz+∆ρC gzˆ (2.7)
whereηis the viscosity andε˙ is the strain-rate tensor.
The steady-state conservation of heat is given by:
ρCP∂T
∂t +ρCP(u· ∇)T= ∇(k· ∇T)+H (2.8)
whereCP is the specific heat capacity, T is temperature, t is time, k is the thermal con- ductivity andH is the heat production.
Non-dimensional conservation equations
It is common in numerical modelling to use non-dimensional forms of the conservation equa- tions in order to avoid the wide range of scales for different quantities that emanates when working with these kind of problems. The fundamental quantities in the conservation equa- tions are time (t), length (l), mass (m), temperature (T) and thermal expansivity (α). The non-dimensional scaling factors for these quantities are, respectively:
t=h2 κ0
t′ (2.9)
l=hl′ (2.10)
m=η0h3 κ0
m′ (2.11)
T=∆T T′+T0 (2.12)
α=α0α′ (2.13)
where h is the mantle thickness (m), κ0 is the reference thermal diffusivity (m2s−1), η0 is the reference viscosity (Pa s), ∆T is the temperature drop across the mantle (K), T0 is the reference surface temperature (K) and α0 is the reference thermal expansivity (K−1). The non-dimensional parameters are indicated as the primed variables, and the corresponding variables that are not primed indicate the dimensional parameters. By applying the scaling on the dimensional conservation equations the non-dimensional conservation equations for mass, momentum and heat become:
∇ ·u=0 (2.14)
− ∇P+ ∇ ·(ηǫ˙)=(RaT−RbC) ˆz (2.15)
∂T
∂t +(u· ∇)T= ∇(k· ∇T) (2.16)
All variables are now non-dimensional, and the notation with primed variables for indic- ating dimensionless values are not needed, and here by dropped. The variables Ra and Rb represent thermal and the chemical Rayleigh numbers respectively, and are defined as:
Ra=ρ0gα0∆T h3 κ0η0
(2.17) Rb=∆ρg h3
κ0η0
(2.18)
These are dimensionless numbers which are associated with buoyancy driven flow. In the chemical Rayleigh number, equation (2.18), ∆ρ represent the density difference between two materials, representing the ambient mantle and a denser component, and is therefore only present in the thermochemical calculations. In this study, the denser component is initialised as a 255 km thick layer above the CMB.
The non-dimensional conservation equations of mass (eq. 2.14), momentum (eq. 2.15) and energy (eq. 2.16) will be solved using the 3D spherical finite-element convection codeCitcomS
(Zhong et al., 2000) for isochemical and thermochemical time-dependent problems of mantle convection.
2.1.2 Boundary and initial conditions
The governing equations listed in Section 2.1.1 describe the flow of a highly viscous fluid, and can be used to solve not only global mantle convection problems but any problem concerning the flow of viscous fluids in the mantle. Thus, it is the boundary and initial conditions that prescribe the particular solution to be obtained from the equations.
The boundary conditions describe how the boundaries of the computational domain be- have. In this work, the boundaries are spherical shells corresponding to the CMB and the surface. Free-slip kinematic boundary conditions, meaning non-zero tangential velocities at the boundary and no mass flow in or out of the boundary, are applied at the CMB, while Dirichlet kinematic boundary conditions are applied at the surface. When using Dirichlet boundary conditions a specific value is assigned to each boundary node, which in this study is the plate velocity extracted from the plate reconstructions. Both the top (surface) and bottom (CMB) boundary in the models have Dirichlet thermal boundary conditions, where the top boundary nodes are assigned non-dimensional temperature values ofT=0and the bottom boundary nodes are assigned non-dimensional temperature values ofT=1.
Additionally, in order to solve partial differential equations, initial conditions have to be given together with the boundary conditions. Initial conditions specify the values for the entire computational mesh at a specific initial time, corresponding to 230 Ma (or model time step 0) in this study. For mantle convection problems the partial differential equations are of first order, which means only one initial condition, for each evolving variable, is needed. Both a thermal initial condition and a compositional initial condition are needed for thermochemical calculations, while only thermal initial conditions are needed in the isochemical calculations.
The horizontally averaged initial temperature profile and initial thermal mantle structures are shown in Figure 2.7.
2.2 Constructing a tectonic plate reconstruction model
In this study, plate reconstructions will be imposed as kinematic boundary conditions at the surface in terms of both plate velocities and subduction parameters. In order to achieve this, the plate reconstructions must have fully closed plate polygons, meaning that all oceanic and continental lithosphere has to be accounted for throughout the desired time period and plate boundaries must intersect.
2.2.1 Absolute and relative reference frames
When describing the motion of an object, it has to be described relative to another object or feature. In the case of global plate reconstructions the motion of a plate can be relative to another plate, referred to as relative plate motion, or a fixed reference frame, referred to as absolute plate motion.
Relative plate reconstructions for the past∼175 Myr (e.g., Torsvik et al., 2008) can be cre- ated from marine geophysical data, by identifying ocean floor magnetic anomalies and fracture zones, and linking them to the time of creation. The magnetic anomalies correspond to re- versals of the polarity in Earth’s magnetic field (Section 1.2.2), thus the width of each stripe is dependent on how fast the plate moves away from the spreading ridge and the time period between the polarity reversals. The timing of the reversals is acquired from a geomagnetic polarity timescale and the width of the anomalies can be extracted from geophysical meas- urements, hence half-spreading rates can be determined. Relative plate reconstructions, for geological recent times, may also be constructed by restoring conjugate margins and fitting continent-ocean boundaries together (e.g. fitting the western coast of South America with the eastern coast of Africa), which may be strengthened by evidence of matching palaeofauna and geophysical features. However, these additional constraints are used to complement the identification of magnetic anomalies and fracture zones which are the primary constraints on relative plate motions for current ocean basins.
However, when reconstructing relict ocean basins, or for spreading systems where one or both of the flanks have been subducted, additional onshore datasets such as geological re- cords (i.e. mapping of sutures, ancient or active volcanic arcs and terrane boundaries) and palaeomagnetics, together with simplifications such as assuming symmetric seafloor spread- ing, are needed.
Absolute reference frames describe plate motions relative to a fixed reference system such as the Earth’s spin axis, and may be based on hotspot tracks, assumed to be fixed or moving, or palaeomagnetic data, with or without corrections of true polar wander. The absolute ref- erence frame is a crucial part of the plate reconstruction model, as use of different reference frames will affect the resulting mantle heterogeneity structure (Shephard et al., 2012).
Dating of hotspot tracks, such as the Hawaiian seamount chain, provide latitude and longitude constraints on plate motions, and serve as absolute reference frames with respect to the mantle. In absolute reference frames based on fixed hotspots, the location of the hotspot with respect to the mantle is assumed to be the same through time. While, in absolute reference frames based on moving hotspots, the location of the hotspot, with respect to the mantle, is assumed to be affected by mantle convection. Hotspot tracks cannot be used to derive absolute plate motions prior to∼130 Ma (e.g., Torsvik et al., 2008), as this is the oldest available data.
Palaeomagnetic data (i.e. data from remnant magnetisation in rocks) can be used through- out the geological history, as long as there are available rocks, and provide latitude and ro- tation constraints on plates. Longitude is not constrained using palaeomagnetic data alone, however, a recent methodology linking kimberlite and LIP eruption sites with the long-term stability of LLSVPs has been suggested to provide longitudinal control (Torsvik et al., 2014).
Absolute reference frames based on palaeomagnetic data are with respect to Earth’s spin axis.
The motion of Earth’s plates (relative to Earth’s spin axis) may either be due to the actual mo- tion of individual plates, or because of true polar wander (TPW). TPW is related to drifting of Earth’s solid exterior as a result of an imbalance of weight relative to Earth’s spin axis (Stein- berger and Torsvik, 2008), which will cause rotation of all plates relative to the spin axis over geologic time scales. The palaeomagnetic data will be affected by TPW, and one has to decide if this should be corrected for in the reference frame or not.
In this study a hybrid reference frame is used, which merges a moving hotspot frame (Torsvik et al., 2008) for times between 100 and 0 Ma with a palaeomagnetically derived TPW corrected reference frame (Steinberger and Torsvik, 2008) for older times. The absolute reference frame is linked to the African plate, which is the plate all the major tectonic plates are linked to in the global plate circuit. This choice is due to the inference that the African plate has moved the least in longitude for the past 170 Ma (around the breakup of Pangaea) (Torsvik et al., 2008), and thus by describing all plate motions relative to Africa the uncertainty from lack of longitude in the palaeomagnetic data is minimized.
2.2.2 GPlates
GPlates is an open-source computer program which enables both interactive visualization and the modification of plate tectonic reconstructions (http://www.gplates.org; Boyden et al., 2011). Functionality that ensures continuously closed plate polygons in the plate reconstruc- tions for all times (Gurnis et al., 2012) is built into GPlates, which is crucial for the use of plate reconstructions in computational geodynamics. Thus, CitcomS formatted velocity files can be extracted from the plate reconstructions using GPlates, and used as time-dependent kinematic boundary conditions in the models.
GPlates require a rotation file (.rot) and vector file(s), in order to visualize and modify the plate reconstructions. A detailed description how to use GPlates can be found in the tutorials and manual on the website (http://www.gplates.org).
The rotation file
The rotation file contains information on how the geometries and features move and rotate through time in the reconstruction. The rotation file is composed of six columns plus optional comments at the end of the lines. The first and last (6t h) column contains the “moving” and
Figure 2.1: Part of the plate hierarchy based on Seton et al. (2012), showing how the top plate, the North-East African plate (NEA), is linked to the Mongol-Okhotsk Ocean (MOK).
The numbers represent the plate ID’s, with ac- ronyms in brackets. The North-East African plate is linked to the absolute reference frame, with Plate ID 001, which in turn is linked to Earth’s spin axis at the top. NEA: North-Easth Africa, AFR: Africa, ANT: Antarctica, SAM:
South America, NWA: North-West Africa, NAM:
North America, EUR: Eurasia, SIB: Siberia, MNG: Mongol, NCH: North China and MOK:
Mongol-Okhotsk Ocean.
“fixed” Plate ID’s, respectively. The Plate ID’s are non-negative integer values assigned to tectonic elements (e.g. plates or island arcs), and is used in the rotation file to define the movement of one tectonic element (“moving” Plate ID) relative to another tectonic element (“fixed” Plate ID).
The relative motion between two plates (tectonic elements) can be described using Euler’s theorem, which states that any point or line on the surface of a sphere (e.g. the Earth’s surface) can be shifted to any other desired position on the surface by rotation about the appropriate axis. This is known as Euler rotation, which in the rotation file is defined by the position of the rotation pole (the axis point of intersect with Earth’s surface) in latitude and longitude, given in the third and fourth column, and the angle of rotation, in the fifth column.
The second column specifies the time the rotation takes place in the plate reconstruction (e.g.
at 230.0 Ma). To describe the motion of one specific plate through time, several different rotation poles are needed (finite rotations). The motion of a plate is described relative to a
“fixed” Plate ID, though the “fixed” plate can change through time, and can be derived for any plate pair via a plate hierarchy.
The rotations of the plates in the rotation file relate to each other through a plate hierarchy or circuit (Figure 2.1). This allows for all the plate motions to be linked back to one Plate ID, corresponding to the absolute reference frame, at the top of the hierarchy. In the plate reconstruction used in this study, all plates are linked to the anchored Plate ID through the African plate, as this gives a longitude control for the palaeomagnetic data (Figure 2.1).
Figure 2.2:A screenshot of GPlates showing present-day velocities and plate boundaries. There are9×9 velocity domain points (yellow dots) per cap in this picture.
Plate boundary files
Together with the rotation file, GPlates requires additional files (e.g. shapefiles or GPML files [GPlates Markup Language]) that link the Plate IDs to defined tectonic features such as plate boundaries, terrane features and other vector data. A vector file, such as that containing plate boundaries and polygons, contains information including line geometries, assigned plate ID’s, the age span of the features and what type of feature it is (e.g. plate polygon or spreading ridge). Information concerning other line attributes, such as subduction dip, age and depth, can also be included in the files, and may be used for assimilating the thermal structure of slabs in numerical models of mantle convection (see Section 2.3).
To modify the base R1 plate reconstruction (Seton et al., 2012) to R2 plate reconstruction (Section 1.4.2) several changes had to be implemented. These included the removal of the northward dipping subduction zone by making it inactive, and addition of a southern dipping subduction for the Mongol-Okhotsk subduction zone. The southern subduction was added by digitizing a polyline along the northern continental boundary of the Amuria terrane, which represents the active subduction zone for Mongol-Okhotsk in R2, for the desired time period.
Then the involved plate polygons (Mongol-Okhotsk, Amuria and Eurasia) were reclosed based on the new plate boundaries between 230 and 150 Ma. From these new files, and using a similar process, the age grids describing the age of the oceanic lithosphere have been created from synthetic isochrons (line of equal age), and had to be generated specific to both R1 and R2.
The reconstructed and resolved plate boundary data, and age grids, were exported as topology line data files in 1 Myr increments in GPlates, and later used in the creation of CitcomS input files for the data assimilation (Section 2.3).
Velocity files
As described in Gurnis et al. (2012), the combination of a file defining time-dependent plate polygons, a rotation file and a loaded grid mesh will generate the plate velocities for the mantle convection models (Figure 2.2). The generation of the mesh and exporting of the velocity files is included within the GPlates functionality. The mesh must be generated according to the lateral mesh used in numerical models (CitcomS), which in this study is 129x129 equally spaced nodes for each cap (12 caps). The velocities were exported in 1 Myr increments, for both R1 and R2.
2.3 Numerical modelling of subduction using data assimilation
Numerical models of mantle convection with plate velocities imposed as kinematic boundary conditions at the surface often generate downwellings with symmetric dip (i.e., double-sided subduction). This is not consistent with the asymmetrical subduction observed in seismic tomography. In order to create more Earth-like asymmetrical subduction, which honours geophysical and geological data near the surface, Bower et al. (2015) developed a method to include asymmetrical slabs in numerical models of mantle convection, and modified CitcomS (Section 2.1) to solve the conservation equations (Section 2.1.1) with data assimilation. The method allows for constraints on the subducting slab in the upper mantle according to a plate reconstruction model.
2.3.1 A prioridata files
The first step to run numerical models with data assimilation in CitcomS, is to createa priori data files. Three different datasets are needed; 1) a set of initial condition, “ic”, files which define the initial temperature structure for the entire computational domain at the first time step; 2) a set of “age” files which define the thermal age of the lithosphere at each surface node for each time step; and 3) a set of “hist” files which define the thermal structure of the slabs in the upper mantle for each time step.
All three datasets are created with output from GPlates, and written as CitcomS-formatted files which are used to perform the data assimilation in CitcomS. The datasets are dependent on the plate reconstruction, hence separate datasets have to be constructed for each plate reconstruction.