Search

Automatic differentiation (AD) is a valuable computing technique that can automatically calculate the derivative of a function. Using the chain rule and algebraic manipulations, AD can save significant human effort by eliminating the need for theoretical derivations, coding, and debugging. Most importantly, it guarantees accurate derivatives, making it a popular choice for many non-linear optimization problems. However, its use in the geophysical inversion has been limited due to difficulties in differentiating the linear-equations solution, which cannot be explicitly defined as an elementary function. To address this issue, we employ an improved AD scheme using implicit differentiation (ADID) that creates a new AD operator that customizes the standard AD scheme to function more efficiently. We demonstrate the effectiveness and validity of ADID using a toy example and compare it with the widely used adjoint equation (AE) approach in a synthetic 2D magnetotelluric (MT) problem. ADID is highly versatile and compatible and can be easily implemented for similar geophysical problems. Finally, we show how ADID can be integrated into 3D MT and 3D direct current resistivity (DC) inversions. • To guarantee the accuracy of the sensitivity and gradient for the geophysical inversion using automatic differentiation. • To extend the functional capabilities of AD in the geophysical inversion using implicit differentiation. • To reduce human effort required for theoretical derivations, coding, and debugging. [ABSTRACT FROM AUTHOR]

Geophysical methods are widely used to gather information about the subsurface as they are non-intrusive and comparably cheap. However, the solution to the geophysical inverse problem is inherently non-unique, which introduces considerable uncertainties. As a partial remedy to this problem, independently acquired geophysical data sets can be jointly inverted to reduce ambiguities in the resulting multi-physical subsurface images. A novel cooperative inversion approach with joint minimum entropy constraints is used to create more consistent multi-physical images with sharper boundaries with respect to the single-method inversions. Here, this approach is implemented in an open-source software and its applicability on electrical resistivity tomography (ERT), seismic refraction tomography (SRT), and magnetic data is investigated. A synthetic 2D ERT and SRT data study is used to demonstrate the approach and to investigate the influence of the governing parameters. The findings showcase the advantage of the joint minimum entropy (JME) stabilizer over separate, conventional smoothness-constrained inversions. The method is then used to analyze field data from Rockeskyller Kopf, Germany. 3D ERT and magnetic data are combined and the results confirm the expected volcanic diatreme structure with improved details. The multi-physical images of both methods are consistent in some regions, as similar boundaries are produced in the resulting models. Because of its sensitivity to hydrologic conditions in the subsurface, observations suggest that the ERT method senses different structures than the magnetic method. These structures in the ERT result do not seem to be enforced on the magnetic susceptibility distribution, showcasing the flexibility of the approach. Both investigations outline the importance of a suitable parameter and reference model selection for the performance of the approach and suggest careful parameter tests prior to the joint inversion. With proper settings, the JME inversion is a promising tool for geophysical imaging, however, this work also identifies some objectives for future studies and additional research to explore and optimize the method. • Detailed review of the theory of joint minimum entropy constraints with additional considerations. • Detailed demonstration of how (joint) minimum entropy constraints influence the (joint) inversion of geophysical data. • Demonstration of approach on joint ERT-SRT inversion including a sensitivity study of important parameters. • Joint minimum entropy constrained inversion of magnetic and electrical field data to image a volcanic structure. • Discussion of factors that have significant influence on success of joint minimum entropy inversion. [ABSTRACT FROM AUTHOR]

Abstract Cross-hole radar is an effective method for geophysical exploration of underground geological conditions. Besides cross-hole radar, some geological information can be obtained by other methods, such as preliminary geological survey, geologic drilling and long-distance geophysical prospecting. This information can work as priori information to help interpret cross-hole radar detection data. Full waveform inversion (FWI) is an advanced inversion method for cross-hole radar, having higher resolution than the ray tomography method. However, traditional FWI for cross-hole radar takes less advantage of the prior information. Therefore, the introduction of the priori information into FWI is worthy of in-depth study to improve inversion quality. Here we use a special penalty function to impose inequality constraint on the objective function, and then calculate a new update step length value different from the traditional FWI. We use this method to introduce the priori information and improve the accuracy and stability of FWI. The reasonable value range of weight coefficients is discussed, which can control the strength of inequality constraint. Several inversion examples are designed to validate the effect of this improved FWI method. The inversion examples for synthetic and field data show that FWI based on inequality constraint for cross-hole radar has better inversion effect than the traditional FWI. To a certain extent, the introduction of inequality constraint can constrain inversion process, improve the inversion stability and provide higher precision inversion results. Highlights • Introduce priori information into full waveform inversion to improve its accuracy and stability. • Apply the inequality constraint representing the parameter variation range into full waveform inversion. • The inversion examples for synthetic and field data has good inversion results. [ABSTRACT FROM AUTHOR]

Abstract The main objective of the full waveform inversion (FWI) is to deliver a velocity model which produces accurate and high resolution depth images. An integral part of FWI is the iterative optimization process aiming at minimizing an objective function describing the misfit between the observed and modelled data with respect to the velocity model parameters. However, FWI commonly suffers from local minima due to errors in data fitting and the starting velocity model. In this study, we develop a methodology to enhance convergence and mitigate the risk of having solutions falling into local minima. We employ optical flow to decompose the source and residual receiver wavefields and map them into the subsurface angle domain. Having evaluated the propagation vectors of the respective cleaner wavefields, subsurface reflection and azimuth angles can be robustly computed. Using such information we can apply angle-dependent preconditioning within the inversion algorithm. During a conventional FWI workflow, larger angles (i.e., transmission wavepaths) are mainly responsible for long wavelength velocity anomalies whereas narrow angles (i.e., reflection wavepaths) influence short wavelength velocity anomalies. In the proposed FWI methodology, we apply an optimized weighting function based on these estimated subsurface reflection angles in order to weight certain angle directions needed for gradient estimation using in parameter updating. Following a hierarchical approach, allowing for larger angles at the early stage of the FWI application, we initially solve for the smooth and longer wavelength velocity updates. Relaxing the opening angle from large to mid to narrow, the gradient direction will first estimate long wavelength velocity model updates (transmission-based wavepaths), following by short wavelength velocity model updates (reflection-based wavepaths). Moreover, angle-dependent-based gradient preconditioning can reduce cycle skips and improve convergence during FWI iterations. Highlights • Decomposing wavefields by optical flow during full waveform inversion is proposed. • A gradient direction is pre-conditioned based on the subsurface reflection angle. • Proposed pre-conditioning could decompose the wavenumber components. • Pre-conditioning provides smoother and less noisy updates in velocity model building. [ABSTRACT FROM AUTHOR]

Inversion of receiver functions is commonly used to recover the S-wave velocity structure beneath seismic stations. Traditional approaches are based on deconvolved waveforms, where the horizontal component of P-wave seismograms is deconvolved by the vertical component. Deconvolution of noisy seismograms is a numerically unstable process that needs to be stabilized by regularization parameters. This biases noise statistics, making it difficult to estimate uncertainties in observed receiver functions for Bayesian inference. This study proposes a method to directly invert observed radial waveforms and to better account for data noise in a Bayesian formulation. We illustrate its feasibility with two synthetic tests having different types of noises added to seismograms. Then, a real site application is performed to obtain the 1-D S-wave velocity structure beneath a seismic station located in the Tengchong volcanic area, Southwestern China. Surface wave dispersion measurements spanning periods from 8 to 65 s are jointly inverted with P waveforms. The results show a complex S-wave velocity structure, as two low velocity zones are observed in the crust and uppermost mantle, suggesting the existence of magma chambers, or zones of partial melt. The upper magma chambers may be the heart source that cause the thermal activity on the surface. [ABSTRACT FROM AUTHOR]

This study is related to the application of sparse signal reconstruction in estimation of Rayleigh-wave modal attenuation coefficients from near-field multichannel seismic data. A two-step approach was proposed. At the first step, the Rayleigh-wave modes are resolved and their phase velocities are estimated by a reconstruction of signal that is represented to be sparse in terms of wave velocity. At the second step, this time the signal is represented to be sparse in terms of attenuation parameter and then recovered optimally from the measurement equation that includes the modal phase velocities estimated above. The attenuation values corresponding to the dominant peak components of the estimated sparse signal vector are chosen as modal attenuation coefficients. Results from synthetic and real-field data show the high performance of proposed approach overcoming noise and the adverse near-field effects. [ABSTRACT FROM AUTHOR]

Macro-velocity models are important for imaging the subsurface structure. However, the conventional objective functions of full waveform inversion in the time and the frequency domain have a limited ability to recover the macro-velocity model because of the absence of low-frequency information. In this study, we propose new objective functions that can recover the macro-velocity model by minimizing the difference between the zero-frequency components of the square of seismic traces. Instead of the seismic trace itself, we use the square of the trace, which contains low-frequency information. We apply several time windows to the trace and obtain zero-frequency information of the squared trace for each time window. The shape of the new objective functions shows that they are suitable for local optimization methods. Since we use the acoustic wave equation in this study, this method can be used for deep-sea marine data, in which elastic effects can be ignored. We show that the zero-frequency components of the square of the seismic traces can be used to recover macro-velocities from synthetic and field data. [ABSTRACT FROM AUTHOR]

Due to the complexity of underground reservoir distribution and the non-uniqueness of geophysical inverse problems, there is still a lack of practical and effective cross-well electromagnetic inversion methods. Our goal is to develop an efficient method to reduce the non-uniqueness of the physical property model recovered in the inversion. Based on a defined framework, we propose a 3-D crosswell electromagnetic sparse inversion based on iterative reweighted least squares (IRLS) method. We perform the inversion using a hybrid norm formulation, where the l p -norm (0 ≤ p ≤ 1) sparse optimization problem can be transformed into a series of weighted l 2 -norm optimization problems. Compared with the traditional norm inversion method, the sparse inversion method can make more effective use of the known physical information and obtain more accurate abnormal body position and shape. Finally, the effectiveness of sparse inversion method is verified by model test and inversion of measured data in mining area. • Through different norm inversion, we can select and obtain the valuable information we need. • The inversion and recovery models of this method have distinct boundary characteristics. • It is of great significance for the effective simulation of complex 3D underground geological body by electromagnetic technique between wells. • Our method is applied to field data examples to obtain accurate information about subsurface anomalies. [ABSTRACT FROM AUTHOR]

We have developed a data-space multiple cross-gradient joint inversion algorithm, and validated it through synthetic tests and applied it to magnetotelluric (MT), gravity and magnetic datasets acquired along a 95 km profile in Benxi-Ji'an area of northeastern China. To begin, we discuss a generalized cross-gradient joint inversion for multiple datasets and model parameters sets, and formulate it in data space. The Lagrange multiplier required for the structural coupling in the data-space method is determined using an iterative solver to avoid calculation of the inverse matrix in solving the large system of equations. Next, using model-space and data-space methods, we inverted the synthetic data and field data. Based on our result, the joint inversion in data-space not only delineates geological bodies more clearly than the separate inversion, but also yields nearly equal results with the one in model-space while consuming much less memory. [ABSTRACT FROM AUTHOR]

In urban and industrialised areas, magnetotelluric (MT) recordings often contain a significant amount of high-power near-field noise which, constituting non-Gaussian noise for most of the recording length, may lead to the failure of the common transfer functions (TF) estimation methods. The most evident consequence of this is the appearance of sharp points in the TF curves versus frequency, deteriorating their inherent smoothness. The TF's smoothness is in fact one of the main criteria adopted by MT practitioners to evaluate the data quality and effectiveness of the estimation procedure; when not reached by the estimation, it may be a-posteriori introduced by splines or filtering procedures, but these approaches lack physical soundness. We propose a new heuristic algorithm aimed at obtaining the maximum TF versus frequency smoothness, through an optimization scheme rejecting EM events in frequency-domain. The algorithm searches for frequency-dependent power thresholds to be applied to the events, rejecting those making the TF sharp. The TF's smoothness is the objective function to be minimized, and the model space is constituted by the infinite set of threshold vectors for each analysed frequency. After the process, the Ordinaly Least Squares (OLS) residual distribution is much closer to a Rayleigh one, allowing for an unbiased estimation. We applied the procedure both to synthetic and real MT data. Measured data has been collected for a geothermal exploration project in the East Africa Rift, where strong coherent noise affected the recordings. The combination of the technique with the remote-reference technique effectively reduced the bias due to the artificial signal. The physical consistency of the resulting TF has been tested by dispersion relationships. • Artificial noise hampers magnetotellurics recordings, and the estimated impedances can be severely biased. • We propose a new estimation algorithm, allowing for smooth and unbiased impedances in case of persistent noise. • The algorithm implements an inverse approach for automatic rejection of noisy portions of the MT data. • We successfully applied the method to synthetic data and real ones from geothermal exploration in East Africa. [ABSTRACT FROM AUTHOR]

Amplitude variation with offset (AVO) inversion, particularly for more than two model parameters, is a highly ill-posed problem and, hence, regularization is indispensable. For common regularization such as norm-regularization, learning can be forged ahead via gradient descent. Alternatively, learning can be achieved by performing gradient descent on the objective and projects the result back into the space of the regularizer; this allows the use of atypical regularization functions. Here, we propose a regularized inverse problem to mitigate the Ill-posedness of the amplitude inversion. The regularization is added to measure the difference in information between the a priori probability density function and the predicted probability density of the inverted parameters, and therefore enhances the consistency between the inversion results and the a priori geological information. Information theory provides a collection of contrast functions which quantify the divergence from one probability distribution to another, such as the relative entropy. The a priori density is approximated by a Gaussian mixture model, obtained from well logs and rock physics data. The mixture model is a density estimator, providing the statistical properties of the model parameters of interest. The likelihood of the data and the divergence are combined in an augmented Lagrangian scheme, to obtain a stable solution that best generate the recorded seismic data and satisfy the geological constraints conveyed by the a priori probability density function. The proposed inversion scheme is applied to the anisotropy AVO inversion, for estimating the elastic and seismic anisotropy parameters of shale formations. Compared to the unconstrained minimization, the P- and S-wave velocity, and ε are better recovered, moreover, density and Thomsen's δ are well-constrained. • Develop a deterministic AVO inversion under a Bayesian framework. • To promote the use of high-dimensional AVO models. • To alleviate ill-posedness of the seismic amplitude inversion. [ABSTRACT FROM AUTHOR]

Wavelet methods have been used in potential fields study to estimate source properties such as depth or structural index, through the analysis of Wavelet Transform Modulus Maxima Lines (WTMML) intersections and slopes at high scales. Little has been done on the study of maximum points of the wavelet diagram, that we call here Maximum Wavelet Coefficient Scales (MWCS). Previous works have shown interesting correlations between MWCS and source depths, depending on the wavelet used in regards to the source nature and the data derivative order. In this paper, we introduce an empirical law involving spectral parameters that have not been studied so far, which allows analytical calculation of the MWCS, knowing the source characteristics and using certain wavelets. In return, the study of MWCS allows recovering source characteristics from the use of a single wavelet, without prior knowledge on the source. We demonstrate through synthetic models that the new capability of predicting the source type and depth according to the wavelet coefficient behaviour allows new ways of potential fields' sources characterization and identification. We show an application of the formula on a real case example in the Uinta Mountains (Utah, USA). [ABSTRACT FROM AUTHOR]

We present a new method to estimate subsurface geomodels using a multi-objective stochastic search technique that allows a variety of direct and indirect measurements to simultaneously constrain the earth model. Inherent uncertainties and noise in real data measurements may result in conflicting geological and geophysical datasets for a given area; a realistic earth model can then only be produced by combining the datasets in a defined optimal manner. One approach to solving this problem is by joint inversion of the various geological and/or geophysical datasets, and estimating an optimal model by optimizing a weighted linear combination of several separate objective functions which compare simulated and observed datasets. In the present work, we consider the joint inversion of multiple datasets for geomodel estimation, as a multi-objective optimization problem in which separate objective functions for each subset of the observed data are defined, followed by an unweighted simultaneous stochastic optimization to find the set of best compromise model solutions that fits the defined objectives, along the so-called “Pareto front”. We demonstrate that geostatistically constrained initializations of the algorithm improves convergence speed and produces superior geomodel solutions. We apply our method to a 3D reservoir lithofacies model estimation problem which is constrained by a set of geological and geophysical data measurements and attributes, and assess the sensitivity of the resulting geomodels to changes in the parameters of the stochastic optimization algorithm and the presence of realistic seismic noise conditions. [ABSTRACT FROM AUTHOR]

We present a spectral-element implementation of full seismic waveform inversion for large heterogeneous HPC systems. In this we address the optimal parallelisation configurations of individual simulations, the large I/O requirements of adjoint simulations, and the scheduling of large numbers of forward and adjoint solves, typical for realistic inversions. Using GPU accelerators allows us to achieve a 3.5–4 times performance improvement over the best known homogeneous implementation. We achieve GPU memory throughput varying from 60 to 80% of the peak bandwidth, thus providing good utilisation of hardware resources. We demonstrate the practical applicability of our developments in a real-data application in the western Mediterranean. With the help of GPU accelerators, we are able to model and invert seismic wave propagation in a frequency band that is broader than permitted by the use of CPUs alone, which helps bridging the traditional gap between crustal and mantle tomography. [ABSTRACT FROM AUTHOR]

A solution to the increasing greenhouse gas emissions is carbon capture and storage. At the Sleipner field, offshore Norway, carbon dioxide is separated from the produced gas and injected into the Utsira formation at approximately 1000 m depth. The Utsira formation is a high porosity saline aquifer that is 200–300 m thick in the injection area. In 1994, two years before the injection started, a seismic survey was acquired. The survey has been followed up by several seismic surveys acquired to monitor the migration of the injected gas. The full waveform inversion (FWI) method is a non-linear inverse method for estimating subsurface elastic parameters. Time-lapse FWI has recently been suggested as a monitoring tool for resolving time-lapse changes directly in the elastic parameter models. This paper presents an application of a three dimensional implementation of isotropic elastic time-lapse FWI applied on time-lapse field data from the Sleipner area. The baseline dataset is the 1994 dataset, whereas the monitor dataset is the survey acquired in 2006 after ten years of injection. Time-lapse FWI is used to estimate separate elastic parameter models for both the baseline and the monitor datasets. The results show that FWI is able to estimate detailed elastic parameter models that produce synthetic data that match the field data. The inverted models are used in a pre-stack depth migration method to yield detailed seismic images of the Sleipner before and after the injection of the gas. The gas cloud after ten years of injection is clearly visible on the monitor seismic images. The baseline and monitor seismic images show events and discontinuities that can explain the migration pathways of the injected gas through the Utsira formation. [ABSTRACT FROM AUTHOR]

Full waveform inversion is currently considered as a promising seismic imaging method to obtain high-resolution and quantitative images of the subsurface. It is a nonlinear ill-posed inverse problem, the main difficulty of which that prevents the full waveform inversion from widespread applying to real data is the sensitivity to incorrect initial models and noisy data. Local optimization theories including Newton's method and gradient method always lead the convergence to local minima, while global optimization algorithms such as simulated annealing are computationally costly. To confront this issue, in this paper we investigate the possibility of applying the trust-region method to the full waveform inversion problem. Different from line search methods, trust-region methods force the new trial step within a certain neighborhood of the current iterate point. Theoretically, the trust-region methods are reliable and robust, and they have very strong convergence properties. The capability of this inversion technique is tested with the synthetic Marmousi velocity model and the SEG/EAGE Salt model. Numerical examples demonstrate that the adaptive subspace trust-region method can provide solutions closer to the global minima compared to the conventional Approximate Hessian approach and the L-BFGS method with a higher convergence rate. In addition, the match between the inverted model and the true model is still excellent even when the initial model deviates far from the true model. Inversion results with noisy data also exhibit the remarkable capability of the adaptive subspace trust-region method for low signal-to-noise data inversions. Promising numerical results suggest this adaptive subspace trust-region method is suitable for full waveform inversion, as it has stronger convergence and higher convergence rate. [ABSTRACT FROM AUTHOR]

We perform a Bayesian inversion of degree-zero spheroidal mode splitting function measurements for radial (1-D) Earth structure, in terms of the Voigt averages of P-wave ( V P ) and S-wave ( V S ) velocities, density, bulk and shear attenuation, using neural networks. The method is flexible and allows us to assess the robustness of features in existing reference models, such as PREM . The Bayesian framework provides a means for quantifying uncertainties in the model parameters and for measuring the information content of the data. The analysis of the information content suggests that the free oscillations constrain most parameters better than body wave travel time data. Our most important findings can be summarised as follows. The data prefer an inner–outer core boundary (ICB) that lies in the depth range 5154.7–5165.7 km, i.e. deeper than in existing reference models; the effect on the travel time of inner-core-sensitive seismic phases is comparable to the estimated noise in such measurements. The density contrast at the ICB (0.73 g cm −3 ) is larger than in PREM (0.60 g cm −3 ) and ak135f (0.56 g cm −3 ), but our range including uncertainties (0.52–0.94 g cm −3 ) encompasses all previous estimates in the literature. The average V P and V S in the D″ region are smaller than in PREM , whereas the mean density is probably larger. The data cannot uniquely determine whether this density excess is restricted to the D″ region or distributed throughout the lower(most) mantle. The data cannot determine with certainty the presence or absence of a discontinuity at 220 km depth for V P , V S and density. If present, the jump in both velocities is likely smaller than in PREM . Shear attenuation parameters in the mantle deviate from PREM in a similar fashion to results from more recent studies. We find a non-homogeneous shear attenuation in the inner core, reinforcing the hypothesis that a distinct ‘innermost inner core’ may exist. The bulk attenuation in the mantle and the outer core is stronger than in PREM . We investigate the influence of radial anisotropy on the inversions and analyse possible trade-offs between (anisotropic) parameters. The largest trade-offs are observed in regions that are believed to be anisotropic, such as the D″ region. This illustrates the need to constrain anisotropy in the (deep) mantle. [ABSTRACT FROM AUTHOR]

Abstract: Recently, inversion of magnetic data to recover a distribution of magnetic susceptibility has been widely used for mineral exploration and other problems. However, the commonly used grid-based techniques have some practical difficulties, namely, degraded resolution with depth, increased computational cost with the size of the problem, and the influence of regional field. Especially, most of inversion techniques employ rectangular grid for division, which is inconvenient to represent complicated magnetic structure and actual topography. We presented the magnetic version of body-growth inversion method with the use of non-rectangular grid. Essentially, this method is a new implementation scheme of 2-D and 3-D magnetic inversion, inherited from the gravity inversion by means of growing bodies, previously developed. For simple modification, we adopt non-rectangular grid to divide the subsurface region into a set of isoparametric finite elements rather than the commonly used rectangular grid of prismatic cells. This allows a better representation of actual topography and complicated magnetic structure in forward modeling. Additionally, the calculations of magnetic field and sensitivity matrix are implemented by the Gauss–Legendre quadrature rather than analytic formulae. We tested the method by using synthetic data for two 2-D and three 3-D models and applied it to field data. Resultantly, we conclude that the method has some advantages such as: better representation of actual topography, automatic imposing of the positive magnetic susceptibility, and possibility to separate the regional field and remove its effect to inversion result. Thus, it seems to be a better alternative to the traditional grid-based techniques when inverting local bodies with complicated shapes under mountainous terrain. [Copyright &y& Elsevier]

Abstract: We propose a Bayesian fusion approach to integrate multiple geophysical datasets with different coverage and sensitivity. The fusion strategy is based on the capability of various geophysical methods to provide enough resolution to identify either subsurface material parameters or subsurface structure, or both. We focus on electrical resistivity as the target material parameter and electrical resistivity tomography (ERT), electromagnetic induction (EMI), and ground penetrating radar (GPR) as the set of geophysical methods. However, extending the approach to different sets of geophysical parameters and methods is straightforward. Different geophysical datasets are entered into a trans-dimensional Markov chain Monte Carlo (McMC) search-based joint inversion algorithm. The trans-dimensional property of the McMC algorithm allows dynamic parameterisation of the model space, which in turn helps to avoid bias of the post-inversion results towards a particular model. Given that we are attempting to develop an approach that has practical potential, we discretize the subsurface into an array of one-dimensional earth-models. Accordingly, the ERT data that are collected by using two-dimensional acquisition geometry are re-casted to a set of equivalent vertical electric soundings. Different data are inverted either individually or jointly to estimate one-dimensional subsurface models at discrete locations. We use Shannon's information measure to quantify the information obtained from the inversion of different combinations of geophysical datasets. Information from multiple methods is brought together via introducing joint likelihood function and/or constraining the prior information. A Bayesian maximum entropy approach is used for spatial fusion of spatially dispersed estimated one-dimensional models and mapping of the target parameter. We illustrate the approach with a synthetic dataset and then apply it to a field dataset. We show that the proposed fusion strategy is successful not only in enhancing the subsurface information but also as a survey design tool to identify the appropriate combination of the geophysical tools and show whether application of an individual method for further investigation of a specific site is beneficial. [Copyright &y& Elsevier]

Abstract: The Jiashian earthquake (M L 6.4) occurred on 4 March 2010. It was the largest inland event in southern Taiwan of 2010. The mainshock location was unexpected since it occurred in an area with relatively low background seismicity. In addition, reports of earthquake focal mechanisms do not fit with any known active fault geometry. In order to understand the origin of this earthquake, especially its rupture process, we perform a joint source inversion by using teleseismic body wave, GPS coseismic displacements and near field ground motion data. In this study, we considered a northwest–southeast trending fault with a northeast dip retrieved from GPS coseismic data and aftershocks distribution. To analyze the detailed slip distribution in space and time, we used near field 3D Green’s functions provided by spectral-element method and a full time–space inversion technique. We find a complex rupture process with several slip patches distributed inside two main asperities. The slip map reveals a mean slip of 12.9cm for a maximum slip of 27.3cm leading to a M w 6.47 for this event. The rupture initiates in the deepest portion of the fault at 20km depth, and propagated upward up to 2km depth to form the two asperities. The source time function of this event revealed two pulses corresponding to the two asperities, for a total duration time of about 16s. Most aftershocks occurred near the upper boundary of the deepest asperity while no aftershocks are located close to the shallowest one. We infer that the locations of these slip patches are related to the surrounding fault systems that may have restricted the rupture propagation during the earthquake. [Copyright &y& Elsevier]

Abstract: This paper describes application of a fast inversion method to recover a 3D susceptibility model from magnetic anomalies. For this purpose, the survey area is divided into a large number of rectangular prisms in a mesh with unknown susceptibilities. Solving the full set of equations is substantially time consuming, and applying an algorithm to solve it approximately can reduce the time significantly. It is shown that the Lanczos bidiagonalization method can be an appropriate algorithm to solve a Tikhonov cost function for this purpose. Running time of the inverse modeling significantly decreases by replacing the forward operator matrix with a matrix of lower dimension. A weighted generalized cross validation method is implemented to choose an optimum value of a regularization parameter. To avoid the natural tendency of magnetic structures to concentrate at shallow depth, a depth weighting is applied. This study assumes that there is no remanent magnetization. The method is applied on a noise-corrupted synthetic data to demonstrate its suitability for 3D inversion. A case study including ground based measurement of magnetic anomalies over a porphyry-Cu deposit located in Kerman providence of Iran, Now Chun deposit, is provided to show the performance of the new algorithm on real data. 3D distribution of Cu concentration is used to evaluate the obtained results. The intermediate susceptibility values in the constructed model coincide with the known location of copper mineralization. [Copyright &y& Elsevier]

Abstract: This study reviews the theory, programming designs and merits of two new Matlab-based routines for the forward and inverse Radon transform. These routines offer users flexible choices of integration path functions to take advantage of improved Radon-domain identification and isolation of seismic phases. Least-squares inversion of frequency components and judicious choices of regularization techniques enables additional noise suppression and signal enhancement in the Radon domain. The forward Radon transform routine has the added benefit of spatial interpolation for irregularly sampled data. The accuracy and applicability of these two new routines are demonstrated using data sets containing long-period SS precursors and high-frequency receiver functions. With minimal modifications these two highly portable, carefully documented Radon-transform routines could be easily adapted for a broad range of applications. [Copyright &y& Elsevier]

Abstract: This paper summarizes part of an ongoing feasibility study that investigates the possible use of the full elastic Born approximation in multipole borehole acoustics. As a first step we exclude the fluid-filled borehole with the motivation that one or two wavelengths away from the fluid-filled borehole, radiating borehole mode amplitudes (e.g., Stoneley wave, formation dipole wave, etc.) are small compared to body wave amplitudes (P-, SV- and SH-waves). Consequently, for scatterers one or two wavelengths away from the fluid-filled borehole, it suffices to only consider their interaction with body waves. In this paper we apply the contrast-source stress-velocity forward scattering (integral equation) formulation for solid configurations in its first order (Born-) approximation (De Hoop, 1995) assuming a multipole force source excitation in a zero-offset configuration. To scrutinize the validity of the Born approximation, we consider the simplest type of scatterer, i.e., one characterized by a (Heaviside) step function change in one or more of the contrast (perturbation) parameters and we derive analytic zero-offset formulas for the scattered wave particle velocity and displacement in both the space–frequency and space–time domain, respectively. We assume the scatterer to be located in the far-field. More complicated layered configurations can easily be derived by superposition of the given solution types. Explicit results are given for the dipole and quadrupole excitation, where the former is allowed to have an arbitrary orientation relative to the scatterer and where the latter one is located in a plane perpendicular to the scatterer. In the time domain it is shown, how the scattered wave field decomposes in a specular and diffuse wave field (two terms borrowed from ‘Optics’), where the former contribution vanishes in the absence of an imaging condition and where the latter is always present. For the dipole case, we subject our results to a sensitivity analysis with respect to the three independent perturbation parameters (i.e., density and two compliance parameters) and we compare these results to a full waveform benchmark code that has implemented the reflectivity method (Kennett, 1983), for a ‘horizontally’ stratified elastic medium. This allowed us to pinpoint the root cause of the observed (small) differences. As it turns out these deviations could be traced back to the inaccurateness of the first order Born scattering coefficients. An additional confirmation of this fact is provided through a comparison between the zero-offset scattering coefficients and the corresponding Zoeppritz reflection coefficients. Most notably, it was found that the PP first order scattering coefficient needs a higher than quadratic correction in two of the three independent perturbation parameters, i.e., the two compliance parameters, δΛ and δM. With respect to the density perturbation parameter (δρ) the PP scattering coefficient correction is quadratic with respect to this perturbation parameter, as is to be expected for a first order approximation. Moreover, also the SS first order scattering coefficient only needs a quadratic correction with respect to its associated perturbation parameters, i.e., δρ and δM. Finally, we give a brief outline on how to numerically implement the Born approximation (employing arbitrary offsets) in a configuration where a source–receiver pair is moving continuously relative to the ‘contrast’ (Geology), as is the case in borehole acoustic applications. [Copyright &y& Elsevier]

Abstract: This paper explores the concept of 3D acoustic modelling and inversion in the Laplace domain. We used the Gaussian quadrature integration method to assemble the element mass matrix of a coarse regular rectangular element containing both water and subsurface parts, approximating the properties of these two parts without refining the element into smaller components. The purpose of this investigation was to model an irregular sea floor using 3D Laplace-domain modelling and inversion without increasing the computational cost incurred when using a small grid interval to precisely delineate the sea floor. The algorithm was verified by comparing the analytical solution to the modelled wavefield obtained using the conventional Laplace-domain modelling method and the modelled wavefield obtained using the Gaussian quadrature integration in a two-layer model (with 50-m, 100-m, and 200-m grid spacing and 5, 10, 20, and 40 Gaussian quadrature integration points) where an oblique water–subsurface interface is given. The results show that, under certain conditions, the Gaussian quadrature integration method can generate relatively accurate and consistent wavefields compared to the conventional modelling method. Moreover, we performed 3D Laplace-domain inversion with this Gaussian quadrature modelling technique to obtain synthetic data with a flat and irregular interface using a 100-m and 200-m grid. The results obtained under the conditions herein show that 3D Laplace-domain waveform inversion with Gaussian quadrature integration can generate a more accurate and cost-efficient velocity model compared to that obtained using conventional Laplace-domain waveform inversion. [Copyright &y& Elsevier]

Abstract: Characterizing the internal structure of landslides is of first importance to assess the hazard. Many geophysical techniques have been used in the recent years to image these structures, and among them is seismic tomography. The objective of this work is to present a high resolution seismic inversion algorithm of first arrival times that minimizes the use of subjective regularization operators. A Quasi-Newton P-wave tomography inversion algorithm has been developed. It is based on a finite frequency assumption for highly heterogeneous media which considers an objective inversion regularization (based on the wave propagation principle) and uses the entire source frequency spectrum to improve the tomography resolution. The Fresnel wavepaths calculated for different source frequencies are used to retropropagate the traveltime residuals, assuming that in highly heterogeneous media, the first arrivals are only affected by velocity anomalies present in the first Fresnel zone. The performance of the algorithm is first evaluated on a synthetic dataset, and further applied on a real dataset acquired at the Super-Sauze landslide which is characterized by a complex bedrock geometry, a layering of different materials and important changes in soil porosity (e.g. surface fissures). The seismic P-wave velocity and the wave attenuation are calculated, and the two tomographies are compared to previous studies on the site. [Copyright &y& Elsevier]

Abstract: Geoelectrical techniques are widely used to monitor groundwater processes, while surprisingly few studies have considered audio (AMT) and radio (RMT) magnetotellurics for such purposes. In this numerical investigation, we analyze to what extent inversion results based on AMT and RMT monitoring data can be improved by (1) time-lapse difference inversion; (2) incorporation of statistical information about the expected model update (i.e., the model regularization is based on a geostatistical model); (3) using alternative model norms to quantify temporal changes (i.e., approximations of l 1 and Cauchy norms using iteratively reweighted least-squares), (4) constraining model updates to predefined ranges (i.e., using Lagrange Multipliers to only allow either increases or decreases of electrical resistivity with respect to background conditions). To do so, we consider a simple illustrative model and a more realistic test case related to seawater intrusion. The results are encouraging and show significant improvements when using time-lapse difference inversion with non l 2 model norms. Artifacts that may arise when imposing compactness of regions with temporal changes can be suppressed through inequality constraints to yield models without oscillations outside the true region of temporal changes. Based on these results, we recommend approximate l 1-norm solutions as they can resolve both sharp and smooth interfaces within the same model. [Copyright &y& Elsevier]

Abstract: Convective flow in the mantle can be thought of (and modeled) as exclusively driven by density heterogeneities in the mantle itself, and the resulting lateral variations in the Earth’s gravity field. With this assumption, and a model of mantle rheology, a theoretical relationship can be found between 3D mantle structure and flow-related quantities that can be measured on the Earth’s surface, like free-air gravity anomalies. This relationship can be used to set up an inverse problem, with 1D mantle viscosity as a solution. In the assumption that seismic velocity anomalies be of purely thermal origin, and related to density anomalies by a simple scaling factor, we invert the large-scale length component of the above-mentioned measurements jointly with seismic observations (waveforms and/or travel times) to derive an accurate 5-layer spherically symmetric model of upper- and lower-mantle viscosity. We attempt to account for non-uniqueness in the inverse problem by exploring the solution space, formed of all possible radial profiles of Earth viscosity, by means of a non-deterministic global optimization method: the evolutionary algorithm (EA). For each sampled point of the solution space, a forward calculation is conducted to determine a map of gravity anomalies, whose similarity to GRACE (gravity recovery and climate experiment) is then measured; the procedure is iterated to convergence, according to EA criteria. The robustness of the inversion is tested by means of synthetic tests, indicating that our gravity data set is able to constrain less than 6 radial layers, each with uniform viscosity. Independently of the tomographic model or the scaling factor adopted to convert seismic velocity into density structure, the EA optimization method finds viscosity profiles characterized by low-viscosity in a depth range corresponding to the transition zone, and relatively uniform elsewhere. [Copyright &y& Elsevier]

Abstract: Will the advent of “petascale” computers be relevant to research in global seismic tomography? We illustrate here in detail two possible consequences of the expected leap in computing capability. First, being able to identify larger sets of differently regularized/parameterized solutions in shorter times will allow to evaluate their relative quality by more accurate statistical criteria than in the past. Second, it will become possible to compile large databases of sensitivity kernels, and update them efficiently in a non-linear inversion while iterating towards an optimal solution. We quantify the expected computational cost of the above endeavors, as a function of model resolution, and of the highest considered seismic-wave frequency. [Copyright &y& Elsevier]

This contribution describes the development of an inverse method to extract fault slip and fault geometry from subsurface horizons derived from seismic reflection datasets. The inversion process uses non-linear minimisation together with a forward model using elastic dislocation (ED) theory. The forward model calculates subsurface vertical displacements from fault parameters such as fault slip, dip, length and width for planar, listric or seismically observed fault geometries. Vertical subsurface displacements are measured at selected stratigraphic horizons from seismic depth sections, and act as observations for the inversion process. The non-linear inversion uses a standard procedure (Powell''s method) to minimise an objective function (misfit parameter). The method has been tested using synthetic data generated by theoretical forward models. The inverse method can successfully recover slip distributions and non-planar, including listric, fault geometries. The method has also been applied to stratigraphic horizons derived from seismic reflection data. Fault slip and fault geometry derived from the inverse method are consistent with direct seismic interpretations, and are reproducible and geologically realistic. This inverse method provides a useful complement to conventional subjective seismic interpretations of fault slip and geometry. ED theory together with the retrieved fault parameters from the application of inverse theory may be used to determine strains for fracture prediction. [Copyright &y& Elsevier]

We addressed the sensitivity of Global Navigation Satellite System (GNSS) radio occultation (RO) measurements to atmospheric horizontal variability based on realistically simulated data. Retrieved parameters from refractivity via pressure, and geopotential height to dry temperature were investigated. The errors in a realistic horizontally variable atmosphere relative to errors in a spherically symmetric atmosphere were quantified based on an ensemble of 60 occultation events. These events have been simulated using ray tracing through a representative European Centre for Medium-Range Weather Forecasts (ECMWF) T213L50 analysis field with and without horizontal variability, respectively. Below ∼7 km height biases and standard deviations of all parameters under spherical symmetry are significantly smaller than corresponding errors in a realistic atmosphere with horizontal variability. The relevance of the geometry of reference profiles was assessed in this context as well. A significant part of the total error below ∼7 km can be attributed to adopting reference profiles vertically at mean tangent point locations instead of extracting them along actual 3D tangent point trajectories. The sensitivity of retrieval products to the angle-of-incidence of occultation rays relative to the boresight direction of the receiving antenna was analyzed for three different azimuth sectors (0–10°, 20–30°, 40–50°) with 20 events in each sector. Below about 7 km, most errors were found to increase with increasing angle of incidence. Dry temperature biases between 7 and 20 km exhibit no relevant increase with increasing angle of incidence, which is favorable regarding the climate monitoring utility of the data. [Copyright &y& Elsevier]

Our understanding of earthquakes is limited by our knowledge, and our description, of the physics of the Earth. When solving for subsurface fault slip, it is common practice to assume minimum complexity for characteristics such as topography, fault geometry and elastic properties. These characteristics are rarely accounted for because our knowledge of them is often partial and they can be difficult to include in simulations. However, topography and bathymetry are known all over the Earth's surface, and recently developed software packages such as SPECFEM-X have simplified the process of including them in calculations. Here, we explore the impact of topography on static slip estimates. We also investigate whether the influence of topography can be accounted for with a zeroth-order correction which accounts for variations in distance between subfaults and the surface of the domain. To this end, we analyze the 2015 M w 7.5 Gorkha, Nepal, and the 2010 M w 8.8 Maule, Chile earthquakes within a Bayesian framework. The regions affected by these events represent different types of topography. Chile, which contains both a deep trench and a major orogen, the Andes, has a greater overall elevation range and steeper gradients than Nepal, where the primary topographic feature is the Himalayan mountain range. Additionally, the slip of the continental Nepal event is well-constrained, whereas observations are less informative in a subduction context. We show that topography has a non-negligible impact on inferred slip models. Our results suggest that the effect of topography on slip estimates increases with limited observational constraints and high elevation gradients. In particular, we find that accounting for topography improves slip estimates where topographic gradients are large. When topography has a significant impact on slip, the zeroth-order correction is not sufficient. • Topography has a non-negligible impact on estimated coseismic slip models. • Impact of topography is greater when data constraints are limited. • Accounting for topography has a significant effect where gradients are high. • Zeroth-order correction is not effective when the topographic effect is significant. [ABSTRACT FROM AUTHOR]

Imaging of the subsurface is central in seismic exploration and a topic of great economic interest. A promising technique for seismic imaging is a wave-equation-based method called full-waveform inversion (FWI). FWI is a data-fitting technique that minimises the difference between the observed data in the seismic records and the simulated data, which is extracted from the solution of the wave equation. Usually, FWI is formulated as an optimisation problem that minimises the least-squares distance. In the perspective of likelihood theory, the minimisation of the least-squares distance assumes a Gaussian distribution of the residual data. In this work, we deal with the q -Gaussian distribution associated with the Tsallis statistics to construct a robust optimisation problem, which we call q -FWI. We tested our method in a typical geophysics velocity model with noisy data. Our results show that q -FWI, based on the q -statistics, is a powerful methodology in noisy environments, especially in the presence of outliers. The long tail of the q -distribution exploits the outliers' information, which helps in the image reconstruction. Furthermore, q -FWI provides better image reconstruction without additional computational cost compared to the traditional approach of using the Gaussian distribution for the residuals. • Full-waveform inversion (FWI) is an intricate problem in seismic imaging. • FWI is an optimisation problem that depends on the difference between observed and simulated data. • We construct a new objective function based on q-Gaussian distribution. • The seismic inversion based on the q-statistics performs better than the usual inversion. • The best results of the q-statistics FWI are related to data strongly contaminated by outliers. [ABSTRACT FROM AUTHOR]

Volcanoes are commonly associated with significant topographic relief; however, most analytical models of volcanoes consider the Earth's surface as flat and use half-space solutions. Finite element method can take into account topography; we applied the equivalent body force method to the Changbaishan volcanic area to account for the influence of topographic effects on the inversion of magma chamber pressurization and to calculate the influence of the magma-induced stresses on the stability of the surrounding faults. The influence of surface topography should not be ignored when calculating the surface deformation caused by the magma chamber pressurization in a volcanic area. The shallower the burial depth, the closer the surface maximum radial displacement u r to the center of the magma chamber. When the slope reaches 30°, the maximum vertical displacement u z is no longer directly above the magma chamber. An ellipsoidal pressurization was used to satisfactorily simulate the GPS and leveling measurements at the Changbaishan volcanic area from 2002 to 2003. The superposition of the magma-induced stresses and the regional tectonic stresses likely cause the present stress directions around the Changbaishan volcano, which are dominated by the EW direction to the west of Tianchi Lake and the NW-NNW direction to the north of Tianchi Lake. The Coulomb failure stress change caused by the magma chamber is conducive to the distribution of the NW-trending earthquake swarms that mainly occur to the southwest and northeast of Changbaishan crater. • Apply the equivalent body force method to the Changbaishan volcanic area to account for the influence of topographic effects. • The influence of topography should not be ignored when calculating the surface deformation in a volcanic area. • The magma-induced Coulomb stress change contributed to the distribution of the NW-trending earthquake swarms near the Changbaishan crater. [ABSTRACT FROM AUTHOR]

We examined the preconditioning effects of model parameterizations using the bulk modulus, logarithmic velocity, slowness, and sloth instead of the P-wave velocity, which is commonly used, for monoparameter acoustic full waveform inversions in the Laplace domain. Introducing the bulk modulus, logarithmic velocity, slowness, and sloth modifies the gradient direction of the velocity parameterization according to the chain rule, and these modifications function as preconditioners for the gradient direction and affect the inversion results. We compared five different parameterization methods for monoparameter Laplace-domain full waveform inversions. We also compared the results obtained using an initial velocity model that varies linearly with those obtained using a constant initial velocity model. The logarithmic velocity, slowness, and sloth parameterization methods provided better inversion results than those from the bulk modulus and velocity parameterization in the numerical examples using the Pluto and SEG/EAGE salt models. • Five parameters are tested for monoparameter Laplace-domain waveform inversions. • Different parameters have different preconditioning effects. • The logarithmic velocity, slowness, and sloth showed good convergence behavior. [ABSTRACT FROM AUTHOR]

Two dimensional surface wave travel time tomography is one of the most practical methods to image the lateral variation of phase/group velocity. In practice, local dispersion curves are reconstructed by performing 2-D inversion of travel times for a series of periods separately. Although 2-D inversion is numerically efficient, it does not consider the dependency between periods and thus can lead to a suboptimal phase/group velocity. To address this issue, in this paper, we propose Simultaneous Tomography of all Periods (STP) method. We consider the model parameters as a 3-D volume of dimensions x, y and T, where the first two parameters are related to spatial location of cells and T denotes period. Then a 3-D inversion is performed which inverts all of the travel times for the model parameters simultaneously. To render the solution unique and stable we employ the Tikhonov regularization in which a 3-D regularization functional is used to couple the periods together and force smooth behavior of the parameters in temporal direction as well as spatial direction. Both l2 and l1 norm regularizations are considered. The former is used to generate a flat model by minimizing the squared l2-norm of the model gradient while the latter is used to force the sparsity of the model in the discrete cosine transform (DCT) domain. The results of numerical experiments using both synthetic and real data show that the proposed STP method using the sparsity regularization performs better than the traditional 2-D tomography method in the sense of retrieving small scale features and also stability against noise. • Over the past few decades, surface wave tomography using dispersion curves, extracted from earthquake waveforms and/or ambient noise, has been widely applied to determine lateral variations of phase or group velocity. • 3-D shear wave velocity model is estimated using dispersion curves either directly or indirectly. Indirect methods are based on reconstruction of 2-D phase/group velocity maps for a series of periods separately. Recently, direct methods have been proposed to invert the dispersion curves directly for shear wave velocity model without intermediate reconstruction of 2-D tomography maps (Boschi & Ekström 2002; An et al. 2009; Feng and An 2010; Gao and Shen 2014; Fang et al. 2015; Zhang et al. 2018). The main disadvantages of these methods are high computational cost and sensitivity of the convergence speed to the accuracy of the initial 3-D model. • In this paper, we propose Simultaneous Tomography of all Periods (STP) method. The results of numerical experiments using both synthetic and real data show that the proposed STP method performs significantly better than the traditional 2-D tomography method in the sense of retrieving small scale features and also stability against noise. [ABSTRACT FROM AUTHOR]

• We can record the electrical field for resistivity and induced polarization surveys. • A new inversion procedure is develop to take advantage of this approach. • The algorithm is extended to 4D using a Kalman filter approach. • Four synthetic case studies demonstrates the usefulness of the approach in hydropgeophysics. Induced polarization is a geophysical method that has gained ground in the last decade in hydrogeophysics. Yet the acquisition of high quality induced polarization data may be challenging using the current technologies mostly because of capacitive coupling effects especially for short durations in the injected current. In addition, making a true 3D induced polarization survey is very tedious and time consuming in field conditions. We discuss the advantages of a new generation of induced polarization equipment composed of individual stations able to measure the two components of the electric field along the ground surface. We show how this approach allows for integrating the whole data avoiding negative apparent resistivity/chargeability data. The use of decentralized recording stations solves the issue of capacitive coupling effects. We present a completely novel induced polarization data inversion methodology based on the measurement of the electric field components. In addition, a robust geostatistical inversion approach is discussed for recovering the conductivity and chargeability fields using the electric field components measured on the ground surface. We also treat the case of time lapse monitoring by using a low-rank Kalman filter approximation, which is computationally very appealing in terms of computational time and storage savings. The effectiveness of the new methodology is demonstrated on several realistic synthetic case studies. These numerical tests show that this electric field-based approach is robust and very promising for applications in hydrogeophysics. [ABSTRACT FROM AUTHOR]

The depth and geometry of the bedrock/regolith interface influence many near-surface processes. Geophysics offer means to image them non-invasively since they correspond to discontinuities in physical properties. For instance, refraction seismic applications can leverage the strong correlation between seismic velocity and porosity, while DC resistivity measurements with their sensitivity to electrical conductivity are indirectly sensitive to porosity and clay content. Combining these two data types is beneficial since their spatial resolution patterns and sensitivities to relevant target properties differ. Here, we propose a probabilistic joint inversion framework that uses DC resistivity and first-arrival seismic travel times to infer a common interface geometry separating two heterogeneous sub-domains. A synthetic test case demonstrates that the method is more accurate in sampling the target interface than inversion results obtained using single datasets. We then applied the joint inversion method to field data from the Calhoun Critical Zone Observatory in South Carolina, where it resolves apparent contradictions in the individual inversion results. It produces a bedrock topography that mirrors surface topography and offers clear indications of a more fractured bedrock below the surface topography lows. • Probabilistic inversion results better describe bedrock properties. • Joint inversion provides more accurate interface estimates. • Joint inversions have an enhanced sensitivity to model errors. [ABSTRACT FROM AUTHOR]