A fully compressible viscoelastic framework for modeling the planetary response to surface loading and tidal forces

Computing the Love numbers for a general spherically symmetric, self-gravitating viscoelastic model is a classical problem in geodynamics, and is still relevant today in geodesy and planetary science. Indeed, the numerical values of the tidal or loading Love numbers are known to be sensitive to the planetary structure and to the internal rheological layering. Hence, measurements of Love numbers for planetary bodies provide constraints on their internal structure, while viscoelastic Love numbers for realistic Earth models are a key ingredient for modeling the response to the redistribution of hydrological and cryospheric loads. Here we present a novel numerical approach to the solution of viscoelastic field equations, aimed at evaluating Love numbers for spherically symmetric ( 1 − D) planetary structures characterised by linear viscoelastic rheology, both in the time and in the frequency domains. The solution scheme is efficiently implemented in the new Python framework pyALMA 4 and benchmarked against reference numerical results based on independently developed numerical tools like LoadDef and ISSM. The pyALMA 4 code extends the range of linear rheological laws available to the user and improves upon previous code ALMA 3 in terms of efficiency, avoiding expensive multi-precision computations. One major refinement is the introduction of compressibility, which improves the accuracy of the computed Love numbers. The versatility of pyALMA 4 is confirmed through two case studies of current interest in planetology, namely the isostatic deformations of the Martian lithosphere forced by fluctuations of the polar caps and the tidal response of the Jovian moon Europa in the presence of a subsurface ocean.