SUMMARY Constitutive theory for viscoelasticity has broad application to solid mantle or ice deformations driven by tides, surface mass variations, and post-seismic flow. Geophysical models using higher order viscoelasticity can better accommodate geodetic observations than lower-order theory, typically provided by tensor versions of Maxwell, 4-parameter Burgers or standard linear (Zener) rheology. We derive, for the first time, a mathematical description of a compressible version of the extended Burgers material (EBM) model paradigm which has a distribution function of relaxation spectra. The latter model is often used for parametrizing high temperature background transient responses in the rock physics and mechanics laboratory setting and have demonstrated application to low frequency seismic wave attenuation. A new generalization of this practical anelastic model is presented and applied to the glacial isostatic adjustment momentum equations, thus providing useful guidance for generating initial-value boundary problem-solving software for quite general coding strategies. The solutions for the vertical motion response to a suddenly imposed surface load reveal a short-term transience of substantial amplitude.