Search

Our understanding of earthquake sources is limited by the availability and the quality of observations and the fidelity of our physical models. Uncertainties in our physical models will naturally bias our inferences of subsurface fault slip. These uncertainties will always persist to some level as we will never have a perfect knowledge of the Earth's interior. The choice of the forward physics is thus ambiguous, with the frequent need to fix the value of several parameters such as crustal properties or fault geometry. Here, we explore the impact of uncertainties related to the choice of both fault geometry and elastic structure, as applied to the 2016  M w 6.2 Amatrice earthquake, central Italy. This event, well instrumented and characterized by a relatively simple fault morphology, allows us to explore the role of uncertainty in basic fault parameters, such as fault dip and position. We show that introducing uncertainties in fault geometry in a static inversion reduces the sensitivity of inferred models to different geometric assumptions. Accounting for uncertainties thus helps infer more realistic and robust slip models. We also show that uncertainties in fault geometry and Earth's elastic structure significantly impact estimated source models, particularly if near-fault observations are available. [ABSTRACT FROM AUTHOR]

The estimation of quality factor, Q, plays an important role in many geophysical problems, including Q-compensated seismic imaging, geophysical interpretation, fluid characterization, ground motion predictions and seismic hazard estimations. One of the most widely used approaches for estimating Q is the spectral ratio method (SRM). However, the spectral division in SRM may not be stable due to the spectral nulls. The shaping regularized inversion that treats the spectral division as a regularized least-squares inversion problem can help solve the spectral-nulls problem and make the spectral division stable. In the case of very noisy seismic data, the time–frequency maps can not be optimally obtained and thus the Q estimation will be strongly affected and unstable even with the regularized inversion method. We propose a multichannel Q estimation approach that takes advantage of the multichannel spatial coherence to constrain the inversion so that the estimated Q is spatially continuous. We use a set of synthetic and real data examples to demonstrate the performance of the multichannel Q estimation method. Results show that the proposed method can obtain accurate and more importantly stable Q estimation result even in the case of strong random noise. [ABSTRACT FROM AUTHOR]

This proof-of-concept study presents a parameter-free, linear Backus–Gilbert inversion scheme, tractable for seismic tomography problems. It leads to efficient computations of unbiased tomographic images, accompanied by meaningful resolution and uncertainty informations. Moreover, as there is no need to parametrize the model space in this parameter-free approach, it enables numerically accurate data sensitivity kernels to be effectively exploited in tomographic inversions. This is a major benefit over discrete tomographic methods, for which data sensitivity kernels are often inaccurate, as they are projected on a given model parametrization prior to be exploited in the inversion, and these parametrizations are usually coarse to limit the number of parameters and keep tractable the problems of model estimation and/or appraisal. Therefore, this new tomographic scheme fuels great hopes on better constraining multiscale seismic heterogeneities in the Earth's interior by exploiting accurate data sensitivity kernels, that is, taking full advantage of known wave-propagation physics, and enabling quantitative appraisals of tomographic features. As a remark, since its computational cost grows as a function of the total number of data squared, it may be better suited to handle moderate-size data sets, typically encountered in regional-scale tomography. Theoretical developments are illustrated within a finite-frequency physical framework. A set of |$27\, 070$| teleseismic S -wave time residuals is inverted, with focus on imaging and appraising shear-wave velocity anomalies lying in the mantle below Southeast Asia, in the 350–1410 km depth range. [ABSTRACT FROM AUTHOR]

Velocity macromodel building is an essential step of the seismic imaging workflow. Indeed, obtaining acceptable results through migration or full waveform inversion is highly dependent on the kinematic accuracy of the background/initial velocity model. Two decades ago, stereotomography was proposed as an alternative to reflection traveltime tomography, the first relying on semi-automatic picking of locally coherent events associated with small reflection or diffraction segments tied to scatterers in depth by a pair of rays, while the latter on interpretive picking of laterally continuous reflections. The flexibility of stereotomography paved the way for many developments that have shown the efficiency of the method whilst emphasizing on the complementary information carried out by traveltimes and slopes of locally coherent events. A recent formulation recast stereotomography under a matrix-free formulation based on eikonal solvers and the adjoint-state method. In the latter, like in the previous works, the scatterer positions and the velocity field are updated jointly to tackle the ill-famed velocity–position coupling in reflection tomography. Following on from this adjoint-state formulation, we propose a new parsimonious formulation of slope tomography that offers the chance to restrain the problem to minimizing the residuals of a single data class being a slope, in search of a sole parameter class being the subsurface velocity field. This parsimonious formulation results from a variable projection, which is implemented by enforcing a consistency between the scatterer coordinates and the velocity macromodel through migration of kinematic attributes. We explain why the resulting reduced-parametrization inversion is more suitable for tomographic problems than the most common joint inversion strategy. We benchmark our method against the complex Marmousi model along with a validation through time domain full waveform inversion and then present the results of a field data case study. [ABSTRACT FROM AUTHOR]

Non-uniqueness in the geophysical inverse problem is well recognized and so too is the ability to obtain solutions with different character by altering the form of the regularization function. Of particular note is the use of ℓ p norms with p ∈ [0, 2] which gives rise to sparse or smooth models. Most algorithms are designed to implement a single ℓ p norm for the entire model domain. This is not adequate when the fundamental character of the model changes throughout the volume of interest. In such cases we require a generalized regularization function where each sub-volume of the model domain has penalties on smallness and roughness and its own suite of ℓ p parameters. Solving the inverse problem using mixed ℓ p norms in the regularization (especially for p < 1) is computationally challenging. We use the Lawson formulation for the ℓ p norm and solve the optimization problem with Iterative Reweighted Least Squares. The algorithm has two stages; we first solve the l 2-norm problem and then we switch to the desired suite of ℓ p norms; there is one value of p for each term in the objective function. To handle the large changes in numerical values of the regularization function when p values are changed, and to ensure that each component of the regularization is contributing to the final solution, we successively rescale the gradients in our Gauss–Newton solution. An indicator function allows us to evaluate our success in finding a solution in which components of the objective function have been equally influential. We use our algorithm to generate an ensemble of solutions with mixed ℓ p norms. This illuminates some of the non-uniqueness in the inverse problem and helps prevent overinterpretation that can occur by having only one solution. In addition, we use this ensemble to estimate the suite of p values that can be used in a final inversion. First, the most common features of our ensemble are extracted using principal component analysis and edge detection procedures; this provides a reference model. A correlation of each member of the ensemble with the reference model, carried out in a windowed domain, then yields a set of p values for each model cell. The efficacy of our technique is illustrated on a synthetic 2-D cross-well example. We then apply our technique to the field example that motivated this research, the 3-D inversion of magnetic data at a kimberlite site in Canada. Since the final regularization terms have different sets of p values in different regions of model space we are able to recover compact regions associated with the kimberlite intrusions, continuous linear features with sharp edges that are associated with dykes and a background that is relatively smooth. The result has a geologic character that would not have been achievable without the use of spatially variable mixed norms. [ABSTRACT FROM AUTHOR]

We present a method for measuring the dispersion of higher mode surface wave phase velocities from a single seismogram using a hierarchical transdimensional Bayesian approach. The 1-D shear velocity profiles down to 800 km depth between sources and seismic stations are regarded as the controlling parameters to tune the phase velocities of fundamental and higher modes. The misfits between synthetics and real waveforms indicate whether the phase velocities are recovered well from the data. We use Monte Carlo Markov chains to approximate the posterior distribution of each model parameters and assess the uncertainties from these probability density functions. These techniques can test models of varying dimensions while being parsimonious, thereby letting the data themselves control the complexity of the solution. Another advantage is that the algorithm can decide how much data noise is needed in order to explain the data without overfitting them. The data noise can be treated as an unknown and different noise levels can be applied to the different time windows considered. The posterior noise distributions can then be used as an indicator of the quality of the waveform fit within each frequency–time window. We considered phase velocities between 50 and 200 s for each mode and performed a reliability analysis to determine which modes and periods are reliably constrained. In this paper, we first present the method and demonstrate its feasibility with synthetic tests, which show that the technique is robust. We then illustrate it with applications to real data. We applied the method to two paths sampling Australia using earthquakes at regional distances and obtained results that agree well with previous studies. The new method can be used in regional and global tomographic studies to obtain phase velocity maps and 3-D models of seismic velocities and anisotropy at depths that are not well resolved by fundamental mode surface waves or body waves. [ABSTRACT FROM AUTHOR]

In coseismic slip distribution inversion, the Laplacian smoothness constraint is used to avoid the rank deficiency of the coefficient Green matrix and ensure the smoothness among the patches of the fault plane. However, by introducing the classic Laplacian smoothness constraint, the inversion of the maximum slip is commonly underestimated. To better solve this problem, here we propose a new method of weighting to the Laplacian second-order smoothness matrices by the slip. The Laplacian smoothness constraint with unequal weights and the classic Laplacian constraint are used to carry out the simulation experiments. The simulation experiments show that the accuracy of the estimated maximum slip from our method is outperforming the classical Laplacian method by 12–19 per cent. Moreover, the estimations of the average slip and moment magnitude have been enhanced, which improved ranging from 4 to 12.5 per cent and 0.4 to 9 per cent over the inversion results of classic Laplacian smoothness constraint, respectively. In order to further validate the general applicability of the proposed method, we conduct the experiments in terms of the inversion of L'Aquila and Taiwan Meinong earthquakes. The maximum slip inversion results show that the inversion of Laplacian smoothness constraint with unequal weights are larger than that of the classic Laplacian smoothness constraint, which is consistent with the conclusion of simulation experiments. In addition, the parameters inversion results of the actual earthquake from the proposed method are in agreement with other studies which also indicate the feasibility and effectiveness of our method. [ABSTRACT FROM AUTHOR]

Seismic inversion is essentially ill-posed because of inaccurate and insufficient data. Existing methods for improving stability and reducing the size of solution space are usually model-driven. These methods often rely on sophisticated mathematical models associated with prior knowledge or expectation of the structure, distribution and correlation of subsurface parameters. In general, model-driven methods are established only for particular usage situations, which are known to have poor adaptivity. Also, because of their abstraction and simplification, these methods often have difficulty achieving satisfactory accuracy and resolution for complex geology. In this paper, we present a data-driven inversion framework using collaborative sparse representation (CSR) for multiparameter simultaneous inversion. We first let the system learn a joint dictionary from well-log data using multicomponent dictionary learning algorithm. Such a learned dictionary captures structural features and correlations of the elastic parameters. Then, in the inversion process, we replace the constraints of mathematical models with that of CSR over the learned dictionary on model parameters. This multiparameter inversion framework incorporates the correlation among different parameters and is data-driven. It only relies upon the widely accepted assumption that subsurface parameters have a certain degree of consistency and lateral continuity within the same geological formation. Tests on a complex synthetic example and field data demonstrate that our algorithm is effective in improving the resolution and accuracy of solutions, especially for the density parameter that is insensitive to amplitude but extremely important for reservoir identification. [ABSTRACT FROM AUTHOR]

In the absence of quantitative relationships, a cross-gradient function can be used to correlate unknown physical properties in a joint inversion of geophysical data sets. It introduces a structural correlation between properties. A commonly used, inexact approach adds a weighted cross-gradient term as a penalty to the cost function being minimized during inversion. This weighting factor needs to be tuned to balance the regularization and cross-gradient terms. In this paper we propose nonlinear weighting for the cross-gradient function which addresses the very different magnitudes of the cross-gradient and regularization terms. This approach also couples the weighting factors for the regularization and correlation terms reducing the number of tuning parameters. The approach is investigated for a synthetic case. Results are also shown for the 3-D joint inversion of high-resolution magnetic and gravity anomaly data from Southern Queensland in Australia with over 30 million cells. [ABSTRACT FROM AUTHOR]

The Ground-source Airborne Time-domain Electromagnetic (GATEM) method is an efficient electromagnetic prospecting technique for engineering, geological surveying and mineral exploration. The electromagnetic response and inversion are obtained via the electrical conductivity of the geological medium, which is conventionally characterized by a homogeneous half-space or a piecewise smooth spatial distribution. However, these approximations hardly reflect the conditions of real geological media because of their rough, self-similar characteristics. This paper primarily considers the GATEM response and inversion with a long grounded wire source in a rough geological medium based on the fractional diffusion. To manage the rough geological medium, the fractional Maxwell equation is solved in the Laplace domain, and then the solution is transformed to the time domain using Guptasarmas digital filter method. According to the mapping relationship between the resistivity and the induced voltage theoretically calculated in a rough geological medium, the neural network method is used to invert the resistivity. The characteristics of the response in a rough medium are analysed. Then, the inversion results are verified for a layered rough geological medium model; these results are closer to those of the real model than the results of the method derived using a smooth spatial distribution. Finally, the method is applied to GATEM field data in Jincheng, Shanxi Province, China. The inversions performed using a rough geological medium result in a reliable conductivity distribution. [ABSTRACT FROM AUTHOR]

Our physical understanding of earthquakes, along with our ability to forecast them, is hampered by limited indications on the current and future state of stress on faults. Integrating indirect observations, laboratory experiments and physics-based numerical modelling to quantitatively estimate this evolution is crucial. However, quantitative integrations are tenuous in light of the scarcity and uncertainty of observations and the difficulty of modelling the physics governing earthquakes. We show that observations and prior physical knowledge, along with their errors, can be efficiently integrated through the statistical framework of ensemble data assimilation (EDA), which is adopted from weather forecasting. To evaluate whether fault stress estimation and forecasting is possible, we perform a perfect model test in a subduction zone setup that mimicks a scaled laboratory experiment. Synthetic noised data on velocities and stresses from one point near the surface are assimilated using an Ensemble Kalman Filter. These data update the velocity and stress states throughout 150 ensemble members, whose dynamics is governed by a seismic cycle model. This visco-elasto-plastic forward model forecasts the system's evolution through solving Navier–Stokes equations with a strongly rate-dependent friction coefficient. The ensemble assimilation of data from a single location provides probabilistic estimates of fault stress and dynamic strength evolution, which capture the true solution exceptionally well. This is possible, because the sampled error covariance matrix contains prior information from the physics that relates velocities, stresses and pressure at the surface to those at the fault. In the analysis step, this covariance allows stress and strength distributions to be reconstructed. In the subsequent forecast step the physical equations are solved to propagate the updated states forward in time. This provides probabilistic information on the likelihood of occurrence of the next earthquake in this synthetic laboratory setting. Throughout the ensemble simulations the forecasting ability for large, quasi-periodic events turns out to be significantly better than that of a periodic recurrence model. For example, it only requires an alarm to sound for 17 per cent instead of 68 per cent of the time to forecast 70 per cent of 21 events. We show that combining our prior knowledge of physical laws with observations through a Bayesian framework provides distinct added value with respect to using observations or numerical models independently. This educational test thus shows vast potential for including physics-based information into probabilistic seismic hazard assessment using EDA. To analyze its real world potential assumptions on an exact representation of the physics in a 2-D simplified system remain to be explored. [ABSTRACT FROM AUTHOR]

A key aspect of geophysical inversion is the ability to model the earth with a low dimensional representation. There exist various approaches to solve the inverse problem. However, most methods do not automatically adapt inverse model complexity or the number of active model parameters as dictated by data noise and sparse receiver coverage, do not quantify inverse model uncertainty or do not work equally well for 1-D, 2-D or 3-D earth models. Low-frequency electromagnetic (EM) inversion, for example, can require for 3-D problems upwards of 106 cells to forward model. Only a small fraction of these cells are effectively resolvable and there are significant trade-offs between them. To address such problems and get around these limitations we present a novel approach to earth model parametrization by using a Gaussian Processes (GP) machine learning (ML) technique, coupled with a parsimonious Bayesian trans-dimensional (trans-D) Markov chain Monte Carlo sampling scheme. One aspect that sets our approach apart from recent spatial dimension agnostic algorithms in the trans-D or ML literature is the ability to specify inversion property priors directly, as opposed to doing so in a transform domain of the property. We develop the theory, describe the effects of specifying different geological priors and apply the trans-D-GP method to a 1-D controlled source EM and 2-D nonlinear regression problem, using actual field data from the Northwest Australian Shelf for the former. The key advantages in using our method are the simplicity of prior specification, parsimonious low dimensional representations and ease of representing large-scale models in 1-D, 2-D or even 3-D with the same parametrization and computer code. [ABSTRACT FROM AUTHOR]

Schematic diagram of new 3D pulse sequence and the simulated echo trains. • A tri-window pulse sequence to collect the full decaying information of porous media is proposed. • A joint inversion algorithm named 'composite-data-processing' to obtain the 3D correlation map is proposed. • A tight oil shale model containing four types of fluid is established and validated the proposed method. Compared with two-dimensional (2D) nuclear magnetic resonance (NMR) technique like correlations among the transversal relaxation time (T 2), the longitudinal relaxation time (T 1), and the diffusion coefficient correlation (D), three-dimensional (3D) NMR technique is superior with the complete measurement of T 2 , T 1 , and D simultaneously. It can solve the problem of overlaps in 2D correlation map and is helpful to characterize relaxation components in unconventional resources such as tight gas and oil shale. However, the existed 3D NMR technique is restricted due to the loss of short relaxation information and the inversion inaccuracy that caused by the incomplete measurement of the diffusion editing window. We developed a tri-window pulse sequence to collect the full decaying information of porous media. In the first window, the inversion-recovery pulse sequence is applied for T 1 encoding. In the second window, D and T 2 are encoded by an adjustable continuous pulse field gradient and echo spacing (TE). In the last window, CPMG with the shortest TE is used to acquire diffusion-free relaxation information. We then proposed a joint inversion algorithm named "composite-data-processing" to obtain the 3D correlation map. The algorithm adopts the dimension reduction technique and the truncated singular value decomposition (TSVD) to speed up the inversion process and enhance the inversion stability. Numerical simulations show that good estimations of the inversion results are obtained at different signal to noise ratios (SNRs). Our results suggest that the novel pulse sequence and inversion algorithm of 3D NMR can be effectively applied to the exploration of unconventional resources. [ABSTRACT FROM AUTHOR]

Geophysical methods provide remotely sensed data that are sensitive to subsurface properties and interfaces. Knowledge about discontinuities is important throughout the Earth sciences: for example, the saltwater/freshwater interface in coastal areas drive mixing processes; the temporal development of the discontinuity between frozen and unfrozen ground is indicative of permafrost development; and the regolith-bedrock interface often plays a predominant role in both landslide and critical-zone investigations. Accurate detection of subsurface boundaries and their geometry is challenging when using common inversion routines that rely on smoothness constraints that smear out any naturally occurring interfaces. Moreover, uncertainty quantification of interface geometry based on such inversions is very difficult. In this paper, we present a probabilistic formulation and solution to the geophysical inverse problem of inferring interfaces in the presence of significant subsurface heterogeneity. We implement an empirical-Bayes-within-Gibbs formulation that separates the interface and physical property updates within a Markov chain Monte Carlo scheme. Both the interface and the physical properties of the two sub-domains are constrained to favour smooth spatial transitions and pre-defined property bounds. Our methodology is demonstrated on synthetic and actual surface-based electrical resistivity tomography data sets, with the aim of inferring regolith-bedrock interfaces. Even if we are unable to achieve formal convergence of the Markov chains for all model parameters, we demonstrate that the proposed algorithm offers distinct advantages compared to manual- or algorithm-based interface detection using deterministic geophysical tomograms. Moreover, we obtain more reliable estimates of bedrock resistivity and its spatial variations. [ABSTRACT FROM AUTHOR]

Debate is ongoing as to which tectonic model is most consistent with the known geology of southeast Australia, formerly part of the eastern margin of Gondwana. In particular, numerous tectonic models have been proposed to explain the enigmatic geological relationship between Tasmania and the mainland, which is separated by Bass Strait. One of the primary reasons for the lack of certainty is the limited exposure of basement rocks, which are masked by the sea and thick Mesozoic–Cenozoic sedimentary and volcanic cover sequences. We use ambient noise tomography recorded across Bass Strait to generate a new shear wave velocity model in order to investigate crustal structure. Fundamental mode Rayleigh wave phase velocity dispersion data extracted from long-term cross-correlation of ambient noise data are inverted using a transdimensional, hierarchical, Bayesian inversion scheme to produce phase velocity maps in the period range 2–30 s. Subsequent inversion for depth-dependent shear wave velocity structure across a dense grid of points allows a composite 3-D shear wave velocity model to be produced. Benefits of the transdimensional scheme include a data-driven parametrization that allows the number and distribution of velocity unknowns to vary, and the data noise to also be treated as an unknown in the inversion. The new shear wave velocity model clearly reveals the primary sedimentary basins in Bass Strait as slow shear velocity zones which extend down to 14 km in depth. These failed rift basins, which formed during the early stages of Australia–Antarctica break-up, appear to be overlying thinned crust, where high velocities of 3.8–4.0 km s−1 occur at depths greater than 20 km. Along the northern margin of Bass Strait, our new model is consistent with major tectonic boundaries mapped at the surface. In particular, we identify an east dipping velocity transition zone in the vicinity of the Moyston Fault, a major tectonic boundary between the Lachlan and Delamerian orogens, which are part of the Phanerozoic accretionary terrane that makes up eastern Australia. A pronounced lineament of high shear wave velocities (∼3.7–3.8 km s−1) in the lower crust of our new model may represent the signature of relict intrusive magmatism from failed rifting in the early stages of Australia–Antarctica break-up along a crustal scale discontinuity in the Selwyn Block microcontinent which joins Tasmania and Victoria. [ABSTRACT FROM AUTHOR]

High-precision observations of the present-day geomagnetic field by ground-based observatories and satellites provide unprecedented conditions for unveiling the dynamics of the Earth's core. Combining geomagnetic observations with dynamo simulations in a data assimilation (DA) framework allows the reconstruction of past and present states of the internal core dynamics. The essential information that couples the internal state to the observations is provided by the statistical correlations from a numerical dynamo model in the form of a model covariance matrix. Here we test a sequential DA framework, working through a succession of forecast and analysis steps, that extracts the correlations from an ensemble of dynamo models. The primary correlations couple variables of the same azimuthal wave number, reflecting the predominant axial symmetry of the magnetic field. Synthetic tests show that the scheme becomes unstable when confronted with high-precision geomagnetic observations. Our study has identified spurious secondary correlations as the origin of the problem. Keeping only the primary correlations by localizing the covariance matrix with respect to the azimuthal wave number suffices to stabilize the assimilation. While the first analysis step is fundamental in constraining the large-scale interior state, further assimilation steps refine the smaller and more dynamical scales. This refinement turns out to be critical for long-term geomagnetic predictions. Increasing the assimilation steps from one to 18 roughly doubles the prediction horizon for the dipole from about  tree to six centuries, and from 30 to about  60 yr for smaller observable scales. This improvement is also reflected on the predictability of surface intensity features such as the South Atlantic Anomaly. Intensity prediction errors are decreased roughly by a half when assimilating long observation sequences. [ABSTRACT FROM AUTHOR]

The gravity gradient tensor has been increasingly used in practical applications. Among them, how to extract information contained in different gravity gradient components is a challenging problem. Gravity gradient joint inversion is one effective method to solve this problem. We integrate different gravity gradient components in a matrix and then apply them in inversion directly. In this paper, we modify the method to get a new gravity gradient joint inversion (NGGJI). The method is based on the reweighted regularized inversion. We choose one component, for example, g zz, and use the other components to build a weighting matrix. Then we apply the weighting matrix in g zz inversion. We present the method to construct the weighting matrix based on a single component and multiple components. We analyse the characteristics of different weighting matrices and the noise effects on weighting matrices. We compare the inversion results obtained from the conventional gravity gradient joint inversion (CGGJI) with the inversion results obtained from the NGGJI. We conclude that the NGGJI's requirement for memory storage is lower and the resolution of the NGGJI inversion results is higher. We apply the method to survey data from Vinton Salt Dome, Louisiana, USA. The results have proved to be consistent with known geological information. [ABSTRACT FROM AUTHOR]

The slip and static stress drop of the 2015 M w 7.8 Gorkha, Nepal earthquake has been studied using its excellent near-source static observations and a newly developed finite fault inversion algorithm with an average stress drop constraint. A series of nonlinear inversions with different target stress drops are conducted to search for the solution that not only most accurately fits the geodetic data, but also has an energy-based stress drop, |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|⁠, matching the target value. Our result reveals that the misfit to the geodetic data gradually decreases when the |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| of the inverted slip model increases from 2 to 7 MPa; but becomes nearly constant when |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| further increases. Hence, only the lower bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$|⁠, that is, |${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$| (∼7 MPa), of the Gorkha earthquake can be constrained with near-fault geodetic measurements, consistent with our previous study using just far-field seismic data. The upper bound of |${\overline {{\rm{\Delta }}\tau } _\mathrm{ E}}$| can be constrained with an extra constraint on the maximum local stress drop. An artificial upper bound can also be reached if using a large subfault size to represent the source. The lower bound of |$\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}}$| leads to the lower bound of the apparent available energy and the upper bound of the radiation efficiency (⁠|${\eta _\mathrm{ R}}$|⁠, 0.09–0.15), though the latter is also sensitive to the determination of the radiated seismic energy. We find that the lower |${\eta _\mathrm{ R}}$| and the reported high rupture velocity (80–93 per cent shear wave speed) can be reconciled by considering the aspect ratio of the dominant slip patch. We recommend using |${\overline {{\rm{\Delta }}{\tau _\mathrm{ E}}} ^{\mathrm{ min}}}$| to replace various slip smoothing constraints to stabilize the finite fault inversion because of its clear physical meaning. [ABSTRACT FROM AUTHOR]

We present a new model of time-dependent flow at low latitudes in the Earth's core between 2000 and 2018 derived from magnetic field measurements made on board the Swarm and CHAMP satellites and at ground magnetic observatories. The model, called CoreFlo-LL.1, consists of a steady background flow without imposed symmetry plus a time-dependent flow that is dominated by geostrophic and quasi-geostrophic components but also allows weak departures from equatorial symmetry. Core flow mode amplitudes are determined at 4-month intervals by a robust least-squares fit to ground and satellite secular variation data. The l 1 norm of the square root of geostrophic and inertial mode enstrophies, and the l 2 norm of the flow acceleration, are minimized during the inversion procedure. We find that the equatorial region beneath the core–mantle boundary is a place of vigorous, localized, fluid motions; time-dependent flow focused at low latitudes close to the core surface is able to reproduce rapid field variations observed at non-polar latitudes at and above Earth's surface. Magnetic field acceleration pulses are produced by alternating bursts of non-zonal azimuthal flow acceleration in this region. Such bursts are prominent in the longitudinal sectors from 80–130°E and 60–100°W throughout the period studied, but are also evident under the equatorial Pacific from 130°E to 150°W after 2012. We find a distinctive interannual alternation in the sign of the non-zonal azimuthal flow acceleration at some locations involving a rapid crossover between flow acceleration convergence and divergence. Such acceleration sign changes can occur within a year or less and, when the structures involved are of large spatial extent, they can give rise to geomagnetic jerks at the Earth's surface. For example, in 2014, we find a change in the sign of the non-zonal azimuthal flow acceleration under the equatorial Pacific as a region of flow acceleration divergence near 130°E changes to a region of flow acceleration convergence. This occurs at a maximum in the amplitude of the time-varying azimuthal flow under the equatorial Pacific and corresponds to a geomagnetic jerk at the Earth's surface. [ABSTRACT FROM AUTHOR]

Obtaining slip distributions for earthquakes results in an ill-posed inverse problem. While this implies that only limited and uncertain information can be recovered from the data, inferences are typically made based only on a single regularized model. Here, we develop an inversion approach that can quantify uncertainties in a Bayesian probabilistic framework for the finite fault inversion (FFI) problem. The approach is suitably efficient for rapid source characterization and includes positivity constraints for model parameters, a common practice in FFI, via coordinate transformation to logarithmic space. The resulting inverse problem is nonlinear and the most probable solution can be obtained by iterative linearization. In addition, model uncertainties are quantified by approximating the posterior probability distribution by a Gaussian distribution in logarithmic space. This procedure is straightforward since an analytic expression for the Hessian of the objective function is obtained. In addition to positivity, we apply smoothness regularization to the model in logarithmic space. Simulations based on surface wave data show that smoothing in logarithmic space penalizes abrupt slip changes less than smoothing in linear space. Even so, the main slip features of models that are smooth in linear space are recovered well with logarithmic smoothing. Our synthetic experiments also show that, for the data set we consider, uncertainty is low at the shallow portion of the fault and increases with depth. In addition, a simulation with a large station azimuthal gap of 180° significantly increases the slip uncertainties. Further, the marginal posterior probabilities obtained from our approximate method are compared with numerical Markov Chain Monte Carlo sampling. We conclude that the Gaussian approximation is reasonable and meaningful inferences can be obtained from it. Finally, we apply the new approach to observed surface wave records from the great Illapel earthquake (Chile, 2015, M w = 8.3). The location and amplitude of our inferred peak slip is consistent with other published solutions but the spatial slip distribution is more compact, likely because of the logarithmic regularization. We also find a minor slip patch downdip, mainly in an oblique direction, which is poorly resolved compared to the main slip patch and may be an artefact. We conclude that quantifying uncertainties of finite slip models is crucial for their meaningful interpretation, and therefore rapid uncertainty quantification can be critical if such models are to be used for emergency response. [ABSTRACT FROM AUTHOR]

The transient electromagnetic (TEM) data acquired in mountainous area are seriously affected by topography, implying that consideration of topographic effect in the interpretation is in urgent need. To deal with this problem, we propose a new method to invert TEM data by simulating the topography. The unstructured tetrahedral grid with flexibility to delineate the rugged topography is adopted to discretize the inverse model. The transmitting source, including magnetic-loop and finite-length electric dipole source, is divided into short electric dipoles to simulate its practical shape, size and attitude. In the numerical simulation, we adopt the edge-based finite-element method (FEM) to discretize the governing equation in space, and the unconditional stable second-order of backward Euler (BE) scheme to discretize in time. The final forward modelling equation is solved iteratively from an initial field for all time channels. The direct solver MUMPS is used to solve the forward problem for efficiency. The calculation of Jacobian matrix of sensitivities is well addressed, especially for the electric source whose initial electric field for a step-off wave is related to the model parameters. Since the full-channel forward modelling matrix is asymmetric, we use a backward recursive method to solve the adjoint forward problem. Considering the limited power of our computational devices, we apply the limited-memory BFGS optimization method to accomplish our 3-D inversion. The inversion algorithm is firstly tested on two synthetic examples, one with a large fixed-loop transmitting source and the other with a long grounded line source as in semi-airborne TEM survey. The inversion results show that the proposed inversion scheme can invert TEM data with sharp topography. The inversion code is further tested on a ground TEM survey data acquired in Shanxi province, China and the results agree well with the geology and other geophysical information. [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]

We introduce a formalism for estimating local spatial averages of the core–mantle boundary (CMB) radial magnetic field and its time derivatives, based on magnetic field observations collected by low-Earth-orbit satellites. This provides a useful alternative to conventional core field modelling based on global spherical harmonic basis functions, where noise in the polar regions maps into all harmonics, and model regularization and spectral truncation are required. A powerful perspective offered by the proposed technique is formal appraisal of the spatial resolution and variance of the resulting field averages. We use the Green's functions for the Neumann boundary value problem to link the satellite observations to the radial magnetic field on the CMB and estimate field averages using a modified Backus–Gilbert inversion approach. Our approach builds on the Subtractive Optimally Localized Averages (SOLA) method developed in helioseismology, that seeks averaging kernels as close as possible to a chosen target kernel. We are able to account for both internal and external field sources and can incorporate data error covariance information, for example describing along-track serial error correlation. As a proof of concept we present a global map collecting local estimates of the radial main field (MF) constructed on a grid at the CMB with one degree spacing in latitude and longitude, derived from 1 month of three component vector magnetic field data collected by the Swarm satellite trio, using data from dark and geomagnetically quiet times. Using sums and differences of the field components taken along track and in the east–west direction we obtain estimates with spatial resolution kernel widths varying between 18° and 54° depending on the latitude, and a standard deviation of approximately 10 μT (i.e. 5  per cent of the mean CMB field amplitude). The morphology of our CMB radial field map agrees well with results from conventional spherical harmonic field models. In a second application, we determine local estimates of the average rate of change, or secular variation (SV), of the radial field at the CMB, initially considering 2-yr time windows, and performing the analysis on data collected by either the Swarm or CHAMP satellites. We obtain stable local estimates of the SV at the CMB, and present maps of estimates with averaging kernel widths of approximately 42°, 33° and 30° on the equator, with corresponding standard derivations of 0.25, 2.5 and 5 μT yr–1. By subtracting SV estimates constructed at different epochs we are able to calculate the local aggregated secular acceleration (SA) and to study its time changes. Differencing SV estimates 2 yr apart, and considering an averaging kernel width of 42° on the equator, we obtain SA maps similar to those found in the CHAOS-6-x7 field model truncated at SH degree 10. Using our approach we are able to directly control the width of the spatial averaging kernel and the length of the time window, enabling us to directly study the robustness of the inferred SA. Pushing to higher resolution in time, considering 1 yr differences of SV estimates constructed using 1 yr windows, we are able to track the evolution of coherent SA structures in time-longitude plots at the equator. At 25°W in mid 2007 we find a distinctive SA 'cross-over' event, with strong, oppositely signed and adjacent, SA features rapidly changing sign within a year. Our method is well suited for studying such spatio-temporally localized SA events at high resolution; there will be further opportunities for such investigations as the time-series of data provided by the Swarm mission lengthens. [ABSTRACT FROM AUTHOR]

The Reykjanes geothermal system is a high-temperature sea water system situated in SW-Iceland. Interferometric analysis of the Sentinel-1 satellite synthetic aperture radar data has been used to determine a time-series of ground deformation induced by geothermal utilization between April 2015 and October 2017. Surface displacements have been estimated at coherent pixels, indicating a steady and linear subsidence within a subcircular bowl centred on the well field at a maximum near-vertical rate of about 25 mm yr−1, together with horizontal contraction. The average line-of-sight displacement rates from ascending and descending tracks are inverted to determine the characteristics of the deformation source at depth, modelling the geothermal reservoir as a body of simple geometry within an elastic half-space. The results indicate a deformation source at about 1-km depth contracting at a rate of (0.7–0.9) × 105 m3yr−1 during the 2015–2017 period. Using pressure and temperature monitoring data at 900-m depth as well as an analysis of the reservoir structure and rock properties, we infer that the recent estimated volume change can be attributed to processes in the steam cap situated in the topmost part of the geothermal reservoir, in the 800–1200 m depth range. Processes involve a combination of compaction under pressure decrease and/or thermal contraction due to cooling of the rocks within or near the steam cap. The steam cap expanded in response to a sudden pressure drop resulting from the increase in extraction of geothermal fluids for a new power plant in 2006. [ABSTRACT FROM AUTHOR]

Reservoir models are numerical representations of the subsurface petrophysical properties such as porosity, volume of minerals and fluid saturations. These are often derived from elastic models inferred from seismic inversion in a two-step approach: first, seismic reflection data are inverted for the elastic properties of interest (such as density, P -wave and S -wave velocities); these are then used as constraining properties to model the subsurface petrophysical variables. The sequential approach does not ensure a proper propagation of uncertainty throughout the entire geo-modelling workflow as it does not describe a direct link between the observed seismic data and the resulting petrophysical models. Rock physics models link the two domains. We propose to integrate seismic and rock physics modelling into an iterative geostatistical seismic inversion methodology. The proposed method allows the direct inference of the porosity, volume of shale and fluid saturations by simultaneously integrating well-logs, seismic reflection data and rock physics model predictions. Stochastic sequential simulation is used as the perturbation technique of the model parameter space, a calibrated facies-dependent rock physics model links the elastic and the petrophysical domains and a global optimizer based on cross-over genetic algorithms ensures the convergence of the methodology from iteration to iteration. The method is applied to a 3-D volume extracted from a real reservoir data set of a North Sea reservoir and compared to a geostatistical seismic AVA. [ABSTRACT FROM AUTHOR]

Following the linearized attenuation tomography from our previous study (Pejić et al., 2017, https://doi.org/10.1002/2016JB013692), we perform hierarchical transdimensional Bayesian tomography of the upper ≈400 km of the inner core, using 398 globally distributed t∗ estimates. The results are in good agreement with the ones obtained through linearized tomography: they show more complex attenuation pattern than the purely hemispherical one, and the noise estimated from the hierarchical inversion is in good agreement with estimates obtained from the Discrepancy Principle in the previous study. The attenuation pattern we observe gives more weight to the geodynamical models that couple the thermal anomalies of the lowermost mantle to the inner core boundary. Key Points: Transdimensional sampling of Voronoi cells on a spherical surface applied to a problem in global tomographyNo hemispherical pattern is observed in attenuation structure of the upper inner coreThe study provides a strong seismological hint for a strongly attenuating inner core [ABSTRACT FROM AUTHOR]

Scanning magnetometers are increasingly used to characterize the magnetization of mineral grains in rock samples. Up-scaling this measurement technique to large numbers of individual particles is hampered by the intrinsic non-uniqueness of potential-field inversion. Here it is shown that this problem can be circumvented by adding tomographic information that determines the location of the possible field sources. Standard potential theory is used to prove a uniqueness theorem that completely characterizes the mathematical background of the corresponding source-localized inversion. It exactly resolves under which conditions a potential-field measurement on a surface can be uniquely decomposed into signals from the different source regions. The intrinsic non-uniqueness of potential-field inversion prevents that the source distribution inside the tomographically outlined regions can be recovered, but the potential field of each region is uniquely defined. For scanning magnetometers in rock magnetism, this result implies that magnetic dipole vectors of large numbers of individual magnetic particles can be reliably reconstructed from surface scans of the magnetic field, if the particle positions are independently determined. This provides an incentive to improve scanning methods for future palaeomagnetic applications. [ABSTRACT FROM AUTHOR]

It is strongly debated whether the interface between the lithosphere and underlying asthenosphere is a temperature-dependent rheological transition, as expected in a thermal convection system, or additionally affected by the presence of melts and/or fluids. Previous surface wave studies of Pacific oceanic lithosphere have found that shear velocity and azimuthal anisotropy vary with seafloor crustal age as expected for a thermal control; however radial anisotropy does not. Various thermomechanical models have been proposed to explain this disparate behaviour. Nonetheless, it is unclear how robust the surface wave constraints are, and this is what we test in this study. We apply a Bayesian model space search approach to three published Pacific surface wave dispersion data sets, two phase-velocity and one combined phase- and group-velocity set, and determine various proxies for the depth of the lithosphere–asthenosphere boundary (LAB) and their uncertainties based on the velocity and radial anisotropy model distributions obtained. In their overall character and pattern with age, the velocity models from different data sets are consistent with each other, although they differ in their values of LAB depths. Uncertainties are substantial (as much as 20 km on LAB depths) and the addition of group-velocity data does not reduce them. Radial anisotropy structures differ even in pattern and display no obvious age dependence. However, given the uncertainties, we cannot exclude that radial anisotropy, azimuthal anisotropy and velocity models actually reflect compatible, age-dependent, LAB depth estimates. The velocity LAB trends are most like those expected for half-space cooling, because velocity differences persist at old ages, below the depth of common plate cooling models. Any direct signature of sub-ridge melt would be too small-scale to be resolved by these data. However, the velocity-increasing effects of dehydration and depletion due to melting below the ridge could explain why LAB proxy depths tend to a minimum of ∼60 km below young ocean floor. [ABSTRACT FROM AUTHOR]

A gravity survey of the Al-Ain region within the United Arab Emirates (UAE) was designated to compare a regional gravity map with surface structural and geophysical studies of the area. The study area is located in the eastern part of the UAE near the Oman border. Gravity data were collected at 452 stations using a Scintrex CG-6 gravimeter. The gravity measurements define the subsurface density distributions in the study area. The Bouguer anomaly map shows values ranging from –21.5 to 10.8 mGal. High Bouguer anomalies were observed in the eastern and western parts of the study area, in good agreement with the location of known thrust faults. Gravity gradient derivative methods (horizontal gradient, analytic signal and tilt angle) were applied to the corrected gravity data to enhance visibility of subsurface geological features. The gradient derivative methods applied to the gravity data outline subsurface structures trending NNE–SSW, NE–SW and E–W. These trends can be correlated with regional tectonic phases and are in agreement with known surface structural trends. The power spectrum of the gravity data shows two distinct geological layers—a deeper layer from 2 to 6 km depth representing Cretaceous formations, and a shallower layer from 0.5 to 2 km depth representing Tertiary formations. The 3-D gravity inversion results show subsurface density distribution characterized by high density contrast of geological bodies located in eastern and western parts of the study area, representing the Semail Ophiolites and Barzaman Formation, respectively. The Hafit Mountain is characterized by a low density contrast body at near surface (–500 kg m–3) beneath which density increases with depth to +50 kg m–3 at 7 km (b.s.l.). This study is the first of its kind to image the regional subsurface and geological features of Al-Ain region with density distribution from ground surface to 15 km (b.s.l.). [ABSTRACT FROM AUTHOR]

Surface wave inversion based on dispersion curves of higher modes is increasingly being applied to estimate near-surface S -wave velocity structures using active or passive sources. Moreover, sensitivity analysis of dispersion curves is important for estimating the confidence of inverted results. In this study, we define a sensitivity function and a sensitivity kernel, give their analytic expressions and then propose a method to estimate the confidence of inverted results and redesign the initial structure through sensitivity analysis. Three kinds of common layered models [velocity-increasing-with-depth (regular), high-velocity layer (HVL), and low-velocity layer (LVL)] are considered. Results show that dispersion curves are more sensitive to S -wave velocity than P -wave velocity and density. The frequency band of the sensitivity function increases as the modes of dispersion curves increase, and the higher modes can penetrate additional depth. For the HVL model, the dispersion curves are particularly sensitive to those of layers above the HVL. By contrast, for the LVL model, the dispersion curves are particularly sensitive to the LVL. Aside from the analysis of sensitivity, the sensitivity functions can also be used for inversion. Thus, we perform an S -wave velocity inversion for synthetic records calculated from the structural model of Shanghai, China. Synthetic noise records are converted into the frequency–velocity domain (c-f panel) by the vector wavenumber transform method. Then, the dispersion curves of multimode are determined from the peaks of the c-f panel. The inverted structure perfectly matches the target structure except the deeper structures that is consistent with the sensitivity analysis of the determined dispersion curves. Finally, we apply the same procedure to the real case in Linyi, Shandong, China. The first three modes of dispersion curves accurately match the data. The synthetic and real examples indicate that the combined use of higher modes offers additional constraint for the velocity structure, which leads to better results of inversion than that of only the fundamental mode. In addition, the sensitivity analysis offers us a tool to estimate the confidence of the inverted structure and redesign the initial structure. [ABSTRACT FROM AUTHOR]

We present the theory for and applications of Hamiltonian Monte Carlo (HMC) solutions of linear and nonlinear tomographic problems. HMC rests on the construction of an artificial Hamiltonian system where a model is treated as a high-dimensional particle moving along a trajectory in an extended model space. Using derivatives of the forward equations, HMC is able to make long-distance moves from the current towards a new independent model, thereby promoting model independence, while maintaining high acceptance rates. Following a brief introduction to HMC using common geophysical terminology, we study linear (tomographic) problems. Though these may not be the main target of Monte Carlo methods, they provide valuable insight into the geometry and the tuning of HMC, including the design of suitable mass matrices and the length of Hamiltonian trajectories. This is complemented by a self-contained proof of the HMC algorithm in Appendix   A. A series of tomographic/imaging examples is intended to illustrate (i) different variants of HMC, such as constrained and tempered sampling, (ii) the independence of samples produced by the HMC algorithm and (iii) the effects of tuning on the number of samples required to achieve practically useful convergence. Most importantly, we demonstrate the combination of HMC with adjoint techniques. This allows us to solve a fully nonlinear, probabilistic traveltime tomography with several thousand unknowns on a standard laptop computer, without any need for supercomputing resources. [ABSTRACT FROM AUTHOR]

Microseismic monitoring has become a tool of choice for the development and optimization of oil and gas production from unconventional reservoirs and is also used to help assess and reduce the risk of induced seismicity. Recently, relative location methods, which can greatly improve the relative location accuracy among nearby events, become popular and attract the attentions of many researchers. However, how to quickly and accurately capture the location of microseismic events remains a big challenge. We introduce a novel relative event localization technique which utilizes the similar spatial distribution characteristics of the propagation paths of the master event and the adjacent microseismic events to simultaneously achieve quick location estimation and reduce the dependency on a velocity model. The ray paths above the suppositional observation plane are reasonably simplified by employing the propagation paths of a master event to act as the equivalent paths of the adjacent events. Moreover, an iterative technique is used to update the suppositional observation points according to the location result of the previous step and then relocate the event to eliminate the location error caused by the previous hypothesis. In addition, a calibration method is applied to correct the location error caused by an inaccurate velocity model. We test the performance of the proposed technique on both synthetic and field data examples. The results suggest the proposed method significantly improves the location accuracy and reduces the dependency on a velocity model by sufficiently exploiting the master event information. Furthermore, the proposed method is excellent in computational efficiency so that can meet the requirements of real-time locating. In addition, the technique can be beneficial in other fields associated with locating, such as aftershock localization of earthquakes, volcanic earthquakes, etc. [ABSTRACT FROM AUTHOR]

A detailed and precise knowledge of ocean bottom topography is essential in many geoscientific and oceanographic applications. Shipborne echo sounding provides the only direct bathymetric method. However, even after decades of applying this technique only a fraction of the global ocean could be covered. Alternatively, gravity data inversion is a feasible method to infer ocean bottom topography since the gravity field correlates with topography at short to medium wavelengths. Gravity field observables are globally provided by dedicated satellite missions like GOCE, GRACE and GRACE-FO and, over the oceans, by satellite altimetry. Regional and local measurements are realised by means of ground-based, shipborne and airborne gravimetry. For the first time in Europe a jet aircraft was used for airborne gravimetry. During the GEOHALO flight campaign over Italy and the Tyrrhenian, Ionian and Adriatic seas in June 2012, the German research aircraft HALO carried an entire suite of geodetic-geophysical instrumentation, including gravity metres. The careful processing of the gravity data acquired by the CHEKAN-AM instrument of the German Research Centre for Geosciences allowed to achieve an accuracy at the mGal level. Subsequently, the Parker–Oldenburg inversion was applied to predict ocean bottom topography along the GEOHALO profiles flown over the ocean. To constrain the parameter space in the inversion we used the Crust1.0 model. Finally, the obtained results were compared to the General Bathymetric Chart of the Oceans in order to estimate the performance of the method. Our study demonstrates that airborne gravimetry aboard a jet aircraft is capable to provide valuable data for regional geoscientific studies. [ABSTRACT FROM AUTHOR]

The subsurface offset is a possible link connecting wave-equation migration methods, such as reverse-time migration, with the angle-domain. However, it describes an artificial reflection geometry that splits the reflection point between incident and scattered waves. The split configuration contradicts basic principles of continuum mechanics, as it represents action-at-a-distance, and has kinematics distinct from those of ordinary physical reflection. The angle-domain image, computed as a Radon transform of the subsurface offset image, implicitly inherits the split configuration, which distinguishes it from conventional angle-domain decomposition. As a consequence, it shows a dissimilar moveout response to erroneous migration velocity, but still carries valuable information about the migration velocity error. Conventional linearized traveltime inversion techniques, applied for velocity optimization, are most likely to fail while inverting angle-domain images produced from subsurface offset split reflections. These images should be labelled differently to highlight the split mechanism, and input to traveltime inversion algorithms only after a proper generalization. We present a modified linearized traveltime inversion formulation, appropriate for split reflections in the image space. The key ingredient in the algorithm is a subsurface offset depended extension of the traveltime change along migrated split reflection rays. A generalized reflection tomography scheme is proposed accordingly to tie imaging errors, in either the subsurface offset domain or its equivalent angle-domain map, to the migration velocity error. [ABSTRACT FROM AUTHOR]

The purpose of this work is to generalize the frozen Gaussian approximation (FGA) theory to solve the 3-D elastic wave equation and use it as the forward modelling tool for seismic tomography with high-frequency data. FGA has been previously developed and verified as an efficient solver for high-frequency acoustic wave propagation (P wave). The main contribution of this paper consists of three aspects: (1) We derive the FGA formulation for the 3-D elastic wave equation. Rather than standard ray-based methods (e.g. geometric optics and Gaussian beam method), the derivation requires to do asymptotic expansion in the week sense (integral form) so that one is able to perform integration by parts. Compared to the FGA theory for acoustic wave equation, the calculations in the derivation are much more technically involved due to the existence of both P and S waves, and the coupling of the polarized directions for SH and SV waves. In particular, we obtain the diabatic coupling terms for SH and SV waves, with the form closely connecting to the concept of Berry phase that is intensively studied in quantum mechanics and topology (Chern number). The accuracy and parallelizability of the FGA algorithm is illustrated by comparing to the spectral element method for 3-D elastic wave equation in homogeneous media. (2) We derive the interface conditions of FGA for 3-D elastic wave equation based on an Eulerian formulation and the Snell's law. We verify these conditions by simulating high-frequency elastic wave propagation in a 1-D layered Earth model. In this example, we also show that it is natural to apply the FGA algorithm to geometries with non-Cartesian coordinates. (3) We apply the developed FGA algorithm for 3-D seismic wave-equation-based traveltime tomography and full waveform inversion, respectively. [ABSTRACT FROM AUTHOR]

We present a method to explore the effective nullspace of nonlinear inverse problems without Monte Carlo sampling. This is based on the construction of an artificial Hamiltonian system where a model is treated as a high‐dimensional particle. Depending on its initial momentum and mass matrix, the particle evolves along a trajectory that traverses the effective nullspace, thereby producing a series of alternative models that are consistent with observations and their uncertainties. Variants of the nullspace shuttle enable hypothesis testing, for example, by adding features or by producing smoother or rougher models. Furthermore, the Hamiltonian nullspace shuttle can serve as a tunable hybrid between deterministic and probabilistic inversion methods: Choosing random initial momenta, it resembles Hamiltonian Monte Carlo; requiring misfits to decrease along a trajectory, it transforms into gradient descent. We illustrate the concept with a low‐dimensional toy example and with high‐dimensional nonlinear inversions of seismic traveltimes and magnetic data, respectively. Plain Language Summary: All geophysical data are inherently limited in their ability to constrain the properties and dynamics of the Earth, such as material parameters deep inside our planet, future atmospheric circulation, or fluid flow in a geothermal reservoir. As a consequence, there are generally many representations (also called models) of the Earth that explain our data to within their uncertainties. The ensemble of these admissible representations is the "nullspace." Here we present an efficient tool to explore the nullspace. This is based on treating an Earth model as an artificial particle under the influence of gravity. Much like a real planet, the particle travels in an orbit. Each position in the orbit corresponds to an admissible model of the Earth. In addition to developing the theory, we demonstrate the efficiency of our method using examples from traveltime tomography and magnetic surveying. Key Points: We propose a method to explore the nullspace of nonlinear inverse problems without Monte Carlo samplingIt has diverse tuning options for systematic hypothesis testing and for the discovery of alternative models with special predefined propertiesIt represents a hybrid between Hamiltonian Monte Carlo and gradient‐based descent methods [ABSTRACT FROM AUTHOR]

The Santos Basin, with an area of about 350.000 km2, is the largest salt basin of the South Atlantic, and due to its high economic hydrocarbon potential, it is a recurrent theme in scientific studies. The salt structures over the region present great importance for hydrocarbon accumulation and the geological/geophysical studies are performed from seismic reflection data, which requires time and efforts for acquisition and data processing. We identify salt structures using a new workflow based on inversion of residual gravity anomalies, where we use the Moho and basement depths obtained from gravity inversion, followed by the calculation of the gravity residual anomaly, assumed to be representative of the salt structures. This workflow is tested for a geological profile in the Santos Basin, and the results are evaluated along a 2D seismic section tied to well markers. The geometry of the stratified salt obtained from gravity inversion correlates with the seismic interpretation, with the advantage of estimating the entire salt package, including halite and stratified salt. With only seismic data, sometimes the stratified salt can be misinterpreted as sediments. The procedure can be applied to identify salt in sedimentary basins where seismic data is unavailable or of low quality. [ABSTRACT FROM AUTHOR]

Injection of CO2 into tight reservoirs produces both gravity change and ground deformation, which provide great opportunities for more accurate coupled inverse modelling. In this study, we incorporate signals generated from several synthetic models to estimate the CO2 distribution in the reservoir. A relationship is found that connects density variations to volumetric changes associated with injected CO2, taking advantage of a common set of model parameters for both gravitational and geo-mechanical inverse modelling. This is achieved by assuming that the injected CO2 increases pressure in the reservoir, which in turn generates extra porosity that is then filled in by the CO2 mass. Tikhonov regularization, supported by the Generalized Cross Validation (GCV) technique for finding the optimized model, is used to solve the ill-posed inverse problems. The results indicate that the uncertainty and ambiguity in gravimetric modelling due to high levels of noise is mitigated by implementing highly accurate ground deformation measurements, which normally have a higher signal to noise ratio. [ABSTRACT FROM AUTHOR]

Receiver function analysis is widely used to make quantitative inferences about the structure below a seismic station. As these observables are mainly sensitive to traveltimes of phases converted and reflected at seismic discontinuities, the resulting inverse problem is highly non-linear, the solution non-unique, and there are strong trade-offs between the depth of discontinuities and absolute velocities. To overcome this difficulty, we propose to measure the misfit between the predicted and observed data with an optimal transport distance instead of the conventional least-squares distance, a strategy that has shown its assets in the context of full waveform inversion. This approach views a seismogram as a distribution of 'mass'. The optimal transport distance between two waveforms is the minimal cost of transporting one waveform onto the other. We test the optimal transport approach on the inversion of a radial P -wave receiver function. We also show how it can be applied to measure the cross-convolution distance between the radial and vertical components, thus avoiding the need for deconvolution associated with the calculation of the receiver function. The resulting misfit function is minimized with a local optimization algorithm to constrain the receiver-side structure. The benefits of this methodology are studied in simple synthetic tests and with real data. In particular, we show that with its increased sensibility to time-shifts, the optimal transport distance reduces the number of local minima in the misfit function, which, in the case of a linearized inversion, significantly reduces the dependency to the starting model and results in a better convergence towards the solution model. A joint inversion of the P -wave receiver function and surface wave dispersion curves is performed at the Hyderabad station in India. [ABSTRACT FROM AUTHOR]

This paper focuses on the combined analysis and interpretation of controlled source electromagnetic (CSEM) and multichannel reflection seismic (MCS) data along one profile in the German North Sea with the goal to reduce ambiguities in interpretation. The shallow water environment of the North Sea is characterized by a complex geological development which includes rifting, several ice age cycles, a propagating shelf margin, mass-transport deposits and salt dome formation. Seismic and electromagnetic methods are sensitive to different physical properties of the seabed and therefore complement each other. We analyse the MCS data with a migration velocity tomography and an amplitude variation with offset analysis and discuss seismic velocities and densities. For true amplitude recovery the amplitude distortions are calibrated with in situ logging data. The CSEM data are analysed in 2-D, for which, for the first time, data were included that were acquired while the instrument was towed on the seafloor in addition to the stationary sites. The CSEM inversions are constrained by seismic horizons. The joint interpretation focuses on two seismic reflectors: One can be interpreted as an unconformity marking a lithological change from fresh water-bearing glacial deposits to compacted sediments below, and the other one as a layer of fine-grained deposits potentially capping patchy shallow gas occurrences. This exemplary case study shows how the combination of both methods can benefit by interpreting complex geology and eliminating ambiguous explanations. [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]

This paper proposes a novel approach to the joint inversion of gravity and magnetic data based on multinary transformation of the model parameters and Gramian constraints. The concept of multinary transformation is a generalization of the binary density inversion for models described by any number of discrete model parameters, while the Gramian constraint enforces the linear relationships between the different model parameters and their attributes or transforms. By combining these two concepts, the joint multinary inversion using Gramian constraints not only makes it possible to explicitly exploit the sharp contrasts of the density and magnetic susceptibility between the host media and anomalous targets in the inversion of gravity and magnetic data, but also provides uniform spatial boundaries of the anomalous targets in the distributions of density and magnetic susceptibility. The novel joint multinary inversion algorithm is demonstrated to be effective in determining the shapes, locations and physical properties of the anomalous targets. We show that this method can be effectively applied to the joint inversion of the full tensor gravity gradiometry data and the total magnetic intensity data, computer simulated for 3-D synthetic models using ternary model transformation. We also tested this joint inversion algorithm in the area of the McFaulds Lake in northwestern Ontario, Canada, and the joint inversion results provide a reasonable geological model with high resolution for the exploration of magmatic chromite deposits. [ABSTRACT FROM AUTHOR]

We investigate the Marchenko redatuming and imaging for carbon sequestration monitoring using synthetic data and observed data from the Frio carbon sequestration experiment. The Marchenko approach has the potential of providing enhanced imaging of the subsurface by fully accounting for multiple scattering. We first present two synthetic examples, with the second using a velocity model derived from the first-arrival traveltimes from the Frio experiment. We next apply the Marchenko redatuming and imaging to the observed data from the Frio experiment and compare this with conventional vertical seismic profile (VSP) corridor stacks. We then apply the Marchenko redatuming and imaging to the near injection zone with the pre- and post-injection data and compare this with conventional processing. In each case, a good agreement is found. By using one-sided illumination and recording, the Marchenko redatuming and imaging provides a potentially cheaper alternative to conventional VSP processing for monitoring CO2injection. [ABSTRACT FROM AUTHOR]

Conventional seismic techniques for detecting the subsurface geologic features are challenged by limited data coverage, computational inefficiency and subjective human factors. We developed a novel data-driven geological feature characterization approach based on pre-stack seismic measurements. Our characterization method employs an efficient and accurate machine learning method to extract useful subsurface geologic features automatically. Specifically, we use kernel ridge regression to account for the nonlinear relationship between seismic data and geological features. We further employ kernel tricks to avoid the explicit nonlinear mapping and infinite dimension of feature space. However, the conventional kernel ridge regression can be computationally prohibitive because of the large volume of seismic measurements. We employ a data reduction technique in combination with the conventional kernel ridge regression method to improve the computational efficiency and reduce memory usage. In particular, we utilize a randomized numerical linear algebra technique, named Nyström method, to effectively reduce the dimensionality of the feature space without compromising the information content required for accurate characterization. We provide thorough computational cost analysis to show the efficiency of our new geological feature characterization methods. We validate the performance of our method in characterizing geologic fault zones because faults play an important role in various subsurface applications. Our numerical examples demonstrate that our new characterization method significantly improves the computational efficiency while maintaining comparable accuracy. Interestingly, we show that our method yields a speedup ratio on the order of ∼102to ∼103in a multicore computational environment. [ABSTRACT FROM AUTHOR]

Ultrafine grains of magnetic minerals provide reliable recordings of both naturally occurring and anthropogenically generated particulate matter in polluted air; magnetic data can be used to understand biogenic iron-cycling in anaerobic environments, as well as pedogenesis and palaeoclimate studies of loess soils. The ultrafine fraction is produced under specific conditions and can be easily recognized by its superparamagnetic (SP) behaviour. Many proxies have been proposed to account for the SP contribution by measuring its susceptibility dependency with frequency (frequency effect) or the magnetization loss after removing an external inducing field. Here we introduce the Superparamagnetic Concentration and Dipole Moment (SPCDM) procedure for quantitative interpretation of SP magnetization. This procedure is well suited to SP carriers with a fast magnetization decay (<1 s), as would be expected for magnetic minerals with a grain size distribution lying below the blocking volume for stable, single-domain (SD) magnetization. SPCDM requires a dedicated experimental procedure to isolate the SP response from the paramagnetic and remnant effects, as observed in samples with mixed contributions. The proposed technique was tested using synthetic, nanoparticles of magnetite and then to characterize the magnetic properties of air particulate matter (PM) sampled at Jânio Quadros tunnel in São Paulo, Brazil. For nano-sized magnetite, SPCDM estimates for dipole moment are invariable with mass concentration and consistent with the published results; estimates for particle concentration are strongly correlated with true mass concentration (R 2= 0.96). For air PM, SPCDM estimates a particle size with a diameter of 7.7 ± 0.1 nm, a kind of ultrafine magnetic material not previously recognized in air pollutants. [ABSTRACT FROM AUTHOR]

One of the main goals of archeomagnetism is to document the secular changes of Earth's magnetic field by laboratory analysis of the magnetization carried by archeological artefacts. Typical techniques for creating a time-dependent model assume a prescribed temporal discretization which, when coupled with sparse data coverage, require strong regularization generally applied over the entire time-series in order to ensure smoothness. Such techniques make it difficult to characterize uncertainty and frequency content, and robustly detect rapid changes. Key to proper modelling (and physical understanding) is a method that places a minimum level of regularization on any fit to the data. Here we apply a transdimensional Bayesian technique based on piecewise linear interpolation to sparse archeointensity data sets, in which the temporal complexity of the model is not set a priori, but is self-selected by the data. The method produces two key outputs: (i) a posterior distribution of intensity as a function of time, a useful tool for archeomagnetic dating, whose statistics are smooth but formally unregularized and (ii) by including the data ages in the model of unknown parameters, the method also produces posterior age statistics of each individual contributing datum. We test the technique using synthetic data sets and confirm agreement of our method with an integrated likelihood approach. We then apply the method to three archeomagnetic data sets all reduced to a single location: one temporally well-sampled within 700 km from Paris (here referred to as Paris700), one that is temporally sparse centred on Hawaii, and a third (from Lübeck, Germany and Paris700) that has additional ordering constraints on age from stratification. Compared with other methods, our average posterior distributions largely agree, however our credible intervals appear to much better reflect the uncertainty during periods of sparse data coverage. Because each ensemble member of the posterior distribution is piecewise linear, we only fit oscillations when required by the data. As an example, we show that an oscillatory signal, associated with temporally localized intensity maxima reported for a sparse Hawaiian data set, is not required by the data. However, we do recover the previously reported oscillation of period 260 yr for the Paris700 data set and compute the probability distribution of the period of oscillation. We further demonstrate that such an oscillation is unresolved when accounting for age uncertainty by using a fixed age and with an artificially inflated error budget on intensity. [ABSTRACT FROM AUTHOR]

We propose a new wavelet-based compressive sensing method to image coseismic radiation sources with teleseismic array data. Compared to the sliding-window-based technique, our method avoids choosing the window length by taking advantages of the adaptive resolution feature of the wavelet transform. Thus artefacts caused by inappropriate window length and waveform truncation can be removed. As waveforms usually become less aligned at late arrivals due to the increasing distance of actual radiation sources with respect to the hypocentre, we introduce a time-correction scheme to realign waveforms with the reference point obtained from beamforming at every time step. Synthetic tests indicate that the new method can better resolve the seismic radiation sources with more reliable locations and source times. To demonstrate its feasibility, we apply the method to the M w9.0 Tohoku earthquake. We find similar spatial rupture and seismic radiation patterns, such as repeating rupture in the region close to the hypocentre and along-dip frequency-dependent radiation, as indicated by previous models from finite-fault inversion and back projection. In addition, our results show structure-controlled rupture characteristics, with the seismic radiation in the first 80 s mainly confined in the off-Miyagi high-wave-speed region. The rupture speed shows clear along-dip variations, with updip rupture slower than that of downdip, which corresponds well with the depth-varying characteristics along the megathrust. [ABSTRACT FROM AUTHOR]

Structural inversion of gravity datasets based on the use of density anomalies to derive robust images of the subsurface (delineating lithologies and their boundaries) constitutes a fundamental non-invasive tool for geological exploration. The use of experimental techniques in geophysics to estimate and interpret di erences in the substructure based on its density properties have proven e cient; however, the inherent non-uniqueness associated with most geophysical datasets make this the ideal scenario for the use of recently developed robust constrained optimization techniques. We present a constrained optimization approach for a least squares inversion problem aimed to characterize 2-Dimensional Earth density structure models based on Bouguer gravity anomalies. The proposed formulation is solved with a Primal-Dual Interior-Point method including equality and inequality physical and structural constraints. We validate our results using synthetic density crustal structure models with varying complexity and illustrate the behavior of the algorithm using di erent initial density structure models and increasing noise levels in the observations. Based on these implementations, we conclude that the algorithm using Primal-Dual Interior-Point methods is robust, and its results always honor the geophysical constraints. Some of the advantages of using this approach for structural inversion of gravity data are the incorporation of a priori information related to the model parameters (coming from actual physical properties of the subsurface) and the reduction of the solution space contingent on these boundary conditions.