Search

In seismic tomography we map the wave speed structure inside the Earth, but we ultimately seek to interpret those images in terms of physical parameters. This is challenging because many parameters can trade‐off with each other to produce a given wave speed. The problem is compounded by the convention of mapping seismic structures as perturbations relative to a 1‐D reference model, rather than absolute wave speeds. Using a full waveform tomography model of Europe as a case study, we quantify the extent to which thermochemical and dynamic properties can be constrained using only S wave speed, expressed in absolute values. The wave speed distributions of this tomography model are compared with 4 million thermochemical models, whose seismic properties are computed via thermodynamic modeling. These models sample the full range of realistic mantle compositions, including variable water and melt contents, and mineral intrinsic anelasticity is taken into account. Intrinsic anelasticity causes waves to travel more slowly at higher temperatures, leading to seismic attenuation, but the sensitivity of the wave speed reduction to temperature is, in turn, controlled by the wave frequency. Global studies of surface waves indicate an anticorrelation between S wave speed and attenuation. We therefore only retain thermochemical models satisfying this anticorrelation. Our study indicates that the frequency dependence of anelasticity, α, depends on temperature or rheology, with α ≈ 0.1 being most appropriate in cold or lithospheric mantle and α ≈ 0.3 in warmer regions (i.e., the asthenosphere). Additionally, the slowest regions require specific compositions and/or a velocity‐weakening mechanism, such as partial melting, elastically accommodated grain boundary sliding, or water., Geochemistry, Geophysics, Geosystems, 19 (10), ISSN:1525-2027

Geophysical Journal International, 194 (1), ISSN:0956-540X, ISSN:1365-246X

Imaging is a field of mathematics and physics that aims to retrieve information about the internal structure of an object that can only be accessed on its boundary. Many imaging methods are based on the following principle: a source outside of the object emits a wave. The wave propagates through the object. Wherever the physical structure of the object changes, scattered waves are induced. These scattered waves are measured by receivers outside of the object, and these scattered data are used to invert for the interior composition of the medium under investigation. The Marchenko integral was originally introduced for one-dimensional inverse scattering problems in the context of quantum mechanics. It can be related to Green's functions and so-called focusing functions - fields that produce a focus when injected into a medium from a single side. About ten years ago, the Marchenko integral was extended to two and three dimensions. This paved the way for, e.g., the elimination of imaging artefacts due to multiple scattering and Green's function retrieval for virtual source locations. However, many questions about the full potential as well as the accuracy of the Marchenko equation in two and three dimensions remain unanswered. In this thesis we present a new derivation for the multidimensional Marchenko integral. Our derivation is based on a generalised framework for wavefield focusing and circumvents the limiting assumptions of the previous extension. As we use partial differential equations rather than integral equations to define focusing functions, it allows for new physical insights. For instance, our approach indicates that it is possible to model Marchenko-type focusing functions with a conventional wave equation. Ultimately, this enables us to study Marchenko-type focusing in different 2D and 3D media and learn about the accuracy of the concept. We present a straightforward modelling approach for 1D as well as a least-squares modelling approach for 2D and 3D. The latter suggests that the Marchenko integral might be inherently approximative in multiple dimensions. We also discuss Green's function retrieval with our newly derived Marchenko integral, i.e. without wavefield decomposition. This method allows for estimating Green's functions for virtual sources inside of the medium. While it requires single-sided scattering data and an estimate of the first arrival of the desired Green's function there is no need to have an actual source or receiver inside of the medium. Our results demonstrate that we can retrieve good estimates of the full-spectrum Green's functions, involving evanescent and refracted waves, which were believed to not be retrievable with the previously derived Marchenko integral. Ultimately, we discuss imaging with these Marchenko-based Green's functions. Being able to include measurements for virtual sources inside of the medium allows for a natural linearisation of the imaging problem. Thus, we use the Marchenko integral to linearise state-of-the-art imaging approaches, similar to full waveform inversion or least-squares reverse time migration, and estimate the scattering potential. Our Marchenko-based linearisation accounts for all orders of scattering and performs slightly better than a single-scattering approximation.

The first study in this thesis is designed to cover both approaches of Marchenko redatuming and Marchenko-based primary estimation in a range of subsampling scenarios, providing insights that can inform acquisition sampling choices as well as processing parameterization and quality control, e.g., to set up appropriate data filters and scaling to accommodate the effects of dipole fields, or to help ensuring that the data interpolation achieves the desired levels of reconstruction quality that minimize subsampling artifacts in Marchenko-derived fields and images. With the aim of handling short-period multiples in laterally-varying media with complex, realistic layering, a updated augmented Marchenko approach is proposed in this thesis providing reliable approximations to the problem. This is achieved through combining an adapted 1.5D approach and additional post-processing with the recently proposed Marchenko-type demultiple formula. The new approach and its limitations are demonstrated with a range of numerical models that capture varying degrees of lateral heterogeneity in the overburden. The benchmark of migrated images after removing short-period multiples suggests promising applications to field data. As a next step in the reflection-data processing chain, the migrated images are nonetheless blurred with uneven amplitudes and low resolution, despite a successful multiples removal procedure. As an important post-processing step, seismic image deblurring aims to enhance both the resolution and amplitude balance. Conventional methods developed for such an inverse problem are usually cumbersome and slowly converging, and require case-dependent user interference, e.g. in the form preconditioning, and the fine-tuning of free parameters for regularization. This thesis proposes instead to address the problem with a physics-based Machine Learning technique --- the (invertible) Recurrent Inference Machine, by explicitly employing the point spread function as the forward operator and thus prior information. With a simple almost-flat synthetic training model and very cheap computational costs, this new approach outperforms the benchmark of the invertible UNet - a commonly-used convolutional neural network architecture - in a series of tests, designed with both complex synthetic models and the field data application. Different from the other methods, the neural networks in this study are trained to map the migrated image into an impedance perturbation model, rather than the reflectivity model, making it closer to the next-level application of the impedance inversion for reservoir characterisation and time-lapse monitoring. In the last part of the thesis, a general and unified framework is proposed, which flexibly incorporates all kinds of interface conditions for any type of wave equations. With the new approach, the implementation of a wave equation solver can be benchmarked by comparing different angles of plane wave simulations to transmission/reflection coefficients from plane-wave analysis with exact boundary conditions. This is illustrated with the case of poroelastic wave equation, which has a wide application in rock physics and reservoir monitoring.

During a megathrust earthquake cycle, the plates accumulate strain in the interseismic stage due to locking of a portion of the megathrust that separates them. This strain is released during the earthquake and following rapid postseismic relaxation. As summarised in Chapter 1, this understanding was built through decades of seismological and geodetic observations and of advances in physics-based models. Models link properties and structures to observable quantities and are crucial to interpreting observed surface deformation in terms of the processes and materials in the inaccessible subsurface. However, different models can approximate the same observations roughly equally well. Furthermore, more sophisticated models do not necessarily reflect better the processes and properties of the subsurface; at the same time, simpler models that can explain some observations do not necessarily reflect the key processes occurring throughout one or multiple earthquake cycles. This thesis uses relatively simple three-dimensional models that still capture the key processes occurring during repeated earthquake cycles, without attempting to reproduce the structure or properties of specific subduction zones, to examine possible explanations of geodetic observations from global navigation satellite systems (GNSSs) at different stages of the cycle. Chapter 2 focuses on the interseismic, horizontal deformation of the overriding plate, which is in apparent contrast with observations of far-reaching coseismic displacement. We estimate the spatial patterns, with uncertainties, of GNSS velocities in South America, Southeast Asia, and northern Japan. Interseismic velocities with respect to the overriding plate generally decrease with distance from the trench with a steep gradient up to a “hurdle”, beyond which the gradient is distinctly lower and velocities are small. The hurdle is located 500-1000 km away from the trench for the trench-perpendicular velocity component, and either at the same distance or closer for the trench-parallel component. The trench-perpendicular hurdle generally follows major tectonic or geological boundaries and seismological contrasts. We formulate and test the hypothesis that both the interseismic hurdle and the coseismic response result from a mechanical contrast in the overriding plate. Our models show that overriding plates with a sufficient contrast respond to locked interseismic convergence similarly to observations. The compliance contrast is probably mainly responsible for the observed hurdle and in turn results from thermal, compositional and thickness contrasts. Chapter 3 is concerned with the increased landward velocities that were recorded onshore after 6 megathrust earthquakes in subduction zone regions adjacent to the ruptured portion. We investigate whether bending can be expected to reproduce this observed enhanced landward motion (ELM). We find that viscous relaxation, but not afterslip, produces ELM when a depth limit is imposed on afterslip. This ELM results primarily from in-plane elastic bending of the overriding plate due to trenchward viscous flow in the mantle wedge near the rupture. Modeled ELM is, however, incompatible with the observations, which are an order of magnitude greater and last longer. This conclusion does not significantly change when varying key model parameters. The observed ELM consequently appears to reflect faster slip deficit accumulation, implying a greater seismic hazard in lateral segments of the subduction zone. In Chapter 4, we study postseismic landward motion observed on the overriding plate in the vicinity of a major megathrust rupture. We argue that relocking of the megathrust, particularly at shallow depths, is needed for postseismic relaxation to produce landward motion on the tip of the overriding plate. We discuss how that this finding is consistent with previous simulations that implicitly relock the megathrust where afterslip is not included. We conclude that the Tohoku megathrust relocked within less than two months of the earthquake. This suggests that the shallow megathrust probably behaves as a true, unstably sliding asperity.

The physical and chemical properties of the Earth’s mantle govern the cause of natural disasters, such as earthquakes and volcanoes. Since we do not have direct access to mantle materials, their properties are often inferred from laboratory measurements and surface observations (e.g. seismic data from earthquake recordings). This thesis addresses some key problems we face while utilising these data to constrain the thermal and chemical properties of the mantle. Firstly, we propose a data-driven approach based on machine learning to explain the laboratory measurements and quantify their uncertainties in the absence of an adequate physical model. Our results show that although conventional approaches based on fitting the measurements to an assumed model may appear better constrained, they could potentially provide biased results. Secondly, we use the data-driven approach to explore which thermochemical parameters can be constrained (and to what extent) with limited seismic observables- wave speeds and density. Our results show that these observables constrain temperature and major chemical parameters (silicon, magnesium, and iron), and they indicate the presence of thermochemical heterogeneities at the lowermost mantle. The dense and slow piles at the bottom of the lower mantle seen in seismic data can be explained by an enrichment in silica and iron content- characteristic feature of enstatite chondrites. The inferred heterogeneities have profound implications for the dynamics of the mantle and outer core. The methodology developed in this thesis is extremely efficient. It can easily incorporate additional observables and thus, has wide applications in the seismology and mineral physics community.

We performed an investigation of the large scale seismic wave speeds and density structure of the Earth’s mantle using free oscillations. Seismic free oscillations, or normal modes, are convenient for analysing low-frequency seismograms in a hetero- geneous Earth. To use these, we must address how to calculate exact seismograms using normal modes, and how to formulate the inverse problem to infer Earth’s 3D structure. The most important findings of this research are: • In order for seismograms to be theoretically exact, full mode coupling calcula- tions must involve an infinite set of modes. In practice, only a finite subset of modes can be used, introducing an error into the seismograms. We found that coupling modes 1-2 mHz above the highest frequency of interest is essential for having sufficiently accurate signals to infer density. • Observations of free oscillations provide important constraints on the heteroge- neous structure of the Earth. This inference problem has usually been addressed by the measurement and interpretation of splitting functions. These can be seen as secondary data extracted from low frequency seismograms. The measurement step necessitates the calculation of synthetic seismograms, but current imple- mentations rely on approximations referred to as self- or group-coupling and do not use fully accurate seismograms. We therefore investigated whether a systematic error might be present in currently published splitting functions. As is well known, the density signal is weak in low-frequency seismograms. Our results suggest this signal is of similar magnitude to the realistic uncertainties associated with currently published splitting functions. Thus, great care must be taken in any attempt to robustly infer details of Earth’s density structure using current splitting functions. • We investigated the problem of inferring density using currently published split- ting functions with properly calibrated uncertainties together with a novel prob- abilistic inversion technique, Hamiltonian Monte Carlo. Models are strongly dependent on damping. We found that shear wave speed models are statisti- cally significant in terms of misfit change, while density and compressional wave speeds are not. Therefore any interpretation of Earth’s mantle density based on splitting functions might be inaccurate. • A promising approach is the direct spectral inversion, which uses spectra di- rectly without the need of splitting functions. We found that misfit changes corresponding to the inferred models are statistically significant even for den- sity and compressional wave speed, but depend on a good starting model. We only used group coupling and relatively low frequency spectra for computa- tional reasons. Full coupling together with high frequencies might solve this long-lasting problem to infer density contrasts in the Earth’s mantle.

My thesis focuses on evolving and (geologically) short-lived plate boundary segments, their segmentation processes and geological imprints in the eastern Mediterranean. In Chapter 2, I investigate the nature/type of the plate boundary between the eastern Aegean region and the Africa plate. The work involves an integrative analysis of geological and geophysical information. I conclude that these surface observations document that the “Pliny-Strabo trench” is a predominantly strike-slip plate boundary. My interpretation is that this plate boundary represents an expression of slab tearing related to an active STEP (c.f. Figure 1.1). The paper represents the first detailed account of surface deformation related to a STEP fault, and constitutes a novel contribution to the understanding of the relation between deep processes and (near-) surface deformation, a key topic in geodynamic research. In Chapter 3, I investigate the location and nature of currently active plate boundaries and other major faults in the southern Anatolia-Aegean region, in the transition region from the Hellenic Arc to the Cyprus Arc. The question is particularly relevant for accessing earthquake hazard. I use mechanical models based on the finite element method. I explore various options for these faults, most of them proposed in the scientific literature, to explore how they would affect the deformation at locations where there are actual observations. The (mis)fit between model predictions and observations allows us to conclude that the active plate boundary is located offshore. The research question that I address in Chapter 4 is what the cause is of deformation in one of the seismically most active fault zones in Europe, the Kefalonia Transform Fault. I present results from a recent full-waveform tomographic model which particularly improves our understanding of the structure of the upper few hundred kilometers of the Earth. The cause of the deformation along the Kefalonia Transform Fault is likely rooted in a fragmented slab that we image for the first time. The geometry of the slab fragment leads me to conclude that it became disconnected from the larger Hellenic slab around 5 Ma, at about the time of opening of the Gulf of Corinth in the overriding plate, which suggests a highly interesting causal relation.

The goal of predicting earthquakes remains elusive despite decades of instrumental observations and research, and a much longer historical record. Even the practice of seismic hazard analysis is a topic of heated debate, in part due to our inability to accurately determine the rate and size distribution of earthquakes that a fault can produce. The main reason for these deficiencies is the lack of a validated, physics-based theory of earthquakes. Heuristic attempts to discover patterns in seismicity based on its phenomenology have produced ambiguous and sometimes contradicting results, due to the relatively short instrumental record of big earthquakes compared to their rate of occurrence. From a geodynamics perspective, earthquakes are bursts of energy release as the lithosphere is loaded due to the motion of tectonic plates at rates of a few centimeters per year. Similar behavior, known as crackling, is observed when shearing granular aggregates. Loosely packed particles behave collectively as a fluid, giving rise to small instabilities only. At a critical packing fraction, the size distribution of the instabilities approaches power law scaling. This suggests that the aggregate is at a phase transition and that long-range correlations are a key characteristic of its macroscopic behavior. Above the critical packing fraction, the collective behavior of the particles is similar to that of a solid. In that solid-like regime, the aggregates alternate between power law distributed event sizes and quasi-periodic stick-slip. A significant number of laboratory studies have employed granular media to explore the dynamics of critical systems in the context of seismicity and fault gouge rheology. These studies have been performed either at low normal stress (< 1 MPa) or to limited shear displacements (< 50 mm), and often under dry conditions. It is not known whether the macroscopic behavior of granular aggregates remains the same under higher normal stress and larger displacements, or in the presence of pressurized water. If not, is it possible to determine what mechanisms are responsible for the change? The rotary shear experiments presented in this thesis expand the envelope of the experimentally tested conditions up to 8 MPa normal stress and 165 mm of shear displacement, \textit{simultaneously}. This enabled us to infer the emergence of correlations under these conditions, through changes in the statistics of granular avalanches. Because the elevated stress conditions do not allow direct visual observation of the glass bead samples, a specially developed AE monitoring system was used to detect and locate the source of crackling. The findings of this thesis highlight the importance of emergent, long range correlations in sheared granular media, as a function of experimental conditions. We infer that the key parameter that determines the scaling of avalanche statistics is the packing fraction, which in turn depends on normal stress, wear rate, and particle size distribution.

Seismology aims at obtaining accurate tomographic images of the Earth’s interior by simulating models to create waveforms that fit recorded seismograms. The resolution of an acquired image greatly depends on the accuracy of the numerical tool used for modelling and the quality of observed data. Using a state-of-the art numerical wave propagation software, I study the structure of global discontinuities. I develop an iterative optimisation methodology for modelling the waveforms by minimising the misfit caused by the existence of topographic structure on discontinuities. Given that the disconti- nuity structure has mainly been studied in a ray theoretical framework, I only use synthetics in order to assess the reliability of conventional methods and to develop a novel approach based on full waveforms and non-linear min- imisation. My study also focuses on the sensitivity of waveforms related to discontinuity structure. Analyses of their exact sensitivity pave the way towards a better comprehension of real data and improvement of the inversion methodologies. To that extent, a new inversion method is proposed which relies on the iterative optimisation of boundary and structural models, with special focus on the structure of global discontinuities using boundary Fréchet derivatives for the first time in an inversion problem. Successive steps for the iterative optimisation of the model are: choose a starting model and select a misfit function to calculate the discrepancies between observed and synthetic data. The objective function is the most important step in the inversion. By computing the derivative of this function, employing time-reversal and adjoint methods, one creates a model update which should minimise the previous misfit. Adjoint methods rely on the inter- action between "forward" and "adjoint" wavefields, propagating from source to receivers and vice versa. This process is iterated until the global misfit value sufficiently reduces. In this thesis, it is shown that this novel approach for imaging discontinuities improves the inference of internal discontinuity structure and provides an integrated method for global seismology. The pro- posed full waveform methodology outperforms ray theory to a great extent and should be used in real data applications.

Nienke Blom studied methods to image variations in density inside the Earth using the technique of waveform tomography, a method based on modeling earthquake vibrations. Her work contributes towards better understanding the interior of the Earth, the dynamics of which drive everything that happens at the Earth's surface. As heavy material in the Earth sinks and light material rises, the Earth's internal structure constantly moves and changes at the speed of the growth of a finger nail. This has its effect on the outside: everything that happens at the Earth's surface -- the shape of continents, the occurrence of earthquakes and volcanism, but also the existence of a breathable atmosphere -- is a consequence of this constant motion. So far, however, it has been very difficult to image the variations in density that cause these motions in an independent and reliable way. During her PhD, Nienke Blom worked on developing ways to image these variations in density using the vibrations that are emitted from earthquakes. By modeling on a supercomputer how earthquake vibrations propagate through the Earth's interior, she deduced how the wavefield is altered by the presence of variations in density. Based on this, she developed a workflow to derive the structure of the Earth's interior from a comparison between measured seismograms and the modeled wavefield, which she is currently applying to the Eastern Mediterranean.

It is very difficult to study mantle convection over periods of millions of years because convection is a non-linear process. In this thesis, I present a new method for studying the underlying statistics in convection patterns. I use neural networks to find patterns in these statistics to make inferences about the mantle and its history. I can make inferences about constant rheological parameters, evolving time dependent parameters, such as the development of LLSVPs, and the compositional, thermal and viscosity structure of the mantle. All of these parameters have important implications for the formation of Earth, evolution of plate tectonics and therefore life, interpretation of geophysical observations and understanding of dynamic processes. They are all current poorly constrained, making any new method potentially powerful. I also use neural networks as a predictive tool. Every inference is made using a Bayesian approach and is therefore fully probabilistic and includes uncertainty estimates. These uncertainty estimates are in themselve novel to geodynamics.

This thesis reports on studies that have been carried out to better image the crustal and upper mantle shear velocity structure around the Gulf of California, Mexico. The Gulf of California forms a part of the plate boundary between the Pacific and North-American plates, where transform motion in the North converts into oceanic spreading in the South. Although currently oblique extension occurs within the gulf, in the past subduction of the oceanic Farallon plate took place along a trench west of what is currently the Baja California Peninsula. Around 12 Ma subduction and spreading ceased and the remaining parts of the Farallon plate were captured by the Pacific plate. After this the Gulf of California started to develop and around 5 Ma the Pacific-North America plate boundary shifted to its present location within the gulf. The study is subdivided into three parts. In the first part we use a spectral-element code to model seismic wave propagation and the adjoint technique to calculate finite-frequency sensitivity kernels for two-station surface wave measurements. We found a dominant sensitivity in the interstation region, close to what is assumed in the two-station method. However, we also show that anomalous structures even far from the assumed propagation path can significantly affect the measurements. To improve regional tomographic inversions using the two-station method, sensitivity kernels should therefore be taken into account. At the regional scale we are interested in this requires high computational costs, which are difficult to reach in practice. In the second part, instead of improving the lateral resolution using the method described above, we improve the vertical resolution by 1-D joint inversions of surface wave and receiver function data recorded by NARS-Baja stations. We use a Neighbourhood Algorithm to obtain shear velocity models beneath the stations. We are able to find sediment and crustal thicknesses below the stations in the area. In the upper mantle we find indications for the presence of subducted plate remnants below the middle and southern parts of the Baja California Peninsula and we identify the top as a sharp interface. The absence of plate remnants below the northern part explains why the middle and southern parts of the Baja California Peninsula move faster to the Northwest with respect to the North-American plate than the northern part. In the third part we investigate the effects of dipping crustal layers on the receiver functions. We find that azimuthal differences is the receiver function data can be explained by a westward dipping Moho below the Baja California Peninsula, and a dipping mid-crustal layer throughout the entire area. For detailed crustal imaging of the Baja California area using receiver functions more complex modelling than a 1-D inversion is recommended.

Numerous problems in the field of seismology require the determination of parameters of a physical model that are compatible with a set of observations and prior assumptions. This type of problem is generally termed inverse problem. While, in many cases, we are able to predict observations, given a particular set of parameters with high accuracy and precision using numerical and analytical methods, often we are unable to directly infer model parameters from observations. This has lead to the development of a large number of inverse methodologies. For many applications it is necessary not only to obtain a set of best fitting parameters, but also an estimate of the attached uncertainties, which arise from the fact that any observation is subject to measurement uncertainties and from other sources of uncertainty that may be present due to numerical limitations or simplifying assumptions. Often, the situation is further complicated by the non-linear nature of many inverse problems, leading to non-Gaussian and possibly non-unique posterior distributions. An important seismological inverse problem is the determination of earthquake source parameters. Source parameters are on one hand needed as an input for subsequent investigations, such as tomographic imaging or finite fault studies. On the other hand source parameter estimates are required in order to provide information to the public as quickly as possible after a new earthquake has initiated - a field known as earthquake early warning. We have developed a methodology for rapid probabilistic moment tensor point source inversions able to cope with a wide variety of data types. Accurate and computationally expensive synthetic forward modelling can be incorporated without significantly altering the time required for the inversions. Our methodology is based on finding an approximation to the conditional posterior probability of source models given observations by smoothly interpolating a set of prior samples. The interpolation is obtained using ensembles of Mixture Density Networks (MDNs) - a class of artificial neural networks, whose output are the parameters of a Gaussian mixture model. Once an ensemble has been constructed, new observations can be inverted within the fraction of a second on a standard desktop computer. The method is therefore well suited for earthquake early warning, where new observations must be inverted routinely and rapidly. We apply the method to two scenarios in Southern California using both synthetic and observed datasets. We invert regional static and dynamic GPS displacement data for the 2010 M 7.2 El Mayor Cucapah earthquake in Baja California to obtain estimates of magnitude, centroid location and depth, and focal mechanism. Moreover, we show how accurate 3-D spectral element simulations can be incorporated into the rapid inversion scheme by means of a model scenario focussing on the 2008 M 5.4 Chino Hills earthquake. We discuss general implications, advantages and limitations of the method and suggest directions for potential future research.

The gross features of seismic observations can be explained by relatively simple spherically symmetric (1-D) models of wave velocities, density and attenuation, which describe the Earth's average(radial) structure. 1-D earth models are often used as a reference for studies on Earth's thermo-chemical structure and dynamics, earthquake location determination and 3-D seismic tomographic models. Therefore, the quality of the latter is intrinsically linked to the robustness of the former. But what is the quality of such 1-D models? Seismic inverse problems are notoriously non-unique; different earth models can explain the data equally well, but may lead to incompatible interpretations of the nature of the Earth's interior and dynamics. Ideally, the assessment of solution quality is an integral part of any inverse method. The main motivation for this thesis was to investigate a means to simultaneously infer Earth structure and quantify the uncertainties in our estimates. A common (Bayesian) approach is to directly sample the posterior model probability density, as is done in Markov Chain Monte Carlo (MCMC) methods via a (guided) random walk. Here, I solve such seismological inference problems using pattern recognition and machine learning techniques. The method developed here, using articial neural networks, is exible and enables me to address specific hypotheses on Earth's structure in a robust, quantitative manner. Using normal mode splitting function coefficients and body wave travel times, I obtain complete statistical descriptions of features of radial Earth structure, in terms of elastic and anelastic structure, anisotropy and depths of major discontinuities. In general, I conclude that a lot can still be learned on 1-D Earth structure from seismic data; ideally, we do so prior to tackling the 3-D tomographic problem. An analysis of the information content suggests that the free oscillations constrain most parameters better than the body wave data. Spheroidal and toroidal mode data constrain the depth extent of the density excess in the lowermost mantle. Furthermore, I show, for the first time, that the average lower mantle is anisotropic below 1900 km depth, challenging the consensus that this part of the mantle is isotropic. It is possible to explain these seismic observations with currently available mineral physics data for lower mantle minerals. Therefore, seismic anisotropy, such as observed here, can provide constraints on mantle flow and deformation mechanisms. However, meaningful geodynamic interpretations require a full 3-D analysis to be made. Finally, I illustrate one pragmatic approach to data dimensionality reduction using autoencoder networks, as a first step towards non-linear seismic waveform inversion using encoded seismograms. The machine learning method adopted in this thesis is a pragmatic approach to solving non-linear Bayesian inverse problems. Their exibility and interpolation capabilities, in combination with the quantitative Bayesian framework, make neural networks well-suited for data sensitivity analysis and testing hypotheses on Earth structure. Rather than a replacement for Monte Carlo methods, I suggest that in the future they are used as a complementary tool, providing an initial assessment of data sensitivity and a lower bound on the information on model parameters that is contained in the data.

Within the last decade, increasing computational power and theoretical advances have initiated the transition from approximation-based tomography to iterative full-waveform tomography in regional and global seismology. Such full-waveform methods allow for the exploitation of complete seismograms to constrain the 3-D structure of the Earth, which potentially leads to much higher resolution compared to approximation-based methods. The method used in this thesis relies on numerical wavefield simulations with the spectral-element method, which honours the full complexity of seismic wave propagation in realistic 3-D media. Gradients are constructed with the adjoint method. Full-waveform models are affected by a number of subjective inversion choices, such as the type of misfit and the regularisation of the gradient. The complexity of seismic wave propagation together with the non-linearity of the problem make it a challenging task to determine optimum inversion strategies. For this thesis, we first performed synthetic experiments to investigate and optimise the capability of full-waveform tomography to constrain small-scale structures in the Earth's mantle. These experiments were focused on the tomographic detection of possible mantle plumes, which are a challenging target for seismic tomography due to their small lateral extent and their extension deep into the lower mantle. Using a misfit based on the time-continuous and amplitude-independent instantaneous phase difference, we succeeded to fully recover an idealised narrow mantle plume. The accurate measurement and inclusion of diffracted waves arriving in the body-wave coda was identified as essential for the recovery in the lower mantle. We further confirmed that simple cross-correlation traveltime misfits do not permit the recovery of plumes in the lower mantle, wavefront healing largely conceals the acquired time delays. We applied these results to the construction of a high-resolution 3-D model of the S-wave speed beneath the North Atlantic region, extending to a depth of 1200 km. Unprecedented details are revealed, giving new insight into the complex tectonics and dynamics of the North Atlantic region. The resolution of particularly interesting features was validated with resolution tests based on Hessian kernels. Striking features of the new model include individual low-velocity conduits beneath the Iceland and Jan Mayen hotspots, extending into the lower mantle. This observation strongly supports the existence of mantle plumes as a feature of mantle convection, which has been debated passionately over the last decades due to a lack of conclusive seismological evidence. Our model is furthermore the first to clearly resolve two separated hotspots associated with Iceland and Jan Mayen. A further highlight is the observation of a wide-spread layer of low-velocity plume material beneath the lithosphere of much of the North Atlantic ocean. Two fingers of this layer extend beneath the continental lithosphere of parts of the British Isles and of southern Norway, providing an explanation for the considerable post-rift uplift that these regions experienced in Neogene times. Gravity-based estimates of present-day dynamic support agree very well with the location of the low-velocity layer.

We present a palaeomagnetic study on 38 lava flows and 20 dykes encompassing the past 1.3 Myr on S. Jorge Island (Azores Archipelago—North Atlantic Ocean). The sections sampled in the southeastern and central/western parts of the island record reversed and normal polarities, respectively. They indicate a mean palaeomagnetic pole (81.3◦N, 160.7◦E, K = 33 and A95 = 3.4◦) with a latitude shallower than that expected from Geocentric Axial Dipole assumption, suggesting an effect of non-dipolar components of the Earth magnetic field. Virtual Geomagnetic Poles of eight flows and two dykes closely follow the contemporaneous records of the Cobb Mountain Subchron (ODP/DSDP programs) and constrain the age transition from reversed to normal polarity at ca. 1.207 ± 0.017 Ma. Volcano flank instabilities, probably related to dyke emplacement along an NNW–SSE direction, led to southwestward tilting of the lava pile towards the sea. Two spatially and temporally distinct dyke systems have been recognized on the island. The eastern is dominated by NNW–SSE trending dykes emplaced before the end of the Matuyama Chron, whereas in the central/western parts the eruptive fissures oriented WNW–ESE controlled the westward growth of the S. Jorge Island during the Brunhes Chron. Both directions are consistent with the present-day regional stress conditions deduced from plate kinematics and tectonomorphology and suggest the emplacement of dykes along pre-existing fractures. The distinct timing and location of each dyke system likely results from a slight shift of the magmatic source.

In this thesis, we investigate the feasability of using higher modes in seismic interferometry. This is done by a theoretical and computation analysis of the approximations involved, development of a method of stacking noise traces, and analysis of a small array in the northern part of the Netherlands. Higher modes are unambiguously present in ambient noise, and can thus be used to improve vertical resolution in the subsequent inversion procedure. The practical issue in their use, remains the problem of distinguishing them from spurious events/noise, and even distinghuising higher modes from the fundamental mode, since at high frequencies they can be in close proximity in a time-frequency diagram. Higher modes can furthermore be lost in the data due to: -their small amplitude compared to the fundamental mode. -intermodal cross-terms that occur with imperfect distributions, but occur even in the case of a perfectly surrounding integration surface. However, if noise is generated in a large band-like distribution (or very far away) the cross-terms should pose no problem. - source-source correlation terms. The signal to noise ratio is much smaller for higher modes. The exact instantaneous phase takes much longer to converge to the true solution than the envelope, implying larger uncertainty for phase velocity measurements. -incoherent noise not associated to wave propagation, or local effects. Directional noise shows much faster convergence for higher modes, again at least for station paths inline in this direction. Also in our data, the noise is strongly directional. It can be argued that the higher mode information would not be retrieved if this were not the case. The higher mode has a different origin from the dominant secondary microseism.The TF-phase weighted stack is a powerful tool to denoise the retrieved seismograms. It is possible to measure the group velocity from the phase stack directly. To retrieve higher modes, the phase stack provides an easy to use, additional tool. Most information remains in relative few TF-stacks from long displacement traces. Verification with additional information was essential for a robust estimate of the 1D S-velocity structure beneath Annerveen, The Netherlands. As an average for the region, the higher mode information gives additional constraint on the deeper sedimentary layers. The tomographic images seem to correlate with the main geological features, but for the higher modes the amount of interstation paths that can be reliably measured is not large enough. This is due to interference with the fundamental mode in the group velocity images and the limited directionality. The choice of a substantial spread in interstation distances led to a broad frequency range of analysis, but also limited the number of measurements per frequency. A recommendation for future studies is to use a strict 1-D configuration along the prevalent noise direction. The advantages are the large possible frequency bandwidth, and the possibility to use and invert the higher modes.

Over the last three decades a vast number of tomographic images has been produced, but the quantitative assessment of their accuracy and uniqueness has only just started. A relatively recent technique in this direction has been made by (Trampert et al., 2004) using probabilistic tomography. It is based on a full model space search, the Neighbourhood Algorithm (NA), with which to solve the inverse problem. This provides a representation of seismic constraints as probability density functions (pdfs) from which one can easily calculate the most likely model and its associated uncertainty. Since a full model space search quickly limits the dimension of the explorable parameter space, an essential condition to apply the NA is to have a small number of unknowns. Probabilistic tomography has been applied to invert surface waves and normal-modes. The aim of this thesis was to incorporate body wave data to increase the depth resolution in the lowermost mantle. The first step, therefore, was to establish a formalism which allows us to treat body waves as surface waves and free oscillations and, thus, separate the inverse problem into its lateral and radial components. In a spherically layered Earth using the path-average approximation, body wave delay times can be inverted in two steps. The first step is the construction of vertical travel-time residual maps as a function of ray parameter. The second step is their depth inversion, less well constrained than the first step but characterized by a small number of parameters and, hence, suitable to be solved by the NA. Such a procedure, called path-average approach, works well for long wavelength structure for the forward and inverse problem. Combing a body wave dataset characterized by a good resolution between 2441 and 2891 km and the advantages of the path-average approach, we provided a new insight into the lowermost mantle based on lateral variations of P-, S-seismic velocities and the corresponding seismologically observed ratio without an explicit depth inversion for seismic wave velocities variation. Following a robust statistical approach, we are able to link these seismic heterogeneities to mineral physics and provide evidence for the existence of post-perovskite near the core-mantle boundary. In the framework of probabilistic tomography, applying the path-average approach we inverted body wave vertical travel-time residual maps, normal-mode splitting functions and surface wave phase velocity maps locally for depth using the NA. We constructed a shear wave, compressional wave, bulk sound speeds and density perturbation model with associated uncertainty. The estimation of the correlation between parameters and the RMS amplitude of seismic variations clearly indicates a thermal and chemical nature of heterogeneity in the lower mantle. Furthermore, our models fully confirm the conclusions achieved by the model of Resovsky & Trampert (2003), produced without body waves. One strength of probabilistic tomography is the possibility to convert easily the pdfs of seismic velocities and density perturbations into pdfs of thermal and compositional anomalies. We constructed a model of temperature, iron and perovskite variations combining mineral physics data with the shear, bulk sound wave-speeds and density variation models described above.

We have reached a stage in seismic tomography where further refinements with classical techniques become very difficult. Advances in numerical methods and computational facilities are providing new opportunities in seismic tomography to enhance the resolution of tomographic mantle images. 3-D numerical simulations of seismic wavefields also allow us to check the reliability of the current mantle images and the classical techniques. In this thesis, using 3-D wave simulations by a spectral element method, we present three specific studies which underline some shortcomings in global mantle tomography. We first investigated crustal corrections for fundamental mode surface waves. We compared the crustal corrections estimated by first order approximations such as great circle approximation, exact ray theory and finite-frequency theory, to those obtained from 3-D wave simulations. We observed that crustal corrections can produce errors larger than those in phase measurements at periods longer than 60 and 80 s for Rayleigh and Love waves, respectively. Rayleigh and Love waves are differently influenced by the crust. Therefore imperfect crustal corrections have potential to produce radial anisotropy which may bias our interpretation of anisotropy in the upper-mantle. Extensions to the great circle approximation do not improve the results because of the highly non-linear behaviour of the crust. Thus either the mantle has to be inverted together with the crust or 3-D background models have to be used for the measurements. We then compared real data with synthetic seismograms computed in various 3-D mantle models. We particularly investigated the effect of damping on global mantle models. In addition, we tried to examine the elastic and anelastic contribution of the models by comparing seismograms computed in 1-D and 3-D crustal and mantle models. Our results show that different 3-D mantle models give statistically similar results although they are different from each other. We observed that a large part of the seismograms remained unexplained, particularly the amplitudes. 3-D velocity models are not enough to explain amplitude information alone. Attenuation may have an important contribution to surface-wave amplitudes whereas body waves are more affected by scattering effects. Lastly, we proposed new misfit functions based on instantaneous phase and cross-correlation measurements for full waveform inversion. In addition, we qualitatively compared them to classical ones such as waveform and travel-time misfits by computing the finite-frequency sensitivity kernels using adjoint methods based upon a spectral element method. We observed that instantaneous phase measurements provide complementary information to the other misfit functions and are promising to exploit more information from a single seismogram.

This thesis is a study on upper mantle shear velocity structure beneath the Gulf of California. Surface wave interstation dispersion data were measured in the Gulf of California area and vicinity to obtain a 3-D shear velocity structure of the upper mantle. This work has particular significance for understanding the tectonic evolution history in this region and the potential processes that control lithospheric rupture and the initiation of seafloor spreading in the Gulf of California. The NARS-Baja project is designed as a passive seismic experiment to explore the mantle beneath the Gulf of California region. This thesis presents the results on the upper mantle structure beneath the Gulf of California from surface wave dispersion analysis. The study contains three aspects: the Rayleigh wave phase velocity structure with azimuthal anisotropy in the Gulf of California region; a shear velocity model of Gulf of California and its vicinity inferred from Rayleigh wave dispersion data; the radially anisotropic shear wave velocity model inferred from Love and Rayleigh wave dispersion . The seismic models of the upper mantle beneath the Gulf of California and surrounding regions show some interesting features. The average shear-velocity structure has a strong low-velocity zone beneath a thin lid. The average phase velocity data are even suggestive of lower velocities than the global average down into the transition zone. The low velocities at asthenospheric depths are indicative of fluids and/or melt. A striking feature of our seismic models is a relatively high shear-velocity anomaly beneath the central-southern gulf extending from a depth of roughly 120 to 160 km. This anomaly is interpreted as a remnant of a subducted Farallon microplate fragment . There is no evidence for such a slab remnant beneath the northern part of the gulf, in agreement with the tectonic interpretation of the region. We therefore interpret this region as having a slab window, similar to the area at the mouth of the gulf. Additional evidence on the upper mantle structure is obtained from anisotropy. An along-rift-axis variation is observed in the azimuthally anisotropic phase velocity maps of the fundamental mode Rayleigh waves, as well as in the radially anisotropic structure . The obtained upper mantle structure sheds light onto the interior of this dynamic plate boundary developing area. Our new tectonic interpretation of evolution the Gulf of California is based upon tectonic reconstructions, volcanic evidence, and our seismic interpretation of a stalled slab remnant. In our model, the development of the slab window beneath the North American lithosphere and pre-existing fracture zones in the subducting plate are essential in the tectonic evolution. The details of the subduction dynamics and the plate boundary relocation can therefore be largely be attributed to pre-existing weak features in the plates. Ridges and fracture zones of the subducted plate can determine tearing and fracturing of the slab, whereas faults and extended back-arc basins in the overriding plate can be essential to the plate boundary relocation.

We present a neural network approach to invert surface wave data for discontinuities and velocity structure in the upper mantle. We show how such a neural network can be trained on a set of random samples to give a continuous approximation to the inverse relation in a compact and computationally efficient form. The trained networks are applied to a real data set consisting of surface wave dispersion measurements. For each inversion, performed on a global grid, we obtain the a posteriori probability density function (pdf) of crustal thickness and independently of vertically averaged crustal shear wave velocity. From this pdfs, any desired statistic such as mean and variance can be computed. The obtained results are compared with current knowledge of crustal structure. Generally our results are in good agreement with other crustal models. However in certain regions such as central Africa and the back arc of the Rocky Mountains we observe a thinner crust than the other models. We also see evidence for thickening of oceanic crust with increasing age. We show that the resulting model allows to compute global crustal corrections for surface wave tomography for periods T>50 s for phase velocities and T>60 s for group velocities. In a second application we invert surface wave overtone measurements for posterior pdfs of the 400- and 660-km discontinuities as well as the transition zone thickness. We find that posterior uncertainties are almost as large as the observed variations of the mean 400- and 660-km depth as well as the mean transition zone thickness. Nevertheless, the model of the mean 660-km topography agrees reasonably well with other global models constrained by SS precursors, giving further evidence for a well defined and sharp 660-km discontinuity. The situation for the 400-km discontinuity however is more complicated. Agreement between four different 400-km topography models is very poor and we find a global average discontinuity depth of 394 km which is significantly lower than 409 and 418 km reported previously. Surface wave overtones and SS precursors should give similar results for a sharp 400-km discontinuity. If the 400-km discontinuity however is broad, it might be that surface wave overtones map an average depth that differs from the bottom of the discontinuity as mapped by SS precursors. Experimental studies suggest that if water concentrations are sufficiently high the 400-km discontinuity might be as broad as 20 km. The observed difference in the mean depth of the 400-km discontinuity between surface wave overtone measurements and SS precursor studies might therefore be indicative of high water concentrations near 400-km depth on a global scale.

In this thesis we present all three stages of the inversion approach proposed by Kennett and Yoshizawa (2002). The three stage inversion approach consists of obtaining fundamental and higher mode Love and Rayleigh wave phase velocity measurements through waveform fitting in the first stage, combining them into multimode phase velocity models using the path average assumption in the second stage and an inversion for local shear wave speed properties to obtain a three-dimensional shear wave velocity model in the third stage. We first present a fully automated Monte Carlo modelspace search approach to measure fundamental and higher mode Love and Rayleigh wave phase velocities. The advantage of this approach is that we do not only obtain the best fitting model but the whole ensemble of models with their corresponding fit. This allows us to obtain mutually consistent uncertainties for the fundamental and higher mode phase velocity measurements. Our phase velocity measurements agree remarkably well with previous studies, but we have been able to enlarge the available dataset dramatically (>350,000 measurements). Surface waves are well suited to observe anisotropy since they carry information about both radial and azimuthal anisotropy. The presence of anisotropy is well established in the Earth's uppermost mantle, but it is more ambiguous in the deeper mantle. Since we use higher mode phase velocities which have a larger sensitivity to deeper structure compared to the fundamental mode, we improve our understanding of anisotropy into the deeper mantle. In the second stage, we thus invert the phase velocity measurements to obtain isotropic and azimuthally anisotropic phase velocity maps. Prior to inversion, we determine the optimum relative weighting for the isotropic and azimuthally anisotropic terms. We found that all fundamental and higher mode measurements require anisotropy. Special care was taken to obtain the posterior model uncertainties, needed to obtain the radial anisotropic shear wave velocity model. The azimuthal anisotropic phase velocity maps give indications of anisotropy at larger depths but a more detailed depth inversion using finite frequency kernels will provide more information on the depth extend of azimuthal anisotropy. Finally, we inverted the local azimuthally averaged phase velocity maps for radial shear wave anisotropy taking the full non-linearity into account using a Monte Carlo modelspace search approach. For the anisotropic model, we decided to compute the total probability of radial anisotropy and likewise that the amplitude of radial anisotropy is above 1% or 2%. We find a lithosphere dominated by fast horizontally polarized shear wave anisotropy (horizontal flow), an asthenosphere dominated by fast vertically polarized shear wave anisotropy (vertical flow) and a transition zone also dominated by fast vertically polarized shear wave anisotropy (vertical flow). The amplitude of anisotropy is likely large (>2%) in the lithosphere and asthenosphere and smaller in the transition zone (1%-2%). The lower mantle is most likely isotropic (or anisotropic with a very small amplitude

Seismic reflection methods are widely used for the detection of hydrocarbons in subsurface structures up to several kilometers depth. However, since most data are acquired at or close to the Earth's surface, it is essential to understand the influence of the near-surface on the acquired data in order that its effects are not interpreted as pertaining to the reservoir. The near-surface effect on seismic data has two main origins: (i) near-surface wave propagation and (ii) wavefield acquisition. Wavefield acquisition comprises both wavefield excitation and wavefield measurement, i.e. source and receiver effects. Since both P- and S-wave velocities are observed to vary rapidly close to the Earth's surface, wave propagation in the near-surface is often very complex. In order to improve our understanding of near-surface wave propagation, we investigated and developed methods to determine near-surface P- and S-wave velocities. For this purpose, we used dense recording geometries including buried geophones, since these material properties cannot be resolved with a conventional acquisition geometry using geophones only at the free surface. One of the methods developed is based on the estimation and inversion of the propagator matrix, and therefore referred to as propagator inversion. The propagator inversion was applied to data which we acquired in Zeist, The Netherlands. We obtained a low near-surface P velocity, namely 270 ± 15 m/s, which is well below the sound velocity in air, and 150 ± 9 m/s for the S velocity. The buried geophone was located at approximately 1.0 m depth, thus the obtained velocities are only representative of the top meter of the near-surface. The second part of this thesis is devoted to the influence of wavefield acquisition. Corrections for source and receiver perturbations are necessary when their behaviour changes within a given survey, and should be performed in the early stages of processing. However, existing techniques, such as surface-consistent deconvolution, require prior processing before even they can be applied. We developed an alternative approach to compensate for source and receiver amplitude perturbations which has the advantage of being purely a raw data preprocessing step. It is applicable to the whole seismic trace, and does not impose additional assumptions on the subsurface. The approach is based on reciprocity of the medium response. This implies that differences between normal and reciprocal traces can be attributed to the source and receiver perturbations. We successfully demonstrated the procedure to compensate for these perturbations on both synthetic and field data. The field data were acquired in Manistee County, Michigan (courtesy of WesternGeco). Along the acquisition line, near-surface conditions change from moist-to-wet sediments to dry sands. The obtained source corrections are strongly correlated to these changing near-surface conditions, whereas the receiver corrections vary more strongly from geophone to geophone. Compensation of the recorded data for the source and receiver perturbations resulted in a significant improvement of the signal-to-noise ratio, both on prestack and poststack data. Finally, the receiver response we found did not agree with the generally accepted damped harmonic oscillator model, implying that this model need to be revised.