<p>Gravity anomalies resulting from a continuous density media are described by a 2D Poisson equation. Current modeling techniques primarily rely on integral solutions to address these anomalies. However, this approach yields 2D gravity anomalies, including vector and tensor of gravitational fields, which inevitably contain integral singularities. To overcome this issue, we present a spectral-element scheme based on Gauss-Lobatto-Chebyshev polynomials to compute the boundary value problem associated with the 2D gravitational potential. By formulating the variational problem related to the 2D gravitational potential, we derive the corresponding spectral-element equation. Then, the resulting linear system of equations from the spectral-element algorithm are solved using a BiCGStab-Jacobi iterative solver. Synthetic 2D models featuring both constant and variable densities are developed to assess the computational accuracy of the proposed algorithm. Additionally, we examine a real-world sedimentary basin model to demonstrate the practical application of the Chebyshev spectral-element forward algorithm. All numerical results indicate that our algorithm can produce gravity responses efficiently and accurately.</p>