In recent years, there are research trends from constant to variable density and low-order to high-order gravitational potential gradients in gravity field modeling. Under the research circumstances, this paper focuses on the variable density model for gravitational curvatures (or gravity curvatures, third-order derivatives of gravitational potential) of a tesseroid and spherical shell in the spatial domain of gravity field modeling. In this contribution, the general formula of the gravitational curvatures of a tesseroid with arbitrary-order polynomial density is derived. The general expressions for gravitational effects up to the gravitational curvatures of a spherical shell with arbitrary-order polynomial density are derived when the computation point is located above, inside, and below the spherical shell. When the computation point is located above the spherical shell, the general expressions for the mass of a spherical shell and the relation between the radial gravitational effects up to arbitrary-order and the mass of a spherical shell with arbitrary-order polynomial density are derived. The influence of the computation point’s height and latitude on gravitational curvatures with the polynomial density up to fourth-order is numerically investigated using tesseroids to discretize a spherical shell. Numerical results reveal that the near-zone problem exists for the fourth-order polynomial density of the gravitational curvatures, i.e., relative errors in $$\log _{10}$$
log
10
scale of gravitational curvatures are large than 0 below the height of about 50 km by a grid size of $$15'\times 15'$$
15
′
×
15
′
. The polar-singularity problem does not occur for the gravitational curvatures with polynomial density up to fourth-order because of the Cartesian integral kernels of the tesseroid. The density variation can be revealed in the absolute errors as the superposition effects of Laplace parameters of gravitational curvatures other than the relative errors. The derived expressions are examples of the high-order gravitational potential gradients of the mass body with variable density in the spatial domain, which will provide the theoretical basis for future applications of gravity field modeling in geodesy and geophysics.