Search

SUMMARY Our ability to forecast earthquakes and slow slip events is hampered by limited information on the current state of stress on faults. Ensemble data assimilation methods permit estimating the state by combining physics-based models and observations, while considering their uncertainties. We use an ensemble Kalman filter (EnKF) to estimate shear stresses, slip rates and the state θ acting on a fault point governed by rate-and-state friction embedded in a 1-D elastic medium. We test the effectiveness of data assimilation by conducting perfect model experiments. We assimilate noised shear-stress and velocity synthetic values acquired at a small distance to the fault. The assimilation of uncertain shear stress observations improves in particular the estimates of shear stress on fault segments hosting slow slip events, while assimilating observations of velocity improves their slip-rate estimation. Both types of observations help equally well to better estimate the state θ. For earthquakes, the shear stress observations improve the estimation of shear stress, slip rates and the state θ, whereas the velocity observations improve in particular the slip-rate estimation. Data assimilation significantly improves the estimates of the temporal occurrence of slow slip events and to a large extent also of earthquakes. Rapid and abrupt changes in velocity and shear stress during earthquakes lead to non-Gaussian priors for subsequent assimilation steps, which breaks the assumption of Gaussian priors of the EnKF. In spite of this, the EnKF still provides estimates that are unexpectedly close to the true evolution. In fact, the forecastability for earthquakes for the same alarm duration is very similar to slow slip events, having a very low miss rate with an alarm duration of just 10 per cent of the recurrence interval of the events. These results confirm that data assimilation is a promising approach for the combination of uncertain physics and indirect, noisy observations for the forecasting of both slow slip events and earthquakes.

SUMMARY The non-uniqueness of the solution of inverse geophysical problem has been recognized for a long-time. Although stochastic inversion methods have been developed, deterministic inversion using subsequent regularization is still more widely applied. This is likely due to their efficiency and robustness, compared to the computationally expensive and sometimes difficult to tune to convergence stochastic methods. Recently, Bayesian evidential learning 1-D imaging has been presented to the community as a viable tool for the efficient stochastic 1-D imaging of the subsurface based on geophysical data. The method has been proven to be as fast, or sometimes even faster, than deterministic solution. However, the method has a significant drawback when dealing with large prior uncertainty as often encountered in geophysical surveys: it tends to overestimate the uncertainty range. In this paper, we provide an efficient way to overcome this limitation through iterative prior resampling (IPR) followed by rejection sampling. IPR adds the posterior distribution calculated at a former iteration to the prior distribution in a subsequent iteration. This allows to sharpen the learning phase of the algorithm and improve the estimation of the final posterior distribution while rejection sampling eliminates models not fitting the data. In this contribution, we demonstrate that this new approach allows BEL1D to converge towards the true posterior distribution. We also analyse the convergence behaviour of the algorithm and derive guidelines for its application. We apply the approach for the interpretation of surface waves dispersion curves but the approach can be generalized to other geophysical methods.

SUMMARY Full-waveform inversion (FWI) of ground-penetrating radar (GPR) data has received particular attention in the past decade because it can provide high-resolution subsurface models of dielectric permittivity and electrical conductivity. In most GPR FWIs, these two parameters are regarded as frequency independent, which may lead to false estimates if they strongly depend on frequency, such as in shallow weathered zones. In this study, we develop frequency-dependent GPR FWI to solve this problem. Using the τ-method introduced in the research of viscoelastic waves, we define the permittivity attenuation parameter to quantify the attenuation resulting from the complex permittivity and to modify time-domain Maxwell’s equations. The new equations are self-adjoint so that we can use the same forward engine to back-propagate the adjoint sources and easily derive model gradients in GPR FWI. Frequency dependence analysis shows that permittivity attenuation acts as a low-pass filter, distorting the waveform and decaying the amplitude of the electromagnetic waves. The 2-D synthetic examples illustrate that permittivity attenuation has low sensitivity to the surface multioffset GPR data but is necessary for a good reconstruction of permittivity and conductivity models in frequency-dependent GPR FWI. As a comparison, frequency-independent GPR FWI produces more model artefacts and hardly reconstructs conductivity models dominated by permittivity attenuation. The 2-D field example shows that both FWIs reveal a triangle permittivity anomaly which proves to be a refilled trench. However, frequency-dependent GPR FWI provides a better fit to the observed data and a more robust conductivity reconstruction in a high permittivity attenuation environment. Our GPR FWI results are consistent with previous GPR and shallow-seismic measurements. This research greatly expands the application of GPR FWI in more complicated media.

SUMMARY Transient electromagnetic (TEM) is an efficient non-invasive method to map electrical conductivity distribution in the subsurface. This paper presents an inversion scheme for 3-D TEM time-lapse (TL) data using a generalized minimum support (MS) norm and its application to monitoring conductivity changes over time. In particular, two challenges for TL TEM applications are addressed: (i) the survey repetition with slightly different acquisition position, that is, because systems are not installed and (ii) non-optimal data coverage above the TL anomalies, for instance, due to the presence of infrastructure that limits the acquisition layout because of coupling. To address these issues, we developed a new TEM TL inversion scheme with the following features: (1) a multimesh approach for model definition and forward computations, which allows for seamless integration of data sets with different acquisition layouts; (2) 3-D sensitivity calculation during the inversion, which allows retrieving conductivity changes in-between TEM soundings and (3) simultaneous inversion of two data sets at once, imposing TL constraints defined in terms of a generalized MS norm, which ensures compact TL changes. We assess the relevance of our implementations through a synthetic example and a field example. In the synthetic example, we study the capability of the inversion scheme to retrieve compact time-lapse changes despite slight changes in the acquisition layout and the effect of data coverage on the retrieval of TL changes. Results from the synthetic tests are used for interpreting field data, which consists of two TEM data sets collected in 2019 and 2020 at the Nesjavellir high-temperature geothermal site (Iceland) within a monitoring project of H2S sequestration. Furthermore, the field example illustrates the effect of the trade-off between data misfit and TL constraints in the inversion objective function, using the tuning settings of the generalized MS norm. Based on the results from both the synthetic and field cases, we show that our implementation of 3-D TL inversion has a robust performance for TEM monitoring.

SUMMARY Airborne systems collecting transient electromagnetic data are able to gather large amounts of data over large areas in a very short time. These data are most often interpreted through 1-D inversions, due to the availability of robust, fast and efficient codes. However, in areas where the subsurface contains complex structures or large conductivity contrasts, 1-D inversions may introduce artefacts into the models, which may prevent correct interpretation of the results. In these cases, 2-D or 3-D inversion should be used. Here, we present a 2.5-D inversion code using 3-D forward modelling combined with a 2-D model grid. A 2.5-D inversion is useful where the flight lines are spaced far apart, in which case a 3-D inversion would not add value in relation to the added computational cost and complexity. By exploiting the symmetry of the transmitter and receiver system we are able to perform forward calculations on a reduced 3-D mesh using only half the domain transecting the centre of the transmitter and receiver system. The forward responses and sensitivities from the reduced 3-D mesh are projected onto a structured 2-D model grid following the flight direction. The difference in forward calculations is within 1.4 per cent using the reduced mesh compared to a full 3-D solution. The inversion code is tested on a synthetic example constructed with complex geology and high conductivity contrasts and the results are compared to a 1-D inversion. We find that the 2.5-D inversion recovers both the conductivity values and shape of the true model with a significantly higher accuracy than the 1-D inversion. Finally, the results are supported by a field case using airborne TEM data from the island of Mayotte. The inverted flight line consisted of 418 soundings, and the inversion spent an average of 6750 s per iteration, converging in 16 iterations with a peak memory usage of 97 GB, using 18 logical processors. In general, the total time of the 2-D inversions compared to a full 3-D inversion is reduced by a factor of 2.5 while the memory consumption was reduced by a factor of 2, reflecting the half-mesh approach.

In seismic exploration methods, imperfect spatial sampling at the surface causes a lack of illumination at the target in the subsurface. The hampered image quality at the target area of interest causes uncertainties in reservoir monitoring and production, which can have a substantial economic impact. Especially in the case of a complex overburden, the impact of surface sampling on target illumination can be significant. The target-oriented acquisition analysis based on wavefield propagation and a known velocity model has been used to provide guidance for optimizing the acquisition parameters. Seismic acquisition design is usually a manual optimization process, with consideration of many aspects. In this study, we develop a methodology that automatically optimizes an irregular receiver geometry when the source geometry is fixed or vice versa. The methodology includes objective functions defined by two criteria: optimizing the image resolution and optimizing the angle-dependent illumination information. We use a two-step parameterization in order to make the problem more linear and, thereby, solve the acquisition design problem by using a gradient descent algorithm. With simple and complex velocity models, we demonstrate that the proposed method is effective, while the involved computational cost is acceptable. Interestingly, the optimization results in our examples show that the conventional uniform geometry already satisfies the resolution requirement, while optimizing for angle coverage can provide a large uplift and is strongly dependent on the velocity model.

The induced seismicity in the Groningen gas field, The Netherlands, presents contrasted spatio-temporal patterns between the central area and the south west area. Understanding the origin of this contrast requires a thorough assessment of two factors: (1) the stress development on the Groningen faults, and (2) the frictional response of the faults to induced stresses. Both factors have large uncertainties that must be honored and then reduced thanks to observational constraints. The present contribution builds upon the previous work of Candela et al. (2019), which was focused on the frictional response, and extends it by taking into account uncertainties in the stress development and by evaluating the relative effect of both factors. Ensembles of induced stress realizations are built by varying the Poisson’s ratio in a poro-elastic model incorporating the 3D complexities of the geometries of the Groningen gas reservoir and its faults, and the historical pore pressure distribution. The a-priori uncertainties in the frictional response are mapped by varying the parameters of a seismicity model based on rate-and-state friction. The uncertainties of each component of this complex physics-based model are honored through an efficient data assimilation algorithm. By assimilating the seismicity data with an Ensemble-Smoother, the prior uncertainties of each model parameter are effectively reduced, and the posterior seismicity rate predictions are consistent with the observations. Our integrated workflow allows us to disentangle the contributions of the main two factors controlling the induced seismicity at Groningen, induced stress development and fault frictional response. Posterior distributions of the model parameters of each modelling component are contrasted between the central and south west area at Groningen. We find that, even after honoring the spatial heterogeneity in stress development across the Groningen gas field, the spatial variability of the observed induced seismicity rate still requires spatial heterogeneity in the fault frictional response. This work is enabled by the unprecedented deployment of an Ensemble-Smoother combined with physics-based modelling over a complex case of reservoir induced seismicity.

SUMMARY We present a solution for long-term direct current resistivity and time-domain induced polarization (DCIP) monitoring, which consists of a monitoring system and the associated software that automates the data collection and processing. This paper describes the acquisition system that is used for remote data collection and then introduces the routines that have been developed for pre-processing of the monitoring data set. The collected data set is pre-processed using digital signal processing algorithms for outlier detection and removal; the resulting data set is then used for the inversion procedure. The suggested processing workflow is tested against a simulated time-lapse experiment and then applied to field data. The results from the simulation show that the suggested approach is very efficient for detecting changes in the subsurface; however, there are some limitations when no a priori information is used. Furthermore, the mean weekly data sets that are generated from the daily collected data can resolve low-frequency changes, making the approach a good option for monitoring experiments where slow changes occur (i.e. leachates in landfills, internal erosion in dams, bioremediation). The workflow is then used to process a large data set containing 20 months of daily monitoring data from a field site where a pilot test of in situ bioremediation is taking place. Based on the time-series analysis of the inverted data sets, we can detect two portions of the ground that show different geophysical properties and that coincide with the locations where the different fluids were injected. The approach that we used in this paper provides consistency in the data processing and has the possibility to be applied to further real-time geophysical monitoring in the future.

SUMMARY Iterative least-squares reverse time migration (LSRTM) is the state-of-the-art linearized waveform inversion method to obtain quantitative subsurface parameters. The main drawback of such an iterative imaging scheme is the significant computational expense of many modelling/adjoint cycles through iterations. In the context of the extended domain, an interesting alternative to LSRTM is the asymptotic linearized direct waveform inversion, providing quantitative results with only a single iteration. This approach was first proposed for constant-density acoustics and recently extended to the variable-density case. The former is based on the application of the asymptotic inverse Born operator, whereas the latter has two more extra steps: building an angle-dependent response of the asymptotic inverse Born operator and then solving a weighted least-squares approach for simultaneous inversion of two acoustic parameters. To examine the importance of accounting for density variations, we compare constant- and variable-density linearized direct waveform inversion techniques applied to a marine real data set from the Eastern Nankai Trough, offshore Japan. The inversion results confirm the efficiency of the asymptotic linearized direct waveform inversion in estimating quantitative parameters within a single iteration. The variable-density direct inversion yields subsurface images that (1) exhibit a superior resolution and (2) better reconstruct the field data than does the constant-density approach, even if the data set does not contain large enough surface offset to fully decompose velocity and density perturbations.

Seismic reservoir characterization is a subfield of geophysics that combines seismic and rock-physics modeling with mathematical inverse theory to predict reservoir variables from the measured seismic data. An open-source comprehensive modeling library that includes the main concepts and tools is still missing. We have developed a Python library called SeReMpy with state-of-the-art of seismic reservoir modeling for reservoir properties characterization using seismic and rock-physics models and Bayesian inverse theory. The most innovative component of the library is the Bayesian seismic and rock-physics inversion to predict the spatial distribution of petrophysical and elastic properties from seismic data. The inversion algorithms include Bayesian analytical solutions of the linear-Gaussian inverse problem and Markov chain Monte Carlo (MCMC) numerical methods for nonlinear problems. The library includes four modules: geostatistics, rock physics, facies, and inversion, as well as several scripts with illustrative examples and applications. We illustrate the use of the functions of the module and develop codes for practical inversion problems using synthetic and real data. The applications include a rock-physics model for the prediction of elastic properties and facies using well-log data, a geostatistical simulation of continuous and discrete properties using well logs, a geostatistical interpolation and simulation of 2D maps of temperature, an elastic inversion of partial stacked seismograms with Bayesian linearized amplitude-variation-with-offset inversion, a rock-physics inversion of partial stacked seismograms with MCMC methods, and a 2D seismic inversion.

SUMMARYThis paper presents a new sparse inversion method based on L1 norm regularization for 3-D magnetic data. In isolation, L1 norm regularization yields model elements which are unconstrained by the input data to be exactly zero, leading to a sparse model with compact and focused structure. Here, we complement the L1 norm with a penalty minimizing total variation, the L1 norm of the model gradients; it is expected that the sharp boundaries of the subsurface structure are not compromised by incorporating this penalty. Although this penalty is widely used in the geophysical inversion studies, it is often replaced by an alternative quadratic penalty to ease solution of the penalized inversion problem; in this study, the original definition of the total variation, that is form of the L1 norm of the model gradients, is used. To solve the problem with this combined penalty of L1 norm and total variation, this study introduces alternative direction method of multipliers, which is a primal-dual optimization algorithm that solves convex penalized problems based on the optimization of an augmented Lagrange function. To improve the computational efficiency of the algorithm to make this method applicable to large-scale magnetic inverse problems, this study applies matrix compression using the wavelet transform and the preconditioned conjugate gradient method. The inversion method is applied to both synthetic tests and real data, the synthetic tests demonstrate that, when subsurface structure is blocky, it can be reproduced almost perfectly.

Since the 1960s, all studies have assumed that a film thickness “h” provides a unique pressure field “p” by resolving the Reynolds equation. However, it is relevant to investigate the film thickness unicity under a given hydrodynamic pressure within the inverse theory. This paper presents a new approach to deduce from an initial film thickness a widespread number of thicknesses providing the same hydrodynamic pressure under a specific condition of gradient pressure. For this purpose, three steps were presented: 1) computing the hydrodynamic pressure from an initial film thickness by resolving the Reynolds equation with Gümbel’s cavitation model, 2) using a new algorithm to generate a second film thickness, 3) comparing and validating the hydrodynamic pressure produced by both thicknesses with the modified Reynolds equation. Throughout three surface finishes: the macro-shaped, micro-textured, and rough surfaces, it has been demonstrated that under a specific hydrodynamic pressure gradient, several film thicknesses could generate the same pressure field with a slight difference by considering cavitation. Besides, this paper confirms also that with different ratios of the averaged film thickness to the root mean square (RMS) similar hydrodynamic pressure could be generated, thereby the deficiency of this ratio to define the lubrication regime as commonly known from Patir and Cheng theory.

SUMMARY Magnetization vector inversion (MVI) has attracted considerable attention in recent years since by this inversion both distribution of the magnitude and direction of the magnetization are obtained; therefore, it is easy to distinguish between the magnetic causative bodies especially when magnetic data are affected by different remanent magnetization. In this research, the compact magnetization vector inversion is presented: a 3-D magnetic modelling is proposed from surface data measurements to obtain compact magnetization distribution. The equations are solved in data-space least squares and the algorithm includes a combination of two weights as depth weighting and compactness weighting in the Cartesian system. The re-weighted compactness weighting matrix handles sparsity constraints imposed on the magnitude of magnetization for varying Lp-norms ($0 \le p \le 2$). The low value of the norm leads to more focused or compact inversion, and using a large value of p obtains a smooth model. The method is validated with two synthetic examples, the first is a cube that has significant remanent magnetization and the second consists of two causative cube bodies with significant different magnetization directions at different depths. The case study is the magnetic data of Galinge iron ore deposit (China) that the apparent susceptibility and magnetization directions are reconstructed. The compact model reveals that the results agree with drilling and geological information.

SUMMARY The separation of regional-residual anomalies plays an important role in the processing of potential field anomalies for obtaining better understandings of the nature of the underground sources. Many methods have been developed to achieve the separation of anomalies that are of distinct wavelengths. On the other hand, fewer studies have addressed the separation of local anomalies from the observed potential field anomalies. In this paper, we introduce a new process for separating localized anomalies from the observations under the Cartesian and spherical coordinates. The separation is achieved using the equivalent source technique and an iterative inversion process which is to refine and finalize the separated local anomalies. Additionally, we introduce an inversion method for determining the equivalent sources that are of varying dimensions, as well as a quantitative measurement to assess the accuracy of the separation process. Verified with synthetic examples, the proposed method could extract arbitrary shaped local anomalies from the rest with low error levels. Subsequently, we apply the method to the construction of a 3-D model of the mantle uplift beneath the Von Kármán crater (VKC) on the Moon. The VKC is the landing site of the Chinese lunar exploration mission Chang'e 4, which lies in the northwestern portion of the South-Pole Aitken (SPA) basin on the far side of the Moon. Multiple generations of mare basalts are identified within the VKC, which indicates a complex geological history of the basin. Insights into the evolutionary history of this region can be obtained by investigating the deep crustal structure of the VKC using topographic and gravity data. Processed with the proposed method, the 3-D structure we obtain provides evidence for separated mantle uplifting events triggered by the two impact events that created the VKC and the Von Kármán M crater, respectively.

SUMMARYThis study examines local geopotential field modelling over a mountainous region in Tanzania using vector airborne gravity data. We use the adjusted spherical cap and rectangular harmonic analyses. Both methods are based on expansion of gravitational potential into a series of orthogonal harmonic basis functions of local support in such a way that the expansion coefficients are determined by gravity observations. All three components of gravity vector are simultaneously inverted to derive the geopotential coefficients. In order to evaluate the accuracy of the local models, independent checkpoints are selected within the study region and around its boundary and the computed gravity vectors are compared with the independent gravity observations. The results show an excellent agreement with root mean square error (RMSE) of

SUMMARYThe injection of CO2 at the Decatur carbon capture and storage site has generated significant microseismic activity, which occurs in distinct spatial clusters up to approximately 2.2 km from the primary injection well. Accurate and precise event locations are vital for the characterization of the microseismicity to help understand the reservoir response to the CO2 injection, whilst enabling the identification of minor faults and fractures below the resolution of conventional active seismic imaging. However, microseismic monitoring of fluid injection sites, such as Decatur, is often performed using a network of borehole sensors often from a single well. While these downhole sensors have excellent detection capabilities, their poor azimuthal coverage limits the ability to precisely determine event locations. We have developed a modified double-difference relocation algorithm suitable for both 1-D and 3-D velocity models, and which incorporates differential back azimuth observations to allow the benefits of the original double-difference algorithm to be applicable to a downhole microseismic monitoring setting. Applying the modified double-difference algorithm to the microseismicity at Decatur, we have successfully relocated 4293 events. The relocation included over 59 million observations for 757 285 event pairs, split across seven geographic regions. Despite the majority of observations being recorded in only two boreholes, with an almost identical azimuthal coverage, the results have shown to be reliable with significantly reduced residuals and low uncertainties associated with the final locations. We have analysed the residuals in terms of their association with each geographic region, data type, station and individual events, to fully appreciate their influence in the inversion and the fit of the data to the final set of event locations. For each region, the relocated seismicity has become less diffuse with improved clustering, and with newly visible linear features often orientated in a NE–SW direction. These results show the potential improvements that can be made to microseismic event locations recorded by a borehole network with a limited and variable azimuthal distribution.

SUMMARYGravity and magnetic data resolve the Earth with variable spatial resolution, and Earth structure exhibits both discontinuous and gradual features. Therefore, model parametrization complexity should be able to address such variability by locally adapting to the resolving power of the data. The reversible-jump Markov chain Monte Carlo (rjMcMC) algorithm provides variable spatial resolution that is consistent with data information. To address the prevalent non-uniqueness in joint inversion of potential field data, we use a novel spatial partitioning with nested Voronoi cells that is explored by rjMcMC sampling. The nested Voronoi parametrization partitions the subsurface in terms of rock types, such as sedimentary, salt and basement rocks. Therefore, meaningful prior information can be specified for each type which reduces non-uniqueness. We apply nonoverlapping prior distributions for density contrast and susceptibility between rock types. In addition, the choice of noise parametrization can lead to significant trade-offs with model resolution and complexity. We adopt an empirical estimation of full data covariance matrices that include theory and observational errors to account for spatially correlated noise. The method is applied to 2-D gravity and magnetic data to study salt and basement structures. We demonstrate that meaningful partitioning of the subsurface into sediment, salt, and basement structures is achieved by these advances without requiring regularization. Multiple simulated- and field-data examples are presented. Simulation results show clear delineation of salt and basement structures while resolving variable length scales. The field data show results that are consistent with observations made in the simulations. In particular, we resolve geologically plausible structures with varying length scales and clearly differentiate salt structure and basement topography.

Conventional seismic source inversion estimates the earthquake rupture process on an assumed fault plane that is determined a priori. It has been a difficult challenge to obtain the fault geometry together with the rupture process by seismic source inversion because of the nonlinearity of the inversion technique. In this study, we propose an inversion method to estimate the fault geometry and the rupture process of an earthquake from teleseismic P waveform data, through an elaboration of our previously published finite-fault inversion analysis (Shimizu et al. 2020). That method differs from conventional methods by representing slip on a fault plane with five basis double-couple components, expressed by potency density tensors, instead of two double-couple components compatible with the fault direction. Because the slip direction obtained from the potency density tensors should be compatible with the fault direction, we can obtain the fault geometry consistent with the rupture process. In practice we rely on an iterative process, first assuming a flat fault plane and then updating the fault geometry by using the information included in the obtained potency density tensors. In constructing a non-flat model-fault plane, we assume for simplicity that the fault direction changes only in either the strike or the dip direction. After checking the validity of the proposed method through synthetic tests, we applied it to the Mw 7.7 2013 Balochistan, Pakistan, and Mw 7.9 2015 Gorkha, Nepal, earthquakes, which occurred along geometrically complex fault systems. The modelled fault for the Balochistan earthquake is a curved strike-slip fault convex to the south-east, which is consistent with the observed surface ruptures. The modelled fault for the Gorkha earthquake is a reverse fault with a ramp-flat-ramp structure, which is also consistent with the fault geometry derived from geodetic and geological data. These results exhibit that the proposed method works well for constraining fault geometry of an earthquake.

SUMMARY We present a method for inverting total-field anomaly data to estimate the geometry of a uniformly magnetized 3-D geological source in the subsurface. The method assumes the total-magnetization direction is known. We approximate the source by an ensemble of vertically juxtaposed right prisms, all of them with the same total-magnetization vector and depth extent. The horizontal cross-section of each prism is a polygon defined by a given number of equi-angularly spaced vertices from 0° to 360°, whose polygon vertices are described by polar coordinates with an origin defined by a horizontal location over the top of each prism. Because our method estimates the radii of each polygon vertex we refer to it as radial inversion. The position of these vertices, the horizontal location of each prism and the depth extent of all prisms are the parameters to be estimated by solving a constrained non-linear inverse problem of minimizing a goal function. We run successive inversions for a range of tentative total-magnetization intensities and depths to the top of the shallowest prism. The estimated models producing the lowest values of the goal function form the set of candidate solutions. To obtain stabilized solutions, we impose the zeroth- and first-order Tikhonov regularizations on the shape of the prisms. The method allows estimating the geometry of both vertical and inclined sources, with a constant direction of magnetization, by using the Tikhonov regularization. Tests with synthetic data show that the method can be of utility in estimating the shape of the magnetic source even in the presence of a strong regional field. Results obtained by inverting airborne total-field anomaly data over the Anitápolis alkaline–carbonatitic complex, in southern Brazil, suggest that the emplacement of the magnetic sources was controlled by NW–SE-trending faults at depth, in accordance with known structural features at the study area.

SUMMARY The multichannel transient electromagnetic method (MTEM), sensitive to both near-surface and deeply buried geoelectric structures, is a useful geophysical tool to map resistive targets. A typical MTEM survey consists of a large number of receiver and transmitter positions. The produced data volume is much larger than a traditional transient electromagnetic method survey, posing computational difficulties for its data inversion and interpretation. We proposed a 3-D parallel inversion algorithm for MTEM data based on nonlinear conjugate gradient method and the moving footprint technique. The electromagnetic response and Fréchet derivative are first calculated in the frequency domain with integral equation method and then cosine transformed into time domain. We not only fully parallelize the forward modelling and Fréchet derivative calculation but also use a moving footprint to speed up the inversion. The footprints for different offsets are carefully evaluated through synthetic models. The effectiveness of the developed algorithm has been demonstrated through synthetic and case studies.

In frequency domain Electromagnetic Induction (EMI) surveys, an image of the electrical conductivity of the subsurface is obtained non-invasively. The electrical conductivity can be related to important subsurface properties such as the porosity, saturation or water conductivity via Archie’s law. The advantage of geophysical EMI surveys is its cost-effectiveness because it is a non-contacting method, one can easily walk with the device or mount in on a vehicle or a helicopter (AEM).The process of finding the conductivity profile from the collected field data is an ill-posed inverse problem. Regularization improves the stability of the inversion and, based on Occam’s razor principle, a smoothing constraint is typically used with a very large number of thin layers. However, the conductivity profiles are not always expected to be smooth. Another alternative is to use a predefined number of layers and to invert for their conductivity and thickness. This can yield sharp contrasts in conductivity. In practice however, the real underground might be either blocky or smooth, or somewhere in between. Those standard constraints are thus not always appropriate.We develop a new minimum-structure inversion scheme in which we transform the model into the wavelet space and impose a sparsity constraint. This sparsity constrained inversion scheme minimizes an objective function with a least-squares data misfit and a sparsity measure of the model in the wavelet domain. With a solid understanding of wavelet theory, a novel and intuitive model misfit term was developed, allowing for both smooth and blocky models, depending on the chosen wavelet basis. A model in the wavelet domain has both temporal (i.e. low and high frequencies) and spatial resolution, and penalizing small-scale coefficients effectively reduces the complexity of the model.Comparing the novel scale-dependent wavelet-based regularization scheme with wavelet-based regularization with no scale-dependence, revealed significantly better results (Figure A and B) w.r.t. the true model. Comparing with standard Tikhonov regularization (Figure C and D) shows that our scheme can recover high amplitude anomalies in combination with globally smooth profiles. Furthermore, the adaptive nature of the inversion method (due to the choice of wavelet) allows for high flexibility because the shape of the wavelet can be exploited to generate multiple representations (smooth, blocky or intermediate) of the inverse model.We have introduced an alternative inversion scheme for EMI surveys that can be extended to any other 1D geophysical method. It involves a new model misfit or regularization term based on the wavelet transform and scale-dependent weighting which can easily be combined with the existing framework of deterministic inversion (gradient-based optimization methods, L-curve criterion for optimal regularization parameter). A challenge remains to select the optimal wavelet, however, the ensemble of inversion results with different wavelets can also be used to qualitatively assess uncertainty.

SUMMARY Time-lapse resistivity tomography bring valuable information on the physical changes occurring inside a geological reservoir. In this study, resistivity monitoring from controlled source electromagnetics (CSEM) data is investigated through synthetic and real data. We present three different schemes currently used to perform time-lapse inversions and compare these three methods: parallel, sequential and double difference. We demonstrate on synthetic tests that double difference scheme is the best way to perform time-lapse inversion when the survey parameters are fixed between the different time-lapse acquisitions. We show that double difference inversion allows to remove the imprint of correlated noise distortions, static shifts and most of the non-linearity of the inversion process including numerical noise and acquisition footprint. It also appears that this approach is robust against the baseline resistivity model quality, and even a rough starting resistivity model built from borehole logs or basic geological knowledge can be sufficient to map the time-lapse changes at their right positions. We perform these comparisons with real land time-lapse CSEM data acquired one year apart over the Reykjanes geothermal field.

SUMMARYElastic full waveform inversion (EFWI) can, theoretically, give high-resolution estimates of the subsurface. However, in practice, the resolution and illumination of EFWI are limited by the bandwidth and aperture of seismic data. The often-present wells in developed fields as well as some exploratory regions can provide complementary information of the target area. We, thus, introduce a regularization technique, which combines the surface seismic and well log data, to help improve the resolution of EFWI. We use deep neural networks to learn the statistical relations between some selected features of the inverted model and the facies interpreted from well logs. The selected features are the means and variances of the inverted velocities defined within Gaussian windows. Using multiple fully connected layers, we train our neural networks to identify the relation between the means and variances at the well location and those from the inverted model. The network is used to map the means and variances extracted from the well to the whole model domain. We then perform another EFWI in which we fit the predicted data to the observed ones as well as fit the model over a Gaussian window to the predicted means and variances. The tests on synthetic and real seismic data demonstrate that the proposed method can effectively improve the resolution and illumination of deep-buried reservoirs, which often encounter poor illumination from seismic data.

In the context of seismic imaging, the full waveform inversion (FWI) is more and more popular. Because of its lower numerical cost, the acoustic approximation is often used, especially at the exploration geophysics scale, both for tests and for real data. Moreover, some research domains such as helioseismology face true acoustic medium for which FWI can be useful. In this work, we show that the general acoustic inverse problem based on limited frequency band data is intrinsically non-unique, making any general acoustic FWI impossible. Our work is based on two tools: particle relabelling and homogenization. On the one hand, the particle relabelling method shows it is possible to deform a true medium based on a smooth mapping into a new one without changing the signal recorded at seismic stations. This is a potentially strong source of non-uniqueness for an inverse problem based a seismic data. Nevertheless, in the elastic case, the deformed medium loses the elastic tensor minor symmetries and, in the acoustic case, it implies density anisotropy. It is therefore not a source of non-uniqueness for elastic or isotropic acoustic inverse problems, but it is for the anisotropic acoustic case. On the other hand, the homogenization method shows that any fine-scale medium can be up-scaled into an effective medium without changing the waveforms in a limited frequency band. The effective media are in general anisotropic, both in the elastic and acoustic cases, even if the true media are isotropic at a fine scale. It implies that anisotropy is in general present and needs to be inverted. Therefore, acoustic anisotropy can not be avoided in general. We conclude, based on a particle relabelling and homogenization arguments, that the acoustic FWI solution is in general non-unique. We show, in 2-D numerical FWI examples based on the Gauss-Newton iterative scheme, the effects of this non-uniqueness in the local optimization context. We numerically confirm that the acoustic FWI is in general non-unique and that finding a physical solution is not possible.

SUMMARY To understand earth processes, geoscientists infer subsurface earth properties such as electromagnetic resistivity or seismic velocity from surface observations of electromagnetic or seismic data. These properties are used to populate an earth model vector, and the spatial variation of properties across this vector sheds light on the underlying earth structure or physical phenomenon of interest, from groundwater aquifers to plate tectonics. However, to infer these properties the spatial characteristics of these properties need to be known in advance. Typically, assumptions are made about the length scales of earth properties, which are encoded a priori in a Bayesian probabilistic setting. In an optimization setting, appeals are made to promote model simplicity together with constraints which keep models close to a preferred model. All of these approaches are valid, though they can lead to unintended features in the resulting inferred geophysical models owing to inappropriate prior assumptions, constraints or even the nature of the solution basis functions. In this work it will be shown that in order to make accurate inferences about earth properties, inferences can first be made about the underlying length scales of these properties in a very general solution basis. From a mathematical point of view, these spatial characteristics of earth properties can be conveniently thought of as ‘properties’ of the earth properties. Thus, the same machinery used to infer earth properties can be used to infer their length scales. This can be thought of as an ‘infer to infer’ paradigm analogous to the ‘learning to learn’ paradigm which is now commonplace in the machine learning literature. However, it must be noted that (geophysical) inference is not the same as (machine) learning, though there are many common elements which allow for cross-pollination of useful ideas from one field to the other, as is shown here. A non-stationary trans-dimensional Gaussian Process (TDGP) is used to parametrize earth properties, and a multichannel stationary TDGP is used to parametrize the length scales associated with the earth property in question. Using non-stationary kernels, that is kernels with spatially variable length scales, models with sharp discontinuities can be represented within this framework. As GPs are multidimensional interpolators, the same theory and computer code can be used to solve geophysical problems in 1-D, 2-D and 3-D. This is demonstrated through a combination of 1-D and 2-D non-linear regression examples and a controlled source electromagnetic field example. The key difference between this and previous work using TDGP is generalized nested inference and the marginalization of prior length scales for better posterior subsurface property characterization.

SUMMARYAn underground nuclear explosion (UNE) couples mechanical energy into crustal rock, which propagates as seismic and acoustic waves. These different physical phenomena transport, by different pathways, to standoff detectors at varying distances. The transport pathways attenuate the original signal but in different ways. Enabled by correct statistical weighting, signal attenuation models can be used to combine these disparate sensor data to estimate the yield of an UNE. Contemporaneous statistical models, used in yield estimation, can be improved with an advanced partition of error for these physical signal propagation models. We present an advanced multivariate approach to error modelling of multiphenomenology physical signatures. In addition to measurement error, our error model represents physical model biases as random with a physics-based covariance structure. To illustrate this proposed framework, we demonstrate the estimation of explosion yield using openly available seismic and acoustic data from chemical single-point explosions.

SUMMARY During the Cenozoic, the Earth experienced multiple first-order geological events that are likely mantle flow related. These include the termination of large-scale marine inundation in North America in the Palaeocene, the late Tertiary rise of Africa relative to other continents and the long-wavelength tilting of Australia since the late Cretaceous, which occurred when the continent approached the southeast Asia subduction systems on its northward passage from Antartica. Here we explore a suite of eight high-resolution, compressible, global mantle flow retrodictions going back to 50 Ma, using an adoint method with $\approx$670 million finite elements. These retrodictions show for the first time that these events emerge jointly as part of global Cenozoic mantle flow histories. Our retrodictions involve the dynamic effects from an upper mantle low-viscosity zone, assimilate a past plate-motion model for the tangential surface velocity field, probe the influence of two different present-day mantle state estimates derived from seismic tomography, and acknowledge the rheological uncertainties of dynamic Earth models by taking in four different realizations for the radial mantle viscosity profile, two of which were published previously. We find the retrodicted mantle flow histories are sensitive to the present-day mantle state estimate and the rheological properties of the Earth model, meaning that this input information is testable with inferences gleaned from the geological record. For a deep mantle viscosity of $1.7\times 10^{22}$ Pa s and a purely thermal interpretation of seismic structure, lower mantle flow velocities exceed 7 cm yr–1 in some regions, meaning they are difficult to reconcile with the existence of a hotspot reference frame. Conversely, a deep mantle viscosity of $10^{23}$ Pa s yields modest flow velocities (

SUMMARY Interfaces are important part of Earth’s layering structure. Here, we developed a new model parametrization and iterative linearized inversion method that determines 1-D crustal velocity structure using surface wave dispersion, teleseismic P-wave receiver functions and Ps and PmP traveltimes. Unlike previous joint inversion methods, the new model parametrization includes interface depths and layer Vp/Vs ratios so that smoothness constraint can be conveniently applied to velocities of individual layers without affecting the velocity discontinuity across the interfaces. It also allows adding interface-related observation such as traveltimes of Ps and PmP in the joint inversion to eliminate the trade-off between interface depth and Vp/Vs ratio and therefore to reduce the uncertainties of results. Numerical tests show that the method is computationally efficient and the inversion results are robust and independent of the initial model. Application of the method to a dense linear array across the Wabash Valley Seismic Zone (WVSZ) produced a high-resolution crustal image in this seismically active region. The results show a 51–55-km-thick crust with a mid-crustal interface at 14–17 km. The crustal Vp/Vs ratio varies from 1.69 to 1.90. There are three pillow-like, ∼100 km apart high-velocity bodies sitting at the base of the crust and directly above each of them are a low-velocity anomaly in the middle crust and a high-velocity anomaly in the upper crust. They are interpreted to be produced by mantle magmatic intrusions and remelting during rifting events in the end of the Precambrian. The current diffuse seismicity in the WVSZ might be rooted in this ancient distributed rifting structure.

SUMMARYSeismic inversion is one of the most commonly used methods in the oil and gas industry for reservoir characterization from observed seismic data. Deep learning (DL) is emerging as a data-driven approach that can effectively solve the inverse problem. However, existing DL-based methods for seismic inversion utilize only seismic data as input, which often leads to poor stability of the inversion results. Besides, it has always been challenging to train a robust network since the real survey has limited labelled data pairs. To partially overcome these issues, we develop a neural network framework with a priori initial model constraint to perform seismic inversion. Our network uses two parts as one input for training. One is the seismic data, and the other is the subsurface background model. The labels for each input are the actual model. The proposed method is performed by log-to-log strategy. The training data set is first generated based on forward modelling. The network is then pre-trained using the synthetic training data set, which is further validated using synthetic data that have not been used in the training step. After obtaining the pre-trained network, we introduce the transfer learning strategy to fine-tune the pre-trained network using labelled data pairs from a real survey to acquire better inversion results in the real survey. The validity of the proposed framework is demonstrated using synthetic 2-D data including both post-stack and pre-stack examples, as well as a real 3-D post-stack seismic data set from the western Canadian sedimentary basin.

SUMMARY In general, the inverse problem of electrical resistivity tomography (ERT) is treated using a deterministic algorithm to find a model of subsurface resistivity that can numerically match the apparent resistivity data acquired at the ground surface and has a smooth distribution that has been introduced as prior information. In this paper, we propose a new deep learning algorithm for processing the 3-D reconstruction of ERT. This approach relies on the approximation of the inverse operator considered as a nonlinear function linking the section of apparent resistivity as input and the underground distribution of electrical resistivity as output. This approximation is performed with a large amount of known data to obtain an accurate generalization of the inverse operator by identifying during the learning process a set of parameters assigned to the neural networks. To train the network, the subsurface resistivity models are theoretically generated by a geostatistical anisotropic Gaussian generator, and their corresponding apparent resistivity by solving numerically 3-D Poisson's equation. These data are formed in a way to have the same size and trained on the convolutional neural networks with SegNet architecture containing a three-level encoder and decoder network ending with a regression layer. The encoders including the convolutional, max-pooling and nonlinear activation operations are sequentially performed to extract the main features of input data in lower resolution maps. On the other side, the decoders are dedicated to upsampling operations in concatenating with feature maps transferred from encoders to compensate the loss of resolution. The tool has been successfully validated on different synthetic cases and with particular attention to how data quality in terms of resolution and noise affects the effectiveness of the approach.

SUMMARY The largest earthquake recorded in South Korea occurred in Gyeongju on 2016 September 12. Since then, understanding the geological structures near the Gyeongju earthquake's epicentre has been demanded. As a preliminary step to understand the geological structures around the epicentre, a 3-D audio-magnetotelluric (AMT) survey was conducted. These 3-D AMT data were imaged through data processing and inversion. However, because the AMT data acquired in Gyeongju were very noisy, conventional imaging procedures did not yield reliable results. To obtain a reliable model inverted from such noisy data, we propose various imaging strategies: an additional data processing technique using the Nyquist diagram after conventional data processing and several inversion strategies related to the selection of data, the weighting of the data, the constraints on the model parameter vector and the Lagrange multiplier used for the regularization. By applying the additional data processing step and several inversion strategies, we were able to successfully invert the noisy field data. The inversion results verify that the data-screening procedure applying the Nyquist diagram remarkably improves the results compared with those obtained using only conventional data processing. The imaging strategies proposed in this case study can be used to image noisy MT/AMT field data for other regions.

SUMMARYUnder favourable circumstances, seismic waveforms corresponding to an ensemble of events related by a common, spatially distributed process collectively exhibit a regular, signal-space geometry. When events in the ensemble have a common, or nearly common, source mechanism, this geometry is a distorted image of the distribution of events in the source region. The signal-space image can be visualized using a relatively simple waveform alignment and projection operation. Ensemble waveform correlation measurements can be inverted to estimate the distribution of the events in the source region, up to an arbitrary rotation, reflection and scaling, with residual distortion. We demonstrate these concepts with synthetic waveforms and with observations of long-wall mining induced seismicity for which substantial ground truth information is available. Our experience with these data has implications for location, correlation detection and machine learning and possible application to studies of repeating events in induced, volcanic and glacial seismicity. Our results place limits on the widely held assumption that waveform correlation is a useful measure of event separation. We suggest that the constraints on event separation need to be evaluated in the context of a population of related events, whose waveforms sample the signal space image of the source region. A better indicator of event separation is the length of the shortest path in signal space along the image.

SUMMARYAmbient noise tomography (ANT) has been considerably used in the last decade in both academic and industrial research. In this work, we propose an innovative technique for ANT based on nonlinear multiscale inversions. Our method relies on a progressive increase in the model parametrization to reduce the nonlinearity of the inverse problem. The developed method is compared with conventional inversion schemes (linear and nonlinear), using different regularization techniques and two different network configurations. The inversion is tested on 22 different synthetic models including classical checkerboard tests. Furthermore, we performed the inversion using real data from a campaign in 2018 at Cumbre Vieja volcano (Canary Islands). The results obtained on both network configurations show an improvement compared to conventional linear and nonlinear inversion schemes, especially when the ray path density is low. This technique does not require expensive computational resources, making it convenient for small-scale industrial applications, especially in the framework of geothermal exploration.

SUMMARY Wave-equation migration velocity analysis (WEMVA) is an image-domain inversion method for velocity model building. Automatic plane-wave WEMVA (PWEMVA) calculates the moveouts of plane-wave common-image gathers (CIGs) by searching a best-fitting parabola with semblance analysis and backprojects residual CIG moveouts into wavefield wave paths with a reflection tomographic kernel. However, 3-D PWEMVA is very computationally expensive because 3-D reflection tomographic inversion requires at least five 3-D reverse-time migrations per iteration and stores two types of source wavefields at model boundaries. We develop a joint inline and crossline PWEMVA method for efficient 3-D velocity model building. We alternatively implement the inline and crossline PWEMVAs with a constraint for each other, in which we iteratively construct the 3-D velocity model update through 1-D spline interpolation of 2-D gradients. The inline and crossline joint inversion is practical since PWEMVA only inverts for low-wavenumber velocity perturbations along wave paths, and the method can take less than 1 per cent of the computational cost of full 3-D PWEMVA. To construct unaliased plane waves for our joint inline and crossline PWEMVA, we develop a 3-D data interpolation method in the frequency–wavenumber (FK) domain to recover regularly and randomly missing traces. The method minimizes the misfit on sufficiently localized data subsets with iterative optimal step lengths and a gradient preconditioner that iteratively selects dominant dips along different azimuths. In numerical experiments, we use a 3-D synthetic seismic data set and a land 3-D field seismic data set acquired at the Farnsworth CO2-EOR (enhanced oil recovery) field to demonstrate the efficacy of our velocity model building and data interpolation methods.

SUMMARYSurface wave tomography is a valuable tool for constraining azimuthal anisotropy at regional scales. However, sparse and uneven coverage of dispersion measurements make meaningful uncertainty estimation challenging, especially when applying subjective model regularization. This paper considers azimuthal anisotropy constrained by measurements of surface wave dispersion data within a Bayesian trans-dimensional (trans-d) tomographic inversion. A recently proposed alternative model parametrization for trans-d inversion is implemented in order to produce more realistic models than previous studies considering trans-d surface wave tomography. The reversible-jump Markov chain Monte Carlo sampling technique is used to numerically estimate the posterior probability density of the model parameters. Isotropic and azimuthally anisotropic components of surface wave group velocity maps (and their associated uncertainties) are estimated while avoiding model regularization and allowing model complexity to be determined by the data information content. Furthermore, data errors are treated as unknown, and solved for within the inversion. The inversion method is applied to measurements of surface wave dispersion from regional earthquakes recorded over northern Cascadia and Haida Gwaii, a region of complex active tectonics but highly heterogeneous station coverage. Results for isotropic group velocity are consistent with previous studies that considered the southern part of the study region over Cascadia. Azimuthal anisotropic fast-axis directions are generally margin-parallel between Vancouver Island and Haida Gwaii, with a small change in direction and magnitude along the margin which may be attributed to the changing tectonic regime (from subduction to transform tectonics). Estimated errors on the dispersion data (solved for within the inversion) reveal a correlation between surface wave period and the dependence of data errors on travel path length. This paper demonstrates the value of considering azimuthal anisotropy within Bayesian tomographic inversions. Furthermore, this work provides structural context for future studies of tectonic structure and dynamics of northern Cascadia and Haida Gwaii, with the aim of improving our understanding of seismic and tsunami hazards.

SUMMARY We revisit the hypocentre–velocity problem, which is of interest in different fields as for example microseismics and seismology. We develop a formulation based on kinematic migration of two picked kinematic attributes in the 2-D case, the traveltime and the slope (horizontal component of the slowness vector), from which we are able to retrieve the location and subsequently the origin time correction and the subsurface parameters mainly velocity. We show how, through a variable projection, the optimization problem boils down to a physically consistent and parsimonious form where the location estimation is projected into the subsurface parameter problem. We present in this study a proof of concept validated by a toy test in two dimensions and a synthetic case study on the Marmousi model. The method presented in this study is extendible to three dimensions by incorporating the crossline slope or the backazimuth as a supplementary attribute.

SUMMARYInstrumental timing and phase errors are a notorious problem in seismic data acquisition and processing. These can be frequency independent, for example due to clock drift, but may also be frequency dependent, for example due to imperfectly known instrument responses. A technique is presented that allows both types of errors to be recovered in a systematic fashion. The methodology relies on the time-symmetry usually inherent in time-averaged cross-correlations of ambient seismic noise: the difference between the arrival time of the direct surface-wave at positive time and the arrival time of the direct surface-wave at negative time is quantified. Doing this for all eligible receiver–receiver pairs of a large-N seismic array, including one or more receivers devoid of instrumental timing errors, the instrumental timing errors of all incorrectly timed receivers can be determined uniquely. Most notably, this is accomplished by means of a weighted least-squares inversion. The weights are based on the receiver–receiver distances and decrease the adverse effect of inhomogeneities in the noise illumination pattern on the recovered instrumental timing errors. Inversion results are furthermore optimized by limiting the inversion to receiver couples that (i) exceed a specific receiver–receiver distance threshold and (ii) whose time-averaged cross-correlations exceed a specific signal-to-noise ratio threshold. Potential frequency dependence of the timing errors is incorporated by means of an iterative, frequency-dependent approach. The proposed methodology is validated using synthetic recordings of ambient seismic surface-wave noise due to an arbitrary non-uniform illumination pattern. The methodology is successfully applied to time-averaged cross-correlations of field recordings of ambient seismic noise on and around the Reykjanes peninsula, SW Iceland.

SUMMARY Probabilistic inversion methods, typically based on Markov chain Monte Carlo, exist that allow exploring the full uncertainty of geophysical inverse problems. The use of such methods is though limited by significant computational demands, and non-trivial analysis of the obtained set of dependent models. Here, a novel approach, for sampling the posterior distribution is suggested based on using pre-calculated lookup tables with the extended rejection sampler. The method is (1) fast, (2) generates independent realizations of the posterior, and (3) does not get stuck in local minima. It can be applied to any inverse problem (and sample an approximate posterior distribution) but is most promising applied to problems with informed prior information and/or localized inverse problems. The method is tested on the inversion of airborne electromagnetic data and shows an increase in the computational efficiency of many orders of magnitude as compared to using the extended Metropolis algorithm.

SUMMARY We use Monte Carlo simulations to explore the effects of earth model uncertainty on the estimation of the seismic source time functions that correspond to the six independent components of the point source seismic moment tensor. Specifically, we invert synthetic data using Green’s functions estimated from a suite of earth models that contain stochastic density and seismic wave-speed heterogeneities. We find that the primary effect of earth model uncertainty on the data is that the amplitude of the first-arriving seismic energy is reduced, and that this amplitude reduction is proportional to the magnitude of the stochastic heterogeneities. Also, we find that the amplitude of the estimated seismic source functions can be under- or overestimated, depending on the stochastic earth model used to create the data. This effect is totally unpredictable, meaning that uncertainty in the earth model can lead to unpredictable biases in the amplitude of the estimated seismic source functions.

We present a new approach to full-waveform inversion (FWI) that enables the assimilation of data sets that expand over time without the need to reinvert all data. This evolutionary inversion rests on a reinterpretation of stochastic Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS), which randomly exploits redundancies to achieve convergence without ever considering the data set as a whole. Specifically for seismological applications, we consider a dynamic mini-batch stochastic L-BFGS, where the size of mini-batches adapts to the number of sources needed to approximate the complete gradient. As an illustration we present an evolutionary FWI for upper-mantle structure beneath Africa. Starting from a 1-D model and data recorded until 1995, we sequentially add contemporary data into an ongoing inversion, showing how (i) new events can be added without compromising convergence, (ii) a consistent measure of misfit can be maintained and (iii) the model evolves over times as a function of data coverage. Though applied retrospectively in this example, our method constitutes a possible approach to the continuous assimilation of seismic data volumes that often tend to grow exponentially., Geophysical Journal International, 224 (1), ISSN:0956-540X, ISSN:1365-246X

Simulation-based probabilistic inversions of 3D magnetotelluric (MT) data are arguably the best option to deal with the non-linearity and non-uniqueness of the MT problem. However, the computational cost associated with the modeling of 3D MT data has so far precluded the community from adopting and/or pursuing full probabilistic inversions of large MT datasets. In this contribution, we present a novel and general inversion framework, driven by Markov chain Monte Carlo (MCMC) algorithms, which combines i) an efficient parallel-in-parallel structure to solve the 3D forward problem, ii) a reduced order technique to create fast and accurate surrogate models of the forward problem, and iii) adaptive strategies for both the MCMC algorithm and the surrogate model. In particular, and contrary to traditional implementations, the adaptation of the surrogate is integrated into the MCMC inversion. This circumvents the need of costly offline stages to build the surrogate and further increases the overall efficiency of the method. We demonstrate the feasibility and performance of our approach to invert for large-scale conductivity structures with two numerical examples using different parameterizations and dimensionalities. In both cases, we report staggering gains in computational efficiency compared to traditional MCMC implementations. Our method finally removes the main bottleneck of probabilistic inversions of 3D MT data and opens up new opportunities for both stand-alone MT inversions and multi-observable joint inversions for the physical state of the Earth’s interior., Facultad de Ciencias Astronómicas y Geofísicas

SUMMARY The ability to quantify structural uncertainty in geological models that incorporate geophysical data is affected by two primary sources of uncertainty: geophysical parameter uncertainty and uncertainty in the relationship between geophysical parameters and geological properties of interest. Here, we introduce an open-source, trans-dimensional Bayesian Markov chain Monte Carlo (McMC) algorithm GeoBIPy—Geophysical Bayesian Inference in Python—for robust uncertainty analysis of time-domain or frequency-domain airborne electromagnetic (AEM) data. The McMC algorithm provides a robust assessment of geophysical parameter uncertainty using a trans-dimensional approach that lets the AEM data inform the level of model complexity necessary by allowing the number of model layers itself to be an unknown parameter. Additional components of the Bayesian algorithm allow the user to solve for parameters such as data errors or corrections to the measured instrument height above ground. Probability distributions for a user-specified number of lithologic classes are developed through posterior clustering of McMC-derived resistivity models. Estimates of geological model structural uncertainty are thus obtained through the joint probability of geophysical parameter uncertainty and the uncertainty in the definition of each class. Examples of the implementation of this algorithm are presented for both time-domain and frequency-domain AEM data acquired in Nebraska, USA.

Author(s): Hetz, G; Datta-Gupta, A; Przybysz-Jarnut, JK; Lopez, JL; Vasco, DW | Abstract: Our limited knowledge of the relationship between changes in the state of an aquifer or reservoir and the corresponding changes in the elastic moduli, that is the rock physics model, hampers the effective use of time-lapse seismic observations for estimating flow properties within the Earth. A central problem is the complicated dependence of the magnitude of time-lapse changes on the saturation, pressure, and temperature changes within an aquifer or reservoir. We describe an inversion methodology for reservoir characterization that uses onset times, the calendar time of the change in seismic attributes, rather than the magnitude of the changes. We find that onset times are much less sensitive than magnitudes to the rock physics model used to relate time-lapse observations to changes in saturation, temperature and fluid pressure. We apply the inversion scheme to observations from daily monitoring of enhanced oil recovery at the Peace River field in Canada. An array of 1492 buried hydrophones record seismic signals from 49 buried sources. Time-shifts for elastic waves traversing the reservoir are extracted from the daily time-lapse cubes. In our analysis 175 images of time-shifts are transformed into a single map of onset times, leading to a substantial reduction in the volume of data. These observations are used in conjunction with bottom hole pressure data to infer the initial conditions prior to the injection, and to update the reservoir permeability model. The combination of a global and local inversion scheme produces a collection of reservoir models that are best described by three clusters. The updated model leads to a nearly 70 percent reduction in seismic data misfit. The final set of solutions successfully predict the observed normalized pressure history during the soak and flow-back into the wells between 82 and 175 days into the cyclic steaming operation.

SUMMARY Bayesian inversion of magnetotelluric (MT) data is a powerful but computationally expensive approach to estimate the subsurface electrical conductivity distribution and associated uncertainty. Approximating the Earth subsurface with 1-D physics considerably speeds-up calculation of the forward problem, making the Bayesian approach tractable, but can lead to biased results when the assumption is violated. We propose a methodology to quantitatively compensate for the bias caused by the 1-D Earth assumption within a 1-D trans-dimensional Markov chain Monte Carlo sampler. Our approach determines site-specific likelihood functions which are calculated using a dimensionality discrepancy error model derived by a machine learning algorithm trained on a set of synthetic 3-D conductivity training images. This is achieved by exploiting known geometrical dimensional properties of the MT phase tensor. A complex synthetic model which mimics a sedimentary basin environment is used to illustrate the ability of our workflow to reliably estimate uncertainty in the inversion results, even in presence of strong 2-D and 3-D effects. Using this dimensionality discrepancy error model we demonstrate that on this synthetic data set the use of our workflow performs better in 80 per cent of the cases compared to the existing practice of using constant errors. Finally, our workflow is benchmarked against real data acquired in Queensland, Australia, and shows its ability to detect the depth to basement accurately.

SUMMARY We present a time-domain matrix-free elastic Gauss–Newton full-waveform inversion (FWI) algorithm. Our algorithm consists of a Gauss–Newton update with a search direction calculated via elastic least-squares reverse time migration (LSRTM). The conjugate gradient least-squares (CGLS) method solves the LSRTM problem with forward and adjoint operators derived via the elastic Born approximation. The Hessian of the Gauss–Newton method is never explicitly formed or saved in memory. In other words, the CGLS algorithm solves for the Gauss–Newton direction via the application of implicit-form forward and adjoint operators which are equivalent to elastic Born modelling and elastic reverse time migration, respectively. We provide numerical examples to test the proposed algorithm where we invert for P- and S-wave velocities simultaneously. The proposed algorithm performs positively on mid-size problems where we report solutions of slight improvement than those computed using the conventional non-linear conjugate gradient method. In spite of the aforementioned limited gain, the theory developed in this paper contributes to a better understanding of time-domain elastic Gauss–Newton FWI.

SUMMARY Seismic attenuation is one of the main factors responsible for degradation of the resolution of seismic data. During seismic wave propagation in attenuation medium, the energy of signal components seriously decreases, especially those with higher frequencies. The seismic attenuation and resolution reduction are generally compensated for with inverse Q filtering in the frequency or time domain. However, the implementation of pre-stack inverse Q filtering is challenging because the traveltime in each layer is not easy to obtain for the pre-stack seismic gather, unless the accurate velocity model is known. In this study, we propose an inverse Q filtering method for the pre-stack seismic gather that uses the local slope and warped mapping to determine the propagation path, and Taylor-expansion-based division is used to stabilize the inversion. The local slope can determine the reflection events with the same ray path, and the inverse warped mapping can transform the attenuation factor from the ${t_0} - p$ (zero-offset traveltime to ray parameter) domain to the $t - x$ (traveltime and offset) domain. The attenuation factor in the ${t_0} - p$ domain is easy to calculate because the traveltimes and Q values in each layer are known. The proposed oriented pre-stack inverse Q filtering method is velocity-independent and suitable for a depth varying Q model. The synthetic and real data examples demonstrated that the method can effectively correct the attenuation and dispersion of seismic waves, and can obtain pre-stack seismic gathers with high resolution.

SUMMARYThe retrieval of high-frequency seismic source–time functions (STFs) of similar earthquakes tends to be an ill-posed problem, causing unstable solutions. This is particularly true when waveforms are complex and band-limited, such as the regional phase Lg. We present a new procedure implementing the multichannel deconvolution (MCD) method to retrieve robust and objective STF solutions. The procedure relies on well-developed geophysical inverse theory to obtain stable STF solutions that jointly minimize the residual misfit, model roughness and data underfitting. MCD is formulated as a least-squares inverse problem with a Tikhonov regularization. The problem is solved using a convex optimization algorithm which rapidly converges to the global minimum while accommodating physical solution constraints including positivity, causality, finiteness and known seismic moments. We construct two L-shaped curves showing how the solution residual and roughness vary with trial solution durations. The optimal damping is chosen when the curves have acceptable levels while exhibiting no oscillations caused by solution instability. The optimal solution duration is chosen to avoid a rapidly decaying segment of the residual curve caused by parameter underfitting. We apply the MCD method to synthetic Lg data constructed by convolving a real Lg waveform with five pairs of simulated STFs. Four pairs consist of single triangular or parabolic pulses. The remaining pair consists of multipulse STFs with a complex, four-spike large STF. Without noise, the larger STFs in all single-pulse cases are well-recovered with Tikhonov regularization. Shape distortions are minor and duration errors are within 5 per cent. The multipulse case is a rare well-posed problem for which the true STFs are recovered without regularization. When a noise of ∼20 per cent is added to the synthetic data, the MCD method retrieves large single-pulse STFs with minor shape distortions and small duration errors (from 0 to 18 per cent). For the multipulse case, the retrieved large STF is overly smeared, losing details in the later portion. The small STF solutions for all cases are less resilient. Finally, we apply the MCD method to Lg data from two pairs of moderate earthquakes in central Asia. The problem becomes more ill-posed owing to lower signal-to-noise ratios (as low as 3) and non-identical Green's functions. A shape constraint of the small STF is needed. For the larger events with M5.7 and 5.8, the retrieved STFs are asymmetric, rising sharply and lasting about 2.0 and 2.5 s. We estimate radiated energies of 2.47 × 1013 and 2.53 × 1013 J and apparent stresses of 1.4 and 1.9 MPa, which are very reasonable. Our results are very consistent with those obtained in a previous study that used a very different, less objective ‘Landweber deconvolution’ method and a pre-fixed small STF duration. Novel improvements made by our new procedure include the application of a convex algorithm rather than a Newton-like method, a procedure for simultaneously optimizing regularization and solution duration parameters, a shape constraint for the smaller STF, and application to the complex Lg wave.

SUMMARY The K-SVD algorithm has been successfully utilized for adaptively learning the sparse dictionary in 2-D seismic denoising. Because of the high computational cost of many singular value decompositions (SVDs) in the K-SVD algorithm, it is not applicable in practical situations, especially in 3-D or 5-D problems. In this paper, I extend the dictionary learning based denoising approach from 2-D to 3-D. To address the computational efficiency problem in K-SVD, I propose a fast dictionary learning approach based on the sequential generalized K-means (SGK) algorithm for denoising multidimensional seismic data. The SGK algorithm updates each dictionary atom by taking an arithmetic average of several training signals instead of calculating an SVD as used in K-SVD algorithm. I summarize the sparse dictionary learning algorithm using K-SVD, and introduce SGK algorithm together with its detailed mathematical implications. 3-D synthetic, 2-D and 3-D field data examples are used to demonstrate the performance of both K-SVD and SGK algorithms. It has been shown that SGK algorithm can significantly increase the computational efficiency while only slightly degrading the denoising performance.

SUMMARY The resistivity imaging method, an effective geophysical technique, has been widely used in environmental, engineering and hydrological fields. The inversion method based on smooth constraint is one of the most commonly used methods. However, this method causes the resistivity to change smoothly and makes it difficult to describe geological boundaries accurately. An accurate description of the target's boundaries often requires a priori information gained with other methods (such as other geophysical methods or geological drilling). To address this issue, a multiscale inversion method is proposed for extracting boundary features and inverting feature parameters from different scales. In this method, a convolution kernel is used to extract the boundary information from the resistivity model. The model parameters are transformed from the spatial domain to the feature domain via a convolutional wavelet transform. The feature parameters of different scales can then be obtained by solving the inversion equation in the feature domain. After that, the resistivity model of the spatial domain is reconverted from the feature domain by deconvolution transform of the inversion result. Numerical simulations and experiments show that the new multiscale resistivity inversion method has the ability to locate and depict boundaries of geological targets with high accuracy.