Fast calculation of gravitational effects using tesseroids with a polynomial density of arbitrary degree in depth

Abstract

Fast and accurate calculation of gravitational effects on a regional or global scale with complex density environment is a critical issue in gravitational forward modelling. Most existing significant developments with tessroid-based modelling are limited to homogeneous density models or polynomial ones of a limited order. Moreover, the total gravitational effects of tesseroids are often calculated by pure summation in these methods, which makes the calculation extremely time-consuming. A new efficient and accurate method based on tesseroids with a polynomial density up to an arbitrary order in depth is developed for 3D large-scale gravitational forward modelling. The method divides the source region into a number of tesseroids, and the density in each tesseroid is assumed to be a polynomial function of arbitrary degree. To guarantee the computational accuracy and efficiency, two key points are involved: (1) the volume Newton’s integral is decomposed into a one-dimensional integral with a polynomial density in the radial direction, for which a simple analytical recursive formula is derived for efficient calculation, and a surface integral over the horizontal directions evaluated by the Gauss–Legendre quadrature (GLQ) combined with a 2D adaptive discretization strategy; (2) a fast and flexible discrete convolution algorithm based on 1D fast Fourier transform (FFT) and a general Toepritz form of weight coefficient matrices is adopted in the longitudinal dimension to speed up the computation of the cumulative contributions from all tesseroids. Numerical examples show that the gravitational fields predicted by the new method have a good agreement with the corresponding analytical solutions for spherical shell models with both polynomial and non-polynomial density variations in depth. Compared with the 3D GLQ methods, the new algorithm is computationally more accurate and efficient. The calculation time is significantly reduced by 3 orders of magnitude as compared with the traditional 3D GLQ methods. Application of the new algorithm in the global crustal CRUST1.0 model further verifies its reliability and practicability in real cases. The proposed method will provide a powerful numerical tool for large-scale gravity modelling and also an efficient forward engine for inversion and continuation problems.