A novel non-singular and numerically stable algorithm for efficient tesseroid gravity forward modeling

Gravity forward modeling, based on Newton’s law of gravitation, describes the relationship between subsurface mass density distribution and observed gravity data. It has broad applications in global geodynamics, resource exploration, and planetary sciences. The tesseroid (spherical prism) is widely used to represent mass density elements on a spherical Earth, but its application is hindered by singularities, numerical instabilities near computation points, and the trade-off between accuracy and efficiency in large-scale modeling. To address these theoretical deficiencies, we propose a novel tesseroid gravity forward modeling algorithm formulated in the spherical harmonic domain. The spherical harmonic coefficients are separated into three parts: (1) an analytical, non-singular expression in the radial direction; (2) an analytical, non-singular expression in the longitudinal direction; and (3) a latitudinal integral, efficiently evaluated without singularities using the Gauss–Legendre quadrature (GLQ). Fast Fourier transform (FFT) technique is combined with spherical harmonic synthesis to further enhance computational efficiency. Validation with spherical shell tests shows that accuracy depends jointly on the spherical harmonic expansion degree and the GLQ order, meaning that higher expansion degrees require more GLQ nodes, while greater computation heights reduce the required order. For a single tesseroid, the proposed method matches the adaptive discretization method in both spatial distribution and magnitude of gravity disturbances, with mean differences below 1 mGal, while reducing computation time by ~ 35%. Application to the Tibet topographic model confirms differences below 2 mGal, with runtime reduced to 45% of the adaptive discretization method. At the global scale, differences remain below 1 mGal, and computation time is reduced to only 0.05% of the adaptive discretization method. These results demonstrate that the proposed algorithm is accurate, robust, and highly efficient, making it well-suited for large-scale gravity forward modeling on a spherical Earth.