# Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data

Isabelle Panet, Sylvain Bonvalot, Clément Narteau, Dominique Remy, Jean-Michel Lemoine

► To cite this version:

Isabelle Panet, Sylvain Bonvalot, Clément Narteau, Dominique Remy, Jean-Michel Lemoine. Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data. Nature Geoscience, 2018, 11, pp.367-373. 10.1038/s41561-018-0099-3. insu-03589357

HAL Id: insu-03589357

https://insu.hal.science/insu-03589357v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Migrating pattern of deformation prior to the Tohoku-Oki earthquake revealed by GRACE data

Isabelle Panet<fnref n="1" />*, Sylvain Bonvalot², Clément Narteau³, Dominique Remy² and Jean-Michel Lemoine⁴

Understanding how and when far-field continuous motions lead to giant subduction earthquakes remains a challenge. An important limitation comes from an incomplete description of aseismic mass fluxes at depth along plate boundaries. Here we analyse Earth's gravity field variations derived from GRACE satellite data in a wide space-time domain surrounding the $M_w$ 9.0 2011 Tohoku-Oki earthquake. We show that this earthquake is the extreme expression of initially silent deformation migrating from depth to the surface across the entire subduction system. Our analysis indeed reveals large-scale gravity and mass changes throughout three tectonic plates and connected slabs, starting a few months before March 2011. Before the Tohoku-Oki earthquake rupture, the gravity variations can be explained by aseismic extension of the Pacific plate slab at mid-upper mantle depth, concomitant with increasing seismicity in the shallower slab. For more than two years after the rupture, the deformation propagated far into the Pacific and Philippine Sea plate interiors, suggesting that subduction accelerated along 2,000 km of the plate boundaries in March 2011. This gravitational image of the earthquake's long-term dynamics provides unique information on deep and crustal processes over intermediate timescales, which could be used in seismic hazard assessment.

Giant earthquakes have mostly been recorded at subduction zones, where an oceanic plate penetrates into the convecting mantle. These events are generally described as an abrupt release of elastic strain accumulated in the upper plate under conditions of high coupling with the downgoing plate. Thus, their occurrence probably depends on patterns of strain accumulation, which are monitored in the upper plate thanks to space geodetic measurements of surface displacements. However, these observations alone are not sufficient to determine how giant earthquakes relate to convective processes and far-field plate motions<fnref n="1" />⁻³. High-accuracy determinations of crustal movements have evidenced slow-slip episodes of the subducting plate downdip of the seismogenic zone⁴,⁵, but it is still not possible to quantify changes in the on-going subduction conditions at greater depths. Yet, occurrences of seismic events with focal depths exceeding 200 km before great shallow earthquakes raise the question of their link with the deeper slab evolution at timescales of months to years⁶. Thus, improvements in our current understanding of giant earthquakes and forecasting capabilities call for a better description of plates and slab movements at depth. The only means to detect these transient mass changes associated with subduction is by observing Earth's gravity. In real time for earthquake early warning, surface gravimeters and seismic networks are now able to detect gravitational field variations before the arrival of seismic waves⁷,⁸. At longer time scales, Earth's mass redistributions are reflected in the monthly time variations of the gravity field monitored by the Gravity Recovery and Climate Experiment (GRACE) satellites since 2002, from near-global to 250-400 km resolution⁹. These data can bridge the gap between local and global patterns of plate dynamics over periods of decades to months, and provide unique information on the processes leading to giant ruptures at subduction zones.

**A gravitational image of the Tohoku-Oki earthquake.** We analysed the long-term dynamics of mass transport before and after

the 11 March 2011, $M_w$ 9.0 Tohoku-Oki earthquake in a 150° × 100° region around Japan. Previous investigations using GRACE and Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite gravity data have focused on the co- and post-seismic variations in the epicentral area and its vicinity¹⁰⁻¹⁶. Here, we search for earthquake-related gravity signals at a regional scale, considering the entire subduction system of the Northwest Pacific, from the Izu-Bonin and Ryukyu arcs to the south, to the Kuril-Kamchatka arc to the north. For that, we develop a space-time analysis of Earth's gravity vector variations reconstructed from the August 2002-June 2014 time series of GRACE Centre National d'Etudes Spatiales/ Groupe de Recherches de Géodésie Spatiale (CNES/GRGS) Release 3v1 monthly Stokes coefficients, with a resolution down to 250 km (ref. ¹⁷; Supplementary Section 1). Our approach differs from previous works in three key ways. First, we investigate large spatial scales, up to 1,600 km, in addition to more localized signals, down to 800 km. Second, we consider the entire variability at timescales of months to years in the data, without using prior models. Third, we enhance the signals following the geometry of the subduction system by analysing time changes of Earth's gravity gradients in relevant, rotated frames (Supplementary Section 2.2). Initially, we reconstruct the gravity gradients tensor every month from the GRACE geoids at different spatial scales, in the local south-east-up spherical frame. By progressively rotating the frame along the local radial axis, we look for gravity signals oriented along the strikes of both plate boundaries and subducted slabs. For each spatial scale and orientation, we map anomalous transient gravity gradients variations near the earthquake time using a wavelet transform of the time series (Supplementary Section 2.3.1). This approach is analogous to searching local velocity peaks in the time series smoothed at different temporal resolutions. The timing and the duration of the peak in a window around March 2011 inform on the predominance of a pre-seismic, co-seismic or post-seismic variation. Successively applied before and after removing a co-seismic step, we can determine

<footnote n="1">¹LASTIG LAREG, IGN, ENSG, Université Paris Diderot, Sorbonne Paris Cité, Paris, France. ²GET, Université de Toulouse, IRD, UMR 5563 CNRS, CNES, Toulouse, France. ³Institut de Physique du Globe de Paris, Université Paris Diderot, UMR 7154 CNRS, Sorbonne Paris Cité, Paris, France. ⁴GET, Université de Toulouse, CNES, UMR 5563 CNRS, IRD, Toulouse, France. *e-mail: isabelle.panet@ign.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1** | Pre-, co- and post-seismic variations of the Tohoku-Oki earthquake gravity signal. **a–d**, Time sequence of the 1,400-km-scale φφ gravity gradients in the local south-east-up spherical frame of unit vectors (e<sub>θ</sub>, e<sub>φ</sub>, e<sub>z</sub>), rotated in the direction Az₁. They correspond to φφ gravity gradients averages for 20–55° clockwise rotations (defining the direction range Az₁) of the frame about the radial axis (Supplementary Section 2.2). Star: 11 March 2011 earthquake epicentre; orange lines: plate boundaries⁴¹; thin violet lines: Pacific slab isodepth contours⁴², every 200 km (a,b) or 100 km (c,d). Tectonic plates: Pacific (Pac), Philippine Sea (PHS), Okhotsk (Ok), Eurasian (Eur); island arcs: Kuril (Ku), Izu-Bonin (IB), Marianna (M), Ryukyu (R), Honshu (H).</figcaption></figure>

whether the different components of a piece-wise linear model of the time series are part of an earthquake-related signal, or not (Supplementary Section 2.3.2).

The obtained gravitational image of the earthquake reveals two striking behaviours (Fig. 1). First, gravity variations start a few months before March 2011. Second, the earthquake signals near the epicentral area are surrounded by broad-scale variations throughout the entire subduction system. Starting in late 2010, regional gravity variations develop on both sides of the Boso triple junction (Fig. 1a). In March 2011, they spread along the subduction boundaries over three tectonic plates, reaching the southern Kuril, the Ryukyu and the Izu-Bonin arcs (Fig. 1b). This large-scale spreading is also found along the strike of the northern Philippine Sea subduction, from the Boso triple junction to south-central Japan (Supplementary Section 2.3.2). It is concomitant with a strong local gravity signal near the epicentre, at the monthly time resolution of GRACE. From April 2011, slower gravity variations progressively concentrate near the epicentral area, first on the trench (Fig. 1c) and then on the oceanic side of the subduction, east of Tohoku (Fig. 1d). In this long-term gravitational picture, the earthquake appears as the most striking event of a wide and silent deformation evolving over months to years.

At locations where the gravity signals start before the earthquake, there is a clear contrast between the pre-seismic behaviour of the time series and their long-term variability (Fig. 2b). West of the Japan subduction zone, we observe a large and uninterrupted gravity gradient increase from December 2010 to March 2011, synchronous with a decrease on the Philippine Sea plate, south of the triple junction (Fig. 3a,b). Over the previous seven years, all the series exhibit a constant and slow trend. For example, when applying our analysis between 2006 and 2009, no anomalous July–February changes are detected around Japan (Fig. 3a,b). These fast variations before the Tohoku-Oki earthquake are large with respect to the usual water signals and noise in this area (Supplementary Sections 3 and 5): on the same order of magnitude as the annual snow cycle over Japan, and at least two times larger than the seasonal ocean mass variations in the surrounding areas (Supplementary Section 3.2). In fact,

these anomalous signals can be detected in the original time series even before correction for the periodic components, especially west of Japan (Fig. 2a). Beyond the seasonal cycle, their unusual amplitude with respect to water signals and noise is evidenced in a statistical test of detection before the earthquake (Supplementary Section 2.4). We decompose the time series into a long-term linear trend and residuals, assuming a Gaussian distribution of independent monthly deviations. Stacking the time series in longitude across the pre-seismic anomalies, the 97.5th percentile level of the August 2002–May 2010 distribution is systematically approached or exceeded starting from December 2010 in the areas where the gravity signals precede the earthquake (Fig. 2c; spatial map in Fig. 3c,d and Supplementary Section 2.4). Furthermore, the residual values remain close to, or above the 90th percentile level until the end of the time series, even after subtracting the co-seismic and post-seismic signals (Fig. 2c and Supplementary Section 4.1). The abrupt gravity changes recorded before the rupture are thus persistent over time and spatially coherent over a wide ocean/islands area, whatever the subsequent seismic signals (Supplementary Section 3.3). Indicating a sudden movement of mass, largely completed before the earthquake, this pattern cannot be explained by water signals and we attribute its origin to a precursory activity at depth.

These results suggest that the Tohoku-Oki earthquake belongs to a broader sequence of large motions starting a few months before the rupture, the initial phase of this sequence being detectable in the gravity data before the event. The space-time distribution of the gravity signals is summarized in Fig. 4. In this diagram, the near-epicentral zone corresponds to the highest and most localized perturbations, in the middle of a larger area undergoing gravity changes over longer timescales. From late 2010, the gravity signals migrate from the deeper side of the Pacific plate subduction zone, where they spatially coincide with the slab at 100–350 km depth, to its shallow and even oceanic side, more than 1,000 km within the Pacific plate interior in March 2011. This migration points to a depth-to-surface dynamics of mass redistributions during the months before and after the earthquake, over a volume that includes and extends far beyond the seismic rupture. The sequence involves commensurate large-scale mass transfers in the
<figure><img src="image_3.png" /><figcaption>**Fig. 2** | Time series of the gravity signals at different stages of their analysis. **a,** Original time series for the 1,400 km scale, Az, φφ gravity gradients (GRACE de-aliasing ocean model restored in the oceanic areas; Supplementary Sections 1 and 4.1). **b,** Time series at the same points, after removing the annual, semi-annual and 161-day cycles. **c,** Time series averaged across the pre-seismic anomaly, for each latitude in the three considered areas, after removing the periodic components, the long-term trend, co-seismic and post-seismic signals. 90th and 97.5th percentiles of the August 2002-May 2010 distribution: violet and pink dashed lines. Blue dots: December 2010 and February 2011 values; pink dots: March 2011 value.</figcaption></figure>

pre-seismic and early post-seismic phases, as indicated by similar amplitudes of the 1,400-km-scale signals for a limited depth range of sources. This regional evolution is further supported by the constant dominant northeast-southwest orientation of the gravity signal dur-

ing the pre-seismic and the co-seismic phases (Supplementary Section 2.3). Consistent with the geometry of the entire subduction system, it can therefore be a broad subduction acceleration releasing accumulated stresses from depth to the surface all around the Boso triple junction.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** | Detection of anomalous gravity variations around Japan. **a,** Cumulated July–February anomalies stacked over 2006–2009, obtained by successively applying the same analysis as in Fig. 1. **b,** Zoom-out of Fig. 1a. The blue lines delineate the studied area. **c,d,** Gravity anomalies outside the 2.5th–97.5th percentile range of the August 2002–May 2010 distributions, resolving the hydrological signals over Asia at the 1,000 km scale (Supplementary Section 3.3; frame rotation 40°). **c,** August 2010, showing intense continental monsoon signals. Note the absence of leakage into the Japan Sea. **d,** February 2011. In **a,b,** anomalies correspond to an abrupt change in trend, which may start before large instantaneous values are reached (**d**).</figcaption></figure>
<figure><img src="image_5.png" /></figure>

**Fig. 4** | Space-time diagram of the 2011 Tohoku-Oki earthquake gravity signals. Cross-section of the pre-, co- and post-seismic anomalies at 1,400 km and 800 km scales. The inset shows the transect line. Each positive or negative lobe of the gravity gradients cross-section is described by a box. Its location and width define the box horizontal dimension. Its time coverage defines the vertical dimension of the box. The colour scale shows the absolute extremum value. Red dashed lines: 800-km-scale boxes; other boxes: 1,400-km-scale signals. Black dots: depth (right vertical scale) of the 10 March 2010–10 March 2011 M2 earthquakes¹³, within a 100 km horizontal distance from the transect. Dashed black line: top of the Pacific slab⁴².

**A deep precursory phase.** In this scenario, we explain the widespread pre-seismic gravity gradient increase west of Japan by a diffuse mass decrease in the slab (that is, a precursory intra-slab extension episode). Such deformation is consistent with the reported acceleration of the seismic release and extensional mechanisms below 50 km depth from January 2011 up to the earthquake time¹⁸. In addition, a significantly larger aseismic motion initiated months before is likely to occur deeper in the slab. Indeed, Fig. 4 shows that the precursory gravity changes are located above a section of the slab where seismicity vanishes. As an equivalent representation, we estimate the amount of quasi-static normal faulting in the slab leading to commensurate gravity signals. We find a slip of 40 cm along a 100 × 1,200 km² fault plane dipping at 60° from 245 to 330 km depth, which corresponds to a Mw 8.4 earthquake over a few months (Supplementary Section 6.2.1). However, at upper mantle depths, such a large-scale increase of the pre-seismic gravity gradient requires a relatively distributed source, more than a localized slip (Supplementary Section 6.1.5), which could be another argument for ductile deformation.

Coherent with a northwest-oriented upper-mantle slab pull downdip of northeastern Japan¹, the gravity signal orientation suggests that the sinking slab extends under the pulling force of its deeper root¹,¹⁹. This transient deformation manifests the long-term slab evolution in a regional configuration where vertically and laterally extensional stresses accumulate at depth in the Pacific plate slab, beneath northeastern Japan. In this zone of transition between two steeply dipping arc systems with different strikes, the Kuril–Kamchatka and Izu–Bonin–Marianna arcs²⁰,²¹, such stresses may result from the coupling with the upper plate³,²², and from the geometry of the subduction system²⁰,²¹. Consistently, evidence of slab thinning and tearing below 250 km depth is found in seismic data (for example, refs ²¹,²³). In addition, the Pacific plate slab beneath Japan appears only partially attached to its subducting plate, which enters into the mantle along a more westerly path than that imposed by the deep pull force¹. South of the Boso triple junction, the deep pull can directly accelerate subduction from depth to the surface, as the slab is fully connected to its oceanic plate and decoupled from

the upper plate along the Izu–Bonin Trench³,²³. In agreement, the negative pre-seismic signal indicates mass influx into the Izu–Bonin subduction zone.

A deep source for the pre-seismic gravity signal is reinforced when investigating the potential impact of shallower sources. First, the northeastern Japanese arc undergoes east–west compressional stresses from the eastward motion of the Amur plate²⁴, ruling out regional extension of the lithosphere in areas of gravity gradient increase. In addition, the observed anomaly does not align with the north–south-trending Amur–Okhotsk plate boundary. Second, the amount of lithospheric mass decrease that would fit our gravity observations corresponds to a 1.8 ± 0.6 cm subsidence over the islands (Supplementary Section 6.2.3)—a small value, but regionally coherent, which could make it detectable in northern Japan. Finally, overriding plate dilation would more probably result from aseismic slip along the Japan subduction; however, the gravity signals would shift eastward, or include unobserved negative anomalies near the trench (Supplementary Section 6.2.2). Thus, the full spatial pattern of the GRACE pre-seismic anomalies is hardly explained in these alternative hypotheses, further supporting deeper motions, faster and wider than the long-term accelerating slip reported along the Japan Trench²⁵⁻²⁷. Interestingly, widespread gravity variations are also detected before the 2015 Mw 7.8 Nepal event²⁸.

**Propagation of deformation within the oceanic plates.** This regional-scale precursory phase explains why, on both oceanic plates around the triple junction, the dimension of the co-seismic signals derived from gravity data far exceeds that predicted from co-seismic or post-seismic slip distributions based on crustal movement and tsunami data²⁹,³⁰, shown in Fig. 5 and Supplementary Section 6.1.5. In this figure, high values of the modelled gravity gradients indicate co-seismic mass decrease from subsidence and/or crustal dilation³¹; low values result from incoming mass flux. The magnitudes of the gravity and crustal deformation-based co-seismic signals agree on the upper plate, but their orientations and spatial extents differ, as already noticed at local scales¹³,¹⁴. The orientation of the crustal movement signals coincides with the alignment
<figure><img src="image_6.png" /><figcaption>**Fig. 5** | Synthetic gravity signals from co-seismic slip, afterslip and mantle visco-elasticity. **a–d**, $Az_1$, 800-km-scale φφ gravity gradients. **a**, GRACE-observed March 2011 anomalies. **b**, Co-seismic slip model<sup>29</sup>. **c,d**, Cross-sections along the transect of Fig. 4 of the φφ gradients from the co-seismic and 279-days afterslip<sup>30</sup> models, with GRACE co-seismic (March 2011) and post-seismic (first year) anomalies. **e**, $Az_1$, 1,400-km-scale φφ gravity gradients variations over 11 March 2012–23 May 2014, from the asthenosphere visco-elastic relaxation triggered by the co-seismic slip<sup>29</sup>. Here we consider a Maxwell asthenospheric rheology with a $5 \times 10^{18}$ Pa s viscosity (Burgers rheology: Supplementary Section 6.1.5).</figcaption></figure>

ment of the Pacific plate subduction east of Honshu, but not with the pre-seismic and co-seismic gravity signals, which align along the regional strike of the subduction. These differences show that the GRACE-based deformation pattern of the plate/slab system at proximity of the giant earthquake is much larger than expected from surface motions and seismic data.

The large-scale spreading of the observed March 2011 gravity anomalies over the ocean (Fig. 1b) suggests that the rupture is followed by a short-term aseismic acceleration of the Pacific and Philippine Sea subduction along 2,000 km of the boundaries and 1,500 km within the plate interiors, implying some lithosphere mobility over the underlying mantle<sup>32,33</sup>. An acceleration of these two plates after the earthquake is also recorded from seismic data<sup>34</sup>,

and explains well the continent/ocean asymmetry of the GRACE anomaly (Supplementary Section 6.3.1). From the amplitude of the signals, we estimate the increase in mass of the oceanic plates near the trench and the decrease in mass in their interiors. We find it equivalent to ~15 cm of horizontal motion towards the trench of a 50-km-thick lithosphere, on the same order of magnitude as the pre-seismic slab deformation (Supplementary Section 6.3). This observation supports regional mass transfers extending as far as the adjacent oceanic plates. As a consequence of these widespread motions, the observed post-seismic gravity variations shift towards the ocean with respect to the signals derived from afterslip or visco-elastic relaxation models using surface data (Fig. 5). Their different spatial structures at small and large scales (Supplementary Section 2.3)
suggest an interplay between continuous slip and viscous relaxation of stress perturbations in the mantle³⁵, already proposed from composite patterns of GRACE post-seismic variations¹⁵. In the endmember case of a pure visco-elastic asthenosphere relaxation (Fig. 5e and Supplementary Section 6.1.5), our results support nonlinear rheologies with low transient viscosities near the trench, and Maxwell viscosities slightly below $5 \times 10^{18}$ Pa s for the slower response of the oceanic asthenosphere (for example, refs ¹⁶,³⁶,³⁷).

Geodesy efficiently highlights a range of transient surface displacements reflecting slow slip near plate boundaries (for example, refs ⁴,⁵). The main outcome of the gravitational image presented here is that episodic mass transfers at intermediate timescales of months are detected within the entire subduction system, from depth to the surface. Some of these deep changes are inaccessible to other observation systems. Data from the Global Navigation Satellite System networks record no surface displacements before the earthquake at the spatial and temporal scales of the pre-seismic gravity variations³⁸. Furthermore, the orientation of both the co-seismic and early post-seismic crustal motions towards the trench east of Honshu³⁹ does not coincide with that of the observed gravity signals (Fig. 5a,b). These deep mass fluxes indeed require highly deformable layers along the slab, which decouple the slab from the surrounding mantle, limiting the effect of its viscosity that would otherwise hinder deformation rates of $10^{-12}$ s⁻¹. Besides, a high stress sensitivity of the strain rate could facilitate large responses to small stress variations, such as the Coulomb stress changes associated with our simplified pre-seismic slip model, on the 1 kPa order of tidal stresses on the earthquake fault plane⁴⁰. At the same time, such a specific rheology at depth may localize strain along the subducted lithosphere and explain the absence of crustal deformations.

Thus, while measurements of strain accumulation in the upper plate near major boundaries point to areas prone to a large seismic rupture, satellite gravity allows us to replace its occurrence in the context of Earth's global dynamics. Our results provide evidence that Earth's convection system includes timescales as short as months as slabs slowly penetrate into the upper mantle. This is where we find the origin of the Tohoku-Oki earthquake, along with detectable precursor signals in the gravity field. Their possible occurrence before the other M9 events in the GRACE era is now to be explored, together with a link between the slow convective motions in the mantle and giant earthquakes.

**Data availability.** The data that support the findings of this study are available from the corresponding author upon request.

Received: 31 August 2017; Accepted: 7 March 2018;
Published online: 9 April 2018

**References**

1. Conrad, C., Bilek, S. & Lithgow-Bertelloni, C. Great earthquakes and slab pull: interaction between seismic coupling and plate-slab coupling. *Earth Planet. Sci. Lett.* **218**, 109–122 (2004).

2. Stein, S. & Okal, E. Ultralong period seismic study of the December 2004 Indian Ocean earthquake and implications for regional tectonics and the subduction process. *Bull. Seis. Soc. Am.* **97**, S279–S295 (2007).

3. Scholz, C. H. & Campos, J. The seismic coupling of subduction zones revisited. *J. Geophys. Res.* **117**, B05310 (2012).

4. Hirose, H., Hirahara, K., Kimata, F., Fujii, N. & Miyazaki, S. A slow thrust slip event following the two 1996 Hyuganada earthquakes beneath the Bungo Channel, southwest Japan. *Geophys. Res. Lett.* **26**, 3237–3240 (1999).

5. Dragerts, H., Wang, K. & James, T. S. A silent slip event on the deeper Cascadia subduction interface. *Science* **292**, 1525–1528 (2001).

6. Mogi, K. Deep seismic activities preceding the three large 'shallow' earthquakes off south-east Hokkaido, Japan—the 2003 Tokachi-oki earthquake, the 1993 Kushiro-oki earthquake and the 1952 Tokachi-oki earthquake. *Earth Planets Space* **56**, 353–357 (2004).

7. Montagner, J.-P. et al. Prompt gravity signal induced by the 2011 Tohoku-Oki earthquake. *Nat. Commun.* **7**, 13349–13355 (2016).

8. Vallée, M. et al. Observations and modeling of the elastogmetry signals preceding direct seismic waves. *Science* **358**, 1164–1168 (2017).

9. Tapley, B., Bettadpur, S., Ries, J., Thompson, P. F. & Watkins, M. GRACE measurements of mass variability in the Earth system. *Science* **305**, 503–505 (2004).

10. Matsuo, K. & Heki, K. Coseismic gravity changes of the 2011 Tohoku-oki earthquake from satellite gravimetry. *Geophys. Res. Lett.* **38**, L00G12 (2011).

11. Wang, L., Shum, C. K., Simons, F. J., Tapley, B. & Dai, C. Coseismic and postseismic deformation of the 2011 Tohoku-oki earthquake constrained by GRACE gravimetry. *Geophys. Res. Lett.* **39**, L07301 (2012).

12. Cambiotti, G. & Sabadini, R. Gravitational seismology retrieving Centroid-Moment-Tensor solution of the 2011 Tohoku earthquake. *J. Geophys. Res.* **118**, 183–194 (2013).

13. Fuchs, M., Bouman, J., Broerse, T., Visser, P. & Vermeersen, B. Observing coseismic gravity change from the Japan Tohoku-Oki 2011 earthquake with GOCE gravity gradiometry. *J. Geophys. Res.* **118**, 5712–5721 (2013).

14. Dai, C. et al. Improved constraints on seismic source parameters of the 2011 Tohoku earthquake from GRACE gravity and gravity gradient changes. *Geophys. Res. Lett.* **41**, 1929–1936 (2014).

15. Tanaka, Y. & Heki, K. Long- and short-term postseismic gravity changes of megathrust earthquakes from satellite gravimetry. *Geophys. Res. Lett.* **41**, 5451–5456 (2014).

16. Han, S.-C., Sauber, J. & Pollitz, F. Broadscale postseismic gravity change following the 2011 Tohoku-oki earthquake and implication for deformation by visco-elastic relaxation and afterslip. *Geophys. Res. Lett.* **41**, 5797–5805 (2014).

17. Bruinsma, S., Lemoine, J.-M., Biancale, R. & Valès, N. CNES/GRGS 10-day gravity field models (release 2) and their evaluation. *Adv. Space Res.* **45**, 587–601 (2010).

18. Bouchon, M. et al. Potential slab deformation and plunge prior to the Tohoku, Iquique and Maule earthquakes. *Nat. Geosci.* **9**, 380–383 (2016).

19. Billen, M. Slab dynamics in the transition zone. *Phys. Earth Planet. Int.* **183**, 296–308 (2010).

20. Yamaoka, K., Fukao, Y. & Kumazawa, M. Spherical shell tectonics: effects of sphericity and inextensibility on the geometry of the descending lithosphere. *Rev. Geophys.* **24**, 27–53 (1986).

21. Obayashi, M., Yoshimitsu, J. & Fukao, Y. Tearing of stagnant slab. *Science* **324**, 1173–1175 (2009).

22. Astiz, L., Lay, T. & Kanamori, H. Large intermediate-depth earthquakes and the subduction process. *Phys. Earth Planet. Inter.* **53**, 80–166 (1988).

23. Fukao, Y. & Obayashi, M. Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity. *J. Geophys. Res.* **118**, 5920–5938 (2013).

24. Taira, A. Tectonic evolution of the Japanese islands arc system. *Annu. Rev. Earth Planet. Sci.* **29**, 109–134 (2001).

25. Ozawa, S. et al. Preceding, coseismic and postseismic slips of the 2011 Tohoku earthquake, Japan. *J. Geophys. Res.* **117**, B07404 (2012).

26. Mavromatis, A., Segall, P. & Johnson, K. A decadal-scale deformation transient prior to the 2011 $M_w$ 9.0 Tohoku-oki earthquake. *Geophys. Res. Lett.* **41**, 4486–4494 (2014).

27. Yokota, Y. & Koketsu, K. A very long-term transient event preceding the 2011 Tohoku earthquake. *Nat. Commun.* **6**, 5934–5938 (2015).

28. Chen, S. et al. Gravity increase before the 2015 $M_w$ 7.8 Nepal earthquake. *Geophys. Res. Lett.* **43**, 111–117 (2016).

29. Minson, S. E. et al. Bayesian inversion for finite fault earthquake source models - II: the 2011 great Tohoku-oki, Japan earthquake. *Geophys. J. Int.* **198**, 922–940 (2014).

30. Perfettini, H. & Avouac, J.-P. The seismic cycle in the area of the 2011 $M_w$ 9.0 Tohoku-Oki earthquake. *J. Geophys. Res.* **119**, 4469–4515 (2014).

31. Han, S.-C., Shum, C.-K., Bevis, M., Ji, C. & Kuo, C.-Y. Crustal dilatation observed by GRACE after the 2004 Sumatra-Andaman earthquake. *Science* **313**, 658–662 (2006).

32. Kawakatsu, H. et al. Seismic evidence for sharp lithosphere-asthenosphere boundaries of oceanic plates. *Science* **324**, 499–502 (2009).

33. Heki, K. & Mitsui, Y. Accelerated Pacific plate subduction following interplate thrust earthquakes at the Japan trench. *Earth Planet. Sci. Lett.* **363**, 44–49 (2013).

34. Uchida, N., Asano, Y. & Hasegawa, A. Acceleration of regional plate subduction beneath Kanto, Japan, after the 2011 Tohoku-oki earthquake. *Geophys. Res. Lett.* **43**, 9002–9008 (2016).

35. Wang, K., Hu, Y. & He, J. Deformation cycles of subduction earthquakes in a visco-elastic Earth. *Nature* **484**, 327–332 (2012).

36. Pollitz, F., Bürgmann, R. & Banerjee, P. Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth. *Geophys. J. Int.* **167**, 397–420 (2006).

37. Sun, T. et al. Prevalence of viscoelastic relaxation after the 2011 Tohoku-oki earthquake. *Nature* **514**, 84–87 (2014).

38. Hasegawa, A. & Yoshida, K. Preceding seismic activity and slow slip events in the source area of the 2011 Mw 9.0 Tohoku-oki earthquake: a review. *Geosci. Lett.* **2**, 6–18 (2015).
39. Ozawa, S. et al. Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki earthquake. *Nature* **475**, 373–376 (2011).

40. Ide, S., Yabe, S. & Tanaka, Y. Earthquake potential revealed by tidal influence on earthquake size-frequency statistics. *Nat. Geosci.* **9**, 834–837 (2016).

41. Bird, P. An updated digital model of plate boundaries. *Geochem. Geophys. Geosyst.* **4**, 1027–1078 (2003).

42. Hayes, G. P., Wald, D. J. & Johnson, R. L. Slab1.0: a three-dimensional model of global subduction zone geometries. *J. Geophys. Res.* **117**, B01302 (2012).

43. *Bulletin of the International Seismological Centre* (2013); http://www.isc.ac.uk/iscbulletin/search/webservices/bulletin/

## Acknowledgements

We thank CNES for financial support through the TOSCA committee. C.N. acknowledges the financial support from the UnivEarthS LabEx programme of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). We thank H. Perfettini for providing co-seismic and post-seismic slip distribution models and their predicted surface displacements. We thank G. Métris for providing the software for the differentiation of spherical harmonics, and B. Romanowicz for important comments on the manuscript and the figures. Valuable reviews, including a review by K. Heki, contributed to significantly improve our work. All figures have been made using the GMT software. This is IPGP contribution number 3926.

## Author contributions

I.P. designed the four-dimensional gravity analysis and performed all data analyses. S.B., D.R. and I.P. conducted the gravity modelling of earthquake-related signals and wrote the corresponding sections of the Supplementary Information. C.N. conceived the pre-seismic statistical data analysis and the space-time diagram. J-M.L. provided information on the GRACE geoid models. All authors discussed the analyses and their results at all stages. I.P. wrote the manuscript and data analysis sections of the Supplementary Information with input from all co-authors.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41561-018-0099-3.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Correspondence and requests for materials** should be addressed to I.P.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.