We have developed a theory based on the direct numerical integration in time for studying the temporal viscoelastic responses of earth models to surface loads. Modelling in the time domain is motivated by the fact that realistic elastically compressible models generate an infinite number of modes and the width of such "continuous spectrum" may cover several orders of magnitude in the Laplacian spectral domain for complicated viscosity stratification, which causes numerical difficulties for the normal-mode method. From our numerical solutions, we have directed our attention on the influences of elastic compressibility, thickness of the lithosphere, the nature of the internal mantle boundaries with density jumps and the viscosity structure near the interface between the lower and the upper mantle. There is a great difference between the responses of compressible and incompressible models mainly for shorter wavelengths and thus incompressible models seem to be inadequate for short-wavelength responses. The sensitivity of the viscoelastic responses to the lithospheric thickness in the presence of a low viscosity asthenosphere is substantial. This points to the need for constructing models with a 3-D viscosity variations, which would yield more realistic lithospheric-asthenospheric structure globally than models with a constant lithospheric thickness. The sensitivity of the Earth's viscoelastic behaviour to the other parameters is weaker but still is noteworthy.