Search

In velocity model building using full waveform inversion, diving waves are a source of the low-frequency information needed to update the low- to intermediate wavenumber components in the velocity model. Because full waveform inversion is a highly nonlinear optimization problem attempting to fit the information from the whole wavefield, it is prone to get stuck in local minimums and experience cycle skipping if the initial model does not contain accurate enough information about the low-wavenumbers to constrain the kinematics. Such initial velocity models will typically be estimated using alternative inversion methods that only apply information from parts of the wavefield in the forward modelling and recorded data. In migration velocity analysis, the unfocusing of the image, caused by inconsistency in the kinematics of reflected waves, is measured and is used to update the velocity model by applying differential semblance optimization. The same concept is applicable when using information from diving waves. For diving waves, little is known about the focusing behaviour, also known as imaging moveout, especially in the presence of azimuthal anisotropy, such as orthorhombic anisotropy, which is a good approximation for volcanic plains and sedimentary basins. Obtaining more knowledge about the focusing behaviour of diving waves for this type of anisotropy may help when choosing parameters to invert during inversion that yield a low trade-off between the parameters. Here, we derive approximate series-based equations describing the imaging moveout of diving qP waves in acoustic factorized orthorhombic anisotropic media. Due to the complexity of the analytical solutions for diving wave properties in such models, we find a vertical transverse isotropic type of simplification that manages to conserve the curvature and horizontal velocity of the diving qP wave. The accuracy of the simplified expressions is compared to the exact form and found to yield small errors. To verify the validity of the derived equations for imaging moveout, we compare the generated surfaces with images containing the shape of the unfocusing obtained using reverse time migration, showing a good fit. Of the anisotropic parameters defining the acoustic factorized orthorhombic model, the |$\epsilon _1$| and |$\epsilon _2$| parameters have the largest influence on the shape of the surface. The anelliptical parameters are shown only to induce minimal changes. Our method to obtain the series expressions describing the imaging moveout in the acoustic factorized orthorhombic model is also applicable to higher symmetry special cases of orthorhombic anisotropy by setting appropriate conditions. Here, this is demonstrated by additionally finding an approximate series expression describing the imaging moveout in acoustic factorized vertical transversely isotropic media, yielding a good match with the image obtained from reverse time migration. [ABSTRACT FROM AUTHOR]

Complex oil and gas reservoirs represented by orthorhombic shale usually have the characteristics of strong anisotropy, wide distribution of non-welded interfaces (i.e. particular boundaries caused by horizontal bedding or fracture development), and ubiquitous initial stress effects. The knowledge of their seismic response characteristics plays an essential role in hydrocarbon exploration in stressed shale reservoirs. However, seismic wave reflection in the orthorhombic media under initial stress remains unclear. To address this issue, we used the theories of acoustoelasticity and elastic anisotropy to derive a stress-dependent effective elastic stiffness tensor. At the same time, we obtain the wave velocity and polarization characteristics directly by solving the Christoffel equation. According to the linear slip theory, we further constructed the exact reflection and transmission coefficient equations for the non-welded interface between isotropic and orthorhombic half-spaces under the initial stress. The effects of the non-welded interface on seismic wave velocity and reflection and transmission coefficients were systematically analysed. Meanwhile, we characterize the law of P -wave reflection coefficient effected by elastic coefficients. Our equations and results potentially lay an equation foundation for orthorhombic reservoirs in high-stress fields and are essential in orthorhombic reservoirs in high-stress geophysical exploration. [ABSTRACT FROM AUTHOR]

The columnar-flow approximation allows the development of computationally efficient numerical models tailored to the study of the rapidly rotating dynamics of Earth's fluid outer core. In this paper, we extend a novel columnar-flow formulation, called Plesio-Geostrophy (PG) by including thermal effects and viscous boundary conditions. The effect of both no-slip and stress-free boundaries, the latter being a novelty for columnar-flow models, are included. We obtain a set of fully 2-D evolution equations for fluid flows and temperature where no assumption is made regarding the geometry of the latter, except in the derivation of an approximate thermal diffusion operator. To test the new PG implementation, we calculated the critical parameters for onset of thermal convection in a spherical domain. We found that the PG model prediction is in better agreement with unapproximated, 3-D calculations in rapidly rotating regimes, compared to another state-of-the-art columnar-flow model. [ABSTRACT FROM AUTHOR]

SUMMARY: The successful recovery of hydrogeophysical parameters through surface-NMR measurements depends on the quality of the signal, which can be significantly degraded by harmonics from multiple noise sources with different fundamental frequencies in urban areas. Accurate estimation of the fundamental frequencies of harmonics is the main step in harmonic noise cancellation-based methods. The existing 1-D and 2-D model-based approaches involve a computationally expensive process that sets limits for processing of large surface-NMR data sets. In addition, the classical Nyman, Gaiser and Saucier estimation (NGSE) algorithm, despite its fast implementation, may not accurately recover harmonic components when there is no prior knowledge of the expected value of the frequency offset between the true fundamental frequencies and their nominal values. This lack of knowledge can make it difficult to accurately estimate the maximum number of harmonics and, consequently, result in an incorrect recovery of the fundamental frequency. To surmount these limitations, we propose an enhanced version of the NGSE approach based on an efficient maximum number of harmonics search approach to process surface-NMR signals corrupted by powerline harmonics with both single and multiple frequency content. We verify the efficiency of our algorithm on a synthetic data set embedded in simulated powerline harmonic signals, and real electromagnetic noise recordings, as well as a real surface-NMR data set. Our numerical experiments confirm that the proposed algorithm can retrieve the multiple fundamental frequencies simultaneously with a significant speedup ranging from 4 to 87 times, depending on whether the signal has single, dual, or triple frequency content, in the overall computation time compared to the model-based methods. [ABSTRACT FROM AUTHOR]

Full-waveform inversion (FWI) is an effective method for imaging subsurface properties using sparsely recorded data. It involves solving a wave propagation problem to estimate model parameters that accurately reproduce the data. Recent trends in FWI have seen a renewed interest in extended methodologies, among which source extension methods leveraging reconstructed wavefields to solve penalty or augmented Lagrangian (AL) formulations have emerged as robust algorithms, even for inaccurate initial models. Despite their demonstrated robustness on synthetic data, challenges remain, such as the lack of a clear physical interpretation and reliance on difficult-to-compute least-squares (LS) wavefields. Moreover, the literature lacks a general and through comparison of these source extension methods with each other and with the standard FWI. This paper is divided into three critical parts. In the first, a novel formulation of these methods is explored within a unified Lagrangian framework. This novel perspective permits the introduction of alternative algorithms that use LS multipliers instead of wavefields. These multiplier-oriented variants appear as regularizations of the standard FWI, are suitable to the time domain, offer tangible physical interpretations, and foster enhanced convergence efficiency. The second part of the paper delves into understanding the underlying mechanisms of these techniques. This is achieved by solving the associated non-linear equations using iterative linearization and inverse scattering methods. The paper provides insight into the role and significance of Lagrange multipliers in enhancing the linearization of the equations. It explains how different methods estimate multipliers or make approximations to increase computing efficiency. Additionally, it presents a new physical understanding of the Lagrange multiplier used in the AL method, highlighting how important it is for improving algorithm performance when compared to penalty methods. In the final section, the paper presents numerical examples that compare different methods within a unified iterative algorithm, utilizing benchmark Marmousi and 2004 BP salt models. [ABSTRACT FROM AUTHOR]

The pseudo-differential operator with symbol | k |α has been widely used in seismic modelling and imaging when involving attenuation, anisotropy and one-way wave equation, which is usually calculated using the pseudo-spectral method. For large-scale problems, applying high-dimensional Fourier transforms to solve the wave equation that includes pseudo-differential operators is much more expensive than finite-difference approaches, and it is not suitable for parallel computing with domain decomposition. To mitigate this difficulty, we present a truncated space-domain convolution method to efficiently compute the pseudo-differential operator |$\sqrt{-\nabla ^2}$|⁠ , and then apply it to viscoacoustic reverse-time migration. Although |$\sqrt{-\nabla ^2}$| is theoretically non-local in the space domain, we take the limited frequency band of seismic data into account, and constrain the approximated convolution stencil to a finite length. The convolution coefficients are computed by solving a least-squares inverse problem in the wavenumber domain. In addition, we exploit the symmetry of the resulting convolution stencil and develop a fast spatial convolution algorithm. The applications of the proposed method in Q -compensated reverse-time migration demonstrate that it is a good alternative to the pseudo-spectral method for computing the pseudo-differential operator |$\sqrt{-\nabla ^2}$|⁠ , with almost the same accuracy but much higher efficiency. [ABSTRACT FROM AUTHOR]

Spectral induced polarization (SIP) and electrical impedance tomography (EIT) are advanced approaches for structural characterisation of soils and rocks based on measurements of the frequency-dependent complex electrical resistivity. Accurate SIP and EIT measurements require the reduction of electrode effects, which have often been attributed to the polarization of metal electrodes in the current pathway. In this study, we demonstrate an alternative electrode effect that influences the accuracy of complex impedance measurements. This alternative electrode effect is caused by a heterogeneous distribution of the contact impedance, which can occur due to ageing or corrosion of the electrode. First, theoretical considerations are given to explain how this electrode effect influences the accuracy of complex impedance measurements. Next, we demonstrate using SIP and EIT measurements that this alternative mechanism better explains the observed electrode effects than previous explanations based only on the presence of metal. The SIP measurements were made on a water-filled container with a metal sphere between two non-polarizing electrodes, which showed that the mere presence of metal in the current pathway only leads to small electrode effects that are typically well below 1 mrad even for small electrode spacings. When the same metal sphere was used as a potential electrode, the observed electrode effect was about 2.5 times larger and varied with rotation. These changes can be explained by a heterogeneous contact impedance resulting in a varying electrode voltage. The EIT measurements were performed with a cylindrical sample holder filled with water. They are used to illustrate how the presented theoretical considerations can be used to investigate the magnitude of the electrode effects for different electrode configurations. We conclude that electrode effects are important to consider for laboratory SIP and EIT measurements and field measurement with closely spaced electrodes. [ABSTRACT FROM AUTHOR]

The standard-linear-solid (SLS) theory works well for viscoelasticity. However, the coupling of amplitude dissipation and phase dispersion makes it impossible to investigate their effects separately by the discontinuous Galerkin method (DGM). In this paper, we have derived newly viscoacoustic wave equations with amplitude dissipation and phase dispersion separation in isotropic and anisotropic media, based on a Fourier method, which is suitable for using a time–space-domain DGM on unstructured meshes. The basic framework of DGM is constructed and the amplitude-dissipation effect and the phase-dispersion effect in viscoacoustic wave equations are investigated. The original equation is first transformed into the frequency–wavenumber domain, where the amplitude dissipation and phase dispersion are separated effectively, and then the decoupled formulation is converted back to the time–space domain. The new equivalent and approximate equations can be obtained. Compared with the original equation, the newly approximated equation enables us to separate the amplitude-loss and phase-delay terms, respectively, and experiences four kinds of effects, namely acoustic effect, only amplitude-dissipation effect, only phase-dispersion effect and both amplitude-dissipation and phase-dispersion effect. Moreover, the stability condition and numerical dispersion for using DGM to solve the new and old equations are presented. Several numerical examples are used to verify the correctness and effectiveness of the modified approximate equations in viscoacoustic isotropic and anisotropic media. The numerical results in a cave and SEG/EAGE salt models demonstrate that the new equations combined with DGM have performances on viscoacoustic media with complex geological structures. [ABSTRACT FROM AUTHOR]

Deep boreholes represent a source of clean energy. Therefore, effective calculations of potential extraction of heat from boreholes for realistic models of the Earth's crust with variable thermal conductivity and diffusivity are needed. We deal with heat extraction in a quasi-steady state from coaxial boreholes where downward and upward flows of pumped fluid (water) are separated by an inner pipe and connected only at the bottom. We first obtain theoretical estimates of heat extraction for a thermally isolated inner pipe and a model of the ground with constant thermal diffusivity and conductivity. Then, we develop a new analytical matrix method for a general layered ground model that enables us to include depth-dependent ground properties as well as heat exchange between the downward and upward flows of fluid in the borehole. Our straightforward and fast approach is thus suitable for various parametric studies or as a tool for benchmarks of numerical software. A key role in heat extraction from coaxial boreholes is played by the inner-pipe thermal resistance. We apply our method to the parametric study showing the dependence of pumped water temperature and total heat extraction from the borehole on realistic borehole geometries under different amounts of water pumping. The calculations are performed for a 3 km deep borehole as the representative of present deep boreholes used for extraction of geothermal energy and for a 10 km deep borehole. Drilling of such a  superdeep borehole has just started in China and our results demonstrate potential limits of geothermal energy extraction from such great depths. [ABSTRACT FROM AUTHOR]

Based on the first-principles method combined with Special Quasi-random Structures package, we investigated the effects of impurities distributions on the perovskite (Pv)–post-perovskite (PPv) phase transition pressures of Fe- and Al-bearing MgSiO3. 61 different distribution configurations were calculated. The results show that the distribution has marginal effect on the Pv–PPv phase transition pressure of Fe2+-bearing MgSiO3. However, the separate of Al3+Al3+ ions can significantly increase the Pv–PPv phase transition pressure to that higher than the core–mantle boundary (CMB) pressure. In the most disordered structures of Al3+Al3+-bearing case constructed by SQS, the Pv–PPv phase transition pressure is higher than 400 GPa. Based on our previous result, we proposed that only Fe2+-bearing MgSiO3 can induce the discontinuous increase in shear wave velocity at the CMB, while Al3+Al3+-, Fe3+Al3+- and Fe3+Fe3+-bearing MgSiO3 cannot. [ABSTRACT FROM AUTHOR]

Analysis of azimuthally varying non-hyperbolic moveout of P -wave traveltimes in orthotropic media can provide valuable information for reservoir fracture characterization and seismic data processing. Moveout velocities and azimuthal orientation of the two vertical symmetry planes plus anellipticity parameters defined for the three symmetry planes of an orthotropic medium can be estimated from the inversion of P -wave traveltimes measured in long-spread, wide-azimuth seismic surveys. Several orthotropic moveout equations have been developed to this date that can be used in the traveltime data inversion. The accuracy of the inversion results is primarily determined by the complexity of the equation and the assumptions (i.e. elastic or acoustic) made to derive it. Here, we show that often simpler moveout equations derived for transverse isotropy media with vertical axis of symmetry (VTI) can be used to estimate anellipticity parameters of the orthotropic media with high accuracy. We make an accuracy comparison analysis between non-hyperbolic moveout equations that are widely used in modelling P -wave traveltimes in VTI media. Our comparison analysis is based on the accuracy of estimating anellipticity parameters defined for two vertical [ η (1) and η (2)] and one horizontal [ η (3)] symmetry planes of a wide range of orthotropic media. The anellipticity parameters are estimated from the inversion of numerical traveltimes that are generated for a single and horizontal orthotropic layer. We show that for the survey geometry used in this analysis, most of the moveout equations can produce reliable estimates of the anellipticity parameters η (1) and η (2). Despite the availability of traveltime data from wide-azimuths and long source–receiver offsets, none of the moveout equations can estimate η (3) as accurate as η (1) and η (2) estimations. We show that while all the equations produce mostly accurate estimates of η (1) and η (2), generalized moveout approximation (GMA) provides the greatest accuracy. We also assess the validity of using GMA to invert P -wave traveltimes for three anellipticity parameters in terms of anisotropy in the background medium, fracture weakness, source–receiver offset (group angle) and traveltime data coverage in both the vertical (dip) and horizontal (azimuthal) planes. We show that, for a reasonable offset range, GMA can be reliably used to estimate anellipticity parameters η (1) and η (2) for an orthotropic medium that is constructed by intensive fracturing of a host medium with strong VTI anisotropy. [ABSTRACT FROM AUTHOR]

The computation required to simulate surface nuclear magnetic resonance (SNMR) data increases proportionally with the number of sequences and number of pulses in each sequence. This poses a particular challenge to modelling steady-state SNMR, where suites of sequences are acquired, each of which require modelling 10–100 s of pulses. To model such data efficiently, we have developed a reformulation of surface NMR forward model, where the geometry of transmit and receive fields are encapsulated into a vector (or set of vectors), which we call B 1-volume-receive (BVR) curves. Projecting BVR curve(s) along complimentary magnetization solutions for a particular sequence amounts to computing the full SNMR forward model. The formulation has the additional advantage that computations for increased transmitter current amounts to a relative translation between the BVR and magnetization solutions. We generate 1-D kernels using BVR curves and standard integration techniques and find the difference is within 2 per cent. Using BVR curves, a typical suite of steady-state kernels can be computed two orders of magnitude faster than previous approaches. [ABSTRACT FROM AUTHOR]

A hybrid implicit-explicit (IMEX) discontinuous Galerkin spectral element time domain (DG-SETD) algorithm is proposed to simulate 3D elastic wave propagation in inhomogeneous media. In this method, the original problem can be divided into a number of well designed subdomains, and the mesh generation of different subdomains is completely independent, thus allowing flexible spatial discretization of complex geometry. The neighboring subdomains are connected by a Riemann transmission condition (RTC), and spectral elements with different interpolation orders can be used in different subdomains to maximize the computational efficiency of multiscale problems to facilitate parallel computing for different subdomains. In particular, the explicit or implicit time iteration scheme can be appropriately selected for a subdomain according to the size of its mesh elements to increase the global time step increment, thus giving higher computational efficiency: For subdomains with coarse meshes, the explicit time integration scheme is used and the time step increment is limited by the Courant−Friedrichs−Lewy (CFL) stability condition; for subdomain with fine structures and thus fine meshes, an implicit time integration scheme is used so that a large time step increment can be used without affecting the stability. In addition, a jump condition of displacement and velocity is introduced to accurately simulating fractures and faults, including lossless and viscous fractures with plane, curve or cross structures. This avoids the volume modeling of the extremely thin fracture structures, and effectively reduces the number of degrees of freedoms (DoFs) in the discretized system without the loss of accuracy. The accuracy, robustness and efficiency of the DG-SETD algorithm are demonstrated by quantitative comparisons of the waveforms with the commercial finite element software COMSOL. [ABSTRACT FROM AUTHOR]

SUMMARY: When the conventional finite-element method is used to simulate the 3-D direct-current (DC) resistivity response over a conductive earth with large complex structures, it requires finely discretized mesh to accurately represent the underground structures. Directly solving the current-conduction problem on the fine mesh will lead to a huge amount of calculation. In this paper, we develop a fast 3-D forward modelling method for DC resistivity method based on smooth multiscale finite-element algorithm. Instead of using the conventional polynomial basis functions, we construct the multiscale basis functions by solving the local boundary-value problems of partial differential equation in parallel on the multiscale meshes. The multiscale basis functions can capture the small-scale heterogeneous information in coarse cells and reflect it to the large scale by assembling macro matrix of coarse mesh. Thus, it enables us to quickly obtain the accurate solution by solving the original current-conduction problem with complex structures on coarse mesh. We further adopt the oversampling technology to improve the forward modelling accuracy. Besides, by combining with the gradient smoothing technology, we avoid establishing the continuous form of multiscale basis functions and their spatial derivative integral operation to rapidly assemble the macro matrix. Finally, the reliability of the proposed algorithm is verified by applying our code to the 3-D complex models and comparing it with the conventional finite-element method. [ABSTRACT FROM AUTHOR]

Earthquake location is a basic seismological problem and has a key role in many quantitative seismic analyses. The significant anomaly of traveltimes observed in mountainous areas with irregular surfaces has made earthquake location a challenge. To cope with this problem, we develop an effective eikonal equation-based earthquake location method based on unstructured mesh for 2-D/3-D isotropic and anisotropic media with irregular surfaces. First, the location misfit function is established by the reciprocity principle. Then, we use a global search algorithm to find the optimal origin times and hypocentres. Afterward, we apply eikonal equation-based master-event relocation method to relocate earthquakes on lateral boundaries of the region where the location results may be biased. To accurately compute the traveltime in 3-D vertically transversely isotropic and tilted transversely isotropic models with irregular surfaces or interfaces, we also propose a 3-D iterative fast sweeping method for eikonal equation on the unstructured tetrahedral mesh. Finally, we verify the proposed method by performing numerical experiments in 2-D/3-D irregular isotropic and anisotropic models. The numerical tests indicate that the proposed eikonal equation-based earthquake location method provides an effective way to find accurate hypocentre and origin time in 2-D/3-D isotropic and anisotropic media with irregular surfaces, even for inhomogeneous complex structures. In addition, the eikonal equation-based master-event location method has also yielded promising relocation results. [ABSTRACT FROM AUTHOR]

Recent advances in magnetic microscopy have enabled studies of geological samples whose weak and spatially nonuniform magnetizations were previously inaccessible to standard magnetometry techniques. A quantity of central importance is the net magnetic moment, which reflects the mean direction and the intensity of the magnetization states of numerous ferromagnetic crystals within a certain volume. The planar arrangement of typical magnetic microscopy measurements, which originates from measuring the field immediately above the polished surface of a sample to maximize sensitivity and spatial resolution, makes estimating net moments considerably more challenging than with spherically distributed data. In particular, spatially extended and nonuniform magnetization distributions often cannot be adequately approximated by a single magnetic dipole. To address this limitation, we developed a multipole fitting technique that can accurately estimate net moment using spherical harmonic multipole expansions computed from planar data. Given that the optimal location for the origin of such expansions is unknown beforehand and generally unconstrained, regularization of this inverse problem is critical for obtaining accurate moment estimates from noisy experimental magnetic data. We characterized the performance of the technique using synthetic sources under different conditions (noiseless data, data corrupted with simulated white noise, and data corrupted with measured instrument noise). We then validated and demonstrated the technique using superconducting quantum interference device microscopy measurements of impact melt spherules from Lonar crater, India and dusty olivine chondrules from the CO chondrite meteorite Dominion Range 08006. Plain Language Summary: Rocks may acquire records of the ambient magnetic field that can be preserved over geological timescales. This occurs when ferromagnetic minerals in those rocks cool, crystallize, or experience a pressure pulse. Magnetic microscopes are very sensitive instruments that can measure weak magnetic fields produced by rock samples at close proximity. They enable the analysis of rocks that were previously inaccessible to traditional paleomagnetic instrumentation. However, special data analysis techniques are required for obtaining the magnetic record preserved in a rock from the measurements collected by a magnetic microscope. We developed a novel analytical technique that allows us to use those instruments to successfully study small rock samples with nonuniform magnetization. We validate and demonstrate the technique using magnetic microscopy measurements of impact spherules from Lonar crater, India and mineral grains extracted from the meteorite Dominion Range 08006. Key Points: A spherical harmonic multipole expansion model can be used to estimate net magnetic moment from magnetic microscopy dataPlanar geometry and indeterminacy of the origin location for the multipole expansion require regularization of an ill‐posed inverse problemAccurate moment estimates can be obtained for measurements of compact geological samples even with low signal‐to‐noise ratios [ABSTRACT FROM AUTHOR]

The determinations of detailed stress states are of great importance for various environmental and engineering investigations, which makes numerical stress modelling a key issue in many fields. We developed a new stress modelling method governed by elastic wave equations using finite-difference scheme. By introducing an artificial damping factor to the particle velocity in wave modelling, the proposed method is able to solve both the dynamic stress evolution and the static stress state of equilibrium. We validate the proposed method both in body force and surface force benchmarks in different scales. With the proposed method, we are able to substantially improve the modelling accuracy of models in unbounded domain by using the perfectly matched layer as the artificial boundary conditions. A 3-D concrete-faced rockfill dam model is further presented as a numerical example of practical investigation. The consistent results with the finite-element method further illustrate the proposed method's applicability. As a minor modification to wave modelling scheme, the proposed stress modelling method is not only accurate for geological models through different scales, but also physically reasonable and easy to implement for geophysicists. [ABSTRACT FROM AUTHOR]

Because the magnetotelluric (MT) method uses natural sources, the electric and magnetic fields recorded in the field acquisition are not directly used but usually converted into other MT response functions for interpretation such as inversion. Considering that inversion results are dependent on types of input data, it can be helpful to analyse different characteristics of MT response functions for inversion. In this study, we examine sensitivity patterns of MT response functions used commonly in MT inversion, which are the impedance tensor, apparent resistivity, phase, tipper, effective impedance and phase tensor; and investigate how their sensitivity patterns affect inversion results. We first describe overall tendencies of 3-D sensitivity patterns of the MT response functions, and then classify the MT response functions into six groups based on 2-D sensitivity patterns computed at the surface, which are briefly called 'surface-sensitivity patterns' in this study. The 'diagonal components of the impedance' and 'off-diagonal components of the phase tensor', which have four petals-shaped surface-sensitivity patterns along the diagonal directions, belong to Group 1, and contribute to imaging 3-D subsurface structures from receivers installed evenly at the surface. Group 2 contains the ' xy -components of the impedance, apparent resistivity and phase' and ' yy -component of the phase tensor' whose surface-sensitivity patterns are linear in the y -axis. The ' yx -components of the impedance, apparent resistivity and phase' and ' xx -component of the phase tensor' that have strong linear surface-sensitivity patterns along the x -axis are classified into Group 3. The MT response functions of Groups 2 and 3 are useful for inversion of structures close to 2-D, whose strike extends along the y - and x -axes, respectively. Groups 4 and 5 include the ' x - and y -components of tipper' that possess linearly aligned two petals-shaped surface-sensitivity patterns in the x - and y -axes, respectively. The tipper can be helpful in imaging both 2-D and 3-D structures. The 'effective impedance' belongs to Group 6, whose surface-sensitivity patterns appear as a small circle. The surface-sensitivity patterns allow the effective impedance to have an advantage in interpretation of 1-D structures. By using several MT response functions for specific cases of 1-D, 2-D and 3-D interpretation of MT data, we investigate whether characteristics of the sensitivity patterns are reflected in modelling (simulating field data) and inversion results, and then suggest optimal MT response functions for those cases. In doing so, we show how to utilize the characteristics of the sensitivity patterns in inversion, and recommend the input MT response functions for inversion according to MT exploration situations. Our study provides basic information on similarities and differences of major MT response functions for inversion and insights on which MT response functions are suitable to increase the feasibility of MT inversion for different field situations based on the sensitivity patterns. [ABSTRACT FROM AUTHOR]

Abstract Recent advances in magnetic microscopy have enabled studies of geological samples whose weak and spatially nonuniform magnetizations were previously inaccessible to standard magnetometry techniques. A quantity of central importance is the net magnetic moment, which reflects the mean direction and the intensity of the magnetization states of numerous ferromagnetic crystals within a certain volume. The planar arrangement of typical magnetic microscopy measurements, which originates from measuring the field immediately above the polished surface of a sample to maximize sensitivity and spatial resolution, makes estimating net moments considerably more challenging than with spherically distributed data. In particular, spatially extended and nonuniform magnetization distributions often cannot be adequately approximated by a single magnetic dipole. To address this limitation, we developed a multipole fitting technique that can accurately estimate net moment using spherical harmonic multipole expansions computed from planar data. Given that the optimal location for the origin of such expansions is unknown beforehand and generally unconstrained, regularization of this inverse problem is critical for obtaining accurate moment estimates from noisy experimental magnetic data. We characterized the performance of the technique using synthetic sources under different conditions (noiseless data, data corrupted with simulated white noise, and data corrupted with measured instrument noise). We then validated and demonstrated the technique using superconducting quantum interference device microscopy measurements of impact melt spherules from Lonar crater, India and dusty olivine chondrules from the CO chondrite meteorite Dominion Range 08006.

Rupture dynamics along a heterogeneous fault is studied in the framework of the recently developed damage‐breakage rheological model. The model utilizes a fault structure accommodating shear and wear in a heterogeneous weak fault‐core a few centimeters thick, separating two blocks with damage level gradually decreasing toward the host rock. We first demonstrate the similarity of the static stress field around a fault, represented by a three‐body tribosystem, to the one developed around a rough frictional interface. We show that both models predict heterogeneous stress field pattern. We then apply simulations of rupture on heterogeneously damaged faults. We show that increasing initial heterogeneity amplitudes is associated with smaller events with lower slip rates. The simulations further allow us to quantify the amount of accumulated damage correlative with wearing. During the propagating rupture, the strength of the fault‐core evolves, leading to higher wear generation along relatively strong zones or barriers. The total wear production in a given event is strongly dominated by the initial damage heterogeneity. Processing of synthetic seismograms shows excess energy radiation of high‐frequency seismic waves comparing to the expected radiation from planar faults. This radiation is enhanced with the increase in the variability of fault strength heterogeneity. Further calculations of the scaled energy show good fit with previous seismological and laboratory observations and demonstrate the impact of fault heterogeneity on radiated energy during dynamic rupture. Therefore, initial fault heterogeneity, manifested by fault‐core strength or by geometrical irregularity controls many aspects of earthquake rupture, including slip displacement and velocity. Plain Language Summary: Earthquakes are triggered when the stress build‐up on rocks exceeds the resistance to fracture of fault material. When faulting occurs, energy is generated in the fault zone and radiated through the rock body to the surface. The ground motions and their durations strongly depend on the earthquake rupture characteristics such as the amount of displacement, velocity, and directivity. Many studies modeled natural faults as planar frictional interfaces that failed and were displaced during earthquakes. Yet, geological observations and seismological measurements demonstrated that this is probably a naive picture, and fault zones are very heterogeneous, including complicated geometry and variation of strength along and across them. Here we examine the effect of such initial heterogeneities on earthquake and show that the heterogeneity of the rock strength in the fault core itself strongly affects earthquake dynamics. In particular, larger magnitude earthquakes are predicted for the more homogenous fault. We also find that higher fault heterogeneity would result in excess high‐frequency radiation. Considering these effects on ground motion around the fault, we show that seismological records may capture these differences in fault characteristics during earthquakes. Key Points: Initial damage level and its heterogeneity along the fault core control dynamics and seismic moment of simulated rupture eventsWear production is affected by the initial damage distribution. During slip, higher wear accumulation is calculated for barriers along faultsHigh variability of fault damage heterogeneity significantly enhances high‐frequency radiation [ABSTRACT FROM AUTHOR]

Carbon capture and storage is an important technology for greenhouse gas mitigation. Monitoring of CO2 storage should, in addition to locating the plume, provide quantitative information on CO2 saturation. We propose a full waveform inversion (FWI) algorithm for the prediction of the spatial distribution of CO2 saturation from time-lapse seismic data. The methodology is based on the application of a rock-physics parametrized FWI scheme that allows for direct updating of reservoir properties. We derive porosity and lithology parameters from baseline data and use them as input to predict CO2 saturation from monitor data. The method is tested on synthetic time-lapse data generated for the Johansen formation model. Practical issues associated with field data applications, such as acquisition limitations, construction of the initial model, noise and uncertainty in the rock physics model, are taken into account in the simulation. The results demonstrate the robustness of our approach for reconstructing baseline and monitor models. We also illustrate the potential of the approach as compared to conventional two-step inversion algorithms, in which an elastic FWI prediction of velocities and density is followed by rock physics inversion. [ABSTRACT FROM AUTHOR]

The eikonal equation for transversely isotropic media is a complex quartic partial differential equation that causes instability in traveltime inversion. The phenomenon of shear wave splitting in a transversely isotropic model with a horizontal symmetry axis (HTI) provides information for inversion. However, the involvement of shear waves makes the computation of forward modelling and inversion more time-consuming than when only using P waves. Most studies of traveltime approximation based on the eikonal equation are based on an acoustic assumption, which does not work for shear waves. Therefore, based on perturbation theory, we derive approximate solutions of the eikonal equation in a heterogeneous elastic HTI model for P, SV and SH waves to improve the traveltime computation efficiency. The approximate solutions of linear, quadratic and Shanks forms are based on the linearization of the traveltime function in terms of fracture weaknesses in an isotropic background, which are suitable for joint P -, SV - and SH -wave traveltime inversion and tomography for fracture parameters. When applied for a homogenous model, the proposed formulation results in a simple and accurate traveltime approximation. The first- and second-order perturbation coefficients are used to analyse the sensitivity of the traveltime to fracture weakness. To improve the accuracy of our approximation, we also expand the approximation of the P -wave traveltime with respect to the anellipticity parameter η in an elliptic background, which is evaluated by performing an error test. Based on the accuracy and stability of the perturbation, we select a set of linear equations with respect to the fracture weaknesses and a set of more accurate nonlinear equations to serve different needs. [ABSTRACT FROM AUTHOR]

The in situ reservoir status monitoring plays a critical role in natural gas hydrate resource production. Considering the complexity of the field environment, a simulation framework for monitoring gas hydrates with cross-hole electrical resistivity tomography (CHERT) was developed to monitor the hydrate distribution during hydrate formation and dissociation. The simulation study comprised both numerical and physical experiments. The optimal CHERT array was designed through a numerical experiment. The effect of applying CHERT was verified through a physical experiment (a high-resistivity medium and hydrate formation experiment). The results show that improper electrode layouts will lead to varying degrees of low amplitude and blur boundary. An optimal CHERT array of a 100-mm electrode rod spacing, 8-mm electrode ring spacing and 48 electrode rings was obtained. The inversion results obtained using this CHERT array scheme can easily distinguish the distribution of high-resistivity targets and yield satisfactory results in hydrate formation experiments. These findings guarantee data processing and interpretation for applying CHERT in gas hydrate experiments and fields. [ABSTRACT FROM AUTHOR]

The separation of the effects of deep-seated sources of potential fields from those of shallower ones is a frequent requirement when interpreting magnetic or gravity fields. A common procedure is estimating the regional, long wavelength, component of the field by analysing the data over an area larger than that of the local feature of interest. The local components are found by subtracting the estimated regional from the observed data. These approaches may have difficulties in their application, as the dataset over large areas may not be available and other local anomalies, in the enlarged area, may prevent a reliable estimate of the regional field. We present an alternative and simple approach to the regional-residual separation problem not requiring the analysis over large areas and aiming at estimating the local, rather than the regional, component. Our method exploits the natural enhancement of short wavelengths obtainable by computing vertical derivatives of potential fields. An equivalent layer source is computed from the vertical derivative and is used to estimate the local field. The optimal differentiation order can be determined by inspecting the obtained results. This parameter may assume even fractional values, so that the method results a very versatile tool. The application to a complex synthetic case and two real data examples demonstrates the utility of this approach. In summary, our method has some peculiar characteristics making it an interesting alternative to currently used approaches to regional-residual separation: (i) it is a local method, so it can work well even when processing datasets relative to areas of limited extension; (ii) unlike most current methods, estimating a smooth regional component, our method directly produces an estimate of the local field and (iii) it is highly versatile, as the key parameter, that is the fractional differentiation order, can be finely adjusted up to obtain an optimal local field. [ABSTRACT FROM AUTHOR]

Fractures are ubiquitous in Earth's upper crust and strongly affect its physical properties. Their visibility for elastic waves hinges on their stiffness, a measure of the relation between changes in fracture aperture and applied stress, which itself significantly increases with increasing stress according to laboratory data. We constrained normal stiffness (and its stress dependence) of two fractures intersecting a borehole in the research mine 'Reiche Zeche' in Freiberg, Saxony, Germany, from pressure-recovery behaviour observed when the pressure in a double-packer interval enclosing them was briefly released during shut-in sequences. For the evaluation of the pressure recovery, we developed a generic model that accounts for the inherent convolution of mechanical closing of the fracture, fluid flow in the fracture, fluid exchange between fracture and the surrounding rock, here addressed as fracture leakage, and effects associated with finite borehole storage. We address the hydraulic system by considering three different 'storage containers', described by non-dimensionless storage parameters: the wellbore with a fixed storage capacity, the fracture with storage capacities related to fluid compressibility and fracture stiffness, and the surrounding rock constituting a container with unlimited storage capacity. To examine the model's characteristics, we investigated pressure transients predicted for specific cases, such as negligible fracture leakage, and subsequently numerically determined the sensitivity of the model predictions to the involved model parameters. The model's application to the field data yield estimates of properties of the fractures (length, aperture and stiffness) and of the surrounding rock mass (product of permeability and specific storage capacity). The parameters derived from the modelling compare well with independent constraints. An analysis of currently available fracture-stiffness data shows that the lack of knowledge of the stress dependence of stiffness severely hampers firm conclusions on scale dependence and the bias between methods. [ABSTRACT FROM AUTHOR]

SUMMARY: We present a computationally efficient method to approximately propagate uncertainty when linearly inverting seismic data for point source, time variable moment tensor components. The method is based on the assumption that the data residual, given by the difference between the observed seismic data and the data predicated by a linear inversion, contains the effects of both data and model uncertainty. Our method uses a distribution of data residuals, added directly to the data, in a pseudo-Monte Carlo scheme. Using the assumption that the data residual is a stochastic process, we use the well-known Karhunen–Loève (KL) theorem to construct a distribution of data residuals, where the required basis functions are constructed using Fourier series. The Fourier series are scaled by a product of a random variable and the real-valued spectral amplitudes of the original data residual's spectrum. Thus, the Fourier series and spectral amplitudes are eigenfunction-eigenvalue pairs used in the KL-based construction of data residual distribution. Using tests with synthetic data, we show that our method compares closely with a Finite Difference Monte Carlo (FDMC) method that we presented previously. More importantly, the method presented here is computationally several orders of magnitude faster than our previous FDMC method, and requires no a priori assumptions of model and/or data uncertainty. [ABSTRACT FROM AUTHOR]

The finite-difference time-domain scheme is augmented in order to treat the modelling of transient electromagnetic signals containing induced polarization effects from 3-D distributions of polarizable media. Compared to the non-dispersive problem, the discrete dispersive Maxwell system contains costly convolution operators. Key components to our solution for highly digitized model meshes are Debye decomposition and composite memory variables. We revert to the popular Cole-Cole model of dispersion to describe the frequency-dependent behaviour of electrical conductivity. Its inversely Laplace-transformed Debye decomposition results in a series of time convolutions between electric field and exponential decay functions, with the latter reflecting each Debye constituents' individual relaxation time. These function types in the discrete-time convolution allow for their substitution by memory variables, annihilating the otherwise prohibitive computing demands. Numerical examples demonstrate the efficiency and practicality of our algorithm.

Studying the seismic reflection characteristics of fractured strata can provide a theoretical basis for the exploration of unconventional oil and gas reservoirs. When several (more than 1) sets of oblique vertical fractures that are rotationally non-invariant are embedded in an isotropic or vertical transversely isotropic (VTI) background, the stratum exhibits equivalent monoclinic anisotropy. However, the exact algorithm of the reflection and transmission coefficients of monoclinic media is algebraically complicated and difficult to use in theoretical analysis and application. In this paper, under the assumption that fractures are filled with fluid, based on the linear-slip theory, we solve the phase velocities of monoclinic media from the Christoffel equations by using the first-order perturbation method and then derive the recursive equations of the PP -, PS 1-, and PS 2-wave reflection coefficients. The theoretical analysis results show that for an interface of monoclinic media, when setting a series of fracture densities for the lower medium, the relative error of the proposed approximate reflection coefficient of the PS 2 wave is less than 5.7 per cent, and those of the PP and PS 1 waves are generally less than 2.5 per cent for a wide range of anisotropy strength and impedance contrast; when setting a series of fracture densities for the upper medium, the relative error of the proposed approximate reflection coefficients of the PP, PS 1 and PS 2 waves are generally less than 6.4, 5.1 and 4.5 per cent, respectively. When degenerated to the horizontal transversely isotropic model, compared with Rüger's approximation, at incidence angles of 1–30°, the two results are almost coincident, but the proposed approximations are closer to the exact solutions. [ABSTRACT FROM AUTHOR]

We introduce a new finite-element (FE) based computational framework to solve forward and inverse elastic deformation problems for earthquake faulting via the adjoint method. Based on two advanced computational libraries, FEniCS and hIPPYlib for the forward and inverse problems, respectively, this framework is flexible, transparent and easily extensible. We represent a fault discontinuity through a mixed FE elasticity formulation, which approximates the stress with higher order accuracy and exposes the prescribed slip explicitly in the variational form without using conventional split node and decomposition discrete approaches. This also allows the first order optimality condition, that is the vanishing of the gradient, to be expressed in continuous form, which leads to consistent discretizations of all field variables, including the slip. We show comparisons with the standard, pure displacement formulation and a model containing an in-plane mode II crack, whose slip is prescribed via the split node technique. We demonstrate the potential of this new computational framework by performing a linear coseismic slip inversion through adjoint-based optimization methods, without requiring computation of elastic Green's functions. Specifically, we consider a penalized least squares formulation, which in a Bayesian setting—under the assumption of Gaussian noise and prior—reflects the negative log of the posterior distribution. The comparison of the inversion results with a standard, linear inverse theory approach based on Okada's solutions shows analogous results. Preliminary uncertainties are estimated via eigenvalue analysis of the Hessian of the penalized least squares objective function. Our implementation is fully open-source and Jupyter notebooks to reproduce our results are provided. The extension to a fully Bayesian framework for detailed uncertainty quantification and non-linear inversions, including for heterogeneous media earthquake problems, will be analysed in a forthcoming paper. [ABSTRACT FROM AUTHOR]

The b value is a critical parameter that has been used to describe spatial and temporal size distribution characteristics in data from both natural seismicity and laboratory rock acoustic emission (AE) tests. However, unlike the magnitude used for b value estimation in earthquakes which corrects for hypocentral distance, the apparent amplitude of the AE, is most commonly used for b value estimation in the rock AE test. This ignores the effect of attenuation, which may modify the frequency–amplitude distribution and affect the b value results. Here, we investigate the attenuation effect on the estimate of b value from synthetic data generated by random sampling of a doubly truncated frequency–amplitude distribution and a known underlying b value from a statistical perspective. We generated AE source amplitudes that obey an exponential distribution with different theoretical b values and applied an attenuation operator to obtain apparent frequency–amplitude data for comparison with the true underlying distribution. The results show that the b value is unchanged within a specific interval of amplitude distribution after accounting for attenuation under the conditions examined here, and that data truncation is the primary cause of modification of the frequency–amplitude distribution in the presence of attenuation, rather than attenuation alone. In addition, we employed the same synthetic data simulation scheme to estimate the minimum data volume for accurate estimation of the underlying b value from the apparent amplitude distribution. We found that the determination of the minimum data volume is influenced by the theoretical b value of the underlying distribution: the higher the theoretical b value, the larger the minimum data volume required, varying from 200 events from b  < 1.5 to 300 for b  > 1.5. [ABSTRACT FROM AUTHOR]

A detailed knowledge of wave reflectivity parametrization for non-welded interface plays an important role in fracture characterization and reservoir prediction. However, the effect of in situ stress on reflectivity formulae for non-welded contact interface remains unclear. Here based on the linear-slip theory and acoustoelasticity theory, we propose the accurate P -wave reflectivity and transmissivity formulae for planar P wave incident upon the non-welded interface between two elastic half-spaces under the in situ vertical stress. The discontinuous boundary conditions coupled with vertical stress is used to represent the non-welded contact interface for the elastic and viscous cases. Then accurate formulae are derived by solving the boundary conditions. Numerical results show that the proposed formulae are frequency-dependent and the non-welded interface performs as high- and low-pass filters for reflected and transmitted P wave at low stress, respectively. Moreover, initial vertical stress can change the maximum amplitudes of reflection and transmission coefficients and the corresponding critical angles. Particularly, vertical stress weakens the fracture-induced anomalous seismic reflections in zones lack of significant impedance contrasts at low incidence angles. The sensitivity analysis suggests a non-welded interface transforms into a welded interface at normal compliance |${\Delta _N} \le 5 \times {10^{ - 11}}$|  m Pa−1. The phases and amplitudes of the four coefficients reach the maximum values at the specific viscosity |$\eta \approx {10^6}$|  Pa·s m−1 for the case of P wave obliquely incident upon the viscous interface at incidence angle 30°. Finally, we describe a method for simulating seismic reflection at the interface of two elastic media with experimental measurement. [ABSTRACT FROM AUTHOR]

Scattered wave imaging provides a powerful tool for understanding Earth's structure. The development of finite frequency kernels for scattered waves has the potential for improving the resolution of both the structure and magnitude of discontinuities in S‐wave velocity. Here we present a 2‐D analytical expression for teleseismic P‐to‐S and S‐to‐P scattered wave finite‐frequency kernels for a homogeneous medium. We verify the accuracy of the kernels by comparing to a spectral element method kernel calculated using SPECFEM2D. Finally, we demonstrate the ability of the kernels to recover seismic velocity discontinuities with a variety of shapes including a flat discontinuity, a discontinuity with a sharp step, a discontinuity with a smooth bump, and an undulating discontinuity. We compare the recovery using the kernel approach to expected recovery assuming the classical common conversion point (CCP) stacking approach. We find that the P‐to‐S kernel increases recovery of all discontinuity structures in comparison to CCP stacking especially for the shallowest discontinuity in the model. The S‐to‐P kernel is less successful but can be useful for recovering the curvature of shallow discontinuity undulations. Finally, although we observe some variability in the amplitude of the kernels along the discontinuities, the kernels show some potential for recovering the magnitude of the velocity contrast across a discontinuity. Plain Language Summary: Seismic waves convert from compressional waves to shear waves and vice versa at sharp discontinuities in seismic velocity. These conversions are useful for tightly constraining the structure and properties of the seismic velocity discontinuities, and this can be used to better determine the physical and chemical properties of the Earth. Typical approaches that use the converted phases to constrain discontinuities assume that the discontinuities are relatively flat, and migrate the converted phase energy to depth, stacking the converted phase energy on a grid. However, discontinuities may have significant topography, for instance in subduction zones or beneath hotspots, rifts, or ridges. There are several different approaches to account for such topography with varying degrees of success and also computational cost. One approach is to use the theoretical sensitivity kernels to distribute the converted phase energy in space. Here we present analytical sensitivity kernels for converted phases and assess their ability to recover strong discontinuity topography. The approach is very fast and computationally efficient. The compressional to shear kernel successfully retrieves all shallow discontinuity structures, while the shear to compressional kernel may be useful for constraining discontinuity curvature in certain circumstances. Key Points: We present 2‐D analytical P‐to‐S and S‐to‐P finite frequency sensitivity kernels that are rapidly calculated with low computational costWe demonstrate the validity of the kernels in comparison to kernels calculated using SPECFEM2DAnalytical P‐to‐S kernels recover discontinuities with sharp or curved topography, and S‐to‐P kernels may be useful in certain situations [ABSTRACT FROM AUTHOR]

Abstract Scattered wave imaging provides a powerful tool for understanding Earth's structure. The development of finite frequency kernels for scattered waves has the potential for improving the resolution of both the structure and magnitude of discontinuities in S‐wave velocity. Here we present a 2‐D analytical expression for teleseismic P‐to‐S and S‐to‐P scattered wave finite‐frequency kernels for a homogeneous medium. We verify the accuracy of the kernels by comparing to a spectral element method kernel calculated using SPECFEM2D. Finally, we demonstrate the ability of the kernels to recover seismic velocity discontinuities with a variety of shapes including a flat discontinuity, a discontinuity with a sharp step, a discontinuity with a smooth bump, and an undulating discontinuity. We compare the recovery using the kernel approach to expected recovery assuming the classical common conversion point (CCP) stacking approach. We find that the P‐to‐S kernel increases recovery of all discontinuity structures in comparison to CCP stacking especially for the shallowest discontinuity in the model. The S‐to‐P kernel is less successful but can be useful for recovering the curvature of shallow discontinuity undulations. Finally, although we observe some variability in the amplitude of the kernels along the discontinuities, the kernels show some potential for recovering the magnitude of the velocity contrast across a discontinuity.

As early as the Palaeozoic era, the Tianshan Mountains constituted a vast tectonic convergent zone within the Eurasian Continent; by this time, they had already experienced multiple intense crustal movements. The Indian Plate is currently advancing northwards by more than 50 mm yr–1; it affects the Tianshan Mountains even though they are located more than 1000 km away. However, many of the geological issues of the Tianshan Mountains, especially in eastern Tianshan, which remains an active tectonic zone today, are not yet fully understood. There have been few studies into the seismogenic characteristics of this intense seismic zone. Here, the full gravity gradient tensor (GGT) within the eastern Tianshan Mountains was calculated using a stable spatial frequency domain algorithm, based on high-spatial resolution (7.2 arcsec or approximately 200 m) Global Gravity Model plus (GGMplus) surface gravity data. In addition, the features of the geological structures in the eastern Tianshan Mountains were interpreted using different combinations of GGT components. Moreover, the possibilities of using different combinations of GGT components to identify the distributions and strikes of faults was discussed regarding the potential earthquake rupture risks in this study area. The results show that (1) the distributions and strikes of the main fault zones in the study area are highly consistent with the linear geological features within the mountains, as extracted from a combination of different GGT components, including terrain effects. (2) Removing the terrain effects revealed that the different components of GGT (derived from the complete Bouguer gravity anomaly) vary gently, showing few linear geological features that could be extracted. This implies that satellite gravity data are not sufficient to analyse the characteristics of small underground geological structures (spatial scales less than ∼10 km). (3) Comparing the spatial distributions of fault traces depicted at different ages with the linear geological features of the study area revealed that discovering more unknown faults increased the coincidence ratio of linear geological features to fault traces in the area (from 27 to 40 per cent). In conclusion, the findings of this study represent a valuable reference for further revealing the seismogenic characteristics of this area, where there is a lack of detailed surface fault structure measurements due to its inaccessibility. The proposed method could also recognize faults in other similar areas with harsh conditions that are not feasible for ground-based surveys, such as forests and glaciers. [ABSTRACT FROM AUTHOR]

Gravity gradient tensors (GGTs) are used to investigate the density of subsurface structures in the Earth's crust and can reduce ambiguity during data interpretation. However, costs and research area restrictions often prevent their application during surveys, thereby limiting their utility. To address this limitation, a matrix equation based on the Taylor series expansion that uses the gravity vector and its neighbours was formulated to obtain the GGTs. Higher-order derivatives of the gravity vector were utilized to constrain the calculation, which improved the accuracy of the transformation. Synthetic data were used to demonstrate that the proposed approach improved accuracy when the radial component of the gravity vector was transformed into GGTs. This approach was also applied to gravity data from the Otway Basin in Australia. Compared with the measured GGT, the results obtained using the proposed approach had a relative error of 0.46. [ABSTRACT FROM AUTHOR]

Geomagnetic data assimilation merges past and present-day observations of the Earth's magnetic field with numerical geodynamo models and the results are used to initialize forecasts. We present a new 'proxy model' that can be used to test, or rapidly prototype, numerical techniques for geomagnetic data assimilation. The basic idea for constructing a proxy is to capture the conceptual difficulties one encounters when assimilating observations into high-resolution, 3-D geodynamo simulations, but at a much lower computational cost. The framework of using proxy models as 'gate-keepers' for numerical methods that could/should be considered for more extensive testing on operational models has proven useful in numerical weather prediction, where advances in data assimilation and, hence, improved forecast skill, are at least in part enabled by the common use of a wide range of proxy models. We also present a large set of systematic data assimilation experiments with the proxy to reveal the importance of localization and inflation in geomagnetic data assimilation. [ABSTRACT FROM AUTHOR]

Simplified solutions of the Bloch equation can lead to inaccurate estimates of hydrogeological parameters from surface nuclear magnetic resonance measurements. Even for single pulse measurements, using simplified forward models is common practice because of the computational intensity of obtaining the full-Bloch solution. These challenges are exacerbated for multipulse sequences. We show parallelizing the full-Bloch solver on a Graphics Processing Unit reduces the solve time by three orders of magnitude. Further optimizations by numerical, analytical and hybrid solutions yield an additional 3× speed up. We simulate the full-Bloch physics for free-induction decay, spin-echo and pseudo-saturation recovery excitation schemes for an unprecedented range of physical scenarios. We explore the time-dependence and relaxation time sensitivity in these solution spaces. Characterizing the solution spaces with polynomials of the relaxation times, the solutions can be rapidly reproduced; a technique known as fast-mapping. By fitting these higher dimensional solution ensembles with polynomials, the original fast-mapping technique is extended to include T 1 at arbitrary times. Accuracy of the 7th order polynomial is such that a minimum 96 per cent of the models are within a ±3 per cent relative error. This permits the rapid reproduction of full-Bloch solutions with a matrix multiplication and opens up surface NMR to time-series based inversion of single and multipulse data. [ABSTRACT FROM AUTHOR]

This 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 < 1.6 mGal over the study area. On the contrary, the RMSEs of global geopotential models against the checkpoints data are 7 mGal for the models up to the maximum degree of 2190 (a resolution of ∼9.1 km) and 5 mGal to 5399 (∼3.7 km). Our local models are significantly more accurate than the state-of-the-art global models and fully exploit the airborne vector data with the measurement error of ∼1 mGal. We also present the regional models constrained only by radial (vertical) or by lateral (horizontal) gravity observations. Those models are considerably less accurate than the one from 3-D gravity data inversion. Lastly, the regional models are validated against topography data. It is found that the gravity–topography correlation is 0.8–0.9 at 100 km, 0.5 at 20 km and higher than the correlations of the global models at all frequencies. The gravity–topography admittances estimated from our regional models indicate ∼130 mGal km−1 and imply the effective density of 2500 kg m−3 for topographical mass. [ABSTRACT FROM AUTHOR]

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. [ABSTRACT FROM AUTHOR]

In this study, we aim to solve the seismic inversion in the Bayesian framework by generating samples from the posterior distribution. This distribution incorporates the uncertainties in the seismic data, forward model, and prior information about the subsurface model parameters; thus, we obtain more information through sampling than through a point estimate (e.g. maximum a posteriori method). Based on the numerical cost of solving the forward problem and the dimensions of the subsurface model parameters and observed data, sampling with Markov chain Monte Carlo (MCMC) algorithms can be prohibitively expensive. Herein, we consider the promising Langevin dynamics MCMC algorithm. However, this algorithm has two central challenges: (1) the step size requires prior tuning to achieve optimal performance and (2) the Metropolis–Hastings acceptance step is computationally demanding. We approach these challenges by proposing an adaptive step-size rule and considering the suppression of the Metropolis–Hastings acceptance step. We highlight the proposed method's potential through several numerical examples and rigorously validate it via qualitative and quantitative evaluation of the sample quality based on the kernelized Stein discrepancy (KSD) and other MCMC diagnostics such as trace and autocorrelation function plots. We conclude that, by suppressing the Metropolis–Hastings step, the proposed method provides fast sampling at efficient computational costs for large-scale seismic Bayesian inference; however, this inflates the second statistical moment (variance) due to asymptotic bias. Nevertheless, the proposed method reliably recovers important aspects of the posterior, including means, variances, skewness and 1-D and 2-D marginals. With larger computational budget, exact MCMC methods (i.e. with a Metropolis–Hastings step) should be favoured. The results thus obtained can be considered a feasibility study for promoting the approximate Langevin dynamics MCMC method for Bayesian seismic inversion on limited computational resources. [ABSTRACT FROM AUTHOR]

Approaches to seismic modelling using integral methods, notably Born–WKBJ (Wentzel–Kramers–Brillouin–Jeffreys), have seen an increase in applications in geophysical prospecting and near-surface exploration. Moreover, due to its linearity characteristic, which can fast and efficiently simulate full-waveform information, the Born–WKBJ based method performs well in seismic inversion. However, the Born approximation is the linearized wavefield simulation by Taylor expansion of the wavefield and omitting the high-order term, which cannot simulate the wavefield accurately, including the amplitude, phase and waveform information. For the seismic inversion, especially the amplitude variation with angle/offset (AVA/AVO) inversion, the amplitude is the most important element for estimating the parameters. In this paper, a Rytov–WKBJ approximation-based method, which can simulate the seismic amplitude information more accurately, is introduced to pre-stack seismic inversion. Besides, in order to improve the resolution of the inversion results, the ℓ1 − 2-norm regularized basis pursuit inversion is introduced to the inversion algorithm. Then, we demonstrate the superiority of the proposed method with the zero-offset and angle-dependent seismic data simulation. The model tests show that the proposed method performs better than the conventional method significantly both on amplitude and phase of the seismic data. Finally, the inversion tests of synthetic and field seismic data indicate that the proposed method can obtain more accurate results with a higher resolution. [ABSTRACT FROM AUTHOR]

The seismograms for the multiscale crustal model due to dislocations are synthesized by a revised direct stiffness matrix method. By extracting the exponential growth terms related to wavenumber and layer thickness, the fast and accurate wavefield modelling can be achieved for the multiscale system with superficial fine layers (the layer thickness and velocity vary from metre level in the near-surface to kilometre level in deep crustal zones). This method allows relatively high-frequency cases of engineering interest (about 10 Hz) to be tackled without extra computations, linking the geophysics to the geotechnical earthquake engineering. The simulations considering superficial fine layers (5–50  m) show that the horizontal peak ground velocities can be amplified twice with superficial velocity decreasing from 0.4 to 0.15 km s–1. A case study using a realistic fine model in Tokyo metropolis elucidates that the displacements are localized within the epicentre distance about 5 km, predicting the displacement responses by factors up to 6.7, 1.1 and 6.7 for radial, tangential and vertical directions in comparison to the simplified model without superficial fine structures. [ABSTRACT FROM AUTHOR]

We present a new implementation of the shortest path method (SPM) that calculates accurate traveltimes in arbitrarily large model spaces without the requirement of large computational times and large amounts of memory, an inherent problem of the Dijkstra's-like algorithms. The multiple graph realizations method is based upon multiple sampling of the model space, using numerous random graphs. The performance of this new method is compared against the conventional way to improve the accuracy of SPM, which is to use denser grids and connectivity stencils of higher order. Our results suggest that although for relatively small models, single runs of the SPM are more suitable to achieve the desired accuracy, in large models, and after a certain level of desired accuracy, this approach becomes inefficient or even unfeasible, as the requirements in memory and computational time increases dramatically. On the contrary our method can achieve the desired accuracy with linear impact in computational time and negligible impact in required memory. [ABSTRACT FROM AUTHOR]

A popular idea is that accretion of sediment at a subduction zone commonly leads to the formation of a subduction channel , which is envisioned as a narrow zone located above a subducting plate and filled with vigorously circulating accreted sediment and exotic blocks. The circulation can be viewed as a forced convection, with downward flow in the lower part of the channel due to entrainment by the subducting plate, and a 'backflow' in the upper part of the channel. The backflow is often cited as an explanation for the exhumation of high-pressure/low-temperature metamorphic rocks from depths of 30 to 50 km. Previous analyses of this problem have mainly focused on the restricted case where the walls bounding the flow are artificially held fixed and rigid. A key question is if this configuration can be sustained on a geologically relevant timescale. We address this question using a coupled pair of corner flows. The pro-corner accounts for accretion and deformation directly above the subducting plate, and the retro-corner corresponds to a deformable region in the overlying plate. The two corners share a medial boundary , which is fully coupled but is otherwise free to rotate and deform. Our results indicate that the maintenance of a stable circulating flow in a narrow pro-corner (<15°) requires an unusually large viscosity ratio, μretro/μpro > 103. For lower viscosity ratios, the medial boundary would rotate rearwards, converting the initially narrow pro-corner into an obtuse geometry. For a stable narrow corner, we show that the backflow within the corner is caused by downward convergence of the incoming flow and an associated downward increase in dynamic pressure, which reaches a maximum at the corner point. The total pressure is thus expected to be much greater than predicted using a lithostatic gradient, which means that estimates of depth from metamorphic pressure would have to be adjusted accordingly. In addition, we show that the velocity fields associated with a forced corner flow and a buoyancy-assisted channel flow are nearly identical. As such, structural geology studies are not sufficient to distinguish between these two processes. [ABSTRACT FROM AUTHOR]

Evidence from seismology, geology and geodynamic studies suggests that regional-scale lower crustal flow occurs in many tectonic settings. Pressure gradients caused by mantle processes and crustal density heterogeneity can provide driving force for lower crustal flow. Numerically modelling such flow can be computationally expensive. However, by exploiting symmetry in the physical system, it is possible to represent the vertical component of flow in terms of its lateral components, thereby reducing the problem's spatial dimension by one. Here, we present a mathematical formulation for flow in a viscous channel below an elastic upper plate, which is optimized for solution by common numerical methods. Our formulation drastically reduces the computational load required to simulate lower crustal flow over large areas and long timescales. We apply this model to two example problems, with and without an elastic upper plate, identifying combinations of parameters that are capable of generating measurable geologic uplift. [ABSTRACT FROM AUTHOR]

In this article, we validate the neuro‐fuzzy kinematic finite‐fault inversion method by studying the rupture process of the Mw6.2, Aug/24/2016, Amatrice, central Italy, earthquake. We jointly invert three different datasets to infer the spatio‐temporal slip distribution, namely static and high‐rate GNSS data (<=0.06 Hz) and strong‐motion data (0.06−0.5 Hz). Each data set is used to constrain a different frequency range of the source model, depending on the sensitivity of the data set. The inferred slip shows a slow nucleation phase at shallow depths of 3–4 km, followed by a bilateral rupture that forms two asperities, one to the NW (Norcia) and another to the SE (Amatrice) of the hypocenter. Our inferred slip is compared with those previously obtained using well‐established methods. In order to select an adequate number of fuzzy basis functions, we propose two alternative procedures, which yield the same general slip features. The first approach consists of ensuring that the inverse problem is formally over‐determined and uses the same number of basis functions at all frequencies. The second approach is based on a maximum‐likelihood analysis of the model misfit and selects a different number of basis functions for each frequency. The maximum‐likelihood approach allows for more basis functions at high frequencies, where more detail in the spatial slip distribution is needed. The solution obtained with the maximum‐likelihood approach provides a more physically plausible source time function, which shows less back slip artifacts. The accurate prediction of high‐rate GNSS traces not used in the inversion attests the robustness of the inferred slip model. Plain Language Summary: The accuracy of the earthquake's slip imaging method proposed in the companion article by Kheirdast et al. (2021), is tested using data recorded during the 2016 Mw6.2, Amatrice, earthquake. Many previously well‐understood features of the rupture process of this earthquake are well retrieved by the proposed method, including a slow nucleation phase, a bilateral rupture and a high‐amplitude slip patch to the NW of the hypocenter. We present a new strategy to spatially discretize the fault plane, which maximizes the inversion power by taking into account the uncertainty at different frequencies. Key Points: The neuro‐fuzzy kinematic finite‐fault inversion is used to study the Amatrice earthquake and the inferred slip model is compared with previous studiesWe present a new maximum‐likelihood strategy to select an adequate, problem‐specific number of fuzzy basis functionsThe inferred source model is able to correctly predict high‐rate GNSS waveforms not used in the inversion [ABSTRACT FROM AUTHOR]

The stability of partial differential equations determines the properties of their solutions. This study focuses on the stability analysis of the equations describing wave propagation in fluids-saturated porous media. We briefly introduce the stability analysis method for the wave propagation equations and discuss the adverse effects on the solutions. In this way, the first part of this paper is mainly devoted to the analysis of the Tuncay and Corapcioglu's (TC) model, which describes the dynamic behaviour of porous media saturated with two immiscible fluids. It is pointed out that the TC model allows spatially bounded but time-exponentially exploding solutions and may yield unstable numerical results. Based on the deduced unstable factors, we construct a stable equivalent fluid (SEF) model. We rigorously analyse the stability of the SEF model using the energy method. For predicting the influence of saturation on wave velocity, the robustness of this model is preserved due to its consistency with the original TC model. Furthermore, the numerical simulations of the wavefields show that the results of the TC model exponentially increase with time after the initial effective wave signal, which does not occur in the SEF model curves. This indicates the necessity of considering the stability from a mathematical point of view during the construction of physical model. It could be useful to merge the mathematical stability theory with the geophysical wave propagation modelling theory. [ABSTRACT FROM AUTHOR]

The mathematical formulation of the Muskingum method, like that of many numerical schemes used for river routing, requires that all upstream river reaches be updated prior to updating the flow rate of any given reach. Due to this topological constraint, such numerical schemes have traditionally been solved in an upstream-to-downstream manner which imposes inherent limitations on the speedup that can be achieved in a parallel computing environment because each computing core has to wait for completion of all cores addressing upstream subbasins prior to starting its own subbasin. The research presented in this paper quantifies the exact influence among river reaches during the update step of the Muskingum method and shows that the influence decreases with increasing distance between two reaches until it becomes too small to be accounted for by floating-point arithmetic. A formal definition of the minimal distance from which the relative influence becomes numerically inexistent—the radius of influence—is presented. Based on this distance, expressed as a number of river reaches, a new estimate of the maximum theoretical speedup that can be achieved by the Muskingum method or by similar numerical schemes is presented and implies large potential gains in computing time when domains are much larger than the radius of influence. An application to the approximately 180,000 river reaches of the Upper Mississippi River Basin at a 15 min time step over 2004 shows a radius of influence on the order of 150 river reaches. The speedup obtained for this application is much higher than previously thought possible but also much lower than could be attained, suggesting that further investigations are necessary.

Geophysical imaging using the inversion procedure is a powerful tool for the exploration of the Earth's subsurface. However, the interpretation of inverted images can sometimes be difficult, due to the inherent limitations of existing inversion algorithms, which produce smoothed sections. In order to improve and automate the processing and interpretation of inverted geophysical models, we propose an approach inspired from data mining. We selected an algorithm known as DBSCAN (Density-Based Spatial Clustering of Applications with Noise) to perform clustering of inverted geophysical sections. The methodology relies on the automatic sorting and clustering of data. DBSCAN detects clusters in the inverted electrical resistivity values, with no prior knowledge of the number of clusters. This algorithm has the advantage of being defined by only two parameters: the neighbourhood of a point in the data space, and the minimum number of data points in this neighbourhood. We propose an objective procedure for the determination of these two parameters. The proof of concept described here is applied to simulated ERT (electrical resistivity tomography) sections, for the following three cases: two layers with a step, two layers with a rebound, and two layers with an anomaly embedded in the upper layer. To validate this approach, sensitivity studies were carried out on both of the above parameters, as well as to assess the influence of noise on the algorithm's performance. Finally, this methodology was tested on real field data. DBSCAN detects clusters in the inverted electrical resistivity models, and the former are then associated with various types of earth materials, thus allowing the structure of the prospected area to be determined. The proposed data-mining algorithm is shown to be effective, and to improve the interpretation of the inverted ERT sections. This new approach has considerable potential, as it can be applied to any geophysical data represented in the form of sections or maps. [ABSTRACT FROM AUTHOR]

A non-iterative first-arrival traveltime inversion method (NFTI) is proposed for building smooth velocity models using seismic diving waves observed on irregular surface. The new ray and traveltime equations of diving waves propagating in smooth media with undulant observation surface are deduced. According to the proposed ray and traveltime equations, an analytical formula for determining the location of the diving-wave turning points is then derived. Taking the influence of rough topography on first-arrival traveltimes into account, the new equations for calculating the velocities at turning points are established. Based on these equations, a method is proposed to construct subsurface velocity models from the observation surface downward to the bottom using the first-arrival traveltimes in common offset gathers. Tests on smooth velocity models with rugged topography verify the validity of the established equations, and the superiority of the proposed NFTI. The limitation of the proposed method is shown by an abruptly-varying velocity model example. Finally, the NFTI is applied to solve the static correction problem of the field seismic data acquired in a mountain area in the western China. The results confirm the effectivity of the proposed NFTI. [ABSTRACT FROM AUTHOR]