Spectral harmonic analysis and synthesis of Earth’s crust gravity field

We developed and applied a novel numerical scheme for a gravimetric forward modelling of the Earth’s crustal density structures based entirely on methods for a spherical analysis and synthesis of the gravitational field. This numerical scheme utilises expressions for the gravitational potentials and their radial derivatives generated by the homogeneous or laterally varying mass density layers with a variable height/depth and thickness given in terms of spherical harmonics. We used these expressions to compute globally the complete crust-corrected Earth’s gravity field and its contribution generated by the Earth’s crust. The gravimetric forward modelling of large known mass density structures within the Earth’s crust is realised by using global models of the Earth’s gravity field (EGM2008), topography/bathymetry (DTM2006.0), continental ice-thickness (ICE-5G), and crustal density structures (CRUST2.0). The crust-corrected gravity field is obtained after modelling and subtracting the gravitational contribution of the Earth’s crust from the EGM2008 gravity data. These refined gravity data mainly comprise information on the Moho interface and mantle lithosphere. Numerical results also reveal that the gravitational contribution of the Earth’s crust varies globally from 1,843 to 12,010 mGal. This gravitational signal is strongly correlated with the crustal thickness with its maxima in mountainous regions (Himalayas, Tibetan Plateau and Andes) with the presence of large isostatic compensation. The corresponding minima over the open oceans are due to the thin and heavier oceanic crust.     

On the local multiscale determination of the Earth’s disturbing potential from discrete deflections of the vertical

As a first approximation, the Earth is a sphere; as a second approximation, it may be considered an ellipsoid of revolution. The deviations of the actual Earth’s gravity field from the ellipsoidal “normal” field are so small that they can be understood to be linear. The splitting of the Earth’s gravity field into a “normal” and a remaining small “disturbing” field considerably simplifies the problem of its determination. Under the assumption of an ellipsoidal Earth model, high observational accuracy is achievable only if the deviation (deflection of the vertical) of the physical plumb line, to which measurements refer, from the ellipsoidal normal is not ignored. Hence, the determination of the disturbing potential from known deflections of the vertical is a central problem of physical geodesy. In this paper, we propose a new, well-promising method for modelling the disturbing potential locally from the deflections of the vertical. Essential tools are integral formulae on the sphere based on Green’s function with respect to the Beltrami operator. The determination of the disturbing potential from deflections of the vertical is formulated as a multiscale procedure involving scale-dependent regularized versions of the surface gradient of the Green function. The modelling process is based on a multiscale framework by use of locally supported surface curl-free vector wavelets.      

Lunar ranging instrument for Chandrayaan-1

Lunar Laser Ranging Instrument (LLRI) proposed for the first Indian lunar mission Chandrayaan-1 is aimed to study the topography of the Moon’s surface and its gravitational field by precisely measuring the altitude from a polar orbit around the Moon. Altimetry data close to the poles of the Moon would also be available from the instrument, which was not covered by earlier missions. This instrument supplements the terrain mapping camera and hyperspectral imager payloads on Chandrayaan-1. The instrument consists of a diode pumped Nd:YAG pulsed laser transmitter having 10 nsec pulse width and a receiver system. The receiver system features 17 cm diameter Ritchey—Chrétien collecting optics, Si Avalanche Photo Detector (APD), preamplifiers, constant fraction discriminators, time-of-flight measurement unit and spacecraft interface. Altimeter resolution of better than 5 m is targeted. The received signal strength of LLRI depends on laser pulse backscatter from the Moon’s surface. Moon’s surface being a poor reflector, the choice of receiver size and its type and the selection of detector play an important role in getting a good signal-to-noise ratio and in turn achieving the target resolution. At the same time, the spacecraft puts a limitation on payload size and weight. This paper discusses the proposed LLRI system for Chandrayaan-1 and signal-to-noise ratio estimation.     

A seismic landslide hazard analysis with topographic effect,a case study in the 99 Peaks region,Central Taiwan

It has been known that ground motion amplitude will be amplified at mountaintops; however, such topographic effects are not included in conventional landslide hazard models. In this study, a modified procedure that considers the topographic effects is proposed to analyze the seismic landslide hazard. The topographic effect is estimated by back analysis. First, a 3D dynamic numerical model with irregular topography is constructed. The theoretical topographic amplification factors are derived from the dynamic numerical model. The ground motion record is regarded as the reference motion in the plane area. By combining the topographic amplification factors with the reference motions, the amplified acceleration time history and amplified seismic intensity parameters are obtained. Newmark’s displacement model is chosen to perform the seismic landslide hazard analysis. By combining the regression equation and the seismic parameter of peak ground acceleration and Arias intensity, the Newmark’s displacement distribution is generated. Subsequently, the calculated Newmark’s displacement maps are transformed to the hazard maps. The landslide hazard maps of the 99 Peaks region, Central Taiwan are evaluated. The actual landslide inventory maps triggered by the 21 September 1999, Chi-Chi earthquake are compared with the calculated hazard maps. Relative to the conventional procedure, the results show that the proposed procedures, which include the topographic effect can obtain a better result for seismic landslide hazard analysis. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Unraveling reservoir compaction parameters through the inversion of surface subsidence observations

In an attempt to derive more information on the parameters driving compaction, this paper explores the feasibility of a method utilizing data on compaction-induced subsidence. We commence by using a Bayesian inversion scheme to infer the reservoir compaction from subsidence observations. The method’s strength is that it incorporates all the spatial and temporal correlations imposed by the geology and reservoir data. Subsequently, the contributions of the driving parameters are unravelled. We apply the approach to a synthetic model of an upscaled gas field in the northern Netherlands. We find that the inversion procedure leads to coupling between the driving parameters, as it does not discriminate between the individual contributions to the compaction. The provisional assessment of the parameter values shows that, in order to identify adequate estimate ranges for the driving parameters, a proper parameter estimation procedure (Markov Chain Monte Carlo, data assimilation) is necessary.     

Benford’s Law Applied to Hydrology Data—Results and Relevance to Other Geophysical Data

Benford’s Law gives the expected frequencies of the digits in tabulated data and asserts that the lower digits (1, 2, and 3) are expected to occur more frequently than the higher digits. This study tested whether the law applied to two large earth science data sets. The first test analyzed streamflow statistics and the finding was a close conformity to Benford’s Law. The second test analyzed the sizes of lakes and wetlands, and the finding was that the data did not conform to Benford’s Law. Further analysis showed that the lake and wetland data followed a power law. The expected digit frequencies for data following a power law were derived, and the lake data had a close fit to these expected digit frequencies. The use of Benford’s Law could serve as a quality check for streamflow data subsets, perhaps related to time or geographical area. Also, with the importance of lakes as essential components of the water cycle, either Benford’s Law or the expected digit frequencies of data following a power law could be used as an authenticity and validity check on future databases dealing with water bodies. We give several applications and avenues for future research, including an assessment of whether the digit frequencies of data could be used to derive the power law exponent, and whether the digit frequencies could be used to verify the range over which a power law applies. Our results indicate that data related to water bodies should conform to Benford’s Law and that nonconformity could be indicators of (a) an incomplete data set, (b) the sample not being representative of the population, (c) excessive rounding of the data, (d) data errors, inconsistencies, or anomalies, and/or (e) conformity to a power law with a large exponent.     

Increase in thermal groundwater due to a flowing well near the Songshan hot spring in Beijing,China

The Songshan hot spring occurs in granite in Yanqing County in northwestern Beijing, China. TDS of the hot water ranges from 459 to 475 mg/L and pH varies between 8.6 and 9.13. The water is of Na–SO₄ type. Isotopic analyses indicate that the hot spring is meteoric in origin and receives recharge from precipitation in the northern and northwestern granite mountain with elevation of about 1,600–1,800 m. The depth of circulation of the thermal groundwater is estimated to be 2,240 m below the spring’s threshold and the temperature of the geothermal reservoir, 76°C. The residence time of the thermal groundwater is estimated to be about 52 years. A flowing well near the spring has chemical compositions and formation conditions similar to the spring. The discharge of the flowing well is approximately eighteen times larger than that of the spring and the residence time of the former (about 15.4 years) is about three times smaller than that of the latter. Although the well and spring are close to each other, the well’s larger flow rates, indicated residence time and high hydraulic head suggest that the well taps a separate, but genetically similar flow system.     

Flexural modeling of the neoproterozoic Araguaia belt,central Brazil

Flexural modeling of bending of the southern and southeastern borders of the Amazon lithospheric plate under the western border of the Goiás Massif and western Parnaı́ba basin was constrained by 1070 gravity stations between 5°–14°S and 46°–52.5°W. Topography and aeromagnetic data were also used to estimate the loads of the Araguaia thrust belt. A sequence of Bouguer gravity anomaly lows (−80 to −40 mGal) is located over the Araguaia thrust belt and Cenozoic sediments of the Ilha do Bananal basin. Bouguer anomalies over the Amazon craton, to the west of the thrust belt, are higher than −20 mGal. Towards the east, over the Goiás Massif, the São Francisco craton and the Paleozoic to Mesozoic Parnaı́ba basin, anomalies range from −70 to −20 mGal. Comparison between topography and gravity along profiles perpendicular to the cratonic borders and across the Araguaia thrust belt shows that the long-wavelength gravity anomalies are best explained by bending of the Amazon plate caused by loads such as the observed topography, the thrust-sheets of the Araguaia belt and the remnants of ancient island-arc system in the Goiás massif. The thickness of the Araguaia thrust belt together with the Cenozoic sediments was estimated using aeromagnetic data and it ranges from 6 to 8 km. This load was used to calculate the minimum effective elastic thickness T_e for the Amazon plate. T_e=80 km was estimated by comparing the observed Bouguer anomalies with the gravity anomalies caused by bending of the crust-mantle interface of a broken elastic plate model. These results support the proposition that the Araguaia belt formed during the collision and suture of the Amazon and the São Francisco lithospheric plates, in late Proterozoic times.     

Benford’s Law in Time Series Analysis of Seismic Clusters

Benford’s analysis is applied to the recurrence times of approximately 17,000 seismic events in different geological contexts of Italy over the last 6 years, including the Mt. Etna volcanic area and the seismic series associated with the destructive M _w 6.3, 2009 L’Aquila earthquake. A close conformity to Benford’s law and a power-law probability distribution for the recurrence times of consecutive events is found, as typical of random multiplicative processes. The application of Benford’s law to the recurrence event times in seismic series of specific seismogenic regions represents a novel approach, which enlarges the occurrence and relevance of Benford-like asymmetries, with implications on the physics of natural systems approaching a power law behaviour. Moreover, we propose that the shift from a close conformity of Benford’s law to Brownian dynamics, observed for time separations among non-consecutive events in the study seismic series, may be ruled by a periodical noise factor, such as the effects of Earth tides on seismicity tuning.     

Using Newton's law and geophysical bounds on mass density contrast to ensure consistency between gravity and height data

In this paper we advocate the use of Newton's law of gravitational attraction to ensure perfect consistency between gravity and height data. Starting with the absolute gravity on the topography we decompose this signal into a number of quantities associated with physics of the system. To model gravitational attraction from topography we use DTM and Newton's law of gravitational attraction. A residual part of the gravity signal is interpreted as inconsistency between gravity and heights. In the paper we discuss a method by which such inconsistency (at least in principle) can be decomposed into a “gravity error” and a “terrain error”. In practice such separation is not possible because the two types of error are nearly 100% correlated. The inconsistency can be interpreted as a measure of ambiguity of the gravity-terrain models which are consistent with a set of measured/interpolated data. We discuss the influence of such ambiguity on the accuracy of the geoid for the investigated area of Jutland, Denmark.     

Geophysical studies of crustal structure of the Ongul Islands and the Northern Mizuho Plateau, East Antarctica

Explosion seismic experiments, gravity measurements and aeromagnetic surveys were made in the northern Mizuho Plateau including the Ongul Islands, East Antarctica, from 1979 to 1982 by the Japanese Antarctic Research Expeditions. The objective of these field operations was to determine the crustal structure along the 300 km-long oversnow traverse route between Syowa and Mizuho Stations. Three big shots were fired; at sea near Syowa Station, in an ice hole near Mizuho Station and in an ice hole between both stations. Twenty-seven temporal seismic stations were set up along the route. Gravity measurements were carried out at 30 points along this route. Aeromagnetic surveys over the area were made four times.In the seismic experiments, clear refracted waves from the Conrad (estimated depth 30 km) and the Moho (estimated depth 40 km) discontinuities were recorded. No layer with a velocity of less than 6 km/s was found in the Ongul Islands nor beneath the ice sheet in the surveyed area. The P-wave velocity in the upper layer varies with depth from 6.0 km/s on the surface to 6.4 km/s at a depth of 13 km. Comparing the observed record section with synthetic seismograms, it was derived that the Conrad was not associated with a sharp velocity discontinuity, but a linear velocity increase of 0.55 km/s in a transition zone of 2.4 km thick. Velocities of P^* and P_n were determined as 6.95 km/s and 7.93 km/s assuming a flat layered structure.Bouguer gravity anomalies could not be calculated along the whole profile because of a lack of data on bedrock topography, so reduced gravity anomalies were calculated. These anomalies indicate no abrupt changes of the bedrock topography.     

Global earth isostatic model using smoothed Airy-Heiskanen and Vening Meinesz hypotheses

Isostatic hypotheses are used for different purposes in geophysics and geodesy. The Erath crustal thickness modelling is more complicated than the classical isostatic models. In this study we try to modify Airy-Hesiskanen model, utilizing a smoothing factor, to a model with regional or global isostatic model through a modern solution of the gravimetric-isostatic Vening Meinesz model and CRUST.0. In Airy-Hesiskanen’s theory there is no correlation between neighbouring crustal columns, while this must be the case in reality due to the elasticity of the Earth. The idea is to keep the simplicity of the Airy-Heiskanen model, because it needs only the topographic data, and change the model which becomes to a model with regional/global isostatic model. The isostatic assumption for compensating the topographic potential is incomplete, as there are other geophysical phenomena which should be considered. Using the isostatic hypothesis for determining the depth of crust causes some disturbing signals in the gravity anomaly (approximately 285 mGal), which influence the crustal thickness determination. In this paper a simple method use for removing these effects. Spherical harmonic potential coefficients of the topographic compensation masses are used for modifying Airy-Heiskanen’s model in a least-square adjustment procedure by estimating smoothing factor. The numerical analysis shows that below degree 10, the modified Airy-Hesiskanen and Vening Meinesz models are close together. Smoothing factors for modifying the Airy-Hesiskanen model vary from 0.75 to 0.64 between degrees 200 and 2159.     

Modelling the geoid and sea-surface topography in coastal areas

This paper looks at the relation between the time-averaged level of the sea surface and a gravimertic geoid, as determined in coastal areas. Measurements in local regions can now be accurate enough to demonstrate that the geoid and mean sea level are not even parallel to each other, let alone identical. The accuracy and pattern structure of surface gravity data in some shelf seas is comparable with those on land, so that a marine geoid can be derived from surface data without using satellite altimetry. The geodetic objective is then to combine the two to determine sea surface topography. In principle, gravimetric studies provide the absolute datum so that local oceanographic models on the shelf can be combined with sea surface topography models related to the global ocean circulation. In contrast, sea surface topography information near deep ocean coasts must come from external sources and satellite altimetry used to give the gravity data needed to offset the less good coverage by ship-borne gravimetry.Marine Bouguer anomalies enable two specific problems of gravity anomaly patterns near the continent ocean transition to be overcome. The necessary extension of Stokes' condensation reduction is developed and illustrated along a north-south profile from the Mediterranean across the Cote d'Azur. The effect on gravity of deep ocean water introduces a geoid correction in the form of a dipolar ridge whose amplitude at the shore is about 11 cm. In addition to geostrophic currents, a semi-quantitative model for the thermohaline effects on sea surface topography is discussed in relation to sea level differences between the Atlantic and Mediterranean.In considering appropriate algorithms for local geoid computation, Kirby's Iterative Fourier Combination routine for combining altimetry and surface gravity is extended to account for global sea surface topography. The impact of very fast spherical harmonic analysis algorithms is discussed and a simple physical model is given which explains the short coherence lengths found for the global gravity field. This necessary assumption for any local geoid computation was hitherto purely empirical.Finally, the use of land data such as tide gauges, ellipsoidal heights from GPS, and orthometric heights from first order levelling are reviewed as ways of corroborating geodetic estimates of sea surface topography and its relation to levelling datums. Successful examples are given from southern England.     

Energetic analysis of global topography

Diagrams of the various distributions of the potential energy of the continental and oceanic-bottom topography as a quantitative indices of energy consumption by tectonic processes have been calculated and plotted for the first time on the basis of high-resolution ETOPO 2′ digital data on a 2′ × 2 ′ uniform network with height intervals of 100 m. The potential energy of the continents and oceanic bottom is mainly stored in mountainous regions of Asia and South America with an anomalous type of energy distribution and in mid-oceanic ridges, respectively. The topography of the oceanic bottom bears five times more energy than the continental topography. With allowance for the isostatic compensation of denudation, the bulk reserve of the potential energy of the global topography is estimated at 10³³ erg, a value which amounts to only an insignificant fraction of other items in the energetic budget of the Earth. Fifty percent of this energy has been gained over the last 70 Ma. Currently, the rate of accumulation of the potential energy of the topography is 10²⁵ erg/yr, i.e., 1000 times less than the bulk geothermal flow and 10 times less than energy released by all earthquakes during a year. So far, there are no grounds to bind the potential energy of the Earth’s topography with a specific source.     

Old inherited origin for the present near-bimodal topography of Africa

Africa’s landscape is dominated by a manifold of second-order epeirogenic structures superimposed on a first-order bimodal topography. Bivariate regression analysis of Africa’s surface topography shows that this is a complexly folded surface with regionally elevated areas in southern and eastern Africa, and a topographically low northern and western Africa. The apparent spatial relationships between these features are analysed using anomaly correlation between surface topography and free-air gravity anomalies. Occurrences of positively correlated features between gravity and topography in Africa are found to be limited to second-order epeirogenic features. Geophysical modelling and geologic evidence indicate that Africa’s bimodal topography is genetically distinct from these second-order features, and linked to sources as deep as the sublithospheric mantle. The age, measured and modelled elevation of the bimodal topography require that topographic uplift of south-central Africa be episodic. We infer from our findings together with relative sea-level changes, that the near-bimodality of Africa’s topography is an ancient feature inherited at least from upper Paleozoic times. Our reconstructed paleotopography suggests that Africa was largely a low-lying continent dominated by its cratons, and that basement distribution disregards the present-day uplift patterns of Africa.     

Gravity field recovery from airborne gradiometer data using collaction and taking into account correlated errors

Least squares collocation is a very comprehensive method for gravity field modelling, since it may use known noise characteristics of the data. In many earlier applications the errors affecting the data were considered uncorrelated, mainly due to the difficulty in estimating the systematic character of such kind of errors. In this study, error covariance functions of airborne gravity gradiometer data are estimated by comparing model covariance functions with empirical covariance functions of the gravity gradiometer data. The model covariance functions were estimated from accurate surface gravity data and continuated upward to the height of the airborne measurements using the covariance propagation law. The estimated error covariance functions were modeled as finite ones and used as an additional information for the prediction of gravity anomalies from gravity gradiometer data. The assessment of the prediction results was made by comparing the gravity values predicted from the airborne gradient data and showed up to 25% improvement compared to not using correlated errors.     

A Practical Procedure for the Measurement of Fragmentation by Blasting by Image Analysis

Summary. The operation of a digital image analysis system in a limestone quarry is described. The calibration of the system, required in order to obtain moderately reliable fragmentation values, is done from muckpile sieving data by tuning the image analysis software settings so that the fragmentation curve obtained matches as close as possible the sieving. The sieving data have also been used to extend the fragment size distribution curves measured to sizes below the system’s optical resolution and to process the results in terms of fragmented rock, discounting the material coming from a loose overburden (natural fines) that is cast together with the fragmented rock. Automatic and manual operation modes of the image analysis are compared. The total fragmentation of a blast is obtained from the analysis of twenty photographs; a criterion for the elimination of outlier photographs has been adopted using a robust statistic. The limitations of the measurement system due to sampling, image processing and fines corrections are discussed and the errors estimated whenever possible. An analysis of consistency of the results based on the known amount of natural fines is made. Blasts with large differences in the amount of fines require a differentiated treatment, as the fine sizes tend to be the more underestimated in the image analysis as they become more abundant; this has been accomplished by means of a variable fines adjustment factor. Despite of the unavoidable errors and the large dispersion always associated with large-scale rock blasting data, the system is sensitive to relative changes in fragmentation.     

Estimation of Vertical Continuous Stochastic Parameters from Seismic Reflection Data

This paper presents a scheme to invert reflected seismic wave fields for vertical spatial statistics of the Earth’s continuous stochastic velocity distribution. The forward model is based on a single scattering theory, where a seismogram is modeled as a convolution of a source wavelet and stochastic reflectivity, which obeys a von Kármán autocorrelation function. Tests with synthetic data indicate that the method presented here recovers vertical characteristic length from seismic reflection data that is similar to that of the original velocity model; however, the inversion does not reliably recover the Hurst exponent, ν. The input wavelet’s power spectrum is a free parameter in the inversion and is estimated as well, but is band-limited relative to the original source wavelet. Noise degrades the quality of the estimates of the stochastic parameters such that the recovered stochastic parameters are biased according to the power spectrum of the noise.