Developments in the physical modelling of the Earth’s interior and space geodesy make it possible nowadays to exploit the gravity signature and deformation patterns, including their time variations, caused by megathrust earthquakes at subduction zones with moment magnitude Mw higher than 9.0, as the 2004 Sumatran and 2011 Tohoku-Oki ones. In order to achieve these goals, it has been necessary to develop realistic self-gravitating, compressible Earth’s models, stratified in terms of density and rheological properties of the Earth’s interior. This new class of models allows us to comprehend some not yet fully appreciated mathematical aspects of viscous stress relaxation, such as the interplay between discrete and continuous relaxation spectra depending on the style of density stratification of the viscoelastic mantle or the effects on the gravity fields of sometimes used simplified treatments of compressibility. We show how GRACE and GOCE data allow us to invert for the mass redistribution, inner volume variations, gravity perturbations and surface deformation affecting areas and volumes larger than those embedding the gouge of the earthquakes, including the slip distribution over the fault surface. A correct interpretation of the mass redistribution process for megathrust earthquakes is important for the understanding of the physics of the ocean and Solid Earth coupling, causing the tsunami which struck Sumatra and Thailand and the eastern coast of Japan due to the huge amount of water washed out from the epicentral region as seen from GRACE data and physical modelling. In the present study we focus on the physics of the co-seismic and post-seismic gravity changes due to a Mw=7.0scenario normal-fault earthquake, comparable to the 1980 Irpinia earthquake. Our modelling provides the earthquake gravity effects within the perspective of the upcoming Next Generation Gravity Missions (NGGM), designed to detect the gravity anomalies caused by earthquake magnitudes as low as Mw=7 as well as the gravity anomalies due to the active tectonic processes responsible for the earthquakes. It is expected that the time-dependent gravity will be exploited at the GOCE spatial resolution and at GRACE time resolution, or even better, thanks to a new class of payload instrumentation, based on a laser ranging system measuring the distance variation between one or two pairs of satellites flying in formation at altitudes of about 300 km and at about 100 km separation, each satellite differently affected by the gravity changes of our Planet.