AbstractThe geodetic observations of static deformations, including gravity perturbations and displacement fields due to huge earthquakes, are understood and explained using recent dislocation theories. Due to multiple possible mechanisms for the post-seismic phase of earthquakes, the dominant mechanism may change at different spatiotemporal ranges for different earthquake types. Accurate forward and inverse modeling of post-seismic deformations is valuable and needed information for geoscience communities. The existing methods for calculating gravitational viscoelastic relaxation can be improved or simplified to make them more suitable for more realistic Earth models and/or to overcome the poor convergence performance and/or overflow risks during numerical calculations. In this study, a simple and effective method for calculating the post-seismic relaxation deformations is proposed. This method is different from previous methods, such as the normal mode summation and rectangle integration methods. The proposed method consists of a rational functional approximation of the integral kernel and a transformation of the numerical inverse Laplace transform (NILT) into an alternating series summation using the residual theorem. Then the intrinsic oscillation and overflow risks are thoroughly suppressed. The accuracy of the calculated Green’s functions can be easily controlled by choosing a suitable parameter. In addition, the proposed method also has applicability in different Earth models with linear rheological profiles.