Search

We present a novel approach for imaging global mantle discontinuities based on full-waveform inversion (FWI). Over the past decades, extensive research has been done on imaging mantle discontinuities at approximately 400 km and 670 km depth. Accurate knowledge of their topography can put strong constraints on thermal and compositional variations and hence geodynamic modelling. So far, however, there is little consensus on their topography. We present an approach based on adjoint tomography, which has the advantage that Fréchet derivatives for discontinuities and measurements, to be inverted for, are fully consistent. Rather than working with real data, we focus on synthetic tests, where the answer is known in order to be able to evaluate the performance of the developed method. All calculations are based on the community code SPECFEM3D_GLOBE. We generate data in fixed 1-D or 3-D elastic background models of mantle velocity. Our ‘data’ to be inverted contain topography along the 400 km and 670 km mantle discontinuities. To investigate the approach, we perform several tests: (i) In a situation where we know the elastic background model 1-D or 3-D, we recover the target topography fast and accurately, (ii) The exact misfit is not of great importance here, except in terms of convergence speed, similar to a different inverse algorithm, (iii) In a situation where the background model is not known, the convergence is markedly slower, but there is reasonable convergence towards the correct target model of discontinuity topography. It has to be noted that our synthetic test is idealised and in a real data situation, the convergence to and uncertainty of the inferred model is bound to be larger. However, the use of data consistent with Fréchet kernels seems to pay off and might improve our consensus on the nature of mantle discontinuities. Our workflow could be incorporated in future FWI mantle models to adequately infer boundary interface topography., Geophysical Journal International, 226 (3), ISSN:0956-540X, ISSN:1365-246X

Probabilistic seismic hazard analysis (PSHA) is the standard method used for designing earthquake-resistant infrastructure. In recent years, several unexpected and destructive earthquakes have sparked criticism of the PSHA methodology. The seismological part of the problem is the true frequency-magnitude distribution of regional seismicity. Two major models exist, the Gutenberg-Richter (G-R) and the Characteristic Earthquake (CE) model, but it is difficult to choose between them. That is because the instrumental, historical, and paleoseimological data available are limited in many regions of interest. Here, we demonstrate how a friction experiment on aggregates of glass beads can produce both regular (CE equivalent) and irregular (G-R equivalent) stick-slip. Using a new rotary shear apparatus we produced and analyzed large catalogs of acoustic emission (AE) events related to stick-slip. The distributions of AE sizes, interevent times, and interevent distances were found to be sensitive to particle size and the applied normal stress, and, to a lesser degree, the stiffness of the loading apparatus. More importantly, the system spontaneously switched behavior for short periods of time. In the context of PSHA, if faults are able to switch behavior as our experimental system does, then justifying the choice of either the CE or the G-R model is impossible based on existing observations.

To establish a consistent framework for seismic wave propagation that accommodates the effects of stress changes, it is critical to take into account the different definitions of stress and their corresponding effects on seismic quantities (e.g., wave speeds) as dictated by continuum mechanics. Revisiting this fundamental theoretical foundation, we first emphasize the role of stress within various forms of the wave equation resulting from different choices of stress definitions. Subsequently, using this basis, we investigate connections among existing theories that describe the variation of elastic moduli as a function of changes in stress. We find that there is a direct connection between predicting stress-induced elastic changes with the well-known third-order elasticity tensor and the recently proposed adiabatic pressure derivatives of elastic moduli. However, each of these approaches has different merits and drawbacks in terms of experimental validation as well as in their use. In addition, we investigate the connection with another general approach that relies on micromechanical structures (e.g., cracks and pores). Although this can be done algebraically, it remains unclear as to which definition of stress and which corresponding constitutive relationship should be considered in practical scenarios. We support our analysis with validations using previously published benchmark experimental data.

SUMMARY The method of ScS reverberation migration is based on a ‘common reflection point’ analysis of multiple ScS reflections in the mantle transition zone (MTZ). We examine whether ray-theoretical traveltimes, slownesses and reflection points are sufficiently accurate for estimating the thickness H of the MTZ, defined by the distance between the 410- and 660-km phase transitions. First, we analyse ScS reverberations generated by 35 earthquakes and recorded at hundreds of seismic stations from the combined Arrays in China, Hi-NET in Japan and the Global Seismic Network. This analysis suggests that H varies by about 30 km and therefore that dynamic processes have modified the large-scale structure of the MTZ in eastern Asia and the western Pacific region. Second, we apply the same procedure to spectral-element synthetics for PREM and two 3-D models. One 3-D model incorporates degree-20 topography on the 410 and 660 discontinuities, otherwise preserving the PREM velocity model. The other model incorporates the degree-20 velocity heterogeneity of S20RTS and leaves the 410 and 660 flat. To optimize reflection point coverage, our synthetics were computed assuming a homogeneous grid of stations using 16 events, four of which are fictional. The resolved image using PREM synthetics resembles the PREM structure and indicates that the migration approach is correct. However, ScS reverberations are not as strongly sensitive to H as predicted ray-theoretically because the migration of synthetics for a model with degree-20 topography on the 410 and 660: H varies by less than 5 km in the resolved image but 10 km in the original model. In addition, the relatively strong influence of whole-mantle shear-velocity heterogeneity is evident from the migration of synthetics for the S20RTS velocity model and the broad sensitivity kernels of ScS reverberations at a period of 15 s. A ray-theoretical approach to modelling long-period ScS traveltimes appears inaccurate, at least for continental-scale regions with relatively sparse earthquake coverage. Additional modelling and comparisons with SS precursor and receiver function results should rely on 3-D waveform simulations for a variety of structures and ultimately the implementation of full wave theory.

Parameter estimation from reflection moveout analysis represents one of the most fundamental problems in subsurface model building. We have developed an efficient moveout inversion method based on the process of automatic flattening of common-midpoint (CMP) gathers using local slopes. We find that as a by-product of this flattening process, we can also estimate reflection traveltimes corresponding to the flattened CMP gathers. This traveltime information allows us to construct a highly overdetermined system and subsequently invert for moveout parameters including normal-moveout velocities and quartic coefficients related to anisotropy. We use the 3D generalized moveout approximation (GMA), which can accurately capture the effects of complex anisotropy on reflection traveltimes as the basis for our moveout inversion. Due to the cheap forward traveltime computations by GMA, we use a Monte Carlo inversion scheme for improved handling of the nonlinearity between the reflection traveltimes and moveout parameters. This choice also allows us to set up a probabilistic inversion workflow within a Bayesian framework, in which we can obtain the posterior probability distributions that contain valuable statistical information on estimated parameters such as uncertainty and correlations. We use synthetic and real data examples including the data from the SEAM Phase II unconventional reservoir model to demonstrate the performance of our method and discuss insights into the problem of moveout inversion gained from analyzing the posterior probability distributions. Our results suggest that the solutions to the problem of traveltime-only moveout inversion from 2D CMP gathers are relatively well constrained by the data. However, parameter estimation from 3D CMP gathers associated with more moveout parameters and complex anisotropic models are generally nonunique, and there are trade-offs among inverted parameters, especially the quartic coefficients.

Interpretation of information available from seismic data in terms of temperature and composition requires an understanding of the physical properties of minerals, in particular, the elastic properties of candidate Earth minerals at the relevant (here, lower mantle) pressure and temperature. A common practice for the bulk elastic properties is to measure volume at a range of pressures and temperatures using experiments or computational methods. These datasets are then typically fit to a pre-determined functional form, or equation of state to allow computation of elastic properties at any other pressure or temperature. However, errors, both random and systematic, limitations in the number of data and choice of pressure marker and scale, as well as different functional forms of equations of state, all contribute to the uncertainties in mineral seismic properties. In an attempt to present a more comprehensive view of these uncertainties, we use neural-network based techniques to infer the relationship among: pressure, temperature, volume, bulk modulus, and thermal expansivity of MgO. We illustrate our approach on experimental data, but an extension to ab initio data is straightforward. The type of neural network used is called a Mixture Density Network (MDN) which is a combination of a conventional feed-forward neural network and a mixture model that consists of Gaussian functions. MDNs are capable of approximating arbitrary probability density functions, which allows us to compute the uncertainties in the predicted equations of state. Since the networks interpolate locally between input samples, pressure-volume-temperature relations are implicitly learned from data without imposing any explicit thermodynamic assumptions or ad-hoc relationships. We use the partial derivatives of the mapping between inputs (pressure and temperature) and output (volume) to compute the isothermal bulk modulus and thermal expansivity. Flexibility of the MDNs allows us to investigate the uncertainty due to certain data in one region of pressure-temperature space without influencing the posterior probability density everywhere. In general, we find that the elastic properties of MgO are well-constrained by experimental data. However, our study highlights regions in which sparse or inconsistent data lead to poorly constrained elastic properties, namely: at low pressure and high temperature (1500 K), and temperatures above 2700 K. While the former conditions are likely not important for the Earth's lower mantle, they are relevant in other planetary bodies such as the Moon and Mars. Comparison with conventional equation of state forms shows that assuming a certain functional form of the pressure-volume-temperature relationship leads to potential bias in uncertainty quantification, because the uncertainties are then specific to the underlying form. In combination with data sets of other lower mantle minerals, this technique should improve uncertainty quantification in interpretations of seismic data.

When a continuum is subjected to an induced stress, the equations that govern seismic wave propagation are modified in two ways. First, the equation of conservation of linear momentum gains terms related to the induced deviatoric stress, and, second, the elastic constitutive relationship acquires terms linear in the induced stress. This continuum mechanics theory makes testable predictions with regard to stress‐induced changes in the elastic tensor. Specifically, it predicts that induced compression linearly affects the prestressed moduli with a slope determined by their local adiabatic pressure derivatives and that induced deviatoric stress produces anisotropic compressional and shear wave speeds. In this article we successfully compare such predictions against ab initio mineral physics calculations for NaCl and MgO., Key Points We compare ab initio calculations of effects of induced stress on elastic parameters with predictions based on a continuum mechanics theoryWe find that the two methods are in good agreement, without the need of higher‐order theories of elasticityThe theory currently in use for accommodating the effects on nonhydrostatic prestress on seismic wave propagation needs to be modified

SUMMARYThe correct implementation of the continuity conditions between different media is fundamental for the accuracy of any wave equation solver used in applications from seismic exploration to global seismology. Ideally, we would like to benchmark a code against an analytical Green’s function. The latter, however, is rarely available for more complex media. Here, we provide a general framework through which wave equation solvers can be benchmarked by comparing plane wave simulations to transmission/reflection (R/T) coefficients from plane-wave analysis with exact boundary conditions (BCs). We show that this works well for a large range of incidence angles, but requires a lot of computational resources to simulate the plane waves. We further show that the accuracy of a numerical Green’s function resulting from a point-source spherical-wave simulation can also be used for benchmarking. The data processing in that case is more involved than for the plane wave simulations and appears to be sufficiently accurate only below critical angles. Our approach applies to any wave equation solver, but we chose the poroelastic wave equation for illustration, mainly due to the difficulty of benchmarking poroelastic solvers, but also due to the growing interest in imaging in poroelastic media. Although we only use 2-D examples, our exact R/T approach can be extended to 3-D and various cases with different interface configurations in arbitrarily complex media, incorporating, for example, anisotropy, viscoelasticity, double porosities, partial saturation, two-phase fluids, the Biot/squirt flow and so on.

We propose a new approach for generating synthetic earthquakes based on the physics of soft glasses. The continuum approach produces yield-stress materials based on Lattice-Boltzmann simulations. We show that if the material is stimulated below yield stress, plastic events occur, which have strong similarities to seismic events. Based on a suitable definition of displacement in the continuum, we showthat the plastic events obey a Gutenberg-Richter lawwith exponents similar to those for real earthquakes. We also find that the average acceleration, energy release, stress drop and interoccurrence times scale with the same exponent. Furthermore, choosing a suitable definition for aftershocks, we show that they follow Omori's law. Finally, the far field power spectra of elastic waves generated by these plastic events decay as ω-2 similar to those observed for seismic waves. Our approach is fully self-consistent and all quantities can be calculated at all scales without the need of ad hoc friction or statistical assumptions. We herefore suggest that our approach may lead to new insights into the physics connecting the micro- and macroscales of earthquakes.

We present an accelerated full-waveform inversion based on dynamic mini-batch optimization, which naturally exploits redundancies in observed data from different sources. The method rests on the selection of quasi-random subsets (mini-batches) of sources, used to approximate the misfit and the gradient of the complete data set. The size of the mini-batch is dynamically controlled by the desired quality of the gradient approximation. Within each mini-batch, redundancy is minimized by selecting sources with the largest angular differences between their respective gradients, and spatial coverage is maximized by selecting candidate events with Mitchell’s best-candidate algorithm. Information from sources not included in a specific mini-batch is incorporated into each gradient calculation through a quasi-Newton approximation of the Hessian, and a consistent misfit measure is achieved through the inclusion of a control group of sources. By design, the dynamic mini-batch approach has several main advantages: (1) The use of mini-batches with adaptive size ensures that an optimally small number of sources is used in each iteration, thus potentially leading to significant computational savings; (2) curvature information is accumulated and exploited during the inversion, using a randomized quasi-Newton method; (3) new data can be incorporated without the need to re-invert the complete data set, thereby enabling an evolutionary mode of full-waveform inversion. We illustrate our method using synthetic and real-data inversions for upper-mantle structure beneath the African Plate. In these specific examples, the dynamic mini-batch approach requires around 20 per cent of the computational resources in order to achieve data and model misfits that are comparable to those achieved by a standard full-waveform inversion where all sources are used in each iteration., Geophysical Journal International, 221 (2), ISSN:0956-540X, ISSN:1365-246X

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., Key Points Knowing the S wave speed precisely with a loose constraint on mantle attenuation places strong constraints on the temperatureThe frequency dependence of anelasticity α is variable: ~0.1 in cold/lithospheric mantle and ~0.3 in warm/asthenospheric mantleEither partial melting or elastically accommodated grain boundary sliding is required to explain the slowest wave speeds at 130 km

The robust and automated determination of earthquake source parameters on a global and regional scale is important for many applications in seismology. We present a novel probabilistic method to invert a wide variety of (waveform) data for point‐source parameters in real time using pattern recognition. Inferences are made in the form of marginal probability density functions for point‐source parameters and incorporate realistic posterior uncertainty estimates. The neural‐network‐based method is calibrated using samples from the prior distribution, which are synthetic data vectors, and corresponding sources located in a predefined monitoring volume. Once a set of trained neural networks is available, inversions are fast with very moderate demands on computational resources: an inversion takes less than a second on a standard desktop computer. Uncertainties in the layered Earth model are taken into account in the Bayesian framework and increase the robustness of the results with respect to neglected 3D heterogeneities. Moreover, we find that the method is very robust with respect to perturbations such as observational noise and missing data and therefore is potentially well suited for automated and real‐time tasks, such as earthquake monitoring and early warning. We demonstrate the method by means of synthetic tests and by inverting an observed high‐rate Global Positioning System displacement dataset for the 2010 M w 7.2 El Mayor–Cucapah event. Our results are compatible with published point‐source estimates for this event within the respective uncertainty bounds. Online Material: Additional information on the neural network methodology and implementation details, tables on neural network parameters, crustal model and reference double‐couple solution, and figures showing prediction error, crustal models, normalized displacements, and histograms of weighted parameters.

Using spectral-element simulations of wave propagation, we investigated the sensitivity of seismic waveforms, recorded on transverse components, to upper-mantle discontinuities in 1-D and 3-D background models. These sensitivity kernels, or Fréchet derivatives, illustrate the spatial sensitivity to model parameters, of which those for shear wave speed and the surface topography of internal boundaries are discussed in this paper. We focus on the boundaries at 400 and 670 km depth of the mantle transition zone. SS precursors have frequently been used to infer the topography of upper-mantle discontinuities. These seismic phases are underside reflections off these boundaries and are usually analysed in the distance range of 110°–160°. This distance range is chosen to minimize the interference from other waves. We show sensitivity kernels for consecutive time windows at three characteristic epicentral distances within the 110°–160° range. The sensitivity kernels are computed with the adjoint method using synthetic data. From our simulations we can draw three main conclusions: (i) The exact Fréchet derivatives show that in all time windows, and also in those centred on the SS precursors, there is interference from other waves. This explains the difficulty reported in the literature to correct for 3-D shear wave speed perturbations, even if the 3-D structure is perfectly known. (ii) All studies attempting to map the topography of the 400 and 670 km discontinuities to date assume that the traveltimes of SS precursors can be linearly decomposed into a 3-D elastic structure and a topography part. We recently showed that such a linear decomposition is not possible for SS precursors, and the sensitivity kernels presented in this paper explain why. (iii) In agreement with previous work, we show that other parts of the seismograms have greater sensitivity to upper-mantle discontinuities than SS precursors, especially multiply bouncing S waves exploiting the S-wave triplications due to the mantle transition zone. These phases can potentially improve the inference of global topographic variations of the upper-mantle discontinuities in the context of full waveform inversion in a joint inversion for (an)elastic parameters and topography., Geophysical Journal International, 210 (3), ISSN:0956-540X, ISSN:1365-246X

Despite the ever increasing availability of computational power, real-time source inversions based on physical modeling of wave propagation in realistic media remain challenging. We investigate how a nonlinear Bayesian approach based on pattern recognition and synthetic 3-D Green's functions can be used to rapidly invert strong-motion data for point source parameters by means of a case study for a fault system in the Los Angeles Basin. The probabilistic inverse mapping is represented in compact form by a neural network which yields probability distributions over source parameters. It can therefore be evaluated rapidly and with very moderate CPU and memory requirements. We present a simulated real-time inversion of data for the 2008 Mw 5.4 Chino Hills event. Initial estimates of epicentral location and magnitude are available ∼14 s after origin time. The estimate can be refined as more data arrive: by ∼40 s, fault strike and source depth can also be determined with relatively high certainty.

The results of mantle convection simulations are fully determined by the input parameters and boundary conditions used. These input parameters can be for initialisation, such as initial mantle temperature, or can be constant values, such as viscosity exponents. However, knowledge of Earth-like values for many input parameters are very poorly constrained, introducing large uncertainties into the simulation of mantle flow. Convection is highly non-linear, therefore linearised inversion methods cannot be used to recover past configurations over more than very short periods of time, which makes finding both initial and constant simulation input parameters very difficult. In this paper, we demonstrate a new method for making inferences about simulation input parameters from observations of the mantle temperature field after billions of years of convection. The method is fully probabilistic. We use prior sampling to construct probability density functions for convection simulation input parameters, which are represented using neural networks. Assuming smoothness, we need relatively few samples to make inferences, making this approach much more computationally tractable than other probabilistic inversion methods. As a proof of concept, we show that our method can invert the amplitude spectra of temperature fields from 2D convection simulations, to constrain yield stress, surface reference viscosity and the initial thickness of primordial material at the CMB, for our synthetic test cases. The best constrained parameter is yield stress. The reference viscosity and initial thickness of primordial material can also be inferred reasonably well after several billion years of convection., Physics of the Earth and Planetary Interiors, 257, ISSN:0031-9201, ISSN:1872-7395

Seismic free oscillations, or normal modes, provide a convenient tool to calculate low-frequency seismograms in heterogeneous Earth models. A procedure called ‘full mode coupling’ allows the seismic response of the Earth to be computed. However, in order to be theoretically exact, such calculations must involve an infinite set of modes. In practice, only a finite subset of modes can be used, introducing an error into the seismograms. By systematically increasing the number of modes beyond the highest frequency of interest in the seismograms, we investigate the convergence of full-coupling calculations. As a rule-of-thumb, it is necessary to couple modes 1–2 mHz above the highest frequency of interest, although results depend upon the details of the Earth model. This is significantly higher than has previously been assumed. Observations of free oscillations also provide important constraints on the heterogeneous structure of the Earth. Historically, this inference problem has 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 implementations rely on approximations referred to as self- or group-coupling and do not use fully accurate seismograms. We therefore also investigate whether a systematic error might be present in currently published splitting functions. We find no evidence for any systematic bias, but published uncertainties must be doubled to properly account for the errors due to theoretical omissions and regularization in the measurement process. Correspondingly, uncertainties in results derived from splitting functions must also be increased. As is well known, density has only a weak signal in low-frequency seismograms. Our results suggest this signal is of similar scale to the true uncertainties associated with currently published splitting functions. Thus, it seems that great care must be taken in any attempt to robustly infer details of Earth's density structure using current splitting functions.

We demonstrate how effects of induced stress may be incorporated in seismic modelling and inversion. Our approach is motivated by the accommodation of pre-stress in global seismology. Induced stress modifies both the equation of motion and the constitutive relationship. The theory predicts that induced pressure linearly affects the unstressed isotropic moduli with a slope determined by their adiabatic pressure derivatives. The induced deviatoric stress produces anisotropic compressional and shear wave speeds; the latter result in shear wave splitting. For forward modelling purposes, we determine the weak form of the equation of motion under induced stress. In the context of the inverse problem, we determine induced stress sensitivity kernels, which may be used for adjoint tomography. The theory is illustrated by considering 2-D propagation of SH waves and related Frechet derivatives based on a spectral-element method.

This work presents the first comprehensive 3-D model of the crust beneath The Netherlands. To obtain this model, we designed the NARS-Netherlands project, a dense deployment of broadband stations in the area. Rayleigh and Love wave group velocity dispersion was measured from ambient noise cross-correlations. Azimuthally anisotropic group velocity maps were then constructed and the isotropic part was used to determine a shear wave speed model that includes the effects of radial anisotropy. Employing the Neighbourhood Algorithm for the depth inversion, we obtained probabilistic estimates of the radially anisotropic model parameters. We found that the variations in the thickness of the top layer largely match the transition from sediments of Permian age to those of Carboniferous age. Regions of high faulting density such as the West Netherlands Basin and Roer Valley Graben are recognized in our model by their negative radial anisotropy (V SH − V SV < 0). The model has a mid-crustal discontinuity at a depth of around 13 km and the average Moho depth is 33 km, with most of its depth variations within 2 km. Specifically, a localized Moho uplift to a depth of 29 km is found within Roer Valley Graben, in the Campine region in Belgium. Furthermore, our Rayleigh and Love wave group velocity data at periods of around 20 s show evidence for azimuthal anisotropy with a NW-SE fast direction. This anisotropy is likely related to NW-SE rock fabric in the lower crust thought to originate from the Caledonian deformation.