Numerous observations pertaining to the magnitude 9.0 2011 Tohoku-oki earthquake (offshore Japan) have led to new understanding of subduction zone earthquakes. By synthesizing published research results and our own findings, we explore what has been learned about fault behavior and Earth rheology from the observation and modeling of crustal deformation before, during, and after the earthquake. Before the earthquake, megathrust locking models based on land-based geodetic observations correctly outlined the along-strike location of the future rupture zone. Their incorrect definition of the locking pattern in the dip direction demonstrates the need to model the effects of interseismic viscoelastic stress relaxation and stress shadowing. The observation of decade-long accelerated slip downdip of the future rupture zone raises new questions on fault mechanics. During the earthquake, seafloor geodetic measurements revealed huge coseismic displacements (up to 31 m). Modeling of bathymetry difference before and after the earthquake suggests >60 m of coseismic slip of the most seaward 40 km of the fault in the main rupture area, with the slip peaking at the trench. Large differences in shallow slip between published rupture models are due mainly to the near absence of near-trench deformation measurements, but model simplifications in fault and seafloor geometry also bear large responsibility. After the earthquake, seafloor geodetic measurements provided unambiguous evidence for the dominance of viscoelastic relaxation in short-term postseismic deformation. There is little deep afterslip in the fault area where the decade-long pre-earthquake slip acceleration is observed. Investigating the physical processes responsible for the complementary spatial distribution of pre-slip and afterslip calls for new scientific research.


The 2011 moment magnitude (Mw) 9.0 Tohoku-oki earthquake (offshore Japan) is by far the most well recorded great earthquake in history. Massive seismic, tsunami, and geodetic data were recorded by global and local monitoring networks, augmented by other types of data collected through marine surveys and ocean drilling. These data have exposed the anatomy of a megathrust earthquake in unprecedented detail and enabled new research on the geodynamics of earthquake cycles. In this article, we discuss what has been learned from geodetically observed crustal deformation before, during, and after this earthquake. We do not attempt a comprehensive review, but only select what we deem to be the most relevant material from published works and include our own findings. We focus on static and quasi-static deformation, although the scientific value of the lessons learned can be fully appreciated only in the context of multidisciplinary observations and modeling that encompass dynamic rupture processes, tsunami generation, and fault zone geology.

When addressing earthquake deformation, it is important to point out a rather unique aspect of the Tohoku-oki rupture: the concentration of large slip on a relatively small part of the fault. As is widely recognized and will be further discussed in the Coseismic Deformation section of this paper, the maximum fault slip exceeds 50 or 60 m. However, the strike length of the rupture zone is only ∼300 km, typical of earthquakes of Mw ∼8.1 based on the statistics of Wells and Coppersmith (1994) and in sharp contrast with the rupture zone of the Mw 9.2 2004 Sumatra earthquake, Indonesia (Fig. 1). The two events shown in Figure 1 appear to represent two end-member types of giant earthquakes, with the Tohoku-oki type being very compact and the Sumatra type being very spread out. Accordingly, they contrast with each other in rupture duration, ∼2.5 min for Tohoku-oki and 9.0 min for Sumatra (Lay and Kanamori, 2011). Whatever controls the contrasting rupture behavior is an important subject of future research. Some of the lessons learned from the Tohoku-oki earthquake may be generally applicable to all subduction earthquakes, but some may be unique to the compact type.

For the purpose of this study, it is useful to reiterate two important characteristics of the Japan Trench subduction zone. First, the Japan Trench is an end-member cold-slab subduction zone, owing to the old age (ca. 130 Ma) of the subducting Pacific plate and the fast subduction rate (∼8 cm yr–1). As a result, the megathrust produces small interplate earthquakes to depths of as much as 60 km (Nakamura et al., 2016), although large events tend to be confined to shallower depths (e.g., Yamanaka and Kikuchi, 2004). Other consequences of cold-slab subduction include very active arc volcanism, deep extension of intraslab seismicity (Wada and Wang, 2009), and absence or rarity of episodic slow slip events around the mantle wedge corner accompanied with seismic tremor (Gao and Wang, 2017). Second, except for the southernmost and northernmost areas where seamount subduction is taking place, the surface of the subducting plate, and by inference that of the megathrust, is moderately smooth. The moderate smoothness allows the development of seismic rupture as large as in the 2011 Tohoku-oki event (Wang and Bilek, 2014; Gao and Wang, 2014; Scholl et al., 2015), but the presence of low-amplitude roughness, due to the horst-graben structure of the incoming plate and shortage of trench sediment, is proposed to be responsible for observed spatially variable and complex slip behavior exhibiting seismic rupture, repeating earthquakes, and creep events (Wang and Bilek, 2014).


Fault Locking and Interseismic Stress Relaxation

Before the Mw 9 Tohoku-oki earthquake, several models of interseismic locking of the Japan Trench megathrust had been published. They were based on inversion of land-based geodetic observations, primarily those of the ∼400 continuously monitoring Global Navigation Satellite System (GNSS) stations in northern Honshu that had been operational since the mid-1990s. These models did an excellent job in forecasting the along-strike location of the Tohoku-oki earthquake, portraying the future rupture area as a locked patch. But they all did a poor job in the dip direction, in that the reported locked zone was much farther downdip (landward) than the actual rupture zone in 2011. The locked patch outlined in Figure 2 (black dashed line) approximately represents that determined by Nishimura et al. (2004), Suwa et al. (2006), Hashimoto et al. (2009), and Loveless and Meade (2010).

The excellent along-strike performance of these locking models demonstrates that land-based geodetic data can resolve along-strike variations of megathrust locking at wavelengths comparable to the distance from the network to the offshore seismogenic zone. The poor downdip performance has both observational and theoretical reasons. The observational reason is the lack of near-trench resolution of land-based geodetic measurements, and the theoretical reason is the neglect or insufficient account of viscoelastic stress relaxation (Wang et al., 2012). The observational issue is now widely recognized and does not need further explanation; this difficulty can be overcome by conducting seafloor, near-trench geodetic monitoring.

The theoretical issue of viscoelastic stress relaxation is more fundamental and cannot always be resolved with more observations. As explained by Wang et al. (2012) and Wang and Tréhu (2016), the viscoelastic behavior of the asthenospheric mantle causes stress relaxation not only during postseismic deformation but also during interseismic deformation. After an earthquake, stresses induced by the earthquake are relaxed, causing time-dependent crustal deformation (see the Postseismic Deformation and Afterslip section). In the interseismic period, stresses being built up due to fault locking are relaxed at the same time, even when crustal deformation is no longer changing with time. This interseismic relaxation allows the effect of fault locking to reach far inland to cause crustal deformation, provided that the locking is sufficiently extensive along strike. If the viscoelastic effect is neglected or is insufficiently accounted for, the deformation in the inland area has to be incorrectly attributed to fault locking extending to large depths. This problem has also been numerically demonstrated by Li et al. (2015) using observations in Peru and northern Chile. More poorly understood is the expected interseismic deformation in the vertical direction (or crustal tilt) in a viscoelastic Earth (Wang and Tréhu, 2016). We will visit this issue when discussing vertical coseismic deformation in the Dominance of Viscoelastic Stress Relaxation section.

The vast majority of interseismic locking models worldwide assume an elastic Earth and suffer from the same problem. Therefore, the importance of interseismic stress relaxation is well worth reiterating. For a locked zone of short strike length such as in Costa Rica (Protti et al., 2014), this is less of a problem, because the neighboring creeping zones help to prevent locking-induced deformation from reaching very far inland. The longer a locked zone is along strike and the longer it has been locked, the more difficult the problem becomes, to the extent that the value of interseismic geodetic observations in constraining the downdip extent of locking can be questionable (Wang and Tréhu, 2016).

Updip Extent of Megathrust Locking

All of the locking models published prior to the 2011 earthquake as referenced above feature megathrust creep at the trench at the subduction rate, including the area that would eventually host the Tohoku-oki earthquake. It is a feature that also appears in locking models of many other subduction zones. It is not required by the geodetic data in most places, particularly not the Japan Trench before the 2011 earthquake where land-based measurements had no resolving power near the trench (Loveless and Meade, 2011). It appears to have been introduced for inversion convenience: it is simpler to fix slip deficit at zero at the trench. It may also be the result of confusing fault properties with fault motion (Wang and Dixon, 2004): the shallowest part of the megathrust is often thought to exhibit velocity strengthening and thus weak mechanical locking, and this assumed condition is often mistakenly translated into creeping at full speed.

Prior to the 2011 earthquake, almost no seafloor data were collected or processed that could be used to constrain the locking or creeping state of the shallow megathrust near the trench. The motion of two pioneer seafloor GNSS stations (Fig. 2) reported by Matsumoto et al. (2008) was consistent with more complete locking in the future main rupture area (station MYGI) and some creeping to the south (station FUKU). More seafloor GNSS data collected prior to but published after the 2011 earthquake (Sato et al., 2013) would still be inadequate in resolving locking or creeping near the trench because of their relatively large distance from the trench.

Until near-trench monitoring becomes routinely available, the interseismic locking or creeping state of the shallow megathrust depends on the preference of the researcher who does the geodetic inversion. Kinematically, full-speed creep at the trench could not possibly be sustained over the entire interseismic period at the Japan Trench, or there would not have been enough slip deficit to allow the >50 m coseismic slip in 2011. One can also reasonably argue that because of the effect of stress shadowing, the fault segment updip of a fully locked patch is not expected to slip by itself in a sustained fashion regardless of its own frictional behavior (Wang and Dixon, 2004; Hetland and Simons, 2010; Wang and Tréhu, 2016). Therefore, even if seafloor observations indicate a no-slip state of the shallow fault, such as off Peru (Gagnon et al., 2005), we still do not know whether it is mechanically locked or is just in a stress shadow.

Locking of the shallow megathrust was indirectly and implicitly inferred from seafloor pressure-gauge data collected prior to the 2011 earthquake. By analyzing these data, Ito et al. (2013) proposed a month-long slow slip event just before the earthquake in the shallow part of the megathrust within the future main rupture area (Fig. 2), and possibly a similar event earlier. To argue for slow slip events, the necessary assumption is that the fault is otherwise locked to a high degree.

Pre-Slip Behavior

Retroactive examination of pre-earthquake GNSS time series led to an important discovery. Ozawa et al. (2012), Mavrommatis et al. (2014), and Yokota and Koketsu (2015) independently identified a fault area directly downdip of the 2011 main rupture that exhibited faster creep for ∼9 yr leading to the Tohoku-oki earthquake. Ozawa et al. (2012) described it as aseismic slip associated with relatively large (Mw 7.4–7.7) interplate earthquakes in this area, Mavrommatis et al. (2014) described it as an acceleration of slip rate, and Yokota and Koketsu (2015) described it as a long-duration slow slip event. Most of the slip is inferred to have occurred ∼25–50 km depth (Fig. 2, yellow), although that reported by Ozawa et al. (2012) is slightly farther updip. Mavrommatis et al. (2014) and Yokota and Koketsu (2015) also found that, over the same time window, a patch of the fault in the same depth range but located to the north of the 2011 rupture exhibited decelerated slip, or enhanced locking (Fig. 2, blue). As data examples, time series of data from four of the many GNSS stations that were analyzed by Yokota and Koketsu (2015) are reproduced in Figure 3.

This decadal deep pre-slip is one of the many surprises presented to the scientific community by the Tohoku-oki earthquake. It is poorly understood, although attempts have been made to explain it in terms of interseismic shrinking of the locked zone in the context of rate-and-state–dependent friction (Ozawa et al., 2012; Mavrommatis et al., 2014). However, the depth range of the pre-slip and accompanying large interplate events is not surprising. Given the cold thermal field in this subduction zone, the megathrust is expected to exhibit a frictional behavior in this depth range (Gao and Wang, 2014, 2017). One may think that the deep pre-slip “triggered” the Tohoku-oki rupture, but it is also possible that the pre-slip is an integral part of the initiation process of the Tohoku-oki rupture. Of course, given the limited time span of continuous monitoring (Fig. 3), there is some arbitrariness in defining the “normal” linear trend of motion from the earlier portion of the GNSS time series (Fig. 3), and one cannot exclude the possibility that similar accelerated creep had occurred before the era of space geodesy.


GNSS Displacements and Rupture Models

Japan’s GNSS network provided the most spectacular geodetic measurements of crustal deformation caused by a great earthquake in history (Fig. 4, red arrows). For the first time, seafloor displacements right on top of the rupture zone were measured, and the site nearest the trench (50 km) recorded 31 m of horizontal motion (site GJT3 in Fig. 4A) (Kido et al., 2011). The seafloor measurements provided some of the most critical information for understanding this earthquake. For example, rupture models developed by Ozawa et al. (2011) and Simons et al. (2011) before these seafloor data became available put the maximum slip ∼100 km from the trench, very different from later models that used these data.

For a shallowly buried and shallowly dipping thrust fault, horizontal surface displacements in the rupture area provide robust constraints for fault slip. The vertical component is more sensitive to local tilting caused by slip heterogeneity and off-fault permanent deformation. For example, an abrupt change in fault slip magnitude cannot cause opposite motion of neighboring sites in the horizontal direction but can do so in the vertical direction. Thus, precisely fitting sparse measurements of vertical motion may not lead to an accurate fault slip model. The large horizontal motion of GNSS seafloor sites during the Tohoku-oki earthquake (Fig. 4A) can be readily explained with simple models. The vertical measurements, displaying rather different values over short distances (Fig. 4B), are more demanding of model details.

Many coseismic slip models of the 2011 earthquake have been derived by inverting the GNSS and/or other types of data recorded during the earthquake. Tajima et al. (2013) reviewed 26 slip models, all published in the same year as the earthquake, plus some other models that did not determine fault slip vectors. Brown et al. (2015) included 40 published slip models to derive static stress drop from their slip vector distribution. Building on Brown et al. (2015), Sun et al. (2017) compiled 45 published slip models and summarized the types of data used by each model. Some of the 26 models reviewed by Tajima et al. (2013) were subsequently superseded by later versions published by the same author groups; Brown et al. (2015) and Sun et al. (2017) only included the latest version in their compilation. All of the models that used seismic wave or high-rate GNSS data (e.g., Yokota et al., 2011; Shao et al., 2012; Yue and Lay, 2013), and some of the models that used tsunami waveform data (e.g., Satake et al., 2013), are of the finite-fault type describing spatiotemporal evolution of the rupture. For discussing static deformation in this paper, we only consider the net fault slip portrayed by these models and leave aside the valuable scientific information on the evolution.

The mean slip magnitude distribution of 43 of the 45 models compiled by Sun et al. (2017) plus that of Freed et al. (2017) is shown in Figure 4, with the standard deviation shown in Figure 5A and the ratio of the standard deviation to the mean slip in Figure 5B. Two of the 45 models compiled by Sun et al. (2017) are not included in the mean model, because one of them only provides an estimate of near-trench slip and the other consists only of 12 uniform-slip segments. Displacements of GNSS sites that we calculated using the mean slip model are also shown in Figure 4, for comparison with observed displacements. In taking the mean slip magnitude, we ignored differences in fault geometry and slip rake between models. However, to calculate the GNSS site displacements, we mapped the mean slip magnitude to a three-dimensional (3-D) curved fault (contoured in Fig. 4) and assigned a uniform rake of 84° (thrust faulting). The finite element model for this calculation is the same as that of Sun et al. (2014) and Sun and Wang (2015), except here we used a rigidity of 35 MPa for the top 25 km of the upper plate and 64 MPa for the rest of the model.

On average, the coseismic slip featured by the published models is mostly far offshore with peak values within 100 km of the trench. However, there are substantial differences between models, and the differences become larger with distance offshore. Between latitudes 37°N and 39°N, their standard deviations are typically >10 m within 100 km of the trench and approaching 20 m at the trench (Fig. 5A). The large scatter near the trench is also shown in Figure 6A along a corridor through the main rupture area. At a given distance, the shown slip value of each model is an average of three profiles spaced ∼20 km apart. Some of the slip profiles in Figure 6A do not reach the trench, because their constant-strike fault models did not extend to the trench at this latitude; for map views of fault boundaries of most of the models, see the supporting information file of Brown et al. (2015).

Figure 6B shows fault-zone-average static stress drop of each model as a function of the averaging area encompassed by given coseismic slip contours in the same model. The method of stress drop calculation is as described by Brown et al. (2015). Because the calculation is performed in a uniform elastic half-space (rigidity 40 MPa) and involves mapping slip vectors of all of the published models to the same 3-D curved fault surface (Brown et al., 2015), the stress drop values shown here may be slightly different from those that would be calculated directly from the original models, but the difference is negligibly small if compared with the differences between the models. The average stress drop of the high-slip area encompassed by the 30 m slip contour is between 6 and 14 MPa for most models. The average stress drop over a much larger fault area encompassed by the 5 m slip contour is between 1.5 and 4 MPa for most models. The decreasing trend of the stress drop with increasing averaging area reflects the role of stress increase (negative stress drop) in different parts of the heterogeneous fault. The heterogeneity is considered to be of fundamental importance to earthquake mechanics (Brown et al., 2015).

Understanding Differences between Models

General Remarks

Given the limited data coverage over and near the rupture zone, rupture models for this earthquake will always be non-unique, and it is by no means surprising to see different researchers continue to produce very different results, even if they invert essentially the same data. This is true for any earthquake lacking near-field observations. In this regard, the lack of data itself is often less of a problem than the lack of recognition of model limitations due to the lack of data. Various resolution tests accompanying the published models commonly portray better accuracy than reflected by the differences between models and therefore can be misleading. Non-uniqueness due to limited data coverage can be overcome by collecting and/or incorporating critical data. For example, as will be briefly recapped in the Improving Near-Trench Resolution section, repeated bathymetry measurements have helped to resolve trench-breaching rupture in a local area (Fujiwara et al., 2011; Sun et al., 2017). However, some of the differences between models shown in Figures 5 and 6A must be due to different modeling practices.

Different types of data illuminate different aspects of the rupture process, so that the performance of a rupture model cannot be assessed by examining its predicted net slip alone. For example, discoveries about rupture dynamics such as the findings that high-frequency seismic wave energy is radiated mainly from the deeper part of the rupture zone (e.g., Lay et al., 2012) and that the rupture propagated alternately downdip and updip (Ide et al., 2011) need not at all be accompanied by an accurate prediction of the net slip distribution. The compilations of Figures 4, 5, and 6 thus do not describe the general state of knowledge of the rupture process of the Tohoku-oki earthquake, but only that of static coseismic slip. Because we focus on static deformation in this paper, here we only examine the net slip.

Some of the issues often discussed in the literature regarding the accuracy of slip inversion are of relatively minor importance. One example is material heterogeneity. If we use a uniform rigidity such as 40 MPa or any geologically acceptable value to calculate GNSS site displacements, we only need to scale the mean slip magnitude down by 10% in order to reach the same agreement with observations as shown in Figure 4. Differences between published models are much larger in both slip magnitude and slip distribution (Figs. 5 and 6A). The choice of inversion algorithm, especially with regard to Bayesian versus non-Bayesian approaches, is of similarly minor practical importance. All commonly used modern inversion algorithms seem to yield reasonable results, as long as they properly apply parameter constraints such as boundedness and smoothness. Some models fail to constrain rake directions in low-slip areas, but the impact on the main slip area is small.

Fault Geometry and Seafloor Slope

If the rupture zone is relatively deeply buried, the rupture area is small, and/or the slip is not very large, simplifications of fault geometry may not cause serious problems to inferring fault slip from surface geodetic observations. However, the Tohoku-oki earthquake does not meet any of these conditions, and the simplifications do cause problems.

The use of a planar fault, as in some of the early finite-fault models, causes the dip of the fault to be too high in the shallow part and possibly too low in the deep part of the model as compared to the actual fault dip. The wrong dip introduces errors to fault slip inferred from surface deformation or motion. The use of uniform-slip sub-faults of very large size, with dimensions much larger than the depth of the fault, is not merely a resolution problem; sudden and large changes of slip magnitude across the boundaries of these sub-faults give rise to deformation and stress singularities. If the spurious surface deformation produced by these singularities is used to fit near-field observations—either actual seafloor deformation or tsunami waves due to such deformation—large errors occur in the inferred fault slip.

The vast majority of the inversion models assume a uniform half space or layered Earth with a horizontal surface of no topographic relief (Fig. 7), as required by the derivation of fault slip via Green’s functions in these models. These models suffer a potentially serious problem in dealing with seafloor slope and fault dip near the trench. If a flat-top model uses the actual fault geometry and depth below sea level (Fig. 7A), it cannot simulate trench-breaching rupture. If seafloor deformation produced by such a model is mapped to a sloping seafloor to simulate tsunami generation, the tsunami waves will be modeled incorrectly. If the actual fault geometry is used but the depth at trench is set to zero in order to simulate trench-breaching rupture (Fig. 7B), fault depth below seafloor will be systematically underrepresented. The stiffness of the near-trench part of the upper plate in the model will be much too low to yield correct results.

A dip correction to the shallow part of the fault (Fig. 7C) is a better compromise (Wang et al., 2003) because the fault depth in the model approximately represents the actual fault depth below the sloping seafloor. This correction is important in calculating fault stress drop using a flat-top dislocation model (Brown et al., 2015). It leads to an overestimate of coseismic uplift of the seafloor near the trench that happens to compensate for the effect of the sloping seafloor in tsunami generation that is missing in a flat-top model (Fig. 7C). For tsunami generation, it is the vertical motion of the seafloor relative to seawater that is important (Lotto et al., 2017). If the seafloor slope angle is α, and fault dip near the trench is β, the rise of the near-trench seafloor beneath a fixed water column due to a trench-breaching slip s is u = ssin(α + β) / cosα (Fig. 7C). Here for simplicity we only discuss the rigid-body-translation component of the coseismic deformation. With the adjustment shown in Figure 7C, the surface slope becomes zero, the fault dip becomes approximately α + β, and the modeled seafloor rise u′ ≈ ssin(α + β) = ucosα. Because α is usually a small angle (<10°), cosα ≈ 1, so that u′ ≈ u (Fig. 7C). Even if α were as large as 30°, u′ would underrepresent u by only 13%.

The problem shown in Figure 7 is made acute by the relatively steep lower continental slope at the Japan Trench and the extraordinarily large shallow slip of the Tohoku-oki earthquake. In particular, if a buried rupture in some form of Figure 7A is used, it is not possible to handle the large shallow slip. If a numerical (such as finite element) model is used, the actual seafloor and fault geometry can be simulated, as in the forward models used to produce the results shown in Figure 4 and those to be discussed in the Improving Near-Trench Resolution section. However, only three of the 44 inversion models included in Figures 4, 5, and 6 (models 8 [Kyriakopoulos et al., 2013], 15 [Romano et al., 2014], and 46 [Freed et al., 2017] in Fig. 8), all based on the finite element method, incorporated actual seafloor slope. Among flat-top models, Ito et al. (2011b) (not included in the 44 models) used an approach similar to that illustrated in Figure 7C, and model 3 (Iinuma et al., 2012) in Figure 8 used a hybrid of those illustrated in Figures 7B and 7C (corrections made for fault depth but not fault dip). Some authors may have introduced more sophisticated solutions (e.g., Goshima and Miyazaki [2014], not included in the 44 models). It is not clear how most of the other flat-top models dealt with this problem. Errors incurred may have contributed to the large scatter of results near the trench.

Models Constrained by Seafloor GNSS

We select from Figure 6 a subset of models that included seafloor GNSS data in their inversion and show the subset in Figure 8. In this selection, we also require that at least two planar fault segments are used to accommodate the change of fault dip with depth, and therefore do not include models with a single planar fault. Index numbers assigned to these models are the same as in Sun et al. (2017, their supplementary table 1), so that the reader can readily check model attributes from that table. The newly added model by Freed et al. (2017) is numbered 46.

From the coast to ∼70 km from the trench, the scatter of the slip distribution in Figure 8A is much smaller than in Figure 6A, except for model 15 (reason unclear). For example, excluding model 15, the difference between the largest and smallest slip values at 75 km is ∼42 m in Figure 6A, but it is ∼26 m in Figure 8A. The improvement appears to owe mainly to the inclusion of seafloor GNSS observations. Accounting for the fault dip change plays a role, but three of the constant-strike models (blue) still deviate quite far from the rest of the group within ∼80 km of the coastline. Models that use 3-D realistic fault geometry (red) similar to that contoured in Figure 4 produce coherent results for this distance range, except for models 5 and 15. Therefore, along-strike curvature of the fault appears to be important; stations on land are influenced by integrated contribution from a large fault area.

Somewhat disconcerting and contrary to common expectation, the inclusion of tsunami data by some of these models (Fig. 8A, thicker lines) has not resulted in a better agreement between models near the trench. A similar plot not limited to models that included seafloor GNSS is shown in Sun et al. (2017, their supplementary figure 1) and conveys the same message. This message may not reflect negatively on the value of the tsunami data in helping to resolve shallow fault slip, but it may indicate that the data have not been properly used for this purpose. Fault models are likely oversimplified and/or too coarsely discretized for simulating tsunami wave generation. The problem illustrated in Figure 7 must have played a large role. It is also possible that the hydrodynamics of tsunami propagation is not always adequately incorporated in the inversion procedure.

Compared to Figure 6B, the stress drop results in Figure 8B show less scatter. There is some correlation between the fault-zone average behavior (Fig. 8B) and behavior along the central corridor (Fig. 8A). The model (23) that exhibits consistently higher stress drop than most other models did not include static land-based GNSS observations (Fig. 4) and poorly resolves near-coast slip. The two models (13 and 15) that yield the highest stress drop in the high-slip area (30 m contour) show slip distributions quite different from those of other models. Model 10 yields non-physical average stress increase over large fault areas (e.g., encompassed by the 5 m slip contour), which reflects stability issues in the inversion, and it shows larger slip than most other models near the coast.

Improving Near-Trench Resolution

The only way to improve near-trench resolution is to make near-trench observations. However, at the time of the earthquake, the easternmost GNSS site was still 50 km from the trench (Fig. 4). Pressure decrease, indicating uplift, recorded by a seafloor gauge some 20 km from the trench (Fig. 4B) (Ito et al., 2011b; Iinuma et al., 2012) was the only reliable direct near-trench geodetic measurement during this earthquake. The other near-trench measurements described by Ito et al. (2011b) are too uncertain to be quantitatively useful. Tsunami waves recorded nearby may contain valuable information, but extracting such information to constrain fault slip appears to be a challenging task at present, as discussed above.

An unusual set of observations that proves to be useful in constraining near-trench slip is the differential bathymetry obtained from surveys before and after the earthquake. The collection and processing of the observed differential bathymetry (ODB) have been described by Fujiwara et al. (2011), and the derivation of synthetic differential bathymetry (SDB) using a finite element model of realistic seafloor slope and fault geometry has been described by Sun et al. (2017). Given the widely present problem of using a flat-top model for a sloping seafloor, illustrated in Figure 7, and the unsatisfactory performance of tsunami inversion, the SDB modeling is the only work that has quantitatively addressed the trench-breaching slip of the Tohoku-oki earthquake in a self-consistent fashion, albeit only in a very local area.

The main results of Sun et al. (2017) are summarized in Figure 9. The ODB (Fig. 9A, right) is from surveys in 1999 and 2011 along a trench-normal profile covering the most seaward portion of the corridor used for Figures 6A and 8A and crossing the trench. The best fault slip model shown in the left panel of Figure 9B was obtained by minimizing the difference between the SDB (in the right panel of Fig. 9B) and ODB through grid searching in the parameter space as described by Sun et al. (2017). It features an average slip of 62 m over the most seaward 40 km of the fault, with the slip gently increasing toward the trench by 5 m over this distance. The SDB modeling helps to reject scenarios of large increase or decrease of slip toward the trench (Figs. 9C and 9D), which would indicate large net strengthening or weakening, respectively, of the shallow fault during the earthquake. By referring to analyses of core samples from a nearby drill site (Chester et al., 2013; Kirkpatrick et al., 2015), Sun et al. (2017) proposed that the finding is consistent with a model of delayed dynamic weakening.

There must have been significant along-strike variations in the shallow slip, but there is only one additional ODB profile in the main rupture area, located ∼50 km further north but with poorer data quality. The SDB modeling of that profile indicates smaller average slip (42 m) with a large increase toward the trench (by 20 m over 40 km) (Sun et al., 2017). The encouraging results of the SDB modeling demonstrate the value of unconventional seafloor geodetic observations.

On Coseismic Coastal Subsidence

Much of the Pacific coastal area of northern Honshu was observed to be subsiding (or titling seaward) before the 2011 earthquake, especially the area directly landward of the rupture zone (e.g., Suwa et al., 2006; Hashimoto et al., 2009). To some researchers, the observed coastal subsidence during the earthquake (Fig. 4B) was somewhat unexpected because there was an expectation that such a giant earthquake should have undone crustal deformation accrued during the interseismic period. This expectation is based on a misconception introduced by the wide application of the elastic interseismic locking model (see the Fault Locking and Interseismic Stress Relaxation section), which features interseismic deformation as a subdued mirror image of coseismic deformation. Given the viscoelastic Earth rheology, there is no reason why interseismic and coseismic deformation must be mirror images of each other (Wang et al., 2012).

However, deciphering vertical interseismic deformation is one of the outstanding issues in the study of megathrust earthquake processes (Wang and Tréhu, 2016). For active continental margins, processes responsible for straining the crust (horizontal deformation) are generally much better understood than processes responsible for tilting the crust (vertical deformation). As discussed by Wang and Tréhu (2016), non-tectonic and tectonic processes causing vertical deformation on time scales of years to millennia at subduction zones are poorly understood. At present, no theoretical and numerical models can satisfactorily and self-consistently explain vertical crustal deformation in subduction earthquake cycles. Making progress requires multidisciplinary observational and theoretical studies of Earth rheology and the various processes that cause vertical deformation.


Dominance of Viscoelastic Stress Relaxation

Following the Tohoku-oki earthquake, land-based GNSS sites continued to move seaward, in approximately the same direction as their coseismic motion. Shown in Figure 10 are one-year average velocities of these sites based on daily coordinate time series data from the Geospatial Information Authority of Japan following the analysis strategy defined by Nakagawa et al. (2009). The wholesale seaward motion is entirely as expected (Wang et al., 2012), because both viscoelastic relaxation of earthquake induced stresses and afterslip of the megathrust can cause this motion. But the partition of contribution from the two processes was by no means clear at first. Telling these two processes apart using surface deformation data had been a notoriously difficult problem (e.g., Wang, 2007; Bürgmann and Dresen, 2008). For subduction zone earthquakes, there did not seem to be even remotely relevant data to constrain it. However, the situation has now changed owing to one of the most remarkable discoveries after the earthquake—the observation of landward motion of a few seafloor GNSS sites in the main rupture area (Fig. 10A).

Within the first year after the earthquake, seafloor GNSS observations were made by Japan Coast Guard (JCG) and Tohoku University. The one-year average velocities of the seafloor GNSS sites surveyed by JCG shown in Figure 10 (labeled in Fig. 10C) are based on data reported by Watanabe et al. (2014); that of the GNSS site GJT3 surveyed by the Tohoku University (labeled in Fig. 10A) is based on data reported by Sun et al. (2014). Sites that provided the most critical information are KAMS, MYGI, and GJT3. The rather fast landward motion of this contiguous cluster of three sites immediately after the earthquake (Fig. 10A) cannot be explained by the relocking of the megathrust, because the rate of motion in the first year far exceeded the rate of plate convergence here. Neither can it be explained by afterslip, because afterslip would have caused seaward motion of these sites, as is obviously the case for site FUKU (Fig. 10). Therefore, it must be the effect of viscoelastic relaxation. Sun et al. (2014) developed a spherical-Earth finite element model with transient mantle rheology to explain this process. The model was subsequently updated by Sun and Wang (2015). Further updates are shown in Figure 10 and explained in the following subsection.

Sun et al. (2014) explained that the opposing motion of the near-trench area and the rest of the upper plate is the natural consequence of the rupture asymmetry of thrust earthquakes, with the hanging wall undergoing much greater coseismic displacement than the footwall. Sun and Wang (2015) demonstrated that, in a viscoelastic Earth, the opposing motion is a robust feature of postseismic deformation after large subduction earthquakes (Mw > 8). The dividing boundary between the landward and seaward motion roughly overlies the downdip termination of coseismic rupture. As the effect of viscoelastic relaxation diminishes with time, the effect of megathrust locking becomes more dominant and causes the dividing boundary to migrate landward, eventually leading to wholesale landward motion (Wang et al., 2012). Given geologically reasonable parameters and boundary conditions, changing modeling details will not do away with the opposing motion, although they may affect the rates of motion. For example, fitting the early part of the GNSS time series (Sun et al., 2014) instead of fitting only cumulative displacements over a longer time (e.g., Hu et al., 2016; Freed et al., 2017) necessitates the use of initially very low viscosity as represented by the transient rheology. As another example, including a thin low-viscosity layer beneath the subducting plate (model A of Sun et al., 2014) or not including it (model B of Sun et al., 2014) affects predicted rates of vertical seafloor motion.

Following a simple logic, Sun et al. (2014) reached an important conclusion. To explain the fast landward motion of the three seafloor sites, (initially) rather low mantle viscosity is needed. This would then lead to fast seaward motion of the land area, to the extent that much less afterslip is needed downdip of the main rupture to explain the observed motion of the land-based GNSS sites (Fig. 10). Therefore, seafloor geodesy has provided unambiguous evidence for the dominance of viscoelastic relaxation in short-term postseismic deformation. For great subduction earthquakes (Mw > 8), all elastic models overpredict afterslip downdip of the rupture and underpredict afterslip updip of the rupture if afterslip is present (Sun and Wang, 2015).

Modeling Five Years of Postseismic Deformation

The results shown in Figure 10 are obtained with an updated version of the viscoelastic finite element models of Sun et al. (2014) and Sun and Wang (2015). The model structure and rheological parameters are identical to those of the earlier versions. Major differences are as follows. (1) The earlier versions modeled the first three years of postseismic deformation. The results in Figure 10 are for a five-year time window constrained by additional, more recent GNSS observations. (2) The effect of megathrust relocking, very small in the first couple of years, was not included in the earlier versions. The updated model includes the effect of locking using the recipe of Wang et al. (2001) and Hu et al. (2004). (3) The afterslip distribution of Sun et al. (2014) was scaled and modified from that inferred by Ozawa et al. (2012) with an elastic model from 9 months of land-based GNSS observations. Sun and Wang (2015) added a shallow patch of afterslip south of the main rupture zone to explain the observed postseismic motion of site FUKU. In the updated model presented here, afterslip is assigned in accordance with more recent land and marine GNSS observations.

Model-predicted time series are compared with recent seafloor GNSS observations in Figure 11 (Japan Coast Guard, 2016; Tomita et al., 2017) and with observations from eight examples of land sites in Figure 12. One-year average velocity vectors derived from the observed time series for the 2–3 year time window, wherever possible, are shown in Figure 10B. The low signal-to-noise ratio of the data for the 4–5 year time window does not allow reliable derivation of the velocities. Our model is constrained only by the observed horizontal motion of the GNSS sites and compares with the vertical component very poorly. The vertical component is more sensitive to fine details of structural and rheological heterogeneity and may be better explained with a more complex viscosity structure (e.g., Muto et al., 2013; Freed et al., 2017).

Applying the same viscosity values as in the early versions of the model to the longer time window may introduce some uncertainties. As noted by Sun et al. (2014), in the three-year time window considered in their work, the mantle rheology probably was still in a transient phase, such that the low steady-state viscosity used in that model (1.8 × 1018 Pa·s) is probably not the true steady-state value. We expect the viscosity to continue to increase with time after the three years, eventually to the more widely seen values of ∼1019 Pa·s (Wang et al., 2012), but the evolution cannot be fully described by the simple bi-viscous rheology employed by our models. At present, we are not yet in a position to introduce more complex transient rheology. For this reason, uncertainties may become larger in the later part of the five-year model period.

Afterslip distribution undoubtedly evolves with time. If an elastic model is used, it is relatively straightforward to determine afterslip evolution by inverting time-dependent postseismic surface observations. For our forward modeling using a viscoelastic model, we take a simplified approach by assigning afterslip patches of various slip durations to match GNSS observations. Except for the irregularly shaped shallow afterslip patch south of the rupture zone (Fig. 10, patch 4), we use patches of elliptic shape with fixed perimeters. The spatial distribution of the slip magnitude within each elliptic patch is assumed to be symmetric with respect to the center and is assigned using the method of Wang et al. (2013a) with a broadness parameter 0.2. For the southern shallow patch, we obtain the slip distribution by adding slip values of several overlapping elliptical patches. The temporal evolution of slip values in each patch is assumed to obey the function chosen by Hu and Wang (2012). The temporal specifics of the six afterslip patches shown in Figure 10 are given in Table 1. The three small and short-duration patches (2, 5, and 6) are needed to explain relatively fast motion of nearby GNSS sites immediately after the earthquake, as a crude way of reflecting the complex spatiotemporal afterslip evolution in these areas. The model-predicted time series in Figures 11 and 12 include the effects of both viscoelastic relaxation and afterslip.

To explain the motion of seafloor sites MYGI and MYGW (Fig. 11), we assigned some small, very slowly evolving afterslip within the deeper part of the main rupture zone (Fig. 10, patch 3). Because of the newly added locking effect, the predicted dividing boundary between landward and seaward motions with the viscoelastic effect alone is located slightly more westward than in the earlier versions of the model. Patch 3 is introduced so that the dividing boundary is still between MYGI and MYGW. There are some uncertainties in this assignment, because the model behavior within the rupture zone is strongly affected by uncertainties in the coseismic slip model. If this afterslip is real, it serves to make up partially for the smaller coseismic slip in the lower part of the rupture zone.

An important “side product” of this modeling exercise is the finding that the megathrust cannot be uniformly locked along strike after the earthquake. The not fully locked segments in Figure 10 have been assigned 25% locking, i.e., creeping at 75% of the subduction rate. Detailed explanation of this finding and its implications will be presented elsewhere. Here we only highlight a few salient points. Uniform locking will cause model-predicted motion of land GNSS sites north and south of the earthquake area to be too slow. The creeping rate of the not fully locked areas cannot be precisely determined because of tradeoff with other model parameters; e.g., zero or 50% locking are both possible. In the strike direction, the locking pattern shown in Figure 10 is in general agreement with the interseismic locking patterns published prior to the Tohoku-oki earthquake partially illustrated in Figure 2. The creeping behavior south of the 2011 rupture zone is consistent with active seamount subduction in that area (Wang and Bilek, 2014).

Spatial Characteristics of Afterslip Distribution

The determination of deep afterslip for the Tohoku-oki earthquake will never be accurate because of the large distance of surface observations from the relevant fault area (typically >40 km for this earthquake) and tradeoff with other parameters such as mantle viscosity and plate thickness. The determination of near-trench shallow afterslip cannot be accurate at present because of the paucity of near-field observations. However, our numerous trial-and-error tests and comparison with other viscoelastic deformation models for this earthquake suggest the following three robust characteristics of afterslip distribution (Fig. 10). (1) There is no large deep afterslip directly downdip of the main rupture area. (2) There is substantial afterslip northwest of the main rupture (patches 1 and 2). (3) There is large near-trench shallow afterslip south of the main rupture (patches 4 and 5).

Other recent models that invoke viscoelasticity and 3-D slab geometry also converge on the three robust characteristics (Hu et al., 2016; Iinuma et al., 2016; Freed et al., 2017). Note that the model of Hu et al. (2016) differs from other models by using a thin viscoelastic layer to simulate afterslip. They constrained the amount of afterslip shallower than 50 km using repeating earthquakes reported by Uchida and Matsuzawa (2013). Their model also features some afterslip in the 50–100 km depth range, but this is equivalent to the viscous shear in the bottom part of the mantle wedge in the models of Sun et al. (2014), a semantic issue we will not further discuss.

The southern shallow afterslip (characteristic 3 above) must be present if seafloor measurements in that area are considered (Figs. 10, 11). The northern deep afterslip (characteristic 2) appears in all afterslip models for this earthquake, as listed in Table 2, although not as a distinctly standalone patch in Diao et al. (2014). It is needed to explain the fast motion of the land GNSS sites at this latitude. It is also consistent with the occurrence of small repeating earthquakes in this area (Hu et al., 2016). Near-trench seafloor GNSS observations at this latitude do not require the slip to extend very far updip (Figs. 10, 11). The amounts of the northern afterslip in viscoelastic models are systematically smaller than in elastic models if the results are truncated or extrapolated to a similar time window (Table 2). The validity of characteristic 1 summarized above, that is, the absence of significant deep afterslip directly downdip of the main rupture (the central patch in Table 2), is perceived to be controversial and will be discussed in the Deep Afterslip Disagreement section.

It is illuminating to compare these characteristics with the patterns of coseismic rupture and deep pre-slip (Fig. 13). It is not surprising that afterslip tends to be absent in the area of primary coseismic slip (Pritchard and Simons, 2006). Neither is it surprising that large shallow afterslip occurs south of the rupture zone. It has been an area of creep with smaller earthquakes, to some degree associated with seamount subduction (Wang and Bilek, 2014), and the creep is expected to accelerate in response to the giant earthquake. The lack of significant deep afterslip directly downdip of the main rupture (characteristic 1) is surprising, but it is made less surprising by the complimentary spatial distribution of the deep afterslip and deep pre-slip seen in Figure 13. What physical processes cause this type of slow slip to occur before or after a large earthquake invites exciting new scientific research. One cannot help but to further speculate that the large afterslip northwest of the rupture zone (characteristic 2) may have some relation with the enhanced locking before the earthquake (Fig. 13, blue), again inviting new research.

The Deep Afterslip Disagreement

Postseismic deformation models that addressed deep afterslip at 30–50 km depths are listed in Table 2. Although they all feature a northern area of large afterslip, they differ in the afterslip behavior directly downdip of the main rupture, called the central patch in Table 2. In contrast with Figure 10, all of the elastic models show a central patch of large afterslip. These elastic models were all published in 2014 or earlier. Not guided by seafloor observations reported in Watanabe et al. (2014) and Sun et al. (2014) to invoke viscoelasticity, these models had to use large deep afterslip to explain the seaward motion of the land GNSS sites.

What seems intriguing is that the viscoelastic models of Diao et al. (2014) and Yamagiwa et al. (2015) also show large afterslip directly downdip of the main rupture, comparable to the elastic models but different from the other viscoelastic models (Table 2). Because the model of Yamagiwa et al. (2015) can be considered an improved version of that of Diao et al. (2014) owing to a similar methodology and model setup, here we only use the former to explain their disagreement with other viscoelastic models. Different models in Table 2 used different time windows. To ease comparison with other models, we show values inferred by Yamagiwa et al. (2015) in three different time windows. For example, the time window of 1.6–9.6 months after the earthquake is comparable with the time window of 1.4–9 months used by Iinuma et al. (2016). Both used viscoelastic models, but Iinuma et al. (2016) obtained much less afterslip for the central patch.

The key question appears to be whether a subducting slab is present in the model. Yamagiwa et al. (2015) used a layered Earth model, with an elastic plate on top of a Maxwell viscoelastic mantle without a subducting slab. As explained by Miyashita (1987), the slab has a first-order control on mantle flow during viscoelastic relaxation and thus on surface deformation. Using numerical tests, Sun and Wang (2015) demonstrated that a slab acts like an anchor in the deep mantle and retards the landward motion of the trench area following a subduction earthquake. They show that a thinner slab leads to faster landward motion of the trench area. In a layered Earth model (i.e., zero slab thickness), the effect goes extreme and enhances the landward motion of the near-trench area. The rate of the landward motion of the near-trench area predicted by Yamagiwa et al. (2015) fits the seafloor GNSS observations because a mantle viscosity value as high as 0.9 × 1019 Pa·s is used, higher than the transient viscosity values used by Sun et al. (2014), Sun and Wang (2015), Hu et al. (2016), and Iinuma et al. (2016) or the asthenospheric viscosity of the Maxwell model of Freed et al. (2017) by one to two orders of magnitude. Such a high viscosity would lead to very slow seaward motion of land GNSS sites. To compensate for this, large deep afterslip must be introduced.

If realistic subduction structure is used, with a subducting slab, the afterslip disagreement is resolved. All viscoelastic models including a slab, such as the last five listed in Table 2, then need, at least initially, a very low viscosity and predict little deep afterslip directly downdip of the main rupture, regardless of whether they use a bi-viscous or Maxwell rheology and whether they fit the GNSS time series or cumulative displacements. For the same reason, they also predict less deep afterslip to the north, if the results are approximately projected to a similar time window (Table 2).


Upon examining relevant publications and results of our own research on the observation and modeling of static crustal deformation associated with the Mw 9 2011 Tohoku-oki earthquake, we conclude that we have learned the following about subduction earthquakes.

Before the Earthquake

  • In the strike direction, Japan Trench interplate locking models published prior to the Tohoku-oki earthquake correctly portrayed the future rupture zone as a segment of full locking. They demonstrate that land-based geodetic observations are capable of resolving along-strike variations in offshore megathrust locking and creep at wavelengths comparable to distances from the network.

  • In the dip direction, locking models published prior to the earthquake incorrectly portrayed a locked patch that is much deeper than the actual rupture in 2011. Besides limited offshore resolution of the land-based geodetic data, one primary reason is the neglect or inadequate incorporation of the effects of interseismic viscoelastic stress relaxation.

  • Land-based geodetic observations offered no information on the locking state of the shallowest part of the Japan Trench megathrust. Full-speed creep of the near-trench part of the fault portrayed in locking models prior to the earthquake was a result of model assumption. Confusing mechanical states such as strong or weak with kinematic states such as locking or creeping is a common mistake in deriving interseismic locking models from geodetic observations.

  • The discovery of decade-long accelerated fault creep downdip of the future rupture zone raises new questions about how giant earthquakes are incubated and initiated. It is not known whether it was a one-time pre-slip as part of the rupture initiation or one of many similar events that had occurred episodically. The discovery motivates further theoretical and experimental investigation.

During the Earthquake

  • Seven seafloor GNSS sites yielded the most remarkable measurements of coseismic displacements ever seen, up to 31 m. Coseismic slip models that used these measurements as constraints show large slip peaking within 100 km of trench.

  • Different types of data may have illuminated different aspects of the rupture process, but models of net slip that included the seafloor GNSS data and realistic fault geometry tend to show better mutual agreement than those that did not. Differences between models near the trench are unavoidable because of the general absence of near-trench observations of static displacements. Differences near the coast indicate needs for model improvements.

  • Static stress drop, if averaged over the fault area encompassed by the 5 m coseismic slip contour, is between 1.5 and 4 MPa for most rupture models that included seafloor GNSS and accounted for fault dip change with depth. Stress drop averaged over the high-slip area encompassed by the 30 m slip contour is >10 MPa in most of these models.

  • Most rupture models employed a flat top. Not properly accounting for the seafloor slope above the rupture zone causes difficulties in simulating trench-breaching rupture and is inferred to be responsible for some of the near-trench differences between published models.

  • Lack of near-trench agreement between models that included tsunami data as constraints suggest that much effort is needed in properly extracting the valuable information on shallow fault slip in these data. Oversimplification of fault and seafloor geometry is one obvious reason for the divergent results, although there may be other factors that require investigation.

  • Bathymetry differences before and after the earthquake provided the only direct, although unconventional, coseismic deformation measurements across the trench. Modeling these differences along a central profile in the main rupture area indicates >60 m of slip of the most seaward 40 km part of the fault, with only a gentle increase toward the trench.

  • The observed coastal subsidence during the earthquake following decades of interseismic coastal subsidence was unexpected to some researchers, because of the expectation that interseismic and coseismic deformation should be mirror images of each other. Although it is clear that this expectation is invalid in a viscoelastic Earth, correctly modeling vertical deformation throughout earthquake cycles remains an unresolved issue at present.

After the Earthquake

  • Fast landward motion of a cluster of three seafloor GNSS sites located in the main rupture area, opposing the seaward motion of all of the other land-based and seafloor sites, provided unambiguous evidence for the dominance of viscoelastic stress relaxation in short-term postseismic deformation.

  • Models that properly account for viscoelastic relaxation suggest three robust characteristics of afterslip distribution: lack of significant deep afterslip directly downdip of the main rupture, moderate deep afterslip slightly north of the main rupture, and large shallow afterslip south of the main rupture.

  • If results are approximately projected to similar time windows, afterslip models that do not include viscoelastic relaxation and/or do not include a subducting slab predict much larger deep afterslip, especially directly downdip of the main rupture. Lessons learned have important implications to modeling afterslip following great earthquakes (M >8) in other subduction zones.

  • The distributions of deep pre-slip and deep afterslip appear to be spatially complementary. The physical processes responsible for the complimentary distributions deserve new scientific research and may hold keys to important fault mechanics.


We thank the many researchers who sent us digital values of their coseismic slip models for the Tohoku-oki earthquake, and K. Koketsu who sent us the GNSS data shown in Figure 3. The paper benefitted from comments made by Harold Tobin and an anonymous reviewer. T.S. was supported by a University of Victoria Ph.D. Fellowship, an Alexander and Helen Stafford MacCathy Muir Graduate Scholarship, a Bob Wright Graduate Scholarship, and a Discovery Grant to K.W. from Natural Sciences and Engineering Research Council of Canada. R.H. and S.K. acknowledge Japan Society for the Promotion of Science KAKENHI grant JP26000002 supporting marine geophysical research in the Japan Trench. This is Geological Survey of Canada contribution 20170258.

Science Editor: Shanaka de Silva
Guest Associate Editor: David W. Scholl
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.