The Problem of Viscoelastic Relaxation of the Earth Solved by a Matrix Eigenvalue Approach Based on Discretization in Grid Space

Viscoelasticity is a geophysical process which operates over intermediate timescales between a few years to millions of years, depending on the ambient thermal conditions of the self-gravitating spherical planet. Topography undulations with time represent geological signatures to both the internal and external loading processes, such as post-glacial rebound, volcanic eruptions, sedimentary loading, meteoritic impacts and mountain building. The span of relaxation timescales or relaxation spectrum in a viscoelastic spherical body has traditionally been determined by employing the Laplacian transform method and the correspondence principle relating the elastic solution to its viscoelastic counterpart.We have devised a novel approach based on the method of lines in which the equations at each angular order in the spherical harmonic expansion are discretized in the radial co-ordinate. The finite-dimensional space spanned by the discretized points in the radial direction of the planetary model then forms the basis of a matrix Eigenvalue problem. The Eigenvalues can be computed very fast because of the availability of public domain software. We can, for instance, compute the entire range of viscoelastic relaxation with computational times from 0.05 to 50 sec using only 30 to 300 radial grid points. The models can have both realistic density and elastic parameter profiles, derived from seismology. We show results here for complicated viscosity profiles with an asthenosphere in the upper mantle and a viscosity hill in the middle portion of the lower mantle. Because of the rapidity of the code, we may use this new method for exploring non-linear inversion problems by parameter sweeping.