Scaling full seismic waveform inversions

The main goal of this research study is to scale full seismic waveform inversions using the adjoint-state method to the data volumes that are nowadays available in seismology. Practical issues hinder the routine application of this, to a certain extent theoretically well understood, method. To a large part this comes down to outdated or flat out missing tools and ways to automate the highly iterative procedure in a reliable way. This thesis tackles these issues in three successive stages. It first introduces a modern and properly designed data processing framework sitting at the very core of all the consecutive developments. The ObsPy toolkit is a Python library providing a bridge for seismology into the scientific Python ecosystem and bestowing seismologists with effortless I/O and a powerful signal processing library, amongst other things. The following chapter deals with a framework designed to handle the specific data management and organization issues arising in full seismic waveform inversions, the Large-scale Seismic Inversion Framework. It has been created to orchestrate the various pieces of data accruing in the course of an iterative waveform inversion. Then, the Adaptable Seismic Data Format, a new, self-describing, and scalable data format for seismology is introduced along with the rationale why it is needed for full waveform inversions in particular and seismology in general. Finally, these developments are put into service to construct a novel full seismic waveform inversion model for elastic subsurface structure beneath the North American continent and the Northern Atlantic well into Europe. The spectral element method is used for the forward and adjoint simulations coupled with windowed time-frequency phase misfit measurements. Later iterations use 72 events, all happening after the USArray project has commenced, resulting in approximately 150`000 three components recordings that are inverted for. 20 L-BFGS iterations yield a model that can produce complete seismograms at a period range between 30 and 120 seconds while comparing favorably to observed data

Investigation of Ambient Seismic Noise using seismic interferometry in the Midwestern United States

The portion of the North American craton occupied by the Midwestern United States is a cratonic platform, where a veneer of Phanerozoic sedimentary strata buries the Precambrian basement up to 7 km. Due to the sediment cover and low topographic relief, the at depth structure of the region remains poorly understood. This region is of particular interest because over the past half-billion years tectonic forces have resulted in the formation of epeirogenic provinces in a stable cratonic interior. Using the OIINK flexible seismic array and the Earthscope Transportable Array, Ambient Seismic Noise Tomography was applied to investigate the crustal structure and produce high-resolution structural models of the region. For our analysis we used the vertical component of seismograms recorded between January 2011 and December 2014, where spurious events were filtered out to establish the background seismic noise of the region. Seismic observations based on the cross-correlations of seismic noise from 46,665 station pairs were used to obtain phase velocities at periods from 4 to 40 s. From these measurements a high resolution model was formed, improving the our understanding of the complex basement geology of the Midwestern United States

Doctor of Philosophy

dissertationWe first study the inverse problem of recovering a complex Schro ?dinger potential from a discrete set of measurements of the solution to the Schro ?dinger equation using different source terms. We solve this problem by generalizing the inverse Born series method to nonlinear mappings between Banach spaces. In this general setting, we show convergence and stability of inverse Born series follow from a single problem- specific bound. We show this bound for the inverse Schro ?dinger problem, and study numerically an application of this inverse problem to transient hydraulic tomography. Additionally, we develop a family of iterative methods based on truncated inverse Born series that are akin to iterative methods based on truncated Taylor series. Next, we study the inverse problem of imaging scatterers in a homogeneous medium when only intensities of wavefields can be measured. Classic imaging meth- ods, such as Kirchhoff migration, rely on phase information contained in full waveform data and thus cannot be used directly with intensity-only data. In situations where scattered wavefields are small compared to the incident wavefields, we can form and solve a linear least squares problem to recover a projection (on a known subspace) of full waveform data from intensity data. We show that for sufficiently high frequencies, this projection gives a Kirchhoff image asymptotically equivalent to the Kirchhoff image obtained from full waveform data. We also generalize this imaging method to using stochastic incident fields with autocorrelation measurements. Finally, we study a mathematical model of grain growth in polycrystalline mate- rials. We review a simplified 1D grain growth model and an entropy-based theory for the evolution of an important statistic harvested from this model, the GBCD. The theory suggests the GBCD evolves according to a Fokker-Planck equation, which we validate numerically. We derive methods to estimate times from the GBCD, thus fitting it to Fokker-Planck time scales. This allows for direct comparisons of the GBCD with the Fokker-Planck solution, where we find qualitative agreement. We alsofind an energy dissipation identity which Fokker-Planck solutions must satisfy. We verify the GBCD satisfies this identity both qualitatively and quantitatively, further validating the Fokker-Planck model of GBCD evolution

A study on block flexible iterative solvers with applications to Earth imaging problem in geophysics

Les travaux de ce doctorat concernent le développement de méthodes itératives pour la résolution de systèmes linéaires creux de grande taille comportant de nombreux seconds membres. L’application visée est la résolution d’un problème inverse en géophysique visant à reconstruire la vitesse de propagation des ondes dans le sous-sol terrestre. Lorsque de nombreuses sources émettrices sont utilisées, ce problème inverse nécessite la résolution de systèmes linéaires complexes non symétriques non hermitiens comportant des milliers de seconds membres. Dans le cas tridimensionnel ces systèmes linéaires sont reconnus comme difficiles à résoudre plus particulièrement lorsque des fréquences élevées sont considérées. Le principal objectif de cette thèse est donc d’étendre les développements existants concernant les méthodes de Krylov par bloc. Nous étudions plus particulièrement les techniques de déflation dans le cas multiples seconds membres et recyclage de sous-espace dans le cas simple second membre. Des gains substantiels sont obtenus en terme de temps de calcul par rapport aux méthodes existantes sur des applications réalistes dans un environnement parallèle distribué. ABSTRACT : This PhD thesis concerns the development of flexible Krylov subspace iterative solvers for the solution of large sparse linear systems of equations with multiple right-hand sides. Our target application is the solution of the acoustic full waveform inversion problem in geophysics associated with the phenomena of wave propagation through an heterogeneous model simulating the subsurface of Earth. When multiple wave sources are being used, this problem gives raise to large sparse complex non-Hermitian and nonsymmetric linear systems with thousands of right-hand sides. Specially in the three-dimensional case and at high frequencies, this problem is known to be difficult. The purpose of this thesis is to develop a flexible block Krylov iterative method which extends and improves techniques already available in the current literature to the multiple right-hand sides scenario. We exploit the relations between each right-hand side to accelerate the convergence of the overall iterative method. We study both block deflation and single right-hand side subspace recycling techniques obtaining substantial gains in terms of computational time when compared to other strategies published in the literature, on realistic applications performed in a parallel environment

Application of Surface wave methods for seismic site characterization

Surface-wave dispersion analysis is widely used in geophysics to infer a shear wave velocity model of the subsoil for a wide variety of applications. A shear-wave velocity model is obtained from the solution of an inverse problem based on the surface wave dispersive propagation in vertically heterogeneous media. The analysis can be based either on active source measurements or on seismic noise recordings. This paper discusses the most typical choices for collection and interpretation of experimental data, providing a state of the art on the different steps involved in surface wave surveys. In particular, the different strategies for processing experimental data and to solve the inverse problem are presented, along with their advantages and disadvantages. Also, some issues related to the characteristics of passive surface wave data and their use in H/V spectral ratio technique are discussed as additional information to be used independently or in conjunction with dispersion analysis. Finally, some recommendations for the use of surface wave methods are presented, while also outlining future trends in the research of this topic

Contribution to the study of efficient iterative methods for the numerical solution of partial differential equations

Multigrid and domain decomposition methods provide efficient algorithms for the numerical solution of partial differential equations arising in the modelling of many applications in Computational Science and Engineering. This manuscript covers certain aspects of modern iterative solution methods for the solution of large-scale problems issued from the discretization of partial differential equations. More specifically, we focus on geometric multigrid methods, non-overlapping substructuring methods and flexible Krylov subspace methods with a particular emphasis on their combination. Firstly, the combination of multigrid and Krylov subspace methods is investigated on a linear partial differential equation modelling wave propagation in heterogeneous media. Secondly, we focus on non-overlapping domain decomposition methods for a specific finite element discretization known as the hp finite element, where unrefinement/refinement is allowed both by decreasing/increasing the step size h or by decreasing/increasing the polynomial degree p of the approximation on each element. Results on condition number bounds for the domain decomposition preconditioned operators are given and illustrated by numerical results on academic problems in two and three dimensions. Thirdly, we review recent advances related to a class of Krylov subspace methods allowing variable preconditioning. We examine in detail flexible Krylov subspace methods including augmentation and/or spectral deflation, where deflation aims at capturing approximate invariant subspace information. We also present flexible Krylov subspace methods for the solution of linear systems with multiple right-hand sides given simultaneously. The efficiency of the numerical methods is demonstrated on challenging applications in seismics requiring the solution of huge linear systems of equations with multiple right-hand sides on parallel distributed memory computers. Finally, we expose current and future prospectives towards the design of efficient algorithms on extreme scale machines for the solution of problems coming from the discretization of partial differential equations

Fast Iterative Solution of the Optimal Transport Problem on Graphs

In this paper, we address the numerical solution of the Optimal Transport Problem on undirected weighted graphs, taking the shortest path distance as transport cost. The optimal solution is obtained from the long-time limit of the gradient descent dynamics. Among different time stepping procedures for the discretization of this dynamics, a backward Euler time stepping scheme combined with the inexact Newton-Raphson method results in a robust and accurate approach for the solution of the Optimal Transport Problem on graphs. It is found experimentally that the algorithm requires solving between $\mathcal{O}(1)$ and $\mathcal{O}(M^{0.36})$ linear systems involving weighted Laplacian matrices, where $M$ is the number of edges. These linear systems are solved via algebraic multigrid methods, resulting in an efficient solver for the Optimal Transport Problem on graphs

The role of accurate earthquake locations in the mapping of the volcanic plumbing system at Soufrière Hills Volcano, Montserrat

Volcanic seismicity is used to monitor volcanic activity worldwide, with seismic monitoring the main method used to monitor Soufrière Hills Volcano, Montserrat. Soufrière Hills began erupting on the 18th July 1995 and has undergone five phases of activity, which were preceded by an increase in seismicity. Earthquake location is a complex problem, with several unknowns; computed hypocenters represent the optimal solution given the information available. Synthetic earthquakes allow the testing of velocity models and location methods to be compared with known earthquake locations. The current location method used at Montserrat Volcano Observatory produced large hypocenter errors, with poor constraints on hypocenters at particular depths. Comparison with other velocity models and location methods shows the Rowe et al. (2004) velocity model with NonLinLoc produces locations that more accurately represent synthetic earthquake locations. This new location method was used to relocate the seismic catalogue at Montserrat from 1995 to 2018, to understand how changes in hypocenters influence interpretations. Relocations show a migration of hypocenters towards St Georges Hill on the 12th-14th August 1995; this alongside computed focal mechanisms, suggests dyke propagation and inflation, differing from previous interpretations. Prior to a Vulcanian Explosion on the 29th July 2008, relocated hypocenters are located in SE Montserrat. The majority were located using four P phases; this has been shown to produce large hypocenter errors with synthetic testing. Therefore, earthquakes were repicked for additional P and S phases to improve locations during this period. This resulted in reduced hypocenter errors, with the majority of earthquakes relocated beneath Soufrière Hills with minimal SE locations. This project highlights the importance of using a robust location method suitable for the region to ensure that outputted hypocenters are trustworthy and accurate. Use of unsuitable methods can influence earthquake patterns and thus interpretations. This impacts understanding of volcanic systems, and ultimately hazard assessment

Dynamic modelling of post-collisional magmatism

This study addresses the question of post-collisional magmatism and its production mechanisms, addressing especially the mantle processes involved. Numerical experiments are conducted to examine the effects of viscosity weakening by subduction related water content increase in the upper mantle and the resulting sub-lithospheric small-scale convection. The models presented incorporate parameterized and thermodynamic melting models, and take into account variable relationships between mantle water content, mantle strength, water extraction by partial melting and related depletion stiffening. The results demonstrate the possible importance of so called ”hydrous activation” of the lithosphere-asthenosphere boundary: The post-collisional loss of the lithospheric mantle can be initiated and augmented by the elevated upper mantle water contents that enhances the sub-lithospheric small-scale convection, increases heat flow into the lithosphere, and produces localized lithosphere thin- ning. The irregular spatial and temporal melting patterns and the mantle melt volumes correspond to typical post-collisional mantle-derived magmatism. The small-scale convection can be localized into an edge-driven convection by significant lithosphere thickness gradients, e.g. craton edges. This helps to understand the uplift and volcanism observed in intraplate orogenic settings and implies the importance of these processes at other locations of lithosphere thickness gradients, e.g. recent collision zones. The lithospheric thinning produced by small-scale convection can initiate whole lithosphere mantle loss via positive feedback mechanisms: gradual thinning of the lithosphere causes partial melting in the lowermost crust, weakening the crust-mantle boundary and providing a detachment mechanism for the lithospheric mantle, leading to stronger lithosphere thinning and, finally, exposure of the lower crust to the hot asthenosphere. Small-scale convection and processes related to or initiated by it offer new insight and future research possibilities in studies of continental collision magmatism

Infrasound as upper atmospheric monitor

Understanding and specification of the higher altitudes of the atmosphere with global coverage over all local times is hampered by the challenges of obtaining direct measurements in the upper atmosphere. Methods to measure the properties of the atmosphere above the stratopause is an active area of scientific research. In this thesis, we revisit the use of infrasound as a passive remote sensing technique for the upper atmosphere. Signals from the Tungurahua volcano in Ecuador are used to investigate the behavior of the upper atmosphere. Depending on the atmospheric conditions, stratospheric, mesospheric and thermospheric arrivals are observed during intervals of explosive volcanic activity. It is found that the travel times and dominant frequencies of the thermospheric arrivals exhibit a coherent variability with periods equal to those of the tidal harmonics. Theoretical predictions using atmospheric specifications show that the stratospheric arrivals are predicted within 1% of the observed value. For thermospheric arrivals, this error can be as high as 10%. The error in thermospheric celerities is found to be in accord with the typical uncertainty in upper atmospheric winds. Given the observed response of the infrasound celerities to upper atmospheric tidal variability, it is suggested that infrasound observations may be used as an additional source of information to constrain the atmospheric specifications in the upper atmosphere. We present corrected wind profiles that have been obtained by minimizing misfits in traveltime and source location using a Bayesian statistics grid search algorithm. Also, a Levenberg-Marquardt search algorithm is developed. Additionally, a new numerical method has been developed to solve the problem of infrasound propagation in a stratified medium with (high Mach number) background flow, based on a modal expansion. The underlying mathematics is by no means new and has been earlier described. This solution goes beyond the effective sound speed approximation, which is typically used in infrasound propagation modeling for computational efficiency reasons. Using the wide-angle high Mach number modal solution, it is shown that traveltimes and shadow zones are under predicted using the effective sound speed approximation, with increasing grazing angle and Mach number