Surface gravity inversion attempts to recover the density contrast distribution in the 3D Earth model for geological interpretation. Since airborne gravity is characterized by large data volumes, large-scale 3D inversion exceeds the capacity of desktop computing resources, making it difficult to achieve the appropriate depth/lateral resolution for geological interpretation. In addition, gravity data are finite and noisy, and their inversion is ill posed. Especially in the absence of a priori geological information, regularization must be introduced to overcome the difficulty of the non-uniqueness of the solutions to recover the most geologically plausible ones. Because the use of Haar wavelet operators has an edge-preserving property and can preserve the sensitivity matrix structure at each level of the multilevel method to obtain faster solvers, we present a multilevel algorithm for large-scale gravity inversion solved by the re-weighted regularized conjugate gradient (RRCG) algorithm to reduce the inversion computational resources and improve the depth/lateral resolution of the inversion results. The RRCG-based multilevel inversion was then applied to synthetic cases and airborne gravity data from the Quest-South project in British Columbia, Canada. Results from synthetic models and field data show that the RRCG-based multilevel inversion is suitable for obtaining density contrast distributions with appropriate horizontal and vertical resolution, especially for large-scale gravity inversions compared to Occam’s inversion.