<p>The stochastic representation of an uncertain shape model allows the dynamic evaluation of its induced gravity signal. This can be also applied for representing a time variable gravity field to model mass changes. The algorithm for estimating variations in gravitational potential is extended for the case of second-order derivatives. Two different harmonic synthesis formulas are used to derive the sought variations: one expressed in spherical coordinates using the traditional associated Legendre functions (ALF) and their derivatives up to the second order, while the other expressed in Cartesian coordinates by including the derived Legendre functions (DLF). The obtained variations are compared in terms of convergence with gravity signal differences referring to the specific shape changes using the line integral analytical approach for three asteroid shape models. Both approaches provide results that differ from the analytical method at a 1E−1 level, while the differences between them are at the 1E−15 level. The obtained results are highly influenced by the geometry of the examined shape model, with the ALF approach providing variations with closer agreement with the analytical method only for the almost spherical shape. Both harmonic synthesis expressions can be used to derive accurate results, as they differ at a very low level, and one can choose based on the convenience of their algorithmic characteristics.</p>