Satellite Gravity Data and Lithospheric Structure: A study of the North Atlantic lithosphere beneath Iceland

Satellite Gravity Data and Lithospheric Structure

A study of the North Atlantic lithosphere beneath Iceland

Michael Corr

Thesis submitted for the degree of Master of Science in Geophysics

(Geodynamics) 60 credits

Department of Geoscience

Faculty of Mathematics and Natural Sciences

UNIVERSITY OF OSLO

Satellite Gravity Data and Lithospheric Structure

A study of the North Atlantic lithosphere beneath Iceland

Michael Corr

© 2019 Michael Corr

Satellite Gravity Data and Lithospheric Structure

A study of the North Atlantic lithosphere beneath Iceland

Supervisors:

Dr. Alexander Minakov, Dr. Sergei Medvedev & Prof. Carmen Gaina

This work is published digitally through DUO http://www.duo.uio.no/

Printed: Reprosentralen, University of Oslo

Abstract

Potential field analysis is known to suffer from the problem of

”non-uniqueness” in solutions. Gravity data interpretation is a prime example. Anomalous bodies can be detected from the gravity field, but information regarding size, depth, and geometry is not easily discerned.

Recently, gravity anomaly studies are augmented by satellite gravity gradient observations, which have the abilities to increase apparent sensitivity of gravity models. This study focuses on the theory of gravity and gravity gradients to implement them for the interpretation of anomalous bodies and reduce uncertainty by linking gravity studies with seismic models. The first part of the study involves developing techniques and testing them with known synthetic examples of a buried horizontal cylinder and a buried solid sphere of anomalous density. The gravity anomalies and gravity gradients are calculated and interpreted, especially in relation to the applicability and influence of edge effects. The second part is to extend this technique to a natural example, the North Atlantic centered by Iceland and the corresponding strong gravity anomaly. A density distribution of the lithosphere and upper mantle is modelled based on the S-wave velocity tomography model SL2013sv (Schaeffer and Lebedev, 2013) and 1D reference density model AK135 (Kennett et al., 1995) using a simple relation between velocity and density (Karato, 1993). Using this density distribution of the mantle, gravity anomalies and their gradients are calculated. The results are compared to observed gravity anomaly models and gravity gradients measured and calculated by the European Space Agency’s ”Gravity field and steady-state Ocean Circulation Explorer”

(GOCE) satellite, respectively. Inferences are made about the lithosphere and upper mantle structure and the benefits of gravity gradients analysis are discussed.

Acknowledgements

I would like to start with thanks to my principal supervisor Dr. Alexander Minakov, at The Centre for Earth Evolution and Dynamics (CEED), for taking me on as a master student. His initial guidance helped me to begin using Matlab for the purpose of numerical modelling and analysis of gravity data and gravity gradients. Conversations with him proved to be interesting and helped direct my studies of the North Atlantic and inspire me to delve deeper into the geodynamics of the lithosphere.

I would like to thank my co-supervisor Carmen Gaina, at CEED, for our conversations about the project and her guidance that helped me greatly.

Special thanks go toward my co-supervisor Sergei Medvedev, at CEED.

His assistance in implementing my Matlab routines in the study area, as well as his comments on my text were very helpful, especially toward the end of my thesis project.

I would like to thank my parents for their love and support as I moved a great distance to pursue this opportunity. I am grateful to be able to talk with them regularly.

Thanks to my friends Petra and Rebecca that I met through the geophysics masters program, and as fellow colleagues at CEED. We have shared many great times and had great conversations. I am thankful for their friendship and moral support over the last two years.

Lastly, I would like to thank Kristina. I am grateful for the time I have spent with her. Her energy and enthusiasm about so many things has inspired me. I am truly lucky to have met her.

ii

Preface

This master’s thesis presents gravity theory and applies it to synthetic models and to a study area located in the North Atlantic centered on Iceland. A density distribution for the lithosphere and upper mantle is calculated using an S-wave velocity tomography model. Figures of gravity field anomaly components, gravitational gradient tensor components, observed gravity field, and observed gravity gradients are presented.

This study has been completed at the Centre of Earth Evolution and Dynamics (CEED), Department of Geosciences, University of Oslo under the supervision of Dr. Alexander Minakov, Dr. Sergei Medvedev and Prof.

Carmen Gaina.

iv

Table of Contents

Abstract i

Acknowledgements ii

Preface iii

List of Figures vii

1 Introduction 1

2 Basic Theory and Methods 4

2.1 Gravitational forces between two masses . . . 4

2.2 Gravitational potential of a 3D body . . . 6

2.3 Gravitational attraction . . . 6

2.4 2D simplifications of gravity calculations . . . 8

2.5 Gravitational gradients . . . 9

3 Synthetic Experiments 13 3.1 Gravity anomaly . . . 13

3.1.1 Buried horizontal cylinder . . . 13

3.1.2 Buried solid sphere . . . 16

3.2 Gravitational gradients . . . 19

3.2.1 Buried horizontal cylinder . . . 19

3.2.2 Buried solid sphere . . . 23

4 Study Area - North Atlantic & Iceland 26 4.1 Data and models . . . 28

4.1.1 Gravity anomaly . . . 28

4.1.2 Gravitational gradients . . . 30

4.1.3 Seismic tomography and density of the mantle . . . 31

4.1.4 Spherical, Cartesian & LNOF coordinate systems . . . 33

4.2 Results . . . 38

4.2.1 Workflow for study area . . . 38

4.2.2 Gravity anomaly . . . 39

4.2.3 Gravitational gradients . . . 43

4.3 Discussion . . . 50

5 Summary and Conclusions 57

References 59

Appendix A: Additional figures 64

Appendix B: Matlab routines 74

Appendix C: Sub-lithosphere upper mantle results 88

vi

List of Figures

2.1 Gravitational attraction between two masses . . . 4

2.2 Gravitational attraction observed at a point in space outside of a 3D anomalous body . . . 7

2.3 3D body approximated by dividing intoN prisms . . . 8

3.1 2D vertical gravity anomaly of horizontal cylinder with varied observation height . . . 14

3.2 3D vertical gravity anomaly of horizontal cylinder . . . 16

3.3 3D vertical gravity anomaly of buried sphere . . . 18

3.4 2D gradient T_xx with varying height of observation . . . 20

3.5 2D gradient T_xz with varying height of observation . . . 20

3.6 2D gradient T_zz with varying height of observation . . . 21

3.7 3D gradients horizontal cylinder . . . 22

3.8 3D gradients buried sphere . . . 24

3.8 3D gradients buried sphere . . . 25

4.1 North Atlantic model grid space . . . 26

4.2 Study area observation grid space . . . 27

4.3 WGM2012 free-air gravity anomaly - North Atlantic . . . 29

4.4 Free-air gravity anomaly calculated at satellite height 225km - North Atlantic . . . 29

4.5 Density Anomaly N-S Cross Section along 20^◦W . . . 32

4.6 S-wave velocity tomography at depth 149km . . . 33

4.7 Spherical coordinate basis - Matlab . . . 35

4.8 Test plot of gravity vector anomalies after spherical transformation . . . 36

4.9 Test plot of GGT after spherical transformation . . . 37

4.10 Vertical gravity anomaly versus free-air gravity anomaly . . . 41

4.11 Gravity anomaly vector components x & y . . . 42

4.12 ModelledT_xx versus GOCE T_{W W} . . . 44

4.13 ModelledT_xy versus GOCE T_{N W} . . . 45

4.14 ModelledTxz versus GOCETW U . . . 46

4.15 ModelledT_yy versus GOCE T_{N N} . . . 47

4.16 ModelledT_yz versus GOCE T_{N U} . . . 48

4.17 ModelledT_zz versus GOCE T_{U U} . . . 49

4.18 Laplace’s equation applied to calculated gradient tensor. . . . 50

4.19 Tomography slices for mantle lithosphere - North Atlantic . . 52

4.20 Tomography slices for asthenosphere - North Atlantic . . . 53

4.21 Gravitational gradientsT_xz &T_yz . . . 55

4.22 Gravitational gradientsT_zz &T_xx . . . 56

A.1 3D gravity anomaly ∆gx, ∆gy for buried cylinder . . . 65

A.2 3D gravity anomaly ∆g_x, ∆g_y for buried sphere . . . 65

A.3 2D gradients CylinderT_xx, T_xz, T_zz . . . 66

A.4 3D gradients CylinderT_xy & T_yy . . . 67

A.5 GOCE - the 6 gravity gradient tensor components . . . 68

A.6 Tomography model depth slices 20km & 36km . . . 69

A.7 Tomography model depth slices 47km & 63km . . . 70

A.8 Tomography model depth slices 122km & 149km . . . 71

A.9 Tomography model depth slices 171km & 241km . . . 72

A.10 Tomography model depth slices 300km & 402km . . . 73

C.1 Gravity anomaly vector components, upper mantle . . . 89

C.2 Gravity gradient tensor, T, upper mantle . . . 90

viii

Chapter 1 Introduction

Iceland and the hot spot associated with its formation are some of the most spectacular and interesting examples on Earth relating to geology and geophysics. What makes Iceland interesting is that no other ridge-centered or intraplate volcanics have similar large-scale characteristics as those seen in this area.

There has been a long debate about the origin and history of the Iceland Hotspot. Is it sourced from deep within the lower mantle and connected to a Large Low Shear Wave Velocity Province (LLSVP)(e.g.

Burke et al., 2008; Burke and Torsvik, 2004; Torsvik et al., 2006) or is it sourced from within the upper mantle (e.g. Foulger et al., 2001; Foulger, 2005)? Many studies have tried to show the former to be the case, but most studies have failed to image a plume conduit that extends deep into the mantle (due to restrictions relating to seismic tomography resolution).

Observations from available data have been able to image the hotspot well within the crust and upper mantle with variability of interpretations below a depth of 660 km (e.g. Foulger and Anderson, 2005; Rickers et al., 2013;

Yuan and Romanowicz, 2017). To resolve this long-term enigma, the wealth of new data and models are delivered with good frequency (e.g. Schaeffer and Lebedev, 2013). There is a need for not only deep Earth interior study, but also better understanding and quantitative representation of shallower levels such as the lithosphere and upper mantle. A promising approach is to combine several geophysical methods to achieve this. The studies involved, however, may be of different origin and accuracy. For example, there has been a gap in wavelengths of data between land and airborne surveys (short wavelength) relative to satellite height surveys (medium to long wavelength). Satellite graviometry, however, is especially sensitive to the short wavelength part of the geopotential (Drinkwater et al., 2008). The

modern approach to collect and use geological and geophysical data provides a new view on studies of the Earth’s interior, with years of collecting, sorting, and processing of data leading to models that are more precise and useful.

Gravity data is limited in its ability to accurately resolve the size, depth and dimensions of potential anomalies due to non-uniqueness in potential field analysis (e.g. Dorman and Lewis, 1972; Oldenburg, 1974; Sebera et al., 2018). The use of gravity gradient anomalies has become increasingly common to deal with this challenge. The gravity gradient is a second rank tensor based on gravitational potential providing a measurement of the variations of the gravity vector components in space making them more sensitive to structure and directional properties of attracting sources (Panet et al., 2014). The different gravity gradient tensor components have the ability to detect edges, depths and dimensions of anomalies (Dubey and Tiwari, 2016). Gravity gradients have increased sensitivity to lateral source distribution making them more useful than the standard vertical gravity field normally used in interpretation (Ebbing et al., 2018). Several studies of the North Atlantic and the Iceland Hotspot have been conducted, but the newly acquired gravity gradients from the European Space Agency’s (ESA) “Gravity field and steady-state Ocean Circulation Explorer”

(GOCE) satellite offer new information. This provides the opportunity to produce and test models of the lithosphere structure that have not been available before. The hope is to provide usable models that can lead to greater understanding of the lithosphere and upper mantle structure and to allow better understanding of mantle dynamics and potentially the source of the Iceland hotspot. Some applications of GOCE gravity data include estimation of crustal thickness, modeling lithospheric structure, as well as investigation of deeper mantle sources (Ebbing et al., 2018).

To further aid in the interpretation, additional sources of information such as seismic tomography can be used to help constrain models. Seismic tomography is a model based on various seismic waves showing the variation of velocities within the Earth. There are a variety of models that will be compared with gravity gradients calculated from a density

2

distribution derived from the high resolution Global S-wave tomography model, SL2013sv, of Schaeffer and Lebedev (2013). SL2013sv is an isotropic global model of the shear wave velocity extending from the crust/upper mantle to the base of the mantle transition zone (Schaeffer and Lebedev, 2013). Making use of the SL2013sv model, a density distribution can be obtained for the upper mantle.

The list of geophysical and geological studies that can contribute to complex studies of the upper mantle and lithosphere may continue with, for example, magnetic field analysis or investigations of post-glacial rebound, but in my study I limit myself to the two above mentioned methods. The goal is to analyze calculated and observed gravity gradient data in the North Atlantic region centered on Iceland, to better understand the origin of the observed gravity signal and its link to the seismic low-speed anomaly in the upper mantle beneath Iceland.

Chapter 2

Basic Theory and Methods

2.1 Gravitational forces between two masses

The basis of gravimetric measurements is Newton’s law of gravitational attraction. It states: ”The magnitude of the gravitational force between two masses is proportional to each mass and inversely proportional to the square of their separation.” Figure 2.1 illustrates this concept.

Figure 2.1: Illustration of the gravitational force between two masses m and mo separated by distance r. Convention defines the unit vector ˆr being directed from the gravitational source at point Q toward the observation point P (Blakely, 1996).

The force is measured between a mass centered at a point Q= (x⁰, y⁰, z⁰) and a particle of mass at a point P = (x, y, z) (Blakely, 1996). In the Cartesian coordinate system this gives the equation:

F~ = Gmm_o

r² r,ˆ (2.1)

wherer is the separation between the two points:

4

r=p

(x−x⁰)²+ (y−y⁰)²+ (z−z⁰)²

and G is the Gravitational Constant 6.673×10⁻¹¹m³kg⁻²s⁻² and ˆr is a unit vector from Qto P.

ˆ

r = ¹_r[(x−x⁰)ˆi+ (y−y⁰)ˆj+ (z−z⁰)ˆk]

If massmo is treated as a particle with a unit magnitude, then the force F~ divided by mass m_o yields the gravitational attraction produced by mass m on point P. This gravitational attraction ~g(P) is a vector measured in the direction that points from the mass m to observation pointP.

~g(P) =−Gm

r² r,ˆ (2.2)

Since ~g is a force divided by mass the units are acceleration, it is commonly referred to as the gravitational acceleration. Gravitational acceleration is a conservative field meaning it can be represented as the gradient of a scalar potential:

~g(P) = ∇U(P), (2.3)

where

U(P) = Gm

r (2.4)

U is the gravitational potential (also called the Newtonian potential).

The gravitational potential is harmonic for all points outside of the mass when the following is satisfied:

∇²U(P) = ∂²U

∂x² +∂²U

∂y² + ∂²U

∂z² = 0 (2.5)

This relation is the Laplace’s equation, which is a special case of Poisson’s equation valid only in space not occupied by mass (Blakely, 1996).

2.2 Gravitational potential of a 3D body

Gravitational potential follows the “principle of superposition”. This states: “The gravitational potential of a collection of masses is the sum of the gravitational attractions of the individual masses”. Meaning that the net force observed at a point in space is simply the vector sum of the forces due to all masses in space (Blakely, 1996). To apply this principle to find gravitational attraction due to a continuous mass distribution, the mass is divided into many, very small elements of mass dm = ρ(x⁰, y⁰, z⁰)dV. ρ(x⁰, y⁰, z⁰) is the density distribution (similar to the polygon method used in: Talwani and Ewing, 1960; Talwani et al., 1959). This gives the equation:

U(P) = G Z

V

dm

r (2.6)

=G Z

V

ρ(Q)

r dV, (2.7)

V is the volume occupied by the mass that the integration will be calculated over (dV = dx’dy’dz’),P is the observation point, Qis the integration point, r is the distance between P and Q.

2.3 Gravitational attraction

For the purpose of this study, I consider observation points located outside of a mass distribution (Figure 2.2). It is also natural to assume that the density is a continuous and limited function and the integral will converge for all points outside of the mass. Furthermore, the differentiation with respect to x, y and z can be moved inside the integral.

6

Figure 2.2: An illustration of a three-dimensional body with a density ρ(x⁰, y⁰, z⁰) and the direction of the gravitational attraction observed at point P(x, y, z) defined in a Cartesian coordinate system (Blakely, 1996).

Using Equation 2.3 in the Cartesian coordinate system, then under the assumption of interchangeable differentiation and integration, components of the vector attraction can be expressed (Sebera et al., 2018):





 g_x g_y g_z





=−G Z

x⁰

Z

y⁰

Z

z⁰

ρ(x⁰, y⁰, z⁰) r³







x−x⁰ y−y⁰ z−z⁰





dx⁰dy⁰dz⁰, (2.8) wherer is defined as in Equation 2.1.

A strategy to deal with complex geologic structures is to divide the hypothetical gravitational sources into N simpler parts (Blakely, 1996).

This concept is illustrated in Figure 2.3. Applying this concept to Equation 2.8 yields the following:





 g_x g_y gz





=−G

N

X

i=1

ρ_i r_i³







x−x⁰_i y−y_i⁰ z−z_i⁰





dV_i, (2.9)

wherem is the component of the gravity vector at the observation point, ρi is the density in the model space at point Q(x⁰, y⁰, z⁰). dVi = dx⁰_idy⁰_idz_i⁰ is the volume of the element being integrated. With N tending to infinity, Equation 2.9 naturally tends to Equation 2.8.

Figure 2.3: Approximation of a 3D body by dividing it into N simpler volumes (Blakely, 1996).

2.4 2D simplifications of gravity calculations

The above sections have discussed the gravitational potential of a 3D body. Some problems can be simplified to 2D to ease conceptualization, as well as modelling. A good example of a 2D model is the buried horizontal cylinder in which the properties along one direction (e.g. the y axis) do not change or changes are insignificant. The gravitational potential of a 2D body with density ρ(x⁰, z⁰) is calculated using (from Blakely, 1996):

U =−2G Z

S

ρ(S)log1

rdS, (2.10)

8

where integration is over the cross-sectional surface S and similarly to 3D, r is the distance between observation point and an element of the body. This distance is defined similar to Equation 2.1, but eliminates the y-direction.

This gives the gravitational acceleration at point P is:

g_x = ∂U

∂x =−2G Z Z

S

ρ(S)x−x⁰

r² dS (2.11)

g_z = ∂U

∂z =−2G Z Z

S

ρ(S)z−z⁰

r² dS (2.12)

dS = dx⁰dz⁰ is the cross-sectional surface area of the piece being integration.

2.5 Gravitational gradients

The second-order derivative of the gravitational potential U is the gravitational gradient tensor(GGT), T. For 3D gravitational gradients, the tensor has nine components, it is defined by:

T=







T_xx T_xy T_xz T_yx T_yy T_yz T_zx T_zy T_zz





 (2.13)

The formulas for each tensor component in Cartesian coordinates (modified from Sebera et al., 2018; Zhang et al., 2000) are:

T_xx =−G

N

X

i=1

∆ρi

r_i⁵ (r²_i −3(x−x⁰_i)²)dV_i (2.14) T_yy =−G

N

X

i=1

∆ρ_i

r⁵_i (r²_i −3(y−y_i⁰)²)dV_i (2.15) T_zz =−G

N

X

i=1

∆ρi

r⁵_i (r²_i −3(z−z_i⁰)²)dV_i (2.16)

T_xy =−G

N

X

i=1

∆ρ_i

r_i⁵ (−3(x−x⁰_i)(y−y⁰_i))dV_i (2.17) Txz =−G

N

X

i=1

∆ρ_i

r_i⁵ (−3(x−x⁰_i)(z−z_i⁰))dVi (2.18) T_yz=−G

N

X

i=1

∆ρ_i

r_i⁵ (−3(y−y_i⁰)(z−z_i⁰))dV_i (2.19) where r_i =p

(x−x⁰_i)²+ (y−y⁰_i)² + (z−z_i⁰)²

Equations 2.14 - 2.19 show that Tij = Tji and the gradient tensor is symmetric, thus only 6 terms of the tensor are independent.

In contrast with the above sections, gravitational gradients can be related to density variations (∆ρ=ρ_o−ρ), as an infinite space with constant density (ρ_o) results in values of zero for the graviy gradients.

For 2D gravitational gradients, the tensor is defined by:

T= T_xx T_xz T_zx T_zz

!

(2.20) I derived the T_xx and T_xz tensors by differentiating Equation 2.11 with respect to x and z (respectively). T_zz is derived by differentiating Equation 2.12 with respect to z.

T_xx =−G

N

X

i=1

∆ρi

r⁴_i (r²_i −2(x−x⁰_i)²)dS_i (2.21) T_zz =−G

N

X

i=1

∆ρ_i

r⁴_i (r_i²−2(z−z⁰_i)²)dS_i (2.22) T_xz =−G

N

X

i=1

∆ρi

r⁴_i (−2(x−x⁰_i)(z−z_i⁰))dS_i (2.23) where ri =p

(x−x⁰_i)²+ (z−z_i⁰)² and dS =dx⁰dz⁰ is the cross-sectional surface area of the piece i(i= 1, ..., N). By matrix symmetryT_xz =T_zx and is not presented.

10

Traditionally, only the vertical (radial) component of the anomalous gravity vector (simply referred to as ”the gravity anomaly”) has been utilized in geodetic and geophysical modelling studies. Conventional gravity data (such as gravity anomalies) shows the strength of the Earth’s gravity field, but is less sensitive to the edges of bodies and contains no directional information (Dubey and Tiwari, 2016). Gravity gradients directly recover sharp signal over the edges of structures and provide a higher potential to uncover information about the center of mass and geometry of the subsurface structure (Butler, 1995; Dubey and Tiwari, 2016). Butler (1995) also states that the gradients of the gravity field provide greater spatial resolution, better definition of structural/geometric boundaries, and an improvement in depth determination when compared to the measured gravity field.

Studies relating to crustal thickness estimates, lithosphere structure modeling (e.g. Fullea et al., 2015), and investigation of deep mantle sources (e.g. Panet et al., 2014) have made use of GOCE gravity gradient data.

They have demonstrated that gravitational gradients provide additional sensitivity to the lateral sources distribution in comparison to conventional, vertical gravity field (e.g. Ebbing et al., 2018).

Gravity gradients tensors (eq. 2.13) are more sensitive to shallower geologic structures (Zhang et al., 2000). The usefulness of gravitational gradients is described by Dubey and Tiwari (2016), they explain that the vertical tensor component, T_zz, provides an estimate of maximum depth and can predict boundary information directly related to the causative body. The other components of the tensor can provide information about the geometry of the body. T_xx can indicate the eastern and western edges of a feature. T_yy indicates the northern and southern edges. T_xz divides the body into approximately symmetric eastern and western halves that highlight the central axis of the body in the north-south direction.

Conversely, Tyz divides the body into northern and southern halves and define the central axis of the body in the east-west direction. The T_xy can be used to detect the corners of nearly rectangular bodies, as well as identify the centerpoint of symmetrical bodies when they are aligned in the

x and y directions.

Applying Laplace’s equation (eq. 2.5) expressed in terms of the gravity gradient tensor yields:

For 2D:

T_xx+T_zz = 0 (2.24)

For 3D:

T_xx+T_yy+T_zz = 0 (2.25)

Equations 2.24 and 2.25 can provide a first order test of accuracy of gravity gradient data and the numerical calculations of the tensor.

12

Chapter 3

Synthetic Experiments

Before applying the theory and methods discussed to a natural example, it is important to make synthetic models to ensure the calculations are done correctly. In this section I developed a set of Matlab routines to model the effects of a buried horizontal cylinder extending to infinity (or to the calculation domain edge) in one direction, as well as a buried sphere characterized by density excess. Both have the same density and radius. I vary height of observation to see the effects on the gravity anomaly and gravitational gradients calculated. The observation space for both the 2D and 3D setup have the same lateral dimensions to illustrate potential edge effects. Presented are the analytical and numerical solutions of gravity anomalies for both 2D and 3D cases and the numerical solutions of the gravity gradients.

3.1 Gravity anomaly

3.1.1 Buried horizontal cylinder

Turcotte and Schubert (2002) define the analytical solution for the vertical gravity anomaly of a buried horizontal infinite cylinder as:

∆g_z = 2πG∆ρa² d+h

x²+ (d+h)² (3.1)

Where a is the radius of the cylinder, d is the depth of the cylinder’s central axis, xis the observation point along the x-axis,h is the elevation of the observation point, and the cylinder is centered at x= 0.

With the third dimension, y, in which the cylinder extends to infinity being eliminated, the density contrast is only a function of x and z. It can

be seen that by increasing the model resolution, the difference between the analytical and numerical solutions is decreased showing the correctness of the numerical procedures implemented (see Figure 3.1).

Figure 3.1: Two dimensional example of a buried horizontal cylinder with

∆ρ = 600kg/m³, depth 5 km & radius 2 km. The vertical gravity anomaly, ∆g_z, analytical solution (eq. 3.1) compares well with the numerical approximation (eq. 2.12). The effect of increasing observation height is displayed in green for both the analytical and numerical solutions. A reduction in amplitude is seen, as well as a broadened edge of the anomaly.

In these experiments, h is varied between 0 & 2 km.

To calculate the vertical component of the gravity anomaly at observation

14

point P, the following equation is used (modified from Zhou, 2009):

∆g_z = 2G Z Z

S

∆ρ(x⁰, z⁰)(z−z⁰)

r² dx⁰dz⁰ (3.2) which is the general form of the 2D line integral to calculate the gravity anomaly in the z direction at any point P along the x direction produced by a 2D mass with a density contrast varying in both the x and z direction (Zhou, 2009). This is similar to equation 2.12 that calculates the gravitational attraction in thezdirection, but instead calculates the gravitational anomaly by using the density contrast at the point of integration.

The equation for the vertical gravity anomaly at pointP(x, y, z) is given by:

∆g_z =G

Z Z Z

V

∆ρ(x⁰, y⁰, z⁰)(z−z⁰)

r³ dx⁰dy⁰dz⁰ (3.3) Similar to equation 2.8, ∆g_x and ∆g_y can be calculated by substituting

∆ρ=ρ_o−ρ in the equation.

When modelling the infinite horizontal cylinder in three dimensions it becomes clear that edge effects may be present. As the model includes the dimension in which the cylinder is meant to extend to infinity, the model limits introduce a sharp edge effectively cutting the cylinder. Figure 3.2 shows the elongated structure of the cylinder, but toward the model limits of the y-axis (infinite direction) an apparent end of the cylinder is seen. This properly illustrates the importance of defining an observation space that has a reasonable buffer distance to the model space. In this synthetic example, the observation space has the same x and y limits as the model space. To approach a more accurate solution a smaller observation window is required.

Figure 3.2: Example of a buried horizontal cylinder with ∆ρ = 600kg/m³, depth 5 km & radius 2 km. The calculated vertical gravity anomaly, ∆g_z, is displayed.

Comparison of calculations of gravity anomalies caused by the horizontal cylinder using Equations 3.1 - 3.3 show the robustness of the numerical routines developed and that edge effects may visibly interfere with results.

3.1.2 Buried solid sphere

Turcotte and Schubert (2002) define the analytical solution for the vertical gravity anomaly of a buried spherical mass as:

∆g_z = 4πG∆ρa³ 3

d+h

(x²+y²+ (d+h)²)^3/2 (3.4) Where a is the radius of the sphere, d is the depth of centre of sphere,x

16

is the observation point along the x-axis, yis the observation point along the y-axis, his the elevation of the observation point, and the sphere is centered at x= 0 & y= 0.

Figure 3.3 presents the comparison of the analytical solution using Equation 3.4 and numerical calculations using Equation 3.3. The results compare well. One of the main purpose of this case is to constrain gravity anomaly and gravity anomaly gradient patterns which will later be used for testing the complex coordinate transformations for the North Atlantic (Section 4.1.4).

Figure 3.3: Example of a buried sphere with ∆ρ = 600kg/m³, depth 5 km

& radius 2 km. a.) Profile crossing the center of the sphere along the x- axis and the associated vertical gravity anomaly, ∆gz, analytical solution compared with the numerical solution. Increasing model resolution leads to a better fit to the analytical solution. With a mesh resolution of 0.2km, analytical and numerical solutions are indistinguishable. b.) A plot of the calculated vertical gravity anomaly, ∆g_z over the xy observation plane. A strong positive anomaly is calculated over the positive density contrast of the

sphere. 18

3.2 Gravitational gradients

As described in Chapter 2, gravitational gradients can be effective to aid in the interpretation of anomalous bodies/density distributions. The wealth of data collected by the GOCE mission also favours the consideration of gradients in more detail. I present examples of gravitational gradients for both 2D and 3D examples. The gradients are measured in units of eotvos (1E = 10⁻⁹s⁻²). There are no analytical solutions that I found published for these cases.

3.2.1 Buried horizontal cylinder

Figure 3.4 shows the 2D gradient component T_xx of the cylinder. With the positive density constrast, there is an associated negative gradient present over the cross section of the body. The maximum is positioned above the central axis of the cylinder. Figure 3.5 displays the component T_xz for the cylinder. This gradient component offers a good indication of the edges of the cylinder, and the sign of the density anomaly in relation to the direction of observation. With the positive density contrast, moving along the observation line toward the body displays a high, a zero magnitude over the central axis, and a low over the opposite edge. Figure 3.6 shows the T_zz component of the cylinder. It shows a strong positive over the cylinder cross section indicating the central axis in a similar fashion to the T_xx component. The results presented also show the effect of altering observation height. With increased observation distance the gravitational gradients show a decrease in magnitude, and a broadening of the anomaly edges (affecting edge interpretation of the causative body).

Figure 3.4: 2D gradients T_xx calculated at observation height of 0 km and 2 km to show the reduction of amplitude with increased observation height.

Figure 3.5: 2D gradients T_xz calculated at observation height of 0 km and 2 km to show the reduction of amplitude with increased observation height. At h = 0 km, T_xz accurately detects the edges of the cylinder at the maximum and minimum gradient values. With increased height the anomaly broadens and the edges have potential to be misinterpreted indicating a cylinder with a much larger diameter.

20

Figure 3.6: 2D gradient T_zz calculated at observation height of 0 km and 2 km to show the reduction of amplitude with increased observation height.

The main test for the numerical calculations is the sum of the principal tensor components using Laplace’s equation (eq. 2.24).

The gradients for the buried horizontal cylinder in three dimensions (Figure 3.7) can be interpreted based on the method described in Section 2.3. Txx is a good indicator of the eastern and western edges of the cylinder. T_yz displays strong edge effects along the north and south edges where the cylinder cannot be effectively modelled as infinite. T_zz is very effective to define the boundaries of the cylinder. T_xz shows a well defined central axis of the cylinder in the y direction. T_yy (Appendix A1) indicating the north and south edges of the body. T_xy (Appendix A1) shows the corners of the ”ends” of the cylinder at the model limits in the North/South direction and the central axis in the x direction.

If the model space were larger than the observation window these apparent edges would be diminished. The benefit of using gravitational gradients to identify the geometry of an anomalous body is clearly demonstrated in these examples.

a.) b.)

c.) d.)

Figure 3.7: Gravitational gradients of a buried horizontal cylinder. a.) T_xz shows the anomaly divided into 2 halves along the x axis. It gives the central anomaly axis in the y direction. b.) T_yz shows the strong edge effects of the cylinder. c.) T_zz shows the geometry of the body. d.) T_xx shows a clearly defined central axis and orientation of the cylinder.

22

3.2.2 Buried solid sphere

For the buried sphere, similar analysis can be done using the gradients (Figure 3.8). T_xx is a good indicator of the eastern and western edges of the sphere. T_yy displays the north and south edges of the sphere. T_zz is very effective to define the geometry/boundaries of the sphere. T_xy effectively locates the centre of the sphere. T_xz shows a well-defined central axis of the sphere in the y direction. T_yz shows the central axis in the x direction.

T_xx+T_yy+T_zz is effectively zero satisfying Laplace’s equation (eq. 2.25).

a.) b.)

c.) d.)

Figure 3.8: Gravitational gradients of a buried sphere. a.) T_xx approximates the dimensions of the east and west edges of the sphere. b.) T_xy shows the anomaly divided into 4 quadrants centered on the centre of the sphere. c.) T_xz shows the anomaly divided into 2 halves along the x axis. Highlighting the central anomaly axis in the y direction. Strongest negative and positive values concentrated on the east and west edges of the sphere, respectively.

d.) T_yy approximates the dimensions of the north and south edges of the sphere.

24

e.) f.)

Figure 3.8: (cont.) Gravitational gradients of a buried sphere. e.) Tyz

shows the anomaly divided into 2 halves along the y axis. Highlighting the central anomaly axis in the x direction. Strongest negative and positive values concentrated on the north and south edges of the sphere, respectively.

f.) T_zz strong anomaly clearly showing the geometry of the sphere.

Reviewing the results from the synthetic examples provides some insight into the benefit of the gravity field measured over an anomalous body, as well as the usefulness of gradient data to detect the approximate geometry of the body. The observations can be roughly extended to natural examples. For example, the cylinder is similar to what can be observed from elongated geologic bodies such as slabs or plate boundaries/spreading ridges, and the sphere is similar to what can be expected for a vertical body such as a plume/hot spot. With this in mind, there is potential to extend these rudimentary models to a real world example of the North Atlantic/Iceland presented in the next chapter.

Chapter 4

Study Area - North Atlantic &

Iceland

With the set of tools developed in Chapter 3 tested using synthetic examples, they can now be extended to a natural example of the North Atlantic. The location of the study area is presented in Figure 4.1. To reduce the edge effects, the calculation domain (black mesh) is greater than the observation area used to present the results (defined by the yellow border) and displayed in Figure 4.2.

Figure 4.1: The model grid space displayed on the globe (black borders) with observation area inset (yellow border).

26

Figure 4.2: The topography and bathymetry (ETOPO1, Amante and Eakins (2009)) of the observation grid space. The study area is smaller than the model space to reduce possible edge effects.

In the following, data and models are presented first for calculating a density distribution of the lithosphere and upper mantle using a seismic tomography model, and second for calculating gravity anomalies and gravity gradient tensors from the density distribution that can be tested and compared with observed gravity data. The large-scale area considered requires accounting for sphericity of the Earth, and that leads to additional analysis of coordinate system transformations.

4.1 Data and models

4.1.1 Gravity anomaly

WGM2012

The World Gravity Map developed in 2012 (Bonvalot et al., 2012) provides a 1^◦x1^◦ grid and maps of worldwide gravity anomalies including Bouguer, isostatic and surface free-air computed in spherical geometry.

The gravity anomalies in WGM2012 are derived from previously available global gravity models (EGM2008 and DTU10). It also includes 1^◦x1^◦ resolution terrain correction derived from ETOPO1 (Amante and Eakins, 2009). They are tied to the node points of the ETOPO1 equiangular mesh and computed by means of spherical harmonics (Balmino et al., 2012).

WGM2012 presents a strong positive (up to 70 mGal, Figure 4.3) gravity anomaly associated with the Iceland hotspot, mid-ocean ridge, and thickened crust. The synthetic experiments of Chapter 3, however, show that the gravity signal broadens and the amplitude decreases with increased distance to the measurement point. The following numerical experiments are to be presented at satellite observation height, 225km. Thus, the surface gravity data (Figure 4.3) was recalculated for the desired elevation (Minakov, unpublished). Figure 4.4 shows that the amplitude of the anomalies are reduced, but the main features remain visible, such as the Mid-Atlantic Ridge (MAR).

28

Figure 4.3: WGM2012 free-air gravity anomaly over study area in the North Atlantic.

Figure 4.4: Free-air gravity anomaly over the North Atlantic from satellite observation height, 225 km (Minakov, unpublished).

4.1.2 Gravitational gradients

GOCE

This thesis will make use of the gravity gradient data collected from the European Space Agency’s (ESA) “Gravity field and steady-state Ocean Circulation Explorer” (GOCE) satellite that operated between 2009 and 2013. This mission operated at mean altitudes of 255 km (during nominal phase) and 225 km in the extended phase and is the first satellite mission to measure the full gravitational gradient tensor (Panet et al., 2014). GOCE combines GPS tracking and gravity gradiometry to determine Earth’s mean gravity field with global accuracy and spatial resolution down to 80 km.

This resolution makes the gravity gradient data products from GOCE very useful in lithosphere scale modeling (Bouman et al., 2013). The gravity field measurements have a resolution of 1 mGal (Drinkwater et al., 2008).

The gravity gradients are in local satellite reference system and the orbit data in a Earth-fixed reference frame. The six gradients that are measured by GOCE are in the Gradiometer Reference Frame (GRF) which co-rotates with the satellite (Bouman et al., 2016). The measurements were done at an observation height of 255 km during the nominal mission phase, and then at 225 km during the final months of the mission. The T_xx,T_yy,T_zz, and Txz gradients are high accuracy within the measurement bandwidth, whereas the Txy and Tyz gradients have two orders of magnitude worse accuracy (Bouman et al., 2016). The measurement bandwidth is between 40 - 750 km half-wavelength (Bouman et al., 2016). The gradients are available in the Local North Oriented Reference Frame (LNOF) that assigns tensor orientations in North, West, & Upward (NWU) (Refer to Appendix A1 where all 6 components of T are presented in LNOF).

30

4.1.3 Seismic tomography and density of the mantle

S-wave Global Tomography Model

Seismic tomography is a model based on seismic wave data showing the variation of wave velocities within the Earth. Using the assumption that the mantle is chemically homogeneous, tomography data reflects variations of temperature in the mantle that can be directly related to density perturbations (Marquart and Schmeling, 2004). Chemical heterogeneities might be responsible for some of the velocity anomalies observed, but Marquart (2006) discuss that this is less important in the mantle transition zone. Negative density anomalies due to hot, buoyant material cause a negative gravity anomaly, with the buoyancy forces related to hot upwellings driving solid state flow (Marquart et al., 2000). This leads to a deflection of the surface and produces an additional gravity effect (Marquart and Schmeling, 2004). The sum of both of these effects gives the observed gravity potential field anomaly.

Making use of the SL2013sv model (Schaeffer and Lebedev, 2013), a density distribution can be obtained for the lithosphere and upper mantle between depths of 20 and 660 km (See Appendix A3 where 10 layers of the tomography model are presented). To calculate densities for the lithosphere and upper mantle using SL2013sv, a 3D reference density grid is constructed (ρref) using the 1D reference density model AK135 (Kennett et al., 1995).

A simple relation between velocity and density is assumed:

ρ_m =ρ_ref(1 +c×dV_s), (4.1) where dVs is the S-wave velocity perturbation in percentage and c is the velocity to density conversion factor. The value of c = 0.2 is proposed for the S-wave velocities in the lithosphere region (Karato, 1993). The value of 0.2 is also used by Mitrovica and Forte (1997) where they find that different values have a relatively small impact on viscosity profiles, but for a value of 0.4 there is a slight reduction in the viscosity profile between depths of 600

and 1000 km. Since this is primarily below the base of the model (660 km), the value of 0.2 is acceptable. Additionally, based on the thermal origin of wave speed anomalies, 0.2 is a midway value suggested based on previous lithosphere studies (e.g. Simmons et al., 2009; Steinberger and Calderwood, 2006; Warners-Ruckstuhl et al., 2012). The result of this calculation (eq.

4.1) is the modeled density (ρ_m).

For use in the gravity anomaly and gravity gradient tensor calculations outlined in the methods section, ∆ρ is required.

∆ρ=ρ_m−ρ_ref (4.2)

This results in a density anomaly varying between +/−50kg/m³

(Figure 4.5 illustrates a north-south density profile along 20^◦W calculated from the tomography model). Figure 4.6 shows a depth slice of 149 km from the tomography model SL2013sv (Schaeffer and Lebedev, 2013), it shows a strong low velocity anomaly (-5% S-wave perturbation) present beneath Iceland and the Mid-Atlantic Ridge between 55^◦ and 75^◦ N.

Figure 4.5: Density cross section along 20^◦W calculated from SL2013sv (Schaeffer and Lebedev, 2013). Values are the difference between modelled density and the reference density model AK135 (Kennett et al., 1995).

32

Figure 4.6: Seismic tomography at a depth slice of 149 km. Values are the S- wave velocity perturbations (%) (SL2013sv, Schaeffer and Lebedev, 2013). A strong low velocity anomaly is present beneath Iceland and the Mid-Atlantic Ridge between 55^◦ and 75^◦ N.

4.1.4 Spherical, Cartesian & LNOF coordinate systems

The use of a geographical grid in discretization of the calculation domain and general curvature of Earth bring additional complications to the calculation first presented in Chapter 3. First, due to the non-orthogonality of the mesh, integration requires additional constraints.

In calculations of element volumes, there is a big difference between a volume defined by the orthogonal regular mesh and in the non-uniform mesh of this study. As you move along a great circle in a geographic reference frame toward the poles, you get distortion of the longitudinal distances. This must be accounted for in the element volume calculations.

The volume of a rectangular prism:

dV =dx⁰dy⁰dz⁰ (4.3)

becomes the volume of a spherical prism (also referred to as a tesseroid) (from Liang et al., 2014) that is a function of radial position and latitude of the element:

dV =R²×cos(Latitude)×dR×dLat×dLon (4.4) where R is the radial position of the element, dR is the radial dimension of the element, dLat is the latitude dimension of the element, and dLon is the longitude dimension of the element. In contrast to a potentially constant orthogonal mesh case (eq. 4.3), dV is strongly dependent on the position of the elemental volume, radius and latitude.

Second, all of the equations in Chapter 2 are for the Cartesian coordinate system, whereas the data and testing models are presented in coordinates linked to the Earth. The input data for the calculations are available in spherical coordinates, whereas the calculation procedures (Chapter 2) require Cartesian coordinates. Calculated gravity vectors and gravity gradient tensors are to be considered in spherical coordinate system.

For comparison with observed satellite data such as GOCE gravity gradients, a conversion must be done between Cartesian and spherical coordinates. As the calculations have been done in a using a Cartesian reference frame the gravity vector components and gravitational gradient tensor components must be converted.

This requires a set of routines (See Appendix B) for converting coordinates, vectors, and tensors between Cartesian and spherical using the Mathworks format of spherical representation of vectors (Figure 4.7).

34

Figure 4.7: Spherical basis vectors defined as a local set of vectors pointing along the radial and angular directions at any point in space (from MathWorks documentation).

To test the correctness of this transformation, the sphere example that was performed in Chapter 3 is performed in a spherical reference frame.

Figures 4.8 and 4.9 show a good agreement and demonstrate the robustness of the routines.

Figure 4.8: Gravity vector test after transformation of ∆g_x, ∆g_y, & ∆g_z to spherical using a test positive density spherical anomaly centered in the North Atlantic. Confirms equations are correct by comparison to synthetic experiment for the buried sphere conducted in Chapter 3.

36

Figure 4.9: Gravity gradient tensor test after transformation of the gradient tensor,T, to spherical using test positive mass anomaly centered in the North Atlantic. This confirms the equations are correct by comparison to synthetic experiment for the buried sphere conducted in Chapter 3.

An additional step must be considered for data comparison in the next section. As the GOCE data is presented in the LNOF, comparison of modelled results with GOCE requires consideration of the gradient tensor components presented in LNOF. The directions defined in LNOF are marked as NWU. N being positive in the north direction, W being positive

in the westward direction, and U being positive up in the radial direction (Bouman et al., 2013). LNOF defines x as the north direction N, y as positive west direction W, and z as the radial up direction U. In the following results section, I refer to the spherical coordinate system xyz with z being positive downward directed to the Earth centre, x is directed positive eastward along parallels, and y is directed positive northward along the great circles. This is consistent with the synthetic experiments detailed in Chapter 3.

4.2 Results

4.2.1 Workflow for study area

The following workflow details the steps taken to calculate the gravity anomalies and gravity gradient tensors for the study area using Matlab routines developed (Appendix B):

1. Define conversion factors and constants

2. Model Parameters (geographical limits, grid resolution) 3. Define observation height and grid extent

4. Define Regular Mesh and Observation Grid 5. Define Mantle density (20 - 660 km)

Load Tomography Data (SL2013sv) Load 1D reference density model (AK135) Convert velocity to density (eq. 4.1)

6. Calculate gravity anomalies (eq. 2.9 using ∆ρ in place of ρ)

7. Calculate gravity gradient tensor (eqs. 2.14, 2.15, 2.16, 2.17, 2.18, &

2.19)

8. Conversion of gravity vectors from Cartesian to spherical coordinates (using Matlab routine B.14 from Appendix B)

9. Conversion of gravity gradient tensors from Cartesian to spherical coordinates (using Matlab routine B.15 from Appendix B)

10. Validate gradient calculations using Laplace’s equation (eq. 2.25) 38

The calculations were performed using the depth interval of SL2013sv (Schaeffer and Lebedev, 2013) between 20 and 660 km. The dimensions of 120x100x80 grid points (resolution of approximately 5 km x 0.5^◦ of latitude x 1^◦ of longitude) and an observation grid with dimensions of 40x41 (resolution of approximately 0.5^◦ of latitude x 1^◦ of longitude) were used in the study area.

It should be noted when comparing the calculated results with observed data that the scales of amplitudes are of different minimum and maximum values. The difference in amplitudes could be attributed to many factors, with shallow sources, such as the crust, likely influencing the observed data.

Other studies have noted the amplitude discrepancy between the calculated and observed field (e.g. Ebbing et al., 2014). Ebbing et al. (2014) found gravity gradients calculated to be almost twice the amplitude of the observed gradients. I do not have crust accounted for in my models, so this could possibly be the reason amplitudes of nearly three to four times that of the observed gradients are calculated. The same can be said for the gravity anomaly, values are nearly three times that which is observed in the free-air gravity anomaly from satellite observation height.

4.2.2 Gravity anomaly

The calculated gravity field anomalies in the model are converted from Cartesian to spherical coordinates. In spherical coordinates they can be compared to the observed gravity anomalies and calculated gravity anomalies.

There is a strong negative (approx. -100mGal) anomaly present in the vertical gravity anomaly (Figure 4.10) centered on Iceland. This is due to the low velocity zone associated with mantle upwelling and the Iceland hotspot transporting hot, less dense material to the surface. The free-air gravity anomaly shows a strong positive (approx. 40mGal) anomaly speculated to be due to a crustal influence in the upper 20 km not modelled.

Figure 4.11 shows a strong negative, linear anomaly to the west of the Mid-Atlantic Ridge (MAR) and a positive east of Iceland in the plot of

∆g_x. The gravity anomaly component ∆g_y shows a negative anomaly south of Iceland and a positive northwest toward Greenland. Both of these results are consistent with the observations made in the synthetic experiments conducted in Chapter 3. The MAR can be approximated by a horizontal cylinder, and deep sources like the Iceland hotspot by a sphere.

40

Figure 4.10: Calculated vertical gravity anomaly, ∆g_z, versus free-air gravity anomaly from satellite observation height, 225 km (Minakov, unpublished).

Figure 4.11: Gravity anomalies ∆gx and ∆gy. ∆gx shows a strong east/west transition crossing the Mid-Atlantic Ridge. There is a strong, linear negative anomaly (-100mGal) west of the ridge that extends roughly parallel to the ridge axis, on the east side of the axis the anomaly is positive. This is expected when the observation point moves eastward over a negative density anomaly. ∆g_y shows a weaker (-40mGal) anomaly on the southeast edge of Iceland. With the orientation of the Mid-Atlantic Ridge being roughly North-South, the y component does not easily detect the trend.

42

4.2.3 Gravitational gradients

The calculated gravitational gradient tensor in the model is converted from Cartesian to spherical coordinates. In spherical coordinates it can be compared to the observed gradients measured and calculated by GOCE. The LNOF convention used in GOCE data assigns x to the North direction (N), y to the West direction (W), and z to the radial direction up (U) (Bouman et al., 2013).

T_xx shows a positive anomaly that runs along the MAR corresponding to a negative density anomaly (Figure 4.12). The GOCE T_{W W} shows a similar trend, but is weaker and the opposite sign due to spectulated crustal influences not accounted for in the model. Txy shows a faint, negative linear anomaly beneath Iceland and north toward Jan Mayen that trends along the MAR (Figure 4.13). The corresponding GOCE T_{N W} component shows similar trends, but is weaker by comparison. T_xz shows a strong positive (3-4E) gradient northeast and south of Iceland, and a strong negative (-4E) linear anomaly running parallel west of the MAR (Figure 4.14). GOCET_{W U} show similar trends, but with lower amplitude. T_yy shows positive gradient beneath Iceland and to the southwest along the MAR. A similar trend is seen in GOCE T_{N N} with a negative gradient centered on Iceland. T_yz shows a strong negative (-3E) gradient on the southeastern edge of Iceland close to the approximate hotspot location (Shorttle et al., 2010). The same negative can be seen in GOCET_{N U}, but slightly offset to the west to be situated on the southern tip of Iceland. T_zz shows much better detail of the strong negative gravity anomaly seen in ∆g_z. The negative anomalies exhibit a similar areal extent as deep negative velocity/density anomalies present between depths of 122-241km of the SL2013sv tomography model (Schaeffer and Lebedev, 2013). The corresponding GOCE T_{U U} component shows a positive anomaly following the ridge axis and seems to likely be affected by shallow crustal material.

Figure 4.12: ModelledT_xx versus GOCET_{W W}. T_xx shows a positive anomaly that runs along the MAR, this corresponds to a negative density anomaly.

The GOCE T_{W W} shows a similar trend, but is weaker and the opposite sign due to possible crustal influences not accounted for in the model.

44

Figure 4.13: Modelled T_xy versus GOCE T_{N W}. T_xy shows a faint, negative linear anomaly beneath Iceland and north toward Jan Mayen that roughly follows the trend of the MAR. The corresponding GOCE T_{N W} component shows similar trends, but weaker in comparison.

Figure 4.14: Modelled T_xz versus GOCE T_{W U}. T_xz shows a strong positive (3-4E) gradient northeast and south of Iceland, and a strong negative (-4E) linear anomaly running parallel west of the MAR. GOCE T_{W U} show similar trends, but fainter.

46

Figure 4.15: Modelled T_yy versus GOCE T_{N N}. T_yy shows positive gradient beneath Iceland and to the southwest along the MAR. A similar trend is seen in GOCE T_{N N} with a negative centered on Iceland.

Figure 4.16: Modelled T_yz versus GOCE T_{N U}. T_yz shows a strong negative (-3E) gradient on the southeastern edge of Iceland near the approximate hotspot location (Shorttle et al., 2010). The same negative can be seen in GOCET_{N U}, but slightly offset to the west to be situated on the southern tip of Iceland.

48

Figure 4.17: Modelled T_zz versus GOCE T_{U U}. T_zz shows much better detail of the strong negative gravity anomaly seen in ∆g_z. The negative anomalies present are similar in areal extent as deep negative velocity/density anomalies seen between depths of 122-241km seen in the SL2013sv tomography model. The corresponding GOCET_{U U} component shows a positive anomaly following the ridge axis and seems to be mostly affected by crustal material.

Using Laplace’s equation (Equation 2.25) the solutions of the gradient tensor can be validated as seen in Figure 4.18. The sum is virtually zero, so it shows that the calculations of the tensor are correct (also verified in the synthetic experiments).

Figure 4.18: Laplace’s equation applied to the study area. From equation 2.25, T_xx +T_yy +T_zz = 0. The sum is effectively zero, so the equation is satisfied.

4.3 Discussion

Combining datasets from different sources can help to constrain interpretation and highlight strengths of particular methods. In Section 4.1, Figure 4.3 shows the Free-air gravity anomaly in the North Atlantic from ground and airborne surveys, whereas Figure 4.4 shows the Free-air gravity anomaly from satellite observation height (225 km). At satellite height there is less detail as anomalies are smoothed and reduced in amplitude with increased observation distance, but major features are easily identifiable. Results in Section 4.2 show the gravity anomalies and

50

gravitational gradients calculated from a density distribution calculated from the seismic tomography model SL2013sv (Schaeffer and Lebedev, 2013). Figure 4.19 shows depth slices from SL2013sv representing the upper and lower portions of the mantle lithosphere in the North Atlantic. Figure 4.20 highlights depth slices within the asthenosphere beneath the North Atlantic. Toward the top of the asthenosphere a strong negative velocity zone is visible, whereas toward the base of the asthenosphere it is less apparent. Although the perturbation is low, a faint negative anomaly is present in the deeper portion of the asthenosphere (Figure 4.20) of the Mid-Atlantic Ridge (MAR) offset to the west of its surface representation.

The elongated influence of the Mid-Atlantic Ridge is seen in shallow depth slices (Figure 4.19), but reduces in magnitude faster than the observed round traces of the hot spot (probable plume conduit) under Iceland with increased depth (see additional figures in Appendix A).

Figure 4.19: Tomography slices of the North Atlantic from SL2013sv. a.) 47 km depth b.) 122 km depth. Slice from the top and base of the mantle lithosphere show concentrated areas of low-velocity, hot material associated with the active Mid-Atlantic Ridge (MAR) spreading centre, as well as the hotspot beneath Iceland.

52