Within the normal mode relaxation theory, we thoroughly analyze and compare the exact compressible and incompressible solution for a viscoelastic stratified Earth model with an approximated analytical one, where the ratio between gravity and radial distance from the Earth’s centre is considered constant in each layer of the model. We implement an algorithm, based on the Runge-Kutta scheme, to integrate in the r variable the spheroidal part of the linearized momentum and Poisson equations. This numerical scheme allows us to unravel the impact of such an approximation. We disclose new aspects of the physics underneath the terms entering the system of differential equations, such as the advection and the buoyancy force. These issues are relevant for a wide range of geophysical applications and timescales, from the 1-102 yr related to postseismic deformation, to the 103 yr of postglacial rebound, to the 106 yr of True Polar Wander. We show that such an approximation affects the buoyancy force term, due to the sensitivity of the tangential displacement component in the compressible deformation, and the advection term, for that part containing an effective unstable radial force, dependent on the radial displacement. Relative errors between exact and approximated Love numbers are sizable for low harmonic degrees (13, 62, 17 per cent and 4, 39, 4 per cent for the elastic and viscous radial, tangential and gravitational Love numbers), but can become lower than 1 per cent at high harmonic degrees, for appropriate choices of the constant gravity term entering advection and buoyancy force. Our findings shed light on the role of the latter in the readjustment of the interfaces where the planet is subject to variations in its physical properties, in quasi-static deformation case. We show that the contribution of the denumerably infinite compressional D-modes can be dealt with accurately by normal mode approach and that the D-mode cluster point does not contribute to the deformation.