Three-dimensional gravity inversion for Moho depth at rifted continental margins incorporating a lithosphere thermal gravity anomaly correction

This paper describes a method for determining Moho depth, lithosphere thinning factor (γ= 1 − 1/β) and the location of the ocean–continent transition at rifted continental margins using 3-D gravity inversion which includes a correction for the large negative lithosphere thermal gravity anomaly within continental margin lithosphere. The lateral density changes caused by the elevated geotherm in thinned continental margin and adjacent ocean basin lithosphere produce a significant lithosphere thermal gravity anomaly which may be in excess of −100 mGal, and for which a correction must be made in order to determine Moho depth accurately from gravity inversion. We describe a method of iteratively calculating the lithosphere thermal gravity anomaly using a lithosphere thermal model to give the present-day temperature field from which we calculate the lithosphere thermal density and gravity anomalies. For continental margin lithosphere, the lithosphere thermal perturbation is calculated from the lithosphere thinning factor (γ= 1 − 1/β) obtained from crustal thinning determined by gravity inversion and breakup age for thermal re-equilibration time. For oceanic lithosphere, the lithosphere thermal model used to predict the lithosphere thermal gravity anomaly may be conditioned using ocean isochrons from plate reconstruction models to provide the age and location of oceanic lithosphere. A correction is made for crustal melt addition due to decompression melting during continental breakup and seafloor spreading. We investigate the sensitivity of the lithosphere thermal gravity anomaly and the predicted Moho depth from gravity inversion at continental rifted margins to the methods used to calculate and condition the lithosphere thermal model using both synthetic models and examples from the North Atlantic.     

Two-dimensional non-linear inversion of VLF-R data using simulated annealing

The VLF-R (very low frequency-resistivity) data, i.e. the apparent resistivity ( ρ _a ) and phase ( φ ) data, were inverted individually and jointly using the VFSA (very fast simulated annealing) global inversion approach. Global inversion results for synthetic data without and with various amounts of random and normally distributed Gaussian noise reveal that the inversion of neither the ρ _a nor φ data alone yields the true parameters of the structures. However, the joint inversion of the ρ _a and φ data yields very good estimates of the model parameters. Five models, representing typical subsurface structures in the shield areas, are studied here. Various models achieved after 10 VFSA runs were used to compute the mean model and the corresponding covariance and correlation matrices, which were used to estimate the uncertainties in the mean model parameters and correlations between the model parameters. We observe that these correlations follow the physics associated with the problem. VLF-R field data due to a nearly vertical contact structure and a very thick dyke-like structure were also inverted to demonstrate the efficacy of the approach in the delineation of the parameters of 2-D structures.     

Gravimetric determination of an intrusive complex under the Island of Faial (Azores): some methodological improvements

We present some improvements of a gravity inversion method to determine the geometry of the anomalous bodies for priori density contrasts. The 3-D method is based on an exploratory process applied, not for the global model, but for the steps of a growth approach. The (positive and/or negative) anomalous structure is described by successive aggregation of cells, while its corresponding gravity field remains nearly proportional to the observed one. Moreover, a simple (e.g. linear) regional trend can be simultaneously adjusted. The corresponding program is applied to new gravity data on the volcanic island of Faial (Azores archipelago). The inversion approach shows a subsurface anomalous structure for the island, the main feature being an elongated high-density body. The body is interpreted as a compact sheeted dyke swarm, emplaced along Faial-Pico Fracture Zone, a leaky transform structure that forms the current boundary between Eurasian and African plates in the Azores area. The new results in this paper are (1) a Bouguer gravity anomaly map, (2) several improvements in the inversion process (robust process, optimal balance fitness/model magnitude), (3) a new gravimetric method for estimating the mean terrain density, (4) a 3-D model for subsurface mass anomalies in Faial and (5) some interpretative conclusions about a main intrusive complex detected under the island as a wall-like structure extending from a depth of 0.5 to 6 km b.s.l., with a N100°E trend and corresponding to an early fissural volcanic episode controlled by the regional tectonics.     

Three-dimensional massively parallel electromagnetic inversion—I. Theory

An iterative solution to the non-linear 3-D electromagnetic inverse problem is obtained by successive linearized model updates using the method of conjugate gradients. Full wave equation modelling for controlled sources is employed to compute model sensitivities and predicted data in the frequency domain with an efficient 3-D finite-difference algorithm. Necessity dictates that the inverse be underdetermined, since realistic reconstructions require the solution for tens of thousands of parameters. In addition, large-scale 3-D forward modelling is required and this can easily involve the solution of over several million electric field unknowns per solve. A massively parallel computing platform has therefore been utilized to obtain reasonable execution times, and results are given for the 1840-node Intel Paragon. The solution is demonstrated with a synthetic example with added Gaussian noise, where the data were produced from an integral equation forward-modelling code, and is different from the finite difference code embedded in the inversion algorithm     

Seabed shear modulus profile inversion using surface gravity (water) wave-induced bottom motion

Summary. The motion of a layered seabed induced by propagating gravity water waves is modelled by the coupling matrix for a massless incompressible elastic bed according to Yamamoto. An amplitude inversion scheme is developed to extract the bottom shear modulus profile from the motion of the seabed at a point on the bed surface using the linear inverse theory combined with the Yamamoto theory of wave-seabed interaction. Numerical tests using synthetic data without noise confirmed that the inversion is unique and consistent. The inversion technique is applied to a set of field data from the Mississippi River Delta. A good agreement is obtained between the inverted profile and the direct measurements.     

Non-linear altimetric geoid inversion for lithospheric elastic thickness and crustal density

The inversion of high-resolution geoid anomaly maps derived from satellite altimetry should allow one to retrieve the lithospheric elastic thickness, T _e , and crustal density, _c . Indeed, the bending of a lithospheric plate under the load of a seamount depends on both parameters, and the associated geoid anomaly is correspondingly dependent on the two parameters. The difference between the observed and modelled geoid signatures is estimated by a cost function, J , of the two variables, T _e and _c . We show that this cost function forms a valley structure along which many local minima appear, the global minimum of J corresponding to the true values of the lithospheric parameters. Classical gradient methods fail to find this global minimum because they converge to the first local minimum of J encountered, so that the final parameter estimate strongly depends on the starting pair of values ( T _e ,   _c ). We here implement a non-linear optimization algorithm to recover these two parameters from altimetry data. We demonstrate from the inversion of synthetic data that this approach ensures robust estimates of T _e and _c by activating two search phases alternately: a gradient phase to find a local minimum of J , and a tunnelling phase through high values of the cost function. The accuracy of the solution can be improved by a search in an iteratively restricted parameter subspace. Applying our non-linear inversion to the Great Meteor Seamount geoid data, we further show that the inverse problem is intrinsically ill-posed. As a consequence, minute geoid (or gravity) data errors can induce large changes in any recovery of lithospheric elastic thickness and crustal density.     

Moment tensor inversion of complex earthquakes

Summary. Moment tensor inversion methods can be applied with success in the determination of source properties of simple earthquakes. However, these methods utilize the assumption of a point source, which is inadequate for modelling many complicated, shallow earthquakes. For complex earthquakes, an inversion using finite faulting models is desirable but the number of parameters involved requires that a good starting model be found or that independent constraints be placed on some of the parameters. A method is presented for low-pass filtering both the data and Green's functions, passing only signals with wavelengths greater than the dimension of the entire fault. The filter tends to smooth complications in the waveforms and allows application of the point source moment tensor inversion. This method is applied to body waves from the 1978 Thessaloniki, Greece, earthquake, the 1971 San Fernando earthquake and to a multiple-point source synthetic model of the San Fernando event. For the Thessaloniki event, although a multiple-source mechanism has been suggested, inversion results before and after filtering were essentially identical, indicating that a point source mechanism is sufficient in modelling the long-period, teleseismic body waves. In the case of the San Fernando earthquake, the point source Green's functions were incapable of simultaneously modelling the P - and SH -waves. Inversion of P -waves alone resulted in extreme parameter resolution problems, but allowed constraint in one axis of the moment tensor and suggested an overall source time function. Inversion of a synthetic San Fernando data set yielded similar results, but allowed an investigation of the shortcomings of the method under controlled circumstances. Although the results may require substantial interpretation, the method presented represents a simple first step in the analysis of complex earthquakes.     

Inversion of controlled-source seismic tomography and gravity data with the self-adaptive wavelet parametrization of velocities and interfaces

A self-adaptive automated parametrization approach is suggested for the sequential inversion of controlled-source seismic tomography and gravity data. The velocities and interfaces are parametrized by their Haar wavelet expansion coefficients. Only those coefficients that are well constrained by the data, as measured by the number of rays that cross the corresponding wavelet function support area and their angular coverage, are inverted for, others are set to zero. This approach results in a reasonable distribution of resolution throughout the model even in cases of irregular ray coverage and does overcome the trade-off between different types of model parameters. A modified sequential inversion approach is suggested to join the traveltimes and gravity anomalies inversion. An algorithm is developed that inverts for smooth velocity and density variations inside the seismic layer, the position of its bottom interface as well as for optimal values of the velocity-to-density regression coefficients. The algorithm makes use of direct (diving), reflected and head (critically refracted) wave traveltimes. The algorithm workflow is demonstrated on a synthetic data example.     

Estimating Model Parameters from Self-Potential Anomaly of 2D Inclined Sheet Using Whale Optimization Algorithm: Applications to Mineral Exploration and Tracing Shear Zones

The use of spontaneous potential (SP) anomalies is well known in the geophysical literatures because of its effectiveness and significance in solving many complex problems in mineral exploration. The inverse problem of self-potential data interpretation is generally ill-posed and nonlinear. Methods based on derivative analysis usually fail to reach the optimal solution (global minimum) and trapped in a local minimum. A new simple heuristic solution to SP anomalies due to 2D inclined sheet of infinite horizontal length is investigated in this study to solve these problems. This method is based on utilizing whale optimization algorithm (WOA) as an effective heuristic solution to the inverse problem of self-potential field due to a 2D inclined sheet. In this context, the WOA was applied first to synthetic example, where the effect of the random noise was examined and the method revealed good results using proper MATLAB code. The technique was then applied on several real field profiles from different localities aiming to determine the parameters of mineralized zones or the associated shear zones. The inversion parameters revealed that WOA detected accurately the unknown parameters and showed a good validation when compared with the published inversion methods.

     

P-wave traveltime and polarization tomography of VSP data

The presence of anisotropy requires that tomographic methods be generalized to account for anisotropy. This generalization allows geological structure to be correctly imaged and allows the anisotropic parameters to be estimated. Use of isotropic inversion for imaging anisotropic structures gives systematic trends in the traveltime and polarization residuals. However, due to the limited directional coverage, the traveltimes along may not be sufficient to study the anisotropic properties of the structure. Polarizations can provide independent information on the structure. Traveltime and polarization inversion are applied to synthetic examples simulating VSP experiments. Transverse isotropy and 1-D structure are assumed. Plots of traveltime and polarization residuals are an important tool to detect the anomalies due to the presence of anisotropy. For receivers located in anisotropic layers, polarization residuals display consistent anomalies of several degrees. The synthetic examples show that even the simple 1-D problem is difficult, when using direct arrivals only. Large a posteriori errors in anisotropic parameters are obtained by traveltime inversion in layers where available incidence angles are less than 45°. Resolution of the tomographic image of VSP data is greatly improved by a combination of traveltime and polarization information. In order to obtain accurate inversion results, the measurement error of polarization data should be kept to within a few degrees.     

Sequential integrated inversion of refraction and wide-angle reflection traveltimes and gravity data for two-dimensional velocity structures

A new algorithm is presented for the integrated 2-D inversion of seismic traveltime and gravity data. The algorithm adopts the 'maximum likelihood' regularization scheme. We construct a 'probability density function' which includes three kinds of information: information derived from gravity measurements; information derived from the seismic traveltime inversion procedure applied to the model; and information on the physical correlation among the density and the velocity parameters. We assume a linear relation between density and velocity, which can be node-dependent; that is, we can choose different relationships for different parts of the velocity–density grid. In addition, our procedure allows us to consider a covariance matrix related to the error propagation in linking density to velocity. We use seismic data to estimate starting velocity values and the position of boundary nodes. Subsequently, the sequential integrated inversion (SII) optimizes the layer velocities and densities for our models. The procedure is applicable, as an additional step, to any type of seismic tomographic inversion.
We illustrate the method by comparing the velocity models recovered from a standard seismic traveltime inversion with those retrieved using our algorithm. The inversion of synthetic data calculated for a 2-D isotropic, laterally inhomogeneous model shows the stability and accuracy of this procedure, demonstrates the improvements to the recovery of true velocity anomalies, and proves that this technique can efficiently overcome some of the limitations of both gravity and seismic traveltime inversions, when they are used independently.
An interpretation of field data from the 1994 Vesuvius test experiment is also presented. At depths down to 4.5 km, the model retrieved after a SII shows a more detailed structure than the model obtained from an interpretation of seismic traveltime only, and yields additional information for a further study of the area.     

Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data

A tomographic inversion technique that inverts traveltimes to obtain a model of the subsurface in terms of velocities and interfaces is presented. It uses a combination of refraction, wide-angle reflection and normal-incidence data, it simultaneously inverts for velocities and interface depths, and it is able to quantify the errors and trade-offs in the final model. The technique uses an iterative linearized approach to the non-linear traveltime inversion problem. The subsurface is represented as a set of layers separated by interfaces, across which the velocity may be discontinuous. Within each layer the velocity varies in two dimensions and has a continuous first derivative. Rays are traced in this medium using a technique based on ray perturbation theory, and two-point ray tracing is avoided by interpolating the traveltimes to the receivers from a roughly equidistant fan of rays. The calculated traveltimes are inverted by simultaneously minimizing the misfit between the data and calculated traveltimes, and the roughness of the model. This 'smoothing regularization' stabilizes the solution of the inverse problem. In practice, the first iterations are performed with a high level of smoothing. As the inversion proceeds, the level of smoothing is gradually reduced until the traveltime residual is at the estimated level of noise in the data. At this point, a minimum-feature solution is obtained, which should contain only those features discernible over the noise.
The technique is tested on a synthetic data set, demonstrating its accuracy and stability and also illustrating the desirability of including a large number of different ray types in an inversion.     

Bayesian inversion with Markov chains—I. The magnetotelluricone-dimensional case

We use Monte Carlo Markov chains to solve the Bayesian MT inverse problem in layered situations. The domain under study is divided into homogeneous layers, and the model parameters are the conductivity of each layer. We use an a priori distribution of the parameters which favours smooth models. For each layer, the a priori and a posteriori distributions are digitized over a limited set of conductivity values.
  The Markov chain relies on updating the model parameters during successive scanning of the domain under study. For each step of the scanning, the conductivity is updated in one layer given the actual value of the conductivity in the other layers. Thus we designed an ergodic Markov chain, the invariant distribution of which is the a posteriori distribution of the parameters, provided the forward problem is completely solved at each step.
  We have estimated the a posteriori marginal probability distributions from the simulated successive values of the Markov chain. In addition, we give examples of complex magnetotelluric impedance inversion in tabular situations, for both synthetic models and field situations, and discuss the influence of the smoothing parameter.     

3D Subsurface Characterization of Banded Iron Formation Mineralization using Large-Scale Gravity Data: A Case Study in Parts of Bharatpur,Dausa and Karauli Districts of Rajasthan,India

Our study interprets large-scale gravity data to delineate concealed banded iron formation (BIF) iron mineralization in India's Rajasthan province. The study area belongs to the Bharatpur, Dausa, and Karauli districts of Rajasthan. We measured 1462 gravity readings to understand the rock types, depth and geometry of the different rock formations in the proposed study area. We also collected representative lithologies from more than 100 locations in the study area and calculated their density values. The measured gravity datasets are investigated via qualitative (e.g., Bouguer anomaly, first derivative and second derivative) and quantitative (radially averaged power spectrum, 3D Euler deconvolution, and 3D inversion) approach. The qualitative methods suggest a general NE–SW orientation of the BIFs, controlled by the general trend of the study area's structural setting. The lithological contact between the Bhilwara and Vindhyan Supergroups is demarcated by a NE–SW trending steep gravity gradient zone. In this area, representative lithologies yield high densities (about 3.746 gm/cc), and the samples identified as BIF represent exploration targets for iron ore. We have also developed our own in-house 3D gravity inversion code in this study. A model space inversion algorithm is converted into a data space using the identity relationship. It makes inversion algorithm very user-friendly on conventional desktop computers. The outcomes from the 3D inversion suggest that the concealed iron ore thickens to the west. This interpretation is also in good correlation with Euler 3D deconvolution of the gravity data.

     

3D Geological modelling and gravity inversion of a structurally complex carbonate area: application for karstified massif localization

This paper proposes a new methodology to improve the location of potential karstified areas by gravity inversion of a 3D geological model. A geological 3D model is built from surface observations, 2D seismic reflection profiles and well data. The reliability of this geological 3D model obtained from integration, interpretation and interpolation of such data is first tested against the structural consistency of the model. Its theoretical gravimetric response is compared to gravity field during the forward problem in order to evaluate the validity/robustness of the geological model. The coherency between the gravity field and the gravimetric response is tested. The litho‐inversion modelling quantifies the distribution of rock density in a probabilistic way, taking into account the geology and physical properties of rocks, while respecting the geological structures represented in the 3D model. The result of the inversion process provides a density distribution within carbonate formations that can be discussed in term of karstification distribution. Thus, lower densities correlate with areas that are strongly karstified. Conversely, higher than mean densities are found in carbonate formations mostly located under marly and impervious formations, preserving carbonate from karstification and paleokarstification.     

Finite cube elements method for calculating gravity anomaly and structural index of solid and fractal bodies with defined boundaries

In this paper, we present the finite cube elements method (FCEM); a novel numerical tool for calculating the gravity anomaly g and structural index SI of solid models with defined boundaries and variable density distributions, tilted or in normal position (e.g. blocks, faulted blocks, cylinders, spheres, hemispheres, triaxial ellipsoids). Extending the calculation to fractal objects, such as Menger sponges of different orders and bodies defined by polyhedrons, demonstrates the robustness of FCEM. In addition, approximating the cube element by a sphere of equal volume makes the calculation of gravitation and related derivatives much simpler. In gravity modelling of a sphere, cubes with edges of 100 m and 200 m achieve a good compromise between running time and overall error.
Displaying the distribution of SI of the studied models on contour maps and profiles will have a strong impact on the forward and inverse modelling of potential field data, especially for Euler deconvolution.
For Menger sponges, plots of gravity elements g and its derivatives show similar patterns independent of fractal order. Moreover, both the pattern and magnitude of SI are independent of fractal order, allowing the use of SI as a new invariant measure for fractal objects. However, SI pattern and magnitude strongly depend on the depth to the buried bodies as do other elements
In this study, we also present a new type of plot; the structural index against distance variation diagrams from which we extract the three critical SI ( CSI ) values, one per axis. The inversion of gravity anomaly data at CSI values gives the optimal mean location of the buried body.     

The relationship between gravity and bathymetry in the Pacific Ocean

Summary. Surface-ship and satellite derived data have been compiled in new free-air gravity anomaly, bathymetry and geoid anomaly maps of the Pacific Ocean basin and its margin. The maps are based on smoothed values of the gravity anomaly, bathymetry and geoid interpolated on to a 90 × 90 km grid. Each smoothed value was obtained by Gaussian filtering measurements along individual ship and subsatellite tracks. The resulting maps resolve features in the gravity, bathymetry and geoid with wavelengths that range from a few hundred to a few thousand kilometres. The smoothed values of bathymetry and geoid anomaly have been corrected for age. The resulting maps show the Pacific ocean basin is associated with a number of ENE–WSW-trending geoid anomaly highs with amplitudes of about ± 5 m and wavelengths of about 3000 km. The most prominent of these highs correlate with the Magellan seamounts–Marshall Gilbert Islands–Magellan rise and the Hess rise–Hawaiian ridge regions. The correlation between geoid anomaly and bathymetry cannot be explained by models of static compensation, but is consistent with a model in which the geoid anomaly and bathymetry are supported by some form of dynamic compensation. We suggest that the dynamic compensation, which characterizes oceanic lithosphere older than 80 Myr, is the result of mantle convection on scales that are smaller than the lithospheric plates themselves.     

地面高精度磁测在西藏恰功铁矿勘查中的应用效果

恰功铁矿是近几年新发现的矿山,然而矿体的形态、大小、埋深、位置、产状、边界等几何特征还没有被认识清楚。为了弄清这些问题,本文通过欧拉反褶积方法对磁异常化极数据进行了反演,反演矿体深度为0-120m;在C-6异常中心选择了两条剖面进行了2.5维拟合反演,反演矿体厚度为20~30m,欧拉反褶积反演结果和2.5维拟合反演结果与ZK32钻孔验证结果相吻合。最后,通过建立恰功矽卡岩型铁矿地质—地球物理找矿模型,为该地区寻找隐伏夕卡岩型铁矿提供了思路。     

基于GRACE卫星和GLDAS系统的地下水水位估算模型——以和田地区克里雅河流域为例

地球重力场的变化是导致陆地水储量变化的重要因素之一,利用GRACE(Gravity Recovery and Climate Experiment)重力场恢复与气候实验重力卫星数据,结合GLDAS(Global Land Data Assimilation Systems)全球陆面数据同化系统和实测地下水位数据,反演和田地区克里雅河流域11年间四季和田地区的陆地水储量动态变化,模拟计算地下水等效水高变化趋势,构建了地下水水位估算模型。研究结果表明：和田地区春、夏两季的陆地水储量呈现出增加趋势,而秋、冬两季出现亏损状态;GRACE地球重力卫星所反演的陆地水储量比GLDAS同化系统所模拟的水资源变化更为剧烈,但2类数据的动态变化拟合度很高;GLDAS水资源等效水高二阶微分、GLDAS水资源变化倒数一阶微分、GRACE陆地水储量变化倒数变化、地下水储量变化一阶微分的敏感程度最高,构建的多元逐步回归模型明显优于线性函数,且水位深度越浅,该估算模型的适用性越高。     

Lithotype discrimination within a gridded model by inverting gravity data

In this paper we present a method which allows delineation of geologic structures in a bi-modal lithotype setting. We propose to use gravity data in combination with a priori information about the density contrast between the two lithotypes. The iterative method uses an objective function with five tunable parameters which need to be set. Using an efficient parameter search, suitable ranges of these are investigated to determine their optimal values, respectively, which in turn, ensures good inversion results.
The approach produces structural images of the subsurface, without the need of an a priori density model; the depth to the top of the inhomogeneity is also retrieved.
Besides synthetic simulations, the methodology has also been applied to a small gravity data set, acquired by the industry over a basinal structure. A consistent, bi-modal image of the bedrock depression is obtained from the data, which, in this case, was the goal. Other potential areas of application include delineation of salt structures and ore deposits.