Application of the Post-Widder Laplace inversion algorithm to postseismic rebound models

The computation of global postseismic rebound in a spherically symmetric, stratified, self-gravitating Earth with Maxwell viscoelastic rheology can be carried out semi-analytically with a normal-mode approach. The solution scheme usually involves the application of standard propagator techniques to the equivalent problem in the Laplace domain; to recover the temporal dependence a numerical Laplace anti-transform is required. This step involves the solution of the so-called "secular equation", whose degree increases linearly with the detail of the stratification modeling, and whose coefficients become extremely ill-conditioned for high harmonic orders. As a result, the practically solvable models are limited to a few viscoelastic layers, and are anyway affected by severe numerical instabilities. To overcome these difficulties, alternative approaches have been explored by several authors, ranging from Runge-Kutta purely numerical integration to the evaluation of Laplace antitransform by a numerical discretization of the Bromwich integral. The Post-Widder algorithm allows the estimation of the Laplace antitransform by sampling numerically the transform on the positive real axis. This method, which has been recently applied to the computation of GIA viscoelastic Love numbers, allows to bypass completely the root-finding procedure while preserving at the same time the analytical normal-mode solution form. In this work, we apply the Post-Widder method to the computation of post-seismic rebound models. We perform a series of benchmarks to optimize the algorithm for speed while checking its stability against earlier results.