Spherical harmonic analysis of the CRUST 2.0 global crustal model

The spherical harmonic analysis of the 2×2 density and stratification information contained in the global crustal model CRUST 2.0 is presented from the viewpoint of gravity field recovery and interpretation. Using a standard Airy/Heiskanen (A/H) isostatic hypothesis and a radially distributed compensation mechanism, two models of topographic/isostatic (t/i) potential harmonic coefficients are obtained up to degree and order 90. The CRUST-derived coefficients are compared with the spectrum of uncompensated topography, with the EGM96 observed gravity field, and with the t/i spectrum based on an A/H hypothesis with a constant compensation depth of 30 km. The signal degree variances of both CRUST models decrease quite smoothly towards degree 90, while the seven-layer model approaches the EGM96 spectrum for degrees 80–90. The significant deviation of the CRUST spectra from the A/H combined spectrum may prove of relevance to local and regional applications investigating the validity of current isostatic hypotheses.Acknowledgments. Sincere thanks go to Nikolaos Pavlis and three unknown referees for their thoughtful comments. Figure 1 has been produced using the mapping package m_map by R. Pawlowicz, which is a MATLAB toolbox that can be freely downloaded from http://www2.ocgy.ubc.ca/~rich/map.html     

Convergence and optimal truncation of binomial expansions used in isostatic compensations and terrain corrections

 A binomial expansion is a powerful tool in geodetic research. It is often used in terrain correction and isostatic compensation. The behaviour, convergence and truncation of the binomial expansion are investigated. The relation of the topographic height H (or the compensation depth), spherical harmonic degree n and the binomial series term m is discussed using theoretical and numerical results. According to the relation, a truncation number M is determined for obtaining an accuracy of 1%, i.e. it can be found how many terms (or power numbers of the topography) should be used in practical calculations. Received: 24 February 1999 / Accepted: 28 June 2000     

The exterior Airy/Heiskanen topographic–isostatic gravity potential, anomaly and the effect of analytical continuation in Stokes' formula

This study deals with the external type of topographic–isostatic potential and gravity anomaly and its vertical derivatives, derived from the Airy/Heiskanen model for isostatic compensation. From the first and the second radial derivatives of the gravity anomaly the effect on the geoid is estimated for the downward continuation of gravity to sea level in the application of Stokes' formula. The major and regional effect is shown to be of order H ³ of the topography, and it is estimated to be negligible at sea level and modest for most mountains, but of the order of several metres for the highest and most extended mountain belts. Another, global, effect is of order H but much less significant Received: 3 October 1997 / Accepted: 30 June 1998     

Determination of Potential coefficients to degree 52 from 5° Mean Gravity Anomalies

A set of 38406 1°×1° mean free air anomalies were used to derive a set of 1507 5° equal area anomalies that were supplemented by 147 predicted anomalies to form a global coverage of 1654 anomalies. These anomalies were used to derive potential coefficients to degree 52 using the summation formulae. In these computations, a smoothing operator was introduced and found to significantly effect the results at higher degrees. In addition, the effects of the atmosphere, spherical approximation and terrain were studied. It was found that the atmospheric effects and spherical approximation effects were about 0.3% of the actual coefficients. The terrain correction effects amounted to 10 to 25% of the low degree coefficients depending on a specific terrain correction model chosen; however, the correction terms found from the models did not yield solutions that agreed better with current satellite derived potential coefficient determinations. Anomalies were computed from the derived potential coefficients for comparison to the original anomalies. These comparisons showed that the agreement between the two anomalies became significantly better as the degree of expansion increased to the maximum considered. These comparisons shed some doubt on the rule of thumb that a block of size θ^° can be represented by a spherical harmonic expansion to 180^°/θ^°.     

Gravity field approximation in the area of Greece with emphasis on local characteristics

The accuracy of the gravity field approximation depends on the amount of the available data and their distribution as well as on the variation of the gravity field. The variation of the gravity field in the Greek mainland, which is the test area in this study, is very high (the variance of point free air gravity anomalies is 3191.5mgal ²). Among well known reductions used to smooth the gravity field, the complete isostatic reduction causes the best possible smoothing, however remain strong local anomalies which disturb the homogeneity of the gravity field in this area. The prediction of free air gravity anomalies using least squares collocation and regional covariance function is obtained within a ±4 ... ±19mgal accuracy depending on the local peculiarities of the free air gravity field. By taking into account the topography and its isostatic compensation with the usual remove-restore technique, the accuracy of the prediction mentioned obove was increased by about a factor of 4 and the prediction results become quite insensitive to the covariance function used (local or regional). But when predicting geoidal heights, in spite of using the smoothed field, the prediction results remain still depend on the covariance function used in such a way that differences up to about 50cm/100km result between relative geoidal heights computed with regional or local covariance functions.     

A new isostatic model of the lithosphere and gravity field

Based on the analysis of various factors controlling isostatic gravity anomalies and geoid undulations, it is concluded that it is essential to model the lithospheric density structure as accurately as possible. Otherwise, if computed in the classical way (i.e. based on the surface topography and the simple Airy compensation scheme), isostatic anomalies mostly reflect differences of the real lithosphere structure from the simplified compensation model, and not necessarily the deviations from isostatic equilibrium. Starting with global gravity, topography and crustal density models, isostatic gravity anomalies and geoid undulations have been determined. The initial crust and upper-mantle density structure has been corrected in a least squares adjustment using gravity. To model the long-wavelength (>2000 km) features in the gravity field, the isostatic condition (i.e. equal mass for all columns above the compensation level) is applied in the adjustment to uncover the signals from the deep-Earth interior, including dynamic deformations of the Earths surface. The isostatic gravity anomalies and geoid undulations, rather than the observed fields, then represent the signals from mantle convection and deep density inhomogeneities including remnants of subducted slabs. The long-wavelength non-isostatic (i.e. the dynamic) topography was estimated to range from –0.4 to 0.5 km. For shorter wavelengths (<2000 km), the isostatic condition is not applied in the adjustment in order to obtain the non-isostatic topography due to regional deviations from classical Airy isostasy. The maximum deviations from Airy isostasy (–1.5 to 1 km) occur at currently active plate boundaries. As another result, a new global model of the lithosphere density distribution is generated. The most pronounced negative density anomalies in the upper mantle are found near large plume provinces, such as Iceland and East Africa, and in the vicinity of the mid-ocean ridge axes. Positive density anomalies in the upper mantle under the continents are not correlated with the cold and thick lithosphere of cratons, indicating a compensation mechanism due to thermal and compositional density.     

Undulation and anomaly estimation using Geos-3 altimeter data without precise satellite orbits

The paper describes results obtained from the processing of 53 Geos-3 arcs of altimeter data obtained during the first weeks after the launch of the satellite in April, 1975. The measurement from the satellite to the ocean surface was used to obtain an approximate geoid undulation which was contaminated by long wavelength errors caused primarily by altimeter bias and orbit error. This long wavelength error was reduced by fitting with a low degree polynomial the raw undulation data to the undulations implied by the GEM 7 potential coefficients, in an adjustment process that included conditions on tracks that cross. The root mean square crossover discrepancy before this adjustment was ±12.4 meters while after the adjustment it was ±0.9 m. These adjusted undulations were used to construct a geoid map in the Geos-3 calibration area using a least squares filter to remove remaining noise in the undulations. Comparing these undulations to ones computed from potential coefficients and terrestrial gravity data indicates a mean difference of 0.25 m and a root mean square difference of ±1.92 m. The adjusted undulations were also used to estimate several 5^o, 2^o, and 1^o anomalies using the method of least squares collocation. The resulting predictions agreed well with known values although the 1^o x 1^o anomalies could not be considered as reliably determined.     

The reduction correction in North America

An inverse Poisson integral technique has been used to determine a gravity field on the geoid which, when continued by analytic free space methods to the topographic surface, agrees with the observed field. The computation is performed in three stages, each stage refining the previous solution using data at progressively increasing resolution (1^o×1^o, 5′×5′, 5/8′×5/8′) from a decreasing area of integration. Reduction corrections are computed at 5/8′×5/8′ granularity by differencing the geoidal and surface values, smoothed by low-pass filtering and sub-sampled at 5′ intervals. This paper discusses 1^o×1^o averages of the reduction corrections thus obtained for 172 1^o×1^o squares in western North America. The 1^o×1^o mean reduction corrections are predominantly positive, varying from −3 to +15mgal, with values in excess of 5mgal for 26 squares. Their mean andrms values are +2.4 and 3.6mgal respectively and they correlate well with the mean terrain corrections as predicted byPellinen in 1962. The mean andrms contributions from the three stages of computation are: 1^o×1^o stage +0.15 and 0.7mgal; 5′×5′ stage +1.0 and 1.6mgal; and 5/8′×5/8′ stage +1.3 and 1.8mgal. These results reflect a tendency for the contributions to become larger and more systematically positive as the wavelengths involved become shorter. The results are discussed in terms of two mechanisms; the first is a tendency for the absolute values of both positive and negative anomalies to become larger when continued downwards and, the second, a non-linear rectification, due to the correlation between gravity anomaly and topographic height, which results in the values continued to a level surface being systematically more positive than those on the topography.     

A synthetic Earth for use in geodesy

 A synthetic Earth and its gravity field that can be represented at different resolutions for testing and comparing existing and new methods used for global gravity-field determination are created. Both the boundary and boundary values of the gravity potential can be generated. The approach chosen also allows observables to be generated at aircraft flight height or at satellite altitude. The generation of the synthetic Earth shape (SES) and gravity-field quantities is based upon spherical harmonic expansions of the isostatically compensated equivalent rock topography and the EGM96 global geopotential model. Spherical harmonic models are developed for both the synthetic Earth topography (SET) and the synthetic Earth potential (SEP) up to degree and order 2160 corresponding to a 5′×5′ resolution. Various sets of SET, SES and SEP with boundary geometry and boundary values at different resolutions can be generated using low-pass filters applied to the expansions. The representation is achieved in point sets based upon refined triangulation of a octahedral geometry projected onto the chosen reference ellipsoid. The filter cut-offs relate to the sampling pattern in order to avoid aliasing effects. Examples of the SET and its gravity field are shown for a resolution with a Nyquist sampling rate of 8.27 degrees. Received: 6 August 1999 / Accepted: 26 April 2000     

A preliminary error analysis of the gravity field recovery from a lunar Satellite-to-Satellite mission

Summary A low cost lunar Satellite-to-Satellite radio tracking mission in a low-low configuration could considerably improve the existing knowledge about the lunar gravity field. The impact of various mission parameters that may contribute to the recovery of the gravity field, such as satellite altitude, satellite separation, mission duration, measurement precision and sampling interval were quantified using the Jekeli-Rapp algorithm. Preliminary results indicate that the gravity field resolution up to harmonic degree 40 to 80 is feasible depending on various mission configurations. Radio tracking data from a six-month mission with a precision of 1 mm s^–1 every 10 s and 300 km satellite separation at 150 km altitude will permit the determination of 5^o×5^o mean gravity anomalies with an error of approximately 15 mgals. Consideration of other unaccounted error sources of instrumental, operational as well as environmental nature may lower this resolution.     

Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies

The availability of high-resolution global digital elevation data sets has raised a growing interest in the feasibility of obtaining their spherical harmonic representation at matching resolution, and from there in the modelling of induced gravity perturbations. We have therefore estimated spherical Bouguer and Airy isostatic anomalies whose spherical harmonic models are derived from the Earth’s topography harmonic expansion. These spherical anomalies differ from the classical planar ones and may be used in the context of new applications. We succeeded in meeting a number of challenges to build spherical harmonic models with no theoretical limitation on the resolution. A specific algorithm was developed to enable the computation of associated Legendre functions to any degree and order. It was successfully tested up to degree 32,400. All analyses and syntheses were performed, in 64 bits arithmetic and with semi-empirical control of the significant terms to prevent from calculus underflows and overflows, according to IEEE limitations, also in preserving the speed of a specific regular grid processing scheme. Finally, the continuation from the reference ellipsoid’s surface to the Earth’s surface was performed by high-order Taylor expansion with all grids of required partial derivatives being computed in parallel. The main application was the production of a 1′ × 1′ equiangular global Bouguer anomaly grid which was computed by spherical harmonic analysis of the Earth’s topography–bathymetry ETOPO1 data set up to degree and order 10,800, taking into account the precise boundaries and densities of major lakes and inner seas, with their own altitude, polar caps with bedrock information, and land areas below sea level. The harmonic coefficients for each entity were derived by analyzing the corresponding ETOPO1 part, and free surface data when required, at one arc minute resolution. The following approximations were made: the land, ocean and ice cap gravity spherical harmonic coefficients were computed up to the third degree of the altitude, and the harmonics of the other, smaller parts up to the second degree. Their sum constitutes what we call ETOPG1, the Earth’s TOPography derived Gravity model at 1′ resolution (half-wavelength). The EGM2008 gravity field model and ETOPG1 were then used to rigorously compute 1′ × 1′ point values of surface gravity anomalies and disturbances, respectively, worldwide, at the real Earth’s surface, i.e. at the lower limit of the atmosphere. The disturbance grid is the most interesting product of this study and can be used in various contexts. The surface gravity anomaly grid is an accurate product associated with EGM2008 and ETOPO1, but its gravity information contents are those of EGM2008. Our method was validated by comparison with a direct numerical integration approach applied to a test area in Morocco–South of Spain (Kuhn, private communication 2011) and the agreement was satisfactory. Finally isostatic corrections according to the Airy model, but in spherical geometry, with harmonic coefficients derived from the sets of the ETOPO1 different parts, were computed with a uniform depth of compensation of 30?km. The new world Bouguer and isostatic gravity maps and grids here produced will be made available through the Commission for the Geological Map of the World. Since gravity values are those of the EGM2008 model, geophysical interpretation from these products should not be done for spatial scales below 5 arc minutes (half-wavelength).     

Evaluation of deflections of the vertical using topographic-isostatic and astrogeodetic data for the area of Greece

The evaluation of deflections of the vertical for the area of Greece is attempted using a combination of topographic and astrogeodetic data. Tests carried out in the area bounded by 35°≤ϕ≤42°, 19°≤λ≤27° indicate that an accuracy of ±3″.3 can be obtained in this area for the meridian and prime vertical deflection components when high resolution topographic data in the immediate vicinity of computation points are used, combined with high degree spherical harmonic expansions of the geopotential and isostatic reduction potential. This accuracy is about 25% better than the corresponding topographic-Moho deflection components which are evaluated using topographic and Moho data up to 120 km around each station, without any combination with the spherical harmonic expansion of the geopotential or isostatic reduction potential. The accuracy in both cases is increased to about 2″.6 when the astrogeodetic data available in the area mentioned above are used for the prediction of remaining values. Furthermore the estimation of datum-shift parameters is attempted using least squares collocation.     

Lunar gravity field determination using SELENE same-beam differential VLBI tracking data

A lunar gravity field model up to degree and order 100 in spherical harmonics, named SGM100i, has been determined from SELENE and historical tracking data, with an emphasis on using same-beam S-band differential VLBI data obtained in the SELENE mission between January 2008 and February 2009. Orbit consistency throughout the entire mission period of SELENE as determined from orbit overlaps for the two sub-satellites of SELENE involved in the VLBI tracking improved consistently from several hundreds of metres to several tens of metres by including differential VLBI data. Through orbits that are better determined, the gravity field model is also improved by including these data. Orbit determination performance for the new model shows improvements over earlier 100th degree and order models, especially for edge-on orbits over the deep far side. Lunar Prospector orbit determination shows an improvement of orbit consistency from 1-day predictions for 2-day arcs of 6 m in a total sense, with most improvement in the along and cross-track directions. Data fit for the types and satellites involved is also improved. Formal errors for the lower degrees are smaller, and the new model also shows increased correlations with topography over the far side. The estimated value for the lunar GM for this model equals 4902.80080±0.0009 km³/s² (10 sigma). The lunar degree 2 potential Love number k ₂ was also estimated, and has a value of 0.0255 ± 0.0016 (10 sigma as well).     

Geoid determination using adapted reference field, seismic Moho depths and variable density contrast

 The traditional remove-restore technique for geoid computation suffers from two main drawbacks. The first is the assumption of an isostatic hypothesis to compute the compensation masses. The second is the double consideration of the effect of the topographic–isostatic masses within the data window through removing the reference field and the terrain reduction process. To overcome the first disadvantage, the seismic Moho depths, representing, more or less, the actual compensating masses, have been used with variable density anomalies computed by employing the topographic–isostatic mass balance principle. In order to avoid the double consideration of the effect of the topographic–isostatic masses within the data window, the effect of these masses for the used fixed data window, in terms of potential coefficients, has been subtracted from the reference field, yielding an adapted reference field. This adapted reference field has been used for the remove–restore technique. The necessary harmonic analysis of the topographic–isostatic potential using seismic Moho depths with variable density anomalies is given. A wide comparison among geoids computed by the adapted reference field with both the Airy–Heiskanen isostatic model and seismic Moho depths with variable density anomaly and a geoid computed by the traditional remove–restore technique is made. The results show that using seismic Moho depths with variable density anomaly along with the adapted reference field gives the best relative geoid accuracy compared to the GPS/levelling geoid. Received: 3 October 2001 / Accepted: 20 September 2002 Correspondence to: H.A. Abd-Elmotaal     

Downward continuation of the free-air gravity anomalies to the ellipsoid using the gradient solution and terrain correction—An attempt of global numerical computations

The formulas for the determination of the coefficients of the spherical harmonic expansion of the disturbing potential of the earth are defined for data given on a sphere. In order to determine the spherical harmonic coefficients, the gravity anomalies have to be analytically downward continued from the earth's surface to a sphere—at least to the ellipsoid. The goal of this paper is to continue the gravity anomalies from the earth's surface downward to the ellipsoid using recent elevation models. The basic method for the downward continuation is the gradient solution (theg ₁ term). The terrain correction has also been computed because of the role it can play as a correction term when calculating harmonic coefficients from surface gravity data. Theg ₁ term and the terrain correction were expanded into the spherical harmonics up to180 ^th order. The corrections (theg ₁ term and the terrain correction) have the order of about 2% of theRMS value of degree variance of the disturbing potential per degree. The influences of theg ₁ term and the terrain correction on the geoid take the order of 1 meter (RMS value of corrections of the geoid undulation) and on the deflections of the vertical is of the order 0.1″ (RMS value of correction of the deflections of the vertical).     

CRUST 2.0模型在均衡位模型构建中的应用研究

基于均衡理论构制地球重力场模型中一个关键问题是扰动质量的确定,不同的地壳密度分布将得到不同的扰动质量。本文立足于面凝聚模型和Airy模型两种均衡补偿机制,研究了CRUST 2.0全球地壳模型在构制高分辨率地球重力场模型中的应用,推导了顾及地球物理信息的两种均衡重力场模型构建公式,分别讨论了CRUST 2.0模型和补偿深度在两种补偿机制中构建重力场模型的贡献。数值分析表明,地壳模型中Moho面深度的算术平均值22.97 km不是最佳补偿深度,而40 km相对最优,补偿深度的对模型超高阶部分影响较小;CRUST 2.0模型能够在361-1080频段内较好地改善原模型;相同补偿深度的面凝聚模型和Airy模型在不同频段的优劣性不一致。     

Improvement of the geoid in local areas by satellite-to-satellite tracking

Summary The concept of satellite-to-satellite tracking measuring the relative velocity of two orbiting satellites spaced some hundreds kilometers on a close orbit, provides now possibilities for the investigation of the Earth’s gravity field. In the paper only medium and short wave length effects affecting the measured relative velocity have been considered. Collocation is used in such an analysis of local geoid improvement, because this method allows to combine heterogeneous data in a consistent way. Covariance functions relevant for the particular case of a circular equatorial orbit are given. Two kinds of observation equations have been formulated. The choice of observation equation with regard to satellites configuration is discussed. It is found that it is sufficient to have a limited number of satellite-to-satellite observations in a 7^o×7^o area around the estimation point with distances between profiles of about 1^o.5 and between the two satellites forming the pair of 200+350 km; the altitude of satellite-to-satellite observations should be as low as possible. The accuracy of the geoid determination strongly depends on the degree and order of the reference field used. An accuracy of about ±1 m can be achieved with an assumed reference field of (40,40). The influence of measuring errors is discussed and it is shown that only satellite-to-satellite observations with accuracy better then 0.1 mm/sec will give an improvement of the geoid. Finally, some results on the combination of low-low satellite-to-satellite tracking and terrestrial gravity data are given. The proposed method seems to be especially interesting for unsurveyed areas. Furthermore, it has the practical advantage that only a local coverage data is needed.     

Glacier and glacial lake classification for change detection studies using satellite data: a case study from Baspa basin,western Himalaya

Using high-resolution Google Earth^TM images in conjunction with Landsat images, the glaciers and lakes in the Baspa basin are classified to explore the recent changes. A total number of 109 glaciers (187 ± 3.7 km²) are mapped and subsequently classified as compound valley glaciers, simple valley glaciers, cirques, niches, glacieretes and ice aprons. The compound and simple valley glaciers contribute 67.1 ± 1.3% and 19.8 ± 0.3% to the total glacier cover of the basin. Similarly, a total number of 129 glacial lakes (0.360 ± 0.007 km²) are identified. From 1976 to 2011, the compound valley glaciers have lost a small area of 10.3 ± 0.03% at a rate of 0.41 ± 0.002 km² a^-1, whereas the niche glaciers have lost higher area of 40.1 ± 0.001% at a rate of 0.04 ± 0.0001 km² a^-1. Change detection of two benchmark glacial lakes revealed a progressive expansion during recent decades. The Baspa Bamak proglacial lake has expanded from 0.020 ± 0.0004 km² (2000) to 0.069 ± 0.001 km² (2011). Due to the complete loss of source ice, another glacial lake has expanded from 0.09 ± 0.001 km² (1994) to 0.10 ± 0.002 km² (2011). During the study period, the mean annual temperature that is T_avg, T_min and T_max have increased significantly at the 95% confidence level by 1.5 ^oC (0.070 °C a^-1), 1.8 ^oC (0.076 °C a^-1) and 1.6 ^oC (0.0071 °C a^-1) from 1985 to 2008. However, the precipitation has decreased significantly from 1976 and 1985 to 2008.     

High degree spherical harmonic expansion of gravity data

A spherical harmonic expansion of the gravity field up to degree and order 200 was carried out. Free air anomaly data over Canada (1⁰×1⁰ block averages) with a range of 211.1 mgal were used for testing. A low degree expansion (N=30) produced a map with a range of 63.6 mgal with contour patterns that could hardly be correlated with the original hand contoured map. A high degree expansion (N=200) on the other hand resulted in a map with a range of 199.8 mgal which quite faithfully reproduced the original including its local variations. Test computations verify that by monitoring the RMS values and the range of the expansion it is possible to arrive at an optimum degree of expansion for a given data set. It was also verified by the computations, that, since the computed expansions essentially have a zero value outside the domain of the input, it is possible to combine the results of separate non-overlapping expansions. Contribution of the Earth Physics Branch No. 900. Presented at the 1977 Spring Meeting, AGU, May 30–June 3, Washington, D.C.     

Simulation of regional gravity field recovery from satellite gravity gradiometer data using collocation and FFT

When planning a satellite gravity gradiometer (SGG) mission, it is important to know the quality of the quantities to be recovered at ground level as a function of e.g. satellite altitude, data type and sampling rate, and signal variance and noise. This kind of knowledge may be provided either using the formal error estimates of wanted quantities using least-squares collocation (LSC) or by comparing simulated data at ground level with results computed by methods like LSC or Fast Fourier Transform (FFT). Results of a regional gravity field recovery in a 10^o×20^o area surrounding the Alps using LSC and FFT are reported. Data used as observations in satellite altitude (202 or161 km) and for comparison at ground level were generated using theOSU86F coefficient set, complete to degree 360. These observations are referred to points across simulated orbits. The simulated quantities were computed for a 45 days mission period and 4 s sampling. A covariance function which also included terms above degree 360 was used for prediction and error estimation. This had the effect that the formal error standard deviation for gravity anomalies were considerably larger than the standard deviations of predicted minus simulated quantities. This shows the importance of using data with frequency content above degree 360 in simulation studies. Using data at202 km altitude the standard deviation of the predicted minus simulated data was equal to8.3 mgal for gravity and0.33 m for geoid heights.