The reduction of aliasing in gravity anomalies and geoid heights using digital terrain data

Observations of gravity can be aliased by virtue of the logistics involved in collecting these data in the field. For instance, gravity measurements are often made in more accessible lowland areas where there are roads and tracks, thus omitting areas of higher relief in between. The gravimetric determination of the geoid requires mean terrain-corrected free-air anomalies; however, anomalies based only on the observations in lowland regions are not necessarily representative of the true mean value over the topography. A five-stage approach is taken that uses a digital elevation model, which provides a more accurate representation of the topography than the gravity observation elevations, to reduce the unrepresentative sampling in the gravity observations. When using this approach with the Australian digital elevation model, the terrain-corrected free-air anomalies generated from the Australian gravity data base change by between 77.075 and −84.335 mgal (−0.193 mgal mean and 2.687 mgal standard deviation). Subsequent gravimetric geoid computations are used to illustrate the effect of aliasing in the Australian gravity data upon the geoid. The difference between 'aliased' and 'non-aliased' gravimetric geoid solutions varies by between 0.732 and −1.816 m (−0.058 m mean and 0.122 m standard deviation). Based on these conceptual arguments and numerical results, it is recommended that supplementary digital elevation information be included during the estimation of mean gravity anomalies prior to the computation of a gravimetric geoid model.     

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.     

Geodetic singularities in the gravity field of a non-homogeneous planet

Summary. An inverse geodetic singularity problem is considered for a non-homogeneous spherical planet. The singularity condition is expressed in terms of the density distribution and of the geometrical parameters of the configuration; the condition for the density distribution is deduced which gives rise to singularities of parabolic type in the external gravity field of the planet. The structure of the gravity field in the neighbourhood of the singularities is investigated in detail together with the behaviour of the gradients of the disturbances in the geodetic coordinates.     

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.     

The continuation inverse problem revisited

The non-uniqueness of the continuation of a finite collection of harmonic potential field data to a level surface in the source-free region forces its treatment as an inverse problem. A formalism is proposed for the construction of continuation functions which are extremal by various measures. The problem is cast in such a form that the inverse problem solution is the potential function on the lowest horizontal surface above all sources, serving as the boundary function for the Dirichlet problem in the upper half-plane. The desired continuation, at the higher level of interest, must then be in the range of the upward continuation operator acting on this boundary function, rather than being allowed the full freedom of itself being part of a Dirichlet problem boundary function. Extremal solutions minimize non-linear functionals of the continuation function, which are re-expressed as different functionals of the boundary function. A crux of the method is that there is no essential distinction between the upward and downward continuation inverse problems to levels above or below data locations. Casting the optimization as a Lagrange multiplier problem leads to an integral equation for the boundary function, which is readily solved in the Fourier domain for a certain class of functionals. The desired extremal continuation is then given by upward continuation. It is found that for some functionals, application of the Lagrange multiplier theorem requires a further restriction on the set of allowable boundary functions: bandlimitedness is a natural choice for the continuation problem. With this imposition, the theory is developed in detail for semi-norm functionals penalizing departure from a constant potential, in the 2-norm and Sobelev norm senses, and illustrated by application for a small synthetic Deep Tow magnetic field data set.     

Size distribution of urban blocks in the Tokyo Metropolitan Region: estimation by urban block density and road width on the basis of normative plane tessellation

The size distribution of urban blocks is important for the characterisation of urban block patterns and is known to follow several parametric statistical distributions. However, it has not previously been analysed on the basis of a normative plane tessellation and in terms of urban block density and mean road width. In this article, we formulate the size distribution of Voronoi cells using the gamma distribution estimated by urban block density and mean road width. We found that (1) both log-normal and gamma distributions can be good candidates for the size distribution of urban blocks at the scale of a region that consists of regular urban blocks and that has a uniform road width; and (2) the size distribution of urban blocks depends on the balance between pattern and width improvement effects. Based on one study region in Tokyo, if the pattern improvement effect is more prominent than the width improvement effect, the mode of urban block sizes tends to be larger than if it is not. These findings are expected to provide scientific support for urban planning (e.g. land readjustment projects).     

A geometrical interpretation of ideal body problems

Summary. For linear geophysical inverse problems, the exercise of finding a greatest lower bound on the uniform norms of positive solutions fitting N data, is shown to have a geometrical counterpart in the N- dimensional space of N -tuples of real numbers. By application of the Fenchel Duality Theorem, we demonstrate that the problem is equivalent to the discovery of a particular hyperplane tangent to a convex set in this space. As examples in the case of two data, the new formulation is applied to the problems of recovering density information from planetary mass and moment of inertia, and from two vertical gravity anomalies.     

Characterization of the depocenters and the basement structure,below the central Chile Andean Forearc: A 3D geophysical modelling in Santiago Basin area

Since the last century, several geological and geophysical studies have been developed in the Santiago Basin to understand its morphology and tectonic evolution. However, some uncertainties regarding sedimentary fill properties and possible density anomalies below the sediments/basement boundary remain. Considering that this is an area densely populated with more than 6 million inhabitants in a highly active seismotectonic environment, the physical properties of the Santiago Basin are important to study the geological and structural evolution of the Andean forearc and to characterize its seismic response and related seismic hazard. Two and three‐dimensional gravimetric models were developed, based on a database of 797 compiled and 883 newly acquired gravity stations. To produce a well‐constrained basement elevation model, a review of 499 wells and 30 transient electromagnetic soundings were used, which contribute with basement depth or minimum sedimentary thickness information. For the 2‐D modelling, a total of 49 gravimetric profiles were processed considering a homogeneous density contrast and independent regional trends. A strong positive gravity anomaly was observed in the centre of the basin, which complicated the modelling process but was carefully addressed with the available constrains. The resulting basement elevation models show complex basement geometry with, at least, eight recognizable depocenters with maximum sedimentary infill of ~ 500 m. The 3‐D density models show alignments in the basement that correlates well with important intrusive units of the Cenozoic and Mesozoic. Along with interpreted fault zones westwards and eastwards of the basin, the observations suggest a structural control of Santiago basin geometry, where recent deformation associated with the Andean contractional deformation front and old structures developed during the Cenozoic extension are superimposed to the variability of river erosion/deposition processes.     

The solution of the inverse gravity potential problem for cylindrical bodies

Summary. The inverse gravity potential problem consists in the determination of the form and the density of the body by its exterior gravity potential. We describe two similar classes of bodies for which this problem has a unique constructive solution.
(1) The first class contains the cylindrical bodies with finite length, arbitrary form of section and ρ( R , ø, z) =ρ₁( z )ρ₂( R , ø) density distribution, where z is the cylindrical coordinate; R , ø are the polar coordinates in a section plane. This class is important for prospecting geophysics in that it allows us to determine in a unique and constructive way, the function ρ₁( R , ø), the length, form and orientation of the cylinder if we know the function ρ₁( z ) and the exterior potential. The classical moment problem of functions is the basis for the solution of this problem.
(2) The analogous problem for the class of the spherical cylinders, or bodies bounded by arbitrary similar sections of two different concentric spheres and the radial lateral surface, appears when bodies of planetary size are studied. (An example of these bodies would be the Moon mascons.) The density distribution of these cylinders is ρ(τ, θ, ø) =ρ₁(τ)ρ₂(θ, ø) where τ, θ, ø are the spherical coordinates. The function ρ₁(θ, ø), length and form of spherical sections can be uniquely determined by exterior potential if we know the function ρ₁(τ). We propose a new constructive method for harmonic continuation of the gravity potential into the region containing the perturbing masses for the solution of the problem.     

The Free Air Geoid for Australia

Summary The free air geoid, which is the co-geoid obtained by the use of free air anomalies in Stokes' integral, is computed for Australia from available gravity data. The set of anomalies used to represent the outer zones had been obtained previously using a combined solution from satellite data and terrestrial gravimetry. The solutions so obtained for the free air geoid are compared with the astrogeodetic determination of the geoid on the Australian Geodetic Datum by Fischer and Slutsky and the accuracy of the comparisons is estimated.     

The relationship between seismic velocity and density in the continental crust — a useful constraint?

Summary . Plots of seismic velocity and density of rock samples show that a range of densities is possible for rocks of each seismic velocity and vice versa. although a single linear relationship is often assumed in crustal gravity calculations. Because of the scatter, whenever rocks of known seismic velocity are converted to density using this relationship, a reduction is made to the resolving power of the resulting gravity calculation. If these rocks reach thicknesses of more than a few kilometres, then the uncertainties become significant when compared with the size of commonly observed gravity anomalies. Examples are considered from the North Sea, Mississippi and Carolina Trough. It is concluded that the use of a seismic velocity measurement as the only indication of rock density does not provide a useful constraint when attempting to reproduce observed gravity variations. An appropriate model for isostatic compensation is probably the most important factor for successful predictions of crustal structure on the basis of gravity data.     

Tidal gravity in Britain: tidal loading and the spatial distribution of the marine tide

Summary. Tidal gravity measurements have been made at six sites in Britain with two nulled LaCoste and Romberg Earth tide gravitymeters. The M ₂ observations from these and two further sites are compared with calculations of the tidal loading from the seas around the British Isles and the major oceans. Models of the M ₂ marine tides are convolved with Green's functions for appropriate radially stratified Earth models. The differences between the M ₂ observations and the theoretical calculations are less than 0.6 μ gals and it is shown that these differences contain further information concerning the errors in the marine tide models. The M ₂ marine tides on the north-west European continental shelf are reasonably well known and this allows a useful test of the feasibility of using tidal gravity measurements for the inverse ocean tide problem in areas where the ocean tides are less well known. The differential gravity loading signal between pairs of gravity stations is shown to be important for considerations of the uniqueness and accuracy of the inverse problem. M ₂ tidal gravity loading maps for the British Isles and Europe have been produced which are of use in making corrections to various geodetic measurements.     

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.     

Geomorphometry on the surface of a triaxial ellipsoid: towards the solution of the problem

Existing algorithms of geomorphometry can be applied to digital elevation models (DEMs) given with plane square grids or spheroidal equal angular grids on the surface of an ellipsoid of revolution or a sphere. Computations on spheroidal equal angular grids are trivial for modelling of the Earth, Mars, the Moon, Venus, and Mercury. This is because: (a) forms of these celestial bodies can be described by an ellipsoid of revolution or a sphere and (b) for these surfaces, there are well-developed theory and algorithms to solve the inverse geodetic problem as well as to determine spheroidal trapezoidal areas. It is advisable to apply a triaxial ellipsoid for describing the forms of small moons and asteroids. However, there are no geomorphometric algorithms intended for such a surface. In this article, first, we formulate the problem of geomorphometric modelling on a triaxial ellipsoid surface. Then, we recall definitions and formulae for coordinate systems of a triaxial ellipsoid and their transformation. Next, we present analytical and computational solutions, which provide the basis for geomorphometric modelling on the surface of a triaxial ellipsoid. The Jacobi solution for the inverse geodetic problem has a fundamental mathematical character. The Bespalov solutions for determination of the length of meridian/parallel arcs and the spheroidal trapezoidal areas are computationally efficient. Finally, we describe easy-to-code algorithms for derivation of local and non-local morphometric variables from DEMs based on a spheroidal equal angular grid of a triaxial ellipsoid.     

Preliminary Results on the 3-dimensional Seismic Structure of the Lithosphere under the USGS Central California Seismic Array

About 1500 readings of teleseismic P -time residuals obtained from the US Geological Survey seismograph network in central California have been used to obtain a three-dimensional image of seismic velocity anomalies for this area by the method of Aki, Christoffersson & Husebye We found that the California network is less suitable than the LASA and NORSAR arrays for this kind of studies because of its greater proportion of peripheral blocks in which the resolution is very poor for the stochastic inverse solution and the random error effect is severe for the generalized inverse solution. Nevertheless, the resultant velocity anomalies show a remarkable correlation with the San Andreas fault zone to a depth of 75 km. The anomaly pattern changes drastically as the depth exceeds 75 km, suggesting that the asthenosphere has been reached.     

Calculation of the variance of mean values for blocks in regional resource evaluation studies

Random sampling of a known covariance function can be used during the process of estimating the variance of means or totals for a spatial random variable in blocks of variable size. One advantage of this method is that the precision of any block variance can be determined at the same time as the integral itself. In two-dimensional space this approach yields sufficiently precise results for continuous spatial random variables with exponential, Gaussian, and spherical covariance functions, as well as for point patterns with exponential covariance density or power-law-type, second-order intensity function. Practical examples of application deal with the areal distribution of felsic volcanic rocks and gold deposits in the Abitibi Volcanic Belt, Canadian Shield. The exponential model yields good results in both cases, but, as an overall fit, the fractal (power-law) model performs better in the characterization of the two-dimensional distribution of the gold deposits.     

Detailed 1×1° gravimetric Indian Ocean geoid and comparison with GEOS-3 radar altimeter geoid profiles

Summary. A new set of 1×1° mean free-air anomalies in the Indian Ocean is determined on the basis of previously published free-air anomaly maps (Talwani & Kahle) and the most recent Lamont surface ship gravity measurements. The data are then used to compute a (total) 1×1° gravimetric Indian Ocean geoid. The computation is carried out by combining the Goddard Space Flight Center (GSFC) GEM-6 geoid and a difference geoid that corresponds to the differences between the set of 1×1° surface gravity values and the GEM-6 gravity anomalies. The difference geoid is highest over the Madagascar Ridge (+ 20 m) and lowest over the Timor Trough (-30 m). The total geoid is compared with GEOS-3 radar altimeter derived geoid profiles and geophysical implications are discussed.     

Source ambiguity from an estimation of the scaling exponent of potential field power spectra

An analysis of the field scaling power spectrum yields useful information about the source distribution, but it is uncertain whether deterministic, random, fractal or mixed approaches have to be used for the interpretation. To this end, the scaling properties of potential field spectra are analysed for a number of different source models of geological interest. Besides the models of Naidu (purely random sources) and Spector and Grant (gross block statistical ensembles) we consider other types of density and magnetization distributions with spectral exponents in the fractal range, such as a single homogeneous body with a random white source distribution. Spectral slopes in the fractal range are obtained.
We also study the effects of important natural sources, such as salt domes and sedimentary basins, representing them with simple Gaussians or combinations of Gaussian signals. The same spectral slopes as for gravity signals generated by 3-D fractal source distributions are found for them. Hence the power law decay of the field is not a characteristic only of fractal source models.
If a 3-D fractal source distribution is assumed a priori , a way of verifying the goodness of the model is to examine the whitened field at source level. The probability that the whitened field derives from a random white population is estimated for synthetic and real anomalies by applying the usual statistical tests.     

Estimation of the pole tide gravimetric factor at the chandler period through wavelet filtering

Wavelet analysis for filtering is used to improve estimation of gravity variations induced by Chandler wobble. This method eliminate noise in superconducting gravimeter (SG) records with bandpass filters derived from Daubechies wavelet. The SG records at four European stations (Brussels, Membach, Strasbourg and Vienna) are analysed in this study. First, the earth tidal constituents are removed from the observed data by using synthetic tides, then the gravity residuals are filtered into a narrow period band of 256–512 d by a wavelet bandpass filter. These data are submitted to three regression analysis methods for estimating the gravimetric factor of the Chandler wobble. After processing by wavelet filtering, SG records can provide amplitude factors δ and phase lags κ of the Chandler wobble with much smaller mean square deviation (MSD) than these provided by former studies. It is mainly because the wavelet method can effectively eliminate instrumental drift and provide smoothed data series for the regression analysis.     

A new method for computing synthetic seismograms

Summary. The computation of theoretical seismograms for models in which the elastic parameters and density vary only with depth (in a plane, cylindrical or spherical geometry) reduces to the solution of an ordinary differential equation plus the evaluation of inverse transformations. In principle, the problem is straightforward. In practice, many techniques and approximations can be used at each stage and many combinations and variants are possible. In this paper, we discuss a new method of evaluating the inverse transforms. Any method can be used to solve the differential equation and we only discuss a few analytic approximations to illustrate the new method. The inverse transformations are a frequency and wavenumber integral. Essentially four techniques can be used to evaluate these depending on the order of integration and whether the wavenumber integral is distorted from the real axis. Three of these have been widely used, but the technique of evaluating the frequency integral first and keeping the wavenumber real is new. In this paper, we discuss some of the advantages of this combination.