Abstract: The forward modeling of a 3D mesh model is the basis of gravity data inversion. High precision and high efficiency forward modeling is helpful to the improvement of the quality of inversion interpretation. In order to solve the problem of high precision and high efficiency gravity forward modeling based on a large-scale surface observation area, this paper presents the gravity anomaly forward modeling method and parallel algorithm of a 3D Tesseroid mesh model in the spherical coordinate system. The forward modeling uses the improved Gauss-Legendre Quadrature integration method to realize the high-precision gravity anomaly calculation based on a large-scale surface observation area, and also uses the MATLAB task parallel algorithm based on OpenMP to realize the high-efficiency forward modeling. The test on the 3D theoretical model and the Eastern China lithospheric model has verified the validity of the proposed method. This method can provide technical support for efficient large-scale gravity field simulation and 3D inversion.