A data-driven approach to local gravity field modelling using spherical radial basis functions

We propose a methodology for local gravity field modelling from gravity data using spherical radial basis functions. The methodology comprises two steps: in step 1, gravity data (gravity anomalies and/or gravity disturbances) are used to estimate the disturbing potential using least-squares techniques. The latter is represented as a linear combination of spherical radial basis functions (SRBFs). A data-adaptive strategy is used to select the optimal number, location, and depths of the SRBFs using generalized cross validation. Variance component estimation is used to determine the optimal regularization parameter and to properly weight the different data sets. In the second step, the gravimetric height anomalies are combined with observed differences between global positioning system (GPS) ellipsoidal heights and normal heights. The data combination is written as the solution of a Cauchy boundary-value problem for the Laplace equation. This allows removal of the non-uniqueness of the problem of local gravity field modelling from terrestrial gravity data. At the same time, existing systematic distortions in the gravimetric and geometric height anomalies are also absorbed into the combination. The approach is used to compute a height reference surface for the Netherlands. The solution is compared with NLGEO2004, the official Dutch height reference surface, which has been computed using the same data but a Stokes-based approach with kernel modification and a geometric six-parameter “corrector surface” to fit the gravimetric solution to the GPS-levelling points. A direct comparison of both height reference surfaces shows an RMS difference of 0.6 cm; the maximum difference is 2.1 cm. A test at independent GPS-levelling control points, confirms that our solution is in no way inferior to NLGEO2004.     

Gravity field determination using boundary element methods

The Boundary Element Method (BEM), a numerical technique for solving boundary integral equations, is introduced to determine the earth's gravity field. After a short survey on its main principles, we apply this method to the fixed gravimetric boundary value problem (BVP), i.e. the determination of the earth's gravitational potential from measurements of the intensity of the gravity field in points on the earth's surface. We show how to linearize this nonlinear BVP using an implicit function theorem and how to transform the linearized BVP into a boundary integral equation using the single layer representation. A Galerkin method is used to transform the boundary integral equation using the single layer representation. A Galerkin method is used to transform the boundary integral equation into a linear system of equations. We discuss the major problems of this approach for setting up and solving the linear system. The BVP is numerically solved for a bounded part of the earth's surface using a high resolution reference gravity model, measured gravity values of high density, and a 50 50 m² digital terrain model to describe the earth's surface. We obtain a gravity field resolution of 1 1 km² with an accuracy of the order 10^–3 to 10^–4 in about 1 CPU-hour on a Siemens/Fujitsu SIMD vector pipeline machine using highly sophisticated numerical integration techniques and fast equation solvers. We conclude that BEM is a powerful numerical tool for solving boundary value problems and may be an alternative to classical geodetic techniques.     

Statistically optimal estimation of Greenland Ice Sheet mass variations from GRACE monthly solutions using an improved mascon approach

We present an improved mascon approach to transform monthly spherical harmonic solutions based on GRACE satellite data into mass anomaly estimates in Greenland. The GRACE-based spherical harmonic coefficients are used to synthesize gravity anomalies at satellite altitude, which are then inverted into mass anomalies per mascon. The limited spectral content of the gravity anomalies is properly accounted for by applying a low-pass filter as part of the inversion procedure to make the functional model spectrally consistent with the data. The full error covariance matrices of the monthly GRACE solutions are properly propagated using the law of covariance propagation. Using numerical experiments, we demonstrate the importance of a proper data weighting and of the spectral consistency between functional model and data. The developed methodology is applied to process real GRACE level-2 data (CSR RL05). The obtained mass anomaly estimates are integrated over five drainage systems, as well as over entire Greenland. We find that the statistically optimal data weighting reduces random noise by 35–69%, depending on the drainage system. The obtained mass anomaly time-series are de-trended to eliminate the contribution of ice discharge and are compared with de-trended surface mass balance (SMB) time-series computed with the Regional Atmospheric Climate Model (RACMO 2.3). We show that when using a statistically optimal data weighting in GRACE data processing, the discrepancies between GRACE-based estimates of SMB and modelled SMB are reduced by 24–47%.     

Gravimetric monitoring of the first field‐wide steam injection in a fractured carbonate field in Oman – a feasibility study

Gas‐Oil Gravity Drainage is to be enhanced by steam injection in a highly fractured, low permeability carbonate field in Oman. Following a successful pilot, field‐wide steam injection is being implemented, first of this type in the world. A dedicated monitoring program has been designed to track changes in the reservoir. Various observations are to be acquired, including, surface deformation, temperature measurements, microseismic, well logs, pressure and saturation measurements to monitor the reservoir. Because significant changes in the reservoir density are expected, time‐lapse gravimetry is also being considered. In this paper we investigate the feasibility of gravimetric monitoring of the thermally enhanced gravity drainage process at the carbonate field in Oman. For this purpose, forward gravity modelling is performed. Based on field groundwater measurements, the estimates of the hydrological signal are considered and it is investigated under what conditions the groundwater influences can be minimized. Using regularized inversion of synthetic gravity data, we analyse the achievable accuracy of heat‐front position estimates. In case of large groundwater variations at the field, the gravity observations can be significantly affected and, consequently, the accuracy of heat‐front monitoring can be deteriorated. We show that, by applying gravity corrections based on local observations of groundwater, the hydrological influences can to a large extent be reduced and the accuracy of estimates can be improved. We conclude that gravimetric monitoring of the heat‐front evolution has a great potential.     

A Comparison of Global and Regional GRACE Models for Land Hydrology

When using GRACE as a tool for hydrology, many different gravity field model products are now available to the end user. The traditional spherical harmonics solutions produced from GRACE are typically obtained through an optimization of the gravity field data at the global scale, and are generated by a number of processing centers around the world. Alternatives to this global approach include so-called regional techniques, for which many variants exist, but whose common trait is that they only use the gravity data collected over the area of interest to generate the solution. To determine whether these regional solutions hold any advantage over the global techniques in terms of overall accuracy, a range of comparisons were made using some of the more widely used regional and global methods currently available. The regional techniques tested made use of either spherical radial basis functions or single layer densities (i.e., mascons), with the global solutions having been obtained from the various major processing centers. The solutions were evaluated using a range of computed statistics over a selection of major river basins, which were globally distributed and ranged in size from 1 to 6 million km². For one of the basins tested, the Zambezi, additional validation tests were conducted through comparisons against a custom designed regional hydrology model of the region. We could not prove that current regional models perform better than global ones. Monthly mean water storage variations agree at the level of 0.02 m equivalent water height. The differences in terms of monthly mean water storage variations between regional and global solutions are comparable with the differences among only global or regional solutions. Typically they reach values of 0.02 m equivalent water heights, which seems to be the level of accuracy of current GRACE solutions for river basins above 1 million km². The amplitudes of the seasonal mass variations agree at the sub-centimetre level. Evident from all of the comparisons shown is the importance that the choice of regularization, or spatial filtering, can have on the solution quality. This was found to be true for global as well as regional techniques.     

Fast numerical solution of the linearized Molodensky problem

 When standard boundary element methods (BEM) are used in order to solve the linearized vector Molodensky problem we are confronted with two problems: (1) the absence of O(|x|⁻²) terms in the decay condition is not taken into account, since the single-layer ansatz, which is commonly used as representation of the disturbing potential, is of the order O(|x|⁻¹) as x→∞. This implies that the standard theory of Galerkin BEM is not applicable since the injectivity of the integral operator fails; (2) the N×N stiffness matrix is dense, with N typically of the order 10⁵. Without fast algorithms, which provide suitable approximations to the stiffness matrix by a sparse one with O(N(logN) ^s ), s≥0, non-zero elements, high-resolution global gravity field recovery is not feasible. Solutions to both problems are proposed. (1) A proper variational formulation taking the decay condition into account is based on some closed subspace of co-dimension 3 of the space of square integrable functions on the boundary surface. Instead of imposing the constraints directly on the boundary element trial space, they are incorporated into a variational formulation by penalization with a Lagrange multiplier. The conforming discretization yields an augmented linear system of equations of dimension N+3×N+3. The penalty term guarantees the well-posedness of the problem, and gives precise information about the incompatibility of the data. (2) Since the upper left submatrix of dimension N×N of the augmented system is the stiffness matrix of the standard BEM, the approach allows all techniques to be used to generate sparse approximations to the stiffness matrix, such as wavelets, fast multipole methods, panel clustering etc., without any modification. A combination of panel clustering and fast multipole method is used in order to solve the augmented linear system of equations in O(N) operations. The method is based on an approximation of the kernel function of the integral operator by a degenerate kernel in the far field, which is provided by a multipole expansion of the kernel function. Numerical experiments show that the fast algorithm is superior to the standard BEM algorithm in terms of CPU time by about three orders of magnitude for N=65 538 unknowns. Similar holds for the storage requirements. About 30 iterations are necessary in order to solve the linear system of equations using the generalized minimum residual method (GMRES). The number of iterations is almost independent of the number of unknowns, which indicates good conditioning of the system matrix. Received: 16 October 1999 / Accepted: 28 February 2001     

Exploring gravity field determination from orbit perturbations of the European Gravity Mission GOCE

 A comparison was made between two methods for gravity field recovery from orbit perturbations that can be derived from global positioning system satellite-to-satellite tracking observations of the future European gravity field mission GOCE (Gravity Field and Steady-State Ocean Circulation Explorer). The first method is based on the analytical linear orbit perturbation theory that leads under certain conditions to a block-diagonal normal matrix for the gravity unknowns, significantly reducing the required computation time. The second method makes use of numerical integration to derive the observation equations, leading to a full set of normal equations requiring powerful computer facilities. Simulations were carried out for gravity field recovery experiments up to spherical harmonic degree and order 80 from 10 days of observation. It was found that the first method leads to large approximation errors as soon as the maximum degree surpasses the first resonance orders and great care has to be taken with modeling resonance orbit perturbations, thereby loosing the block-diagonal structure. The second method proved to be successful, provided a proper division of the data period into orbital arcs that are not too long. Received: 28 April 2000 / Accepted: 6 November 2000