Advances in theoretical seismology and computational inverse problems


Advances in theoretical seismology and computational inverse problems
Co-organized by ESSI1/PS7
Convener: Christian BoehmECSECS | Co-convener: Ebru Bozdag
vPICO presentations
| Mon, 26 Apr, 13:30–15:00 (CEST)

vPICO presentations: Mon, 26 Apr

Chairpersons: Christian Boehm, Ebru Bozdag
Nicola Piana Agostinetti, Christina Dahnér-Lindkvist, and Savka Dineva

Rock elasticity in the subsurface can change in response to natural phenomena (e.g. massive precipitation, magmatic processes) and human activities (e.g. water injection in geothermal wells, ore-body exploitation). However, understanding and monitoring the evolution of physical properties of the crust is a challenging due to the limited possibility of reaching such depths and making direct measurements of the state of the rocks. Indirect measurements, like seismic tomography, can give some insights, but are generally biased by the un-even distribution (in space and time) of the information collected from seismic observations (travel-times and/or waveforms). Here we apply a Bayesian approach to overcome such limitations, so that data uncertainties and data distribution are fully accounted in the reconstruction of the posterior probability distribution of the rock elasticity  We compute a full 4D local earthquake tomography based on trans-dimensional Markov chain Monte Carlo sampling of 4D elastic models, where the resolution in space and time is fully data-driven. To test our workflow, we make use of a “controlled laboratory”: we record seismic data during one month of mining activities across a 800x700x600 m volume of Kiruna mine (LKAB, Sweden). During such period, we obtain about 260 000 P-wave and 240 000 S-wave travel-times coming from about 36000  seismic events. We operate a preliminary selection of the well-located events, using a Monte Carlo search. Arrival-times of about 19 000 best-located events (location errors less than 20m) are used as input to the tomography workflow. Preliminary results indicate that: (1) short-term (few hours) evolutions of the elastic field are mainly driven by seismic activation, i.e. the occurrence of a seismic swarm, close to the mine ore-passes. Such phenomena partially mask the effects of explosions; (2) long-term (2-3 days) evolutions of the elastic field closely match the local measurements of the stress field at a colocated stress cell. 

How to cite: Piana Agostinetti, N., Dahnér-Lindkvist, C., and Dineva, S.: Toward a full 4D seismic tomography: a case study of an active mine, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12015,, 2021.

David Harris and Douglas Dodge

Under special circumstances, waveform observations of seismic events related by a common, spatially-distributed source process exhibit a geometric architecture that is a distorted image of event distribution in the source region.  We describe a prescription for visualizing this signal space image and use a machine-learning algorithm, Isomap, and an algorithm due to Menke to invert collections of waveforms directly for relative event location.  We illustrate concepts and methods with well-characterized induced seismicity at a coal mine in the U.S. state of Utah observed by two local seismic instruments, and with synthetics.  We anticipate application of these methods to repetitive volcanic seismicity, icequakes and induced seismicity.

This is Lawrence Livermore National Laboratory contribution LLNL-ABS-818445.

How to cite: Harris, D. and Dodge, D.: The structure of signal space:  a case study of direct mapping between seismic waveforms and event distribution, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-10146,, 2021.

Ranjith Kunnath

A model that explains the anomalies in the Love wave dispersion in the earth is presented. Conventionally, welded contact between the crust and the upper mantle is assumed, leading to Love wave generation when the earth is excited. However, the observations of SH wave dispersion at seismic frequencies is at variance with this model, at least for some crustal plates (Ekström, 2011). When frictional slip occurs at the crust-upper mantle interface, a new type of interfacial elastic wave called the antiplane slip wave can occur (Ranjith, 2017). It is shown that the antiplane slip waves can explain the observed anomalies in the Love wave dispersion. 

How to cite: Kunnath, R.: Love wave or slip wave?, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-273,, 2021.

Eduardo Almeida, Hamzeh Mohammadigheymasi, Maryam Fathi, Paul Crocker, and Graça Silveira

Polarization analysis is a signal processing tool for decomposing multi-component seismic signals to a set of rectilinearly or elliptically polarized elements. Theoretically, time-frequency polarization methods are the most compatible tool to analyze the intrinsically non-stationary seismic signals. They decompose the signal to a superposition of well-defined polarized elements, localized in the time and frequency domains. However, in practice, they suffer from instability and limited resolution for discriminating between interfering seismic phases in time and frequency, as the time-frequency decomposition methods are generally an underdetermined mapping from the time to the time-frequency domain. Our contribution is threefold: Firstly we obtain the frequency-dependent polarization properties in terms of the eigenvalue decomposition of the Fourier spectra of three-components of the signal. Secondly, by extending from the frequency to the time-frequency domain and using the regularized sparsity-based time-frequency decomposition (Portniaguine and Castagna, 2004) we are able to increase resolution and reduce instability in the presence of noise. Finally, by combining directivity, rectilineary, and amplitude attributes in the time-frequency domain, we extend the time-frequency polarization analysis to extract and filter different seismic phases. By applying this method on synthetic and real seismograms we demonstrate the efficacy of the method in discriminating between the interfering seismic phases in time and frequency, including the body, Rayleigh, Love, and coda waves. This research contributes to the FCT-funded SHAZAM (Ref. PTDC/CTA-GEO/31475/2017) project.

Portniaguine, O., and J. Castagna, 2004, Inverse spectral decomposition, in SEG Technical Program Expanded Abstracts 2004: Society of Exploration Geophysicists, 1786–1789.

How to cite: Almeida, E., Mohammadigheymasi, H., Fathi, M., Crocker, P., and Silveira, G.: Eigenvalue decomposition polarization analysis: A regularized sparsity-based approach, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15267,, 2021.

Akash Kharita and Sagarika Mukhopadhyay

The surface wave phase and group velocities are estimated by dividing the epicentral distance by phase and group travel times respectively in all the available methods, this is based on the assumptions that (1) surface waves originate at the epicentre and (2) the travel time of the particular group or phase of the surface wave is equal to its arrival time to the station minus the origin time of the causative earthquake; However, both assumptions are wrong since surface waves generate at some horizontal distance away from the epicentre. We calculated the actual horizontal distance from the focus at which they generate and assessed the errors caused in the estimation of group and phase velocities by the aforementioned assumptions in a simple isotropic single layered homogeneous half space crustal model using the example of the fundamental mode Love wave. We took the receiver locations in the epicentral distance range of 100-1000 km, as used in the regional surface wave analysis, varied the source depth from 0 to 35 Km with a step size of 5 km and did the forward modelling to calculate the arrival time of Love wave phases at each receiver location. The phase and group velocities are then estimated using the above assumptions and are compared with the actual values of the velocities given by Love wave dispersion equation. We observed that the velocities are underestimated and the errors are found to be; decreasing linearly with focal depth, decreasing inversely with the epicentral distance and increasing parabolically with the time period. We also derived empirical formulas using MATLAB curve fitting toolbox that will give percentage errors for any realistic combination of epicentral distance, time period and depths of earthquake and thickness of layer in this model. The errors are found to be more than 5% for all epicentral distances lesser than 500 km, for all focal depths and time periods indicating that it is not safe to do regional surface wave analysis for epicentral distances lesser than 500 km without incurring significant errors. To the best of our knowledge, the study is first of its kind in assessing such errors.

How to cite: Kharita, A. and Mukhopadhyay, S.: Errors Introduced in Estimation of Surface Wave Phase and Group Velocities Due to Wrong Assumptions: An Assessment Using a Simple Model For Love Wave, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-779,, 2021.

Siegfried Rohdewald

We demonstrate improved resolution in P-wave velocity tomograms obtained by inversion of the synthetic SAGEEP 2011 refraction traveltime data (Zelt 2010) using Wavepath-Eikonal Traveltime Inversion (WET; Schuster 1993) and Wavelength-Dependent Velocity Smoothing (WDVS; Zelt and Chen 2016). We use a multiscale inversion approach and a Conjugate-Gradient based search method. Our default starting model is a 1D-gradient model obtained directly from the traveltime first arrivals assuming diving waves (Sheehan, 2005). As a second approach, we map the first breaks to assumed refractors and obtain a layered starting model using the Plus-Minus refraction method (Hagedoorn, 1959). We compare tomograms obtained using WDVS to smooth the current velocity model grid before forward modeling traveltimes vs. tomograms obtained without WDVS. Results show that WET images velocity layer boundaries more sharply when engaging WDVS. We determine the optimum WDVS frequency iteratively by trial-and-error. We observe that the lower the used WDVS frequency, the stronger the imaged velocity contrast at the top-of-basement. Using a WDVS frequency that is too low makes WDVS based WET inversion unstable exhibiting increasing RMS error, too high modeled velocity contrast and too shallow imaged top-of-basement. To speed up WDVS, we regard each nth node only when scanning the velocity along straight scan lines radiating from the current velocity grid node. Scanned velocities are weighted with a Cosine-Squared function as described by (Zelt and Chen, 2016). We observe that activating WDVS allows decreasing WET regularization (smoothing and damping) to a higher degree than without WDVS.


Hagedoorn, J.G., 1959, The Plus-Minus method of interpreting seismic refraction sections, Geophysical Prospecting, Volume 7, 158-182.

Rohdewald, S.R.C., 2021, SAGEEP11 data interpretation,

Schuster, G.T., Quintus-Bosz, A., 1993, Wavepath eikonal traveltime inversion: Theory. Geophysics, Volume 58, 1314-1323.

Sheehan, J.R., Doll, W.E., Mandell, W., 2005, An evaluation of methods and available software for seismic refraction tomography analysis, JEEG, Volume 10(1), 21-34.

Shewchuk, J.R., 1994, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain,

Zelt, C.A., 2010, Seismic refraction shootout: blind test of methods for obtaining velocity models from first-arrival travel times,

Zelt, C.A., Haines, S., Powers, M.H. et al. 2013, Blind Test of Methods for Obtaining 2-D Near-Surface Seismic Velocity Models from First-Arrival Traveltimes, JEEG, Volume 18(3), 183-194.

Zelt, C.A., Chen, J., 2016, Frequency-dependent traveltime tomography for near-surface seismic refraction data, Geophys. J. Int., Volume 207, 72-88.

How to cite: Rohdewald, S.: Improved interpretation of SAGEEP 2011 blind refraction data using Frequency-Dependent Traveltime Tomography, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-4214,, 2021.

Yi Zhang, Luca de Siena, and Alexey Stovas

In waveform inversion, most of the existing adjoint-state methods are based on the second-order elastic wave equations subject to displacement. The implementation of the acoustic-elastic coupling problem and free-surface in this formulation is not explicit, especially for arbitrary boundaries. The formulation of velocity-deviatoric stress-isotropic pressure can tackle the above issue. We firstly review the difference between velocity stress equations and velocity-deviatoric stress-isotropic pressure equations. Then the adjoint state of the velocity-stress equations are derived, decomposing stresses into their deviatoric and isotropic parts. To simulate the unbounded wavefield, perfectly matched layers (PML) are embedded into the system of equations. It is modified for cheap computation, which avoids PML-related memory variables by applying complex coordinate stretch to three Cartesian axes in parallel.

A 3D velocity-deviatoric stress-isotropic stress formulation is implemented with the staggered finite-difference method for several synthetic models (including anisotropic models). And inversions are then performed to reconstruct the model parameters, which is followed by a sensitivity analysis.

This method has the potential to be used with real data, both for active and passive seismics. However, in its current form, since it does not treat fluid/anisotropic solid interfaces correctly, it is limited to fluid or isotropic solid problems.

How to cite: Zhang, Y., de Siena, L., and Stovas, A.: Seismic waveform modelling and inversion with velocity-deviatoric stress-isotropic pressure formulation, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13094,, 2021.

Kassem Asfour, Roland Martin, Ludovic Bodet, Didier El Baz, Bastien Plazolles, and Javier Abreu-Torres
In order to accurately study the properties of partially saturated unconsolidated media at the near surface scale, or be able to image deeper structures through them, accurate 2D and 3D wave propagation numerical modelling tools are required. The rheology/mechanical properties of such media are frequently extremely complex (nonlinear, anisotropic, ... ), even when considered at dry state and of homogeneous mixture. Experimental observations (both at the laboratory and field scales) show that the seismic wave-field in unconsolidated granular materials remains difficult to interpret within standard methodological frameworks. We present here a numerical study aiming at exploring possible alternative forward modelling approaches to better extract information from recorded signals.
We first present a finite volume method (Asfour et al. 2021) in which exact Riemann solvers are introduced. Solutions are compared to high-order finite-differences (Seismic_CPML code) and spectral finite element (SPECFEM code) solutions. A first series of synthetic cases is shown to benchmark the code at the hundred meters scale with a 100-300Hz wavelet source content. Another synthetic and more reallistic case is then presented with a medium affected by a steep nonlinear velocity gradient in depth, typical of an unconsolidated granular medium (as previously considered  at laboratory scale). For this model, a 1500Hz dominant frequency point source wavelet is considered and fluid saturation is also tested by applying a fluid-solid coupling. First arrival times and PSV-wave dispersion obtained from the different codes are compared.
In a second step, and considering the real data recorded at the laboratory, we apply a more realistic source wavelet (obtained through signal spectrum ratio) and we perform parallelized high-order finite difference simulations (UNISOLVER code) to compare 2D and 3D elastic as well as poroelastic solutions on finely discretized meshes. Computed and observed data are compared. The poroelastic rheology provides better amplitudes in the seismograms and better exhibits some PSV modes in the phase velocity dispersion observations. Sensitivity kernels are also shown for the different rheologies. The different results obtained are now paving the way to seismic inversion at the near surface scale and to image shallow fluid/water saturated layers.

How to cite: Asfour, K., Martin, R., Bodet, L., El Baz, D., Plazolles, B., and Abreu-Torres, J.: Numerical tools to model seismic waves in unconsolidated and partially saturated granular media, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9868,, 2021.

Thomas Möller and Wolfgang Friederich

Modeling waveforms of teleseismic body waves requires the solution of the seismic wave equation in the entire Earth. Since fully-numerical 3D simulations on a global scale with periods of a few seconds are far too computationally expensive, we resort to a hybrid approach in which fully-numerical 3D simulations are performed only within the target region and wave propagation through the rest of the Earth is modeled using methods that are much faster but apply only to spherically symmetric Earth models.

We present a hybrid method that uses GEMINI to compute wave fields for a spherically symmetric Earth model up to the boundaries of a regional box. The wavefield is injected at the boundaries, where wave propagation is continued using SPECFEM-Cartesian. Inside the box, local heterogeneities in the velocity distribution are allowed, which can cause scattered and reflected waves. To prevent these waves from reflecting off the edges of the box absorbing boundary conditions are specifically applied to these parts of the wavefields. They are identified as the difference between the wavefield calculated with SPECFEM at the edges and the incident wavefield.

The hybrid method is applied to a target region in and around the Alps as a test case. The region covers an area of 1800 by 1350 km centered at 46.2°N and 10.87°E and includes crust and mantle to a depth of 600 km. We compare seismograms with a period of up to ten seconds calculated with the hybrid method to those calculated using GEMINI only for identical 1D earth models. The comparison of the seismograms shows only very small differences and thus validates the hybrid method. In addition, we demonstrate the potential of the method by calculating seismograms where the 1D velocity model inside the box is replaced by a velocity model generated using P-wave traveltime tomography.

How to cite: Möller, T. and Friederich, W.: A hybrid method to calculate teleseismic body waves in a regional 3D model using GEMINI and SPECFEM., EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14565,, 2021.

Martin van Driel, Johannes Kemper, Christian Boehm, and Amir Khan

The last decade has witnessed an unprecedented increase in high quality long period seismic data. This is because of the occurrence of several very large Earthquakes that were recorded on an exponentially growing number of broadband seismic stations that are installed in very dense networks such as the USarray. While for surface and body waves, tomography based on full waveform observables and 3D numerical simulations has become a standard tool in seismology, this is not the case for the normal modes. The main reason for this discrepancy is the fact that there is no established method to model the full physics of long period seismology in 3D and compute gradients with respect to material properties at reasonable computational cost.

Here, we present a new approach to the inclusion of the full gravitational response to spectral element wave propagation solvers. We leverage the Salvus meshing software to include the external domain using adaptive mesh refinement and high order shape mapping. Together with Neumann boundary conditions based on a multipole expansion of the right hand side this minimizes the number of additional elements needed. Initial conditions for the iterative solution of the Poisson equation based on temporal extrapolation from previous time steps together with a polynomial multigrid method reduce the number of iterations needed for convergence. In summary this approach reduces the extra cost for simulating full gravity to a similar order as the elastic forces.

We demonstrate the efficacy of the proposed method using the displacement from an elastic global wave propagation simulation at 200s period as the right hand side in the Poisson equation.

How to cite: van Driel, M., Kemper, J., Boehm, C., and Khan, A.: On the modelling of self-gravitation for full 3D global seismic wave propagation, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9439,, 2021.

Navid Hedjazian, Thomas Bodin, and Yann Capdeville

Seismic imaging techniques such as elastic full waveform inversion (FWI) have their spatial resolution limited by the maximum frequency present in the observed waveforms. Scales smaller than a fraction of the minimum wavelength cannot be resolved, only a smoothed version of the true underlying medium can be recovered. Application of FWI to media containing small and strong heterogeneities therefore remains problematic. This smooth tomographic image is related to the effective elastic properties, which can be exposed with the homogenization theory of wave propagation. We study how this theory can be used in the FWI context. The seismic imaging problem is broken down in a two-stage multiscale approach. In the first step, called homogenized full waveform inversion (HFWI), observed waveforms are inverted for a macro-scale, fully anisotropic effective medium, smooth at the scale of the shortest wavelength present in the wavefield. The solution being an effective medium, it is difficult to directly interpret it. It requires a second step, called downscaling, where the macro-scale image is used as data, and the goal is to recover micro-scale parameters. All the information contained in the waveforms is extracted in the HFWI step. The solution of the downscaling step is highly non-unique as many fine-scale models may share the same long wavelength effective properties. We therefore rely on the introduction of external a priori information. In this step, the forward theory is the homogenization itself. It is computationally cheap, allowing to consider geological models with more complexity.

In a first approach to downscaling, the ensemble of potential fine-scale models is described with an object-based parametrization, and explored with a MCMC algorithm. We illustrate the method with a synthetic cavity detection problem. In a second approach, the prior information is introduced by the means of a training image, and the fine-scale model is recovered with a multi-point statistics algorithm. We apply this method on a subsurface synthetic problem, where the goal is to recover geological facies.


How to cite: Hedjazian, N., Bodin, T., and Capdeville, Y.: Downscaling seismic tomography models, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-704,, 2021.

Maria Koroni and Andreas Fichtner

This study is a continuation of our efforts to connect adjoint methods and full-waveform inversion to common beamforming techniques, widely used and developed for signal enhancement. Our approach is focusing on seismic waves traveling in the Earth's mantle, which are phases commonly used to image internal boundaries, being however quite difficult to observe in real data. The main goal is to accentuate precursor waves arriving in well-known times before some major phase. These waves generate from interactions with global discontinuities in the mantle, thus being the most sensitive seismic phases and therefore most suitable for better understanding of discontinuity seismic structure. 

Our work is based on spectral-element wave propagation which allows us to compute exact synthetic waveforms and adjoint methods for the calculation of sensitivity kernels. These tools are the core of full-waveform inversion and by our efforts we aim to incorporate more parts of the waveform in such inversion schemes. We have shown that targeted stacking of good quality waveforms arriving from various directions highlights the weak precursor waves. It additionally makes their traveltime finite frequency sensitivity prominent. This shows that we can benefit from using these techniques and exploit rather difficult parts of the seismogram.  It was also shown that wave interference is not easily avoided, but coherent phases arriving before the main phase also stack well and show on the sensitivity kernels. This does not hamper the evaluation of waveforms, as in a misfit measurement process one can exploit more phases on the body wave parts of seismograms.

In this study, we go a step forward and present recent developments of the approach relating to the effects of noise and a real data experiment. Realistic noise is added to synthetic waveforms in order to assess the methodology in a more pragmatic scenario. The addition of noise shows that stacking of coherent seismic phases is still possible and the sensitivity kernels of their traveltimes are not largely distorted, the precursor waves contribute sufficiently to their traveltime finite-frequency sensitivity kernels.
Using a well-located seismic array, we apply the method to real data and try to examine the possibilities of using non-ideal waveforms to perform imaging of the mantle discontinuity structure on the specific areas. In order to make the most out of the dense array configuration, we try subgroups of receivers for the targeted stacking and by moving along the array we aim at creating a cluster of stacks. The main idea is to use the subgroups as single receivers and create an evaluation of seismic discontinuity structure using information from each stack belonging to a subgroup. 
Ideally, we aim at improving the tomographic images of discontinuities of selected regions by exploiting weaker seismic waves, which are nonetheless very informative.

How to cite: Koroni, M. and Fichtner, A.: Targeted stacking of weak seismic phases for improving mantle discontinuity imaging using full-waveform inversion, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14954,, 2021.

Dirk-Philip van Herwaarden, Michael Afanasiev, Sölvi Thrastarson, and Andreas Fichtner

We present a full-waveform inversion (FWI) of the African plate. Starting from the Collaborative Seismic Earth Model, we invert seismograms that are filtered to 35 s and compute gradients using the adjoint state method. This FWI uses a novel approach that we introduced earlier with the name Evolutionary FWI.

In contrast to conventional waveform inversion, our approach uses dynamically changing mini-batches (subsets of the full dataset) that approximate the gradient of the larger dataset at each iteration. This has three major advantages, (1) We exploit redundancies within the dataset, which results in a reduced computational cost for model updates, (2) The size of the complete dataset does not directly impact the computational cost of an iteration, thereby enabling us to work with larger datasets, and (3) The nature of the algorithm makes it trivial to assimilate new data, as the new data can simply be added to the complete dataset from which the mini-batches are sampled.

The aforementioned advantages enable us to extend the boundaries of what was previously possible for a given computational budget. We perform more than 80 mini-batch iterations and invert waveforms from over 400 unique earthquakes. This has the same cost as 8 iterations with all data. Our latest model clearly images tectonic features such as the Afar triple junction as well as slow zones below areas with dynamic topography, such as the Tibesti and Hoggar mountain ranges.

How to cite: van Herwaarden, D.-P., Afanasiev, M., Thrastarson, S., and Fichtner, A.: Full-waveform inversion of the African Plate, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12899,, 2021.

Solvi Thrastarson, Dirk-Philip van Herwaarden, Lion Krischer, Martin van Driel, Christian Boehm, and Andreas Fichtner

As the volume of available seismic waveform data increases, the responsibility to use the data in an effective way emerges. This requires computational efficiency as well as maximizing the exploitation of the information associated with the data.

In this contribution, we present a long-wavelength Earth model, created by using the data recorded from over a thousand earthquakes, starting from a simple one-dimensional background (PREM). The model is constructed with an accelerated full-waveform inversion (FWI) method which can seamlessly include large data volumes with a significantly reduced computational overhead. Although we present a long-wavelength model, the approach has the potential to go to much higher frequencies, while maintaining a reasonable cost.

Our approach combines two novel FWI variants. (1) The dynamic mini-batch approach which uses adaptively defined subsets of the full dataset in each iteration, detaching the direct scaling of inversion cost from the number of earthquakes included. (2) Wavefield-adapted meshes which utilize the azimuthal smoothness of the wavefield to design meshes optimized for each individual source location. Using wavefield adapted meshes can drastically reduce the cost of both forward and adjoint simulations as well as it makes the scaling of the computing cost to modelled frequencies more favourable.

How to cite: Thrastarson, S., van Herwaarden, D.-P., Krischer, L., van Driel, M., Boehm, C., and Fichtner, A.: Long-Wavelength Earth Model via Accelerated Full-Waveform Inversion, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-12687,, 2021.

Andrew Curtis, Xuebin Zhao, and Xin Zhang

The ultimate goal of a geophysical investigation is usually to find answers to scientific (often low-dimensional) questions: how large is a subsurface body? How deeply does lithosphere subduct? Does a certain subsurface feature exist? Background research reviews existing information, an experiment is designed and performed to acquire new data, and the most likely answer is estimated. Typically the answer is interpreted from geophysical inversions, but is usually biased because only one particular forward function (model-data relationship) is considered, one inversion method is used, and because human interpretation is a biased process. Interrogation theory provides a systematic way to answer specific questions. Answers balance information from multiple forward models, inverse methods and model parametrizations probabilistically, and optimal answers are found using decision theory.

Two examples illustrate interrogation of the Earth’s subsurface. In a synthetic test, we estimate the cross-sectional area of a subsurface low velocity anomaly by interrogating Bayesian probabilistic tomographic maps. By combining the results of four different nonlinear inversion algorithms, the optimal answer is very close to the true answer. In a field data application, we evaluate the extent of the Irish Sea Sedimentary Basin based on the uncertainties in velocity structure derived from Love wave tomography. This example shows that the computational expense of estimating uncertainties adds explicit value to answers.

This study demonstrates that interrogation theory answers realistic questions about the Earth’s subsurface. The same theory can be used to solve different types of scientific problem - experimental design, interpreting models, expert elicitation and risk estimation - and can be applied in any field of science. One of its most important contributions is to show that fully nonlinear estimates of uncertainty are critical for decision-making in real-world geoscientific problems, potentially justifying their computational expense.


How to cite: Curtis, A., Zhao, X., and Zhang, X.: Interrogating Tomographic Uncertainties for Subsurface Structural Information, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-1334,, 2021.

Andrea Zunino, Klaus Mosegaard, Christian Boehm, Lars Gebraad, and Andreas Fichtner

The Hamiltonian Monte Carlo method (HMC) is gaining popularity in the geophysical community to fully address nonlinear inverse problems and related uncertainty quantification. We present here an application of HMC to invert seismic data in the acoustic approximation in the context of reflection seismology. We address a 2-D problem, in the form of a vertical cross section where both source and receivers are located near the surface of the model. To solve the forward problem we utilise the finite-difference method with PML absorbing boundary conditions. The observed data are represented by a set of shotgathers.

The crucial aspect for a successful application of the HMC lies in the capability of performing gradient computations in an efficient manner. To this end, we use the adjont state method to compute the gradient of the misfit functional, which has a computational cost of only about twice that of the forward computation, a very efficient strategy. From the collection of samples characterising the posterior distribution obtained with the HMC, we can derive quantities of interest using statistical analysis and assess uncertainties.

We illustrate an application of this methodology on a synthetic test mimicking the setup encountered in exploration problems.

How to cite: Zunino, A., Mosegaard, K., Boehm, C., Gebraad, L., and Fichtner, A.: Hamiltonian Monte Carlo inversion of seismic reflection data in the acoustic approximation, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9615,, 2021.

Xin Zhang and Andrew Curtis

Seismic full waveform inversion (FWI) produces high resolution images of the subsurface directly from seismic waveforms, and has been applied at global, regional and industrial spatial scales. FWI is traditionally solved using optimization methods, by iteratively updating a model of the subsurface so as to minimize the misfit between observed waveforms and those predicted by linearised (approximate) physics. Due to the nonlinearity of the physical relation between model parameters and seismic waveforms, a good starting model (derived from other methods) and hence strong prior information is required in order that the linearised physical relationships are reasonable in the vicinity of the true solution. Such linearised methods cannot provide accurate estimates of uncertainties, which are required if we are to understand and interpret the results appropriately.

To estimate uncertainties more accurately, nonlinear Bayesian methods have been deployed to solve the FWI problem. Monte Carlo sampling is one such algorithm but it is computationally expensive, and all Markov chain Monte Carlo-based methods are difficult to parallelise fully. Variational inference provides an efficient, fully parallelisable alternative methodology. This is a class of methods that optimize an approximation to a probability distribution describing post-inversion parameter uncertainties. Both Monte Carlo and variational full waveform inversion (VFWI) have been applied previously to solve FWI problems, but only with strong prior information about the velocity structure to limit the space of possible models. Unfortunately such strong information is almost never available in practice. In addition, VFWI has only been applied to wavefield transmission problems in which seismic data are recorded on a receiver array that lies above the structure to be imaged, given known, double-couple (earthquake-like) sources located underneath the same structure. In practice, knowledge of such sources is never definitive, and usually depends circularly on the unknown structure itself. In this study, we present the first application of VFWI to seismic reflection data acquired from known (deliberately fired) near-surface sources. We also apply variational inference (specifically, Stein variational gradient descent) to solve FWI problems with realistic prior information. We perform multiple inversions using data from different frequency ranges, and show that VFWI produces high resolution mean and uncertainty models using both low and high frequency data, and given only weak prior information. This is usually impossible using traditional linearised methods. We conclude that VFWI may be a useful method to produce high resolution images and reliable uncertainties, at least in smaller FWI problems.

How to cite: Zhang, X. and Curtis, A.: Bayesian Full-waveform Inversion with Weak Priors, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-1330,, 2021.

Xuebin Zhao, Andrew Curtis, and Xin Zhang

Seismic travel time tomography is used widely to image the Earth's interior structure and to infer subsurface properties. Tomography is an inverse problem, and computationally expensive nonlinear inverse methods are often deployed in order to understand uncertainties in the tomographic results. Monte Carlo sampling methods estimate the posterior probability distribution which describes the solution to Bayesian tomographic problems, but they are computationally expensive and often intractable for high dimensional model spaces and large data sets due to the curse of dimensionality. We therefore introduce a new method of variational inference to solve Bayesian seismic tomography problems using optimization methods, while still providing fully nonlinear, probabilistic results. The new method, known as normalizing flows, warps a simple and known distribution (for example a Uniform or Gaussian distribution) into an optimal approximation to the posterior distribution through a chain of invertible transforms. These transforms are selected from a library of suitable functions, some of which invoke neural networks internally. We test the method using both synthetic and field data. The results show that normalizing flows can produce similar mean and uncertainty maps to those obtained from both Monte Carlo and another variational method (Stein varational gradient descent), at significantly decreased computational cost. In our tomographic tests, normalizing flows improves both accuracy and efficiency, producing maps of UK surface wave speeds and their uncertainties at the finest resolution and the lowest computational cost to-date, allowing results to be interrogated efficiently and quantitatively for subsurface structure.

How to cite: Zhao, X., Curtis, A., and Zhang, X.: Bayesian Variational Seismic Tomography using Normalizing Flows, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-1455,, 2021.

Lars Gebraad, Sölvi Thrastarson, Andrea Zunino, and Andreas Fichtner

Uncertainty quantification is an essential part of many studies in Earth science. It allows us, for example, to assess the quality of tomographic reconstructions, quantify hypotheses and make physics-based risk assessments. In recent years there has been a surge in applications of uncertainty quantification in seismological inverse problems. This is mainly due to increasing computational power and the ‘discovery’ of optimal use cases for many algorithms (e.g., gradient-based Markov Chain Monte Carlo (MCMC). Performing Bayesian inference using these methods allows seismologists to perform advanced uncertainty quantification. However, oftentimes, Bayesian inference is still prohibitively expensive due to large parameter spaces and computationally expensive physics.

Simultaneously, machine learning has found its way into parameter estimation in geosciences. Recent works show that machine learning both allows one to accelerate repetitive inferences [e.g. Shahraeeni & Curtis 2011, Cao et al. 2020] as well as speed up single-instance Monte Carlo algorithms using surrogate networks [Aleardi 2020]. These advances allow seismologists to use machine learning as a tool to bring accurate inference on the subsurface to scale.

In this work, we propose the novel inclusion of adjoint modelling in machine learning accelerated inverse problems. The aforementioned references train machine learning models on observations of the misfit function. This is done with the aim of creating surrogate but accelerated models for the misfit computations, which in turn allows one to compute this function and its gradients much faster. This approach ignores that many physical models have an adjoint state, allowing one to compute gradients using only one additional simulation.

The inclusion of this information within gradient-based sampling creates performance gains in both training the surrogate and the sampling of the true posterior. We show how machine learning models that approximate misfits and gradients specifically trained using adjoint methods accelerate various types of inversions and bring Bayesian inference to scale. Practically, the proposed method simply allows us to utilize information from previous MCMC samples in the algorithm proposal step.

The application of the proposed machinery is in settings where models are extensively and repetitively run. Markov chain Monte Carlo algorithms, which may require millions of evaluations of the forward modelling equations, can be accelerated by off-loading these simulations to neural nets. This approach is also promising for tomographic monitoring, where experiments are repeatedly performed. Lastly, the efficiently trained neural nets can be used to learn a likelihood for a given dataset, to which subsequently different priors can be efficiently applied. We show examples of all these use cases.


Lars Gebraad, Christian Boehm and Andreas Fichtner, 2020: Bayesian Elastic Full‐Waveform Inversion Using Hamiltonian Monte Carlo.

Ruikun Cao, Stephanie Earp, Sjoerd A. L. de Ridder, Andrew Curtis, and Erica Galetti, 2020: Near-real-time near-surface 3D seismic velocity and uncertainty models by wavefield gradiometry and neural network inversion of ambient seismic noise.

Mohammad S. Shahraeeni and Andrew Curtis, 2011: Fast probabilistic nonlinear petrophysical inversion.

Mattia Aleardi, 2020: Combining discrete cosine transform and convolutional neural networks to speed up the Hamiltonian Monte Carlo inversion of pre‐stack seismic data.

How to cite: Gebraad, L., Thrastarson, S., Zunino, A., and Fichtner, A.: Accelerating inverse problems in seismology using adjoint-based machine learning, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16089,, 2021.

Jonathan Smith, Zachary Ross, Kamyar Azizzadenesheli, and Jack Muir

High resolution earthquake hypocentral locations are of critical importance for understanding the regional context driving seismicity. We introduce a scheme to reliably approximate a hypocenter posterior in a continuous domain that relies on recent advances in deep learning.

Our method relies on a differentiable forward model in the form of a deep neural network, which is trained to solve the Eikonal equation (EikoNet). EikoNet can rapidly determine the travel-time between any source-receiver pair for a non-gridded solution. We demonstrate the robustness of these travel-time solutions are for a series of complex velocity models.

For the inverse problem, we utilize Stein Variational Inference, which is a recent approximate inference procedure that iteratively updates a configuration of particles to approximate a target posterior by minimizing the so-called Stein discrepancy. The gradients of this objective function can be rapidly calculated due to the differentiability of the EikoNet. The particle locations are updated until convergence, after which we utilize clustering techniques and kernel density methods to determine the optimal hypocenter and its uncertainty.

The inversion procedure outlined in this work is validated using a series of synthetic tests to determine the parameter optimisation and the validity for large observational datasets, which can locate earthquakes in 439s per event for 2039 observations. In addition, we apply this technique to a case study of seismicity in the Southern California region for earthquakes from 2019.

How to cite: Smith, J., Ross, Z., Azizzadenesheli, K., and Muir, J.: HypoSVI: Hypocenter inversion with Stein Variational Inference and Physics Informed Neural Networks, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-3371,, 2021.