The Method of Lines in Computing Viscoelastic Relaxation of the Earth

The initial-value (IV) approach to the modeling of viscoelastic relaxation of spherical compressible self-gravitating Earth models to loading processes was designed as an alternative to the approach based on the Laplace transform applied to the time variable. The feature of overcoming the burden of inversion of Laplacian spectra of complex models encouraged extensive developments, including generalization for 2-D and 3-D viscosity distribution. Whereas our original IV formulation utilized the Euler scheme for integrating governing partial differential equations in time as an intrinsic part of the theory and a~series of boundary-value problems had to be solved, now we present a formulation, in which the governing equations are discretized in all spatial variables first (a~technique known as the method of lines). The obtained set of ordinary differential equations forms a linear IV problem in time, which can be numerically integrated by optimized library routines. Despite a~high degree of numerical stiffness emerging due to both underlying physics and applied discretization, appropriate integrating routines reveal the time evolution of relaxation curves of complex Earth models with outstanding speed and stability. It is an important feature of the new formulation that the processed IV problem can be recast into the form of a matrix eigenvalue (EV) problem. The resulting eigenspectrum covers the complete physical relaxation spectrum and also allows us to appraise and alleviate the impact of non-physical modes, which originate in spatial discretization. We show outputs of joint employment of the IV and EV strategies to the forward problem of viscoelastic relaxation of Earth models with complex spatial stratification, discuss applicability of various discretization grids and display attained computational times. Two examples: our code based on the EV strategy running on a Pentium III computes the entire viscoelastic spectrum of a model discretized by 30 radial layers within 0.05~sec per angular order, the same time of 0.05~sec is enough for the IV code to integrate one time step for a~model with 100 layers; the isostatic equilibrium can be reached within few tens of adaptive time steps.