The 3D gravity inverse problem is ill-posed, and regularization technique must be introduced to overcome the non-uniqueness and ill-posedness problem of gravity inversion. The L-2-norm based model regularization has been successfully applied to 3D inversion of field gravity data. However, when sharp boundaries exist between anomalous bodies and their surrounding rocks, the 3D gravity inversion with L-2-norm model regularization recovers density models with smoothed boundaries, which bring some difficulty in subsequent geological interpretation. We have developed a 3D gravity inversion method based on reformulated L-p-norm (0 <= p <= 2) model regularization. In this approach, the L-p-norm regularization term in the objective function is reformulated in the form of L-2-norm by introducing a diagonal matrix, and the value of p can vary within a greater range (i.e., 0 <= p <= 2) rather than 0 <= p < 1 or 1 <= p <= 2. The Gauss Newton minimization scheme is adopted to minimize the reformulated objective function in the weighted model parameters and data space. The method implements a 3D gravity sparse inversion when 0 <= p < 2 and, in particular, becomes the conventional 3D gravity smooth inversion when p = 2. The test of synthetic gravity data shows that, as the value of p decreases from 2 to 0, this method could produce density models with clearer boundaries and density value closer to the true value, and obtain the most satisfactory inversion results when p = 0. In addition, we also investigate the effect of the fine discretization of the mesh on inversion result. At last, this approach is applied to the real gravity data from the San Nicolas deposit in Mexico, and the reconstructed sulfide ore body is in good agreement with the known information. (C) 2021 Elsevier B.V. All rights reserved.