<p>The stochastic modelling of a finite mass distribution can provide a new perspective on the dynamic evaluation of time variable gravity fields. The algorithm for estimating variations of spherical harmonic coefficients implied by corresponding shape changes is implemented for the first-order derivatives of the gravitational potential. The described algorithm uses the spherical harmonic synthesis formula expressed in Cartesian coordinates that includes the derived Legendre functions (DLFs). Here, we expand the estimation process by implementing also the traditional spherical harmonic synthesis formula of normalized associated Legendre functions (ALFs) expressed in spherical coordinates. The variations obtained by applying the two approaches are compared with gravity signal differences induced by the modelled shape changes using the line integral analytical approach. The numerical comparisons refer to three asteroid shape models of Eros, Didymos and Dimorphos. The first-order derivative values provided by the DLF expressions and their variations using ALF are closer to the analytical method’s results. The highest calculated differences refer to Δ<i>V</i>_<i>z</i> with their mean value reaching 37% with respect to the other components obtained by all methods. Finally, the respective harmonic series converge to a fixed numerical value at a maximum expansion degree equal to 15 near Brillouin sphere and 5 as the distance of the computation point increases.</p>