Co- and postseismic slip behaviors extracted from decadal seafloor geodesy after the 2011 Tohoku-oki earthquake

Investigations of the co- and postseismic processes of the 2011 Tohoku-oki earthquake provide essential information on the seismic cycle in the Japan Trench. Although almost all of the source region lies beneath the seafloor, recent seafloor geophysical instruments have enabled to detect the near-field signals of both the coseismic rupture and the postseismic stress relaxation phenomena. Annual-scale seafloor geodesy contributed to refining the postseismic deformation models, specifically to the incorporation of viscoelastic effects. However, because of the insufficiency in the spatial coverage and observation period of seafloor geodetic observations, no consensus on crustal deformation models has been reached, especially on the along-strike extent of the main rupture, even for the coseismic process. To decompose the postseismic transient processes in and around the source region, i.e., viscoelastic relaxation and afterslip, long-term postseismic geodetic observations on the seafloor play an essential role. Here, from decadal seafloor geodetic data, we provide empirical evidence for offshore aseismic afterslip on the rupture edges that had almost decayed within 2–3 year. The afterslip regions are considered to have stopped the north–south rupture propagation due to their velocity strengthening frictional properties. In the southern source region (~ 37° N), despite not being resolved by coseismic geodetic data, shallow tsunamigenic slip near the trench is inferred from postseismic seafloor geodesy as a subsequent viscoelastic deformation causing persistent seafloor subsidence at a geodetic site off-Fukushima. After a decade from the earthquake, the long-term viscoelastic relaxation process in the oceanic asthenosphere is currently in progress and is still dominant not only in the rupture area, but also in the off-Fukushima region.


Introduction
In general, a large fault rupture is followed by postseismic relaxation processes such as viscoelastic relaxation in the asthenosphere and aseismic slip on the fault plane, which lead to transient crustal deformation on the solid Earth's surface (e.g., Wang et al. 2012). Postseismic geodetic data following a megathrust earthquake show the sum of these relaxation processes and deformation caused by interplate coupling due to the secular plate subduction. Along the Japan Trench, postseismic processes were triggered by the 2011 Tohoku-oki earthquake (M w 9.0) and continue in this decade. Clarifying the interplate slip behaviors for the co-and postseismic phase will contribute to the understanding of the frictional state of faults, slow earthquake activities, and seismic cycles in this region.
The Tohoku-oki earthquake caused trench-ward seafloor displacements of several tens of meters Kido et al. 2011), reaching about 50 m at the trench ). Using seafloor geodetic data Kido et al. 2011), an extremely large slip was Open Access *Correspondence: s-watanabe@jodc.go.jp 1 Hydrographic and Oceanographic Department, Japan Coast Guard, 3-1-1, Kasumigaseki, Chiyoda-ku, Tokyo 100-8932, Japan Full list of author information is available at the end of the article estimated at the plate interface shallower than the hypocenter ( Fig. 1) (e.g., Ozawa et al. 2012;Iinuma et al. 2012). Except for the north-south extent of the shallow rupture, which could not be resolved by the geodetic data at that time, a consensus has been reached on the north-south rupture propagation at depths near the hypocenter, where almost all of the fault models produce similar results (summarized by Sun et al. 2017;Wang et al. 2018;Lay 2018;Uchida and Bürgmann 2021). This implies that the rupture in 2011 did not progress to the northern region (> 39° N), even though M w 8 earthquakes have historically occurred in this region ( Fig. 1) (Nagai et al. 2001;Yamanaka and Kikuchi 2004). Investigations of the postseismic behaviors, including the occurrence of afterslip, in the northern and southern regions outside the main rupture area provide useful information on how rupture propagation may have been restrained in a compact region in the depths near the hypocenter.
For the detection of transient postseismic crustal deformation to decompose the elementary processes, sufficiently long-term, high-frequency, and well-distributed geodetic data are required because these sources have different decay times and deformation patterns (e.g., Wang et al. 2012). Although the terrestrial Global Navigation Satellite System (GNSS) observation network has an extremely high spatiotemporal resolution, it cannot easily decompose the transient processes because the two processes of interest cause similar trench-ward movements in the onshore regions and thus cannot be distinguished from one another. In contrast, the viscoelastic relaxation in the oceanic asthenosphere and the afterslip are expected to cause displacements in opposite directions on the seafloor above the main rupture, i.e., landward and trench-ward directions, respectively (e.g., Hu et al. 2016). Therefore, seafloor geodetic observations can be used to decompose transient deformation sources despite their lower temporal resolution compared to that of terrestrial observations. Actually, seafloor geodetic technique detected postseismic landward movements larger than the subduction rate (~ 9 cm/year) (Argus et al. 2011) above the main rupture area, whereas terrestrial geodetic sites showed trench-ward movements. This is a conclusive evidence for the dominance of viscoelastic relaxation in the oceanic asthenosphere below the subducting plate on the postseismic deformation in the main rupture area (Watanabe et al. 2014;Sun et al. 2014;Hu et al. 2016;Tomita et al. 2017;Honsho et al. 2019).
The seafloor observation results stimulated researchers to develop postseismic deformation models incorporating the viscoelastic effects. Table 1 summarizes the postseismic models (1) referencing the seafloor geodetic data; (2) incorporating the afterslip models; (3) modeling the subducting cold slab, and (4) calculating the viscoelastic deformation in an area wider than the latitude range of 37-39° N. Every model indicated that the deformation patterns in the main rupture area can be roughly explained as the viscoelastic response, but that the observed deformation outside of main rupture cannot be reproduced only by viscoelastic relaxation. Hu et al. (2016) proposed a model where the afterslip effects are approximated as a viscoelastic flow in a thin layer whose locations are constrained by repeating earthquakes, which was yet insufficient to explain the vertical motions of the seafloor. Some models put additional afterslip patches in the offshore region to reproduce the geodetic data (e.g., Sun and Wang 2015;Freed et al. 2017;Wang et al. 2018), even though these were only tentative models because of the insufficient spatiotemporal resolution and observation period (at most 5 years) of available seafloor data.
On the other hand, there remains the uncertainty in some coseismic dislocation input for the postseismic deformation modeling. Many researchers had adopted coseismic slip distribution models inverted from seafloor and terrestrial coseismic geodetic data with constraints of smooth slip distribution (e.g., Iinuma et al. 2012;Freed et al. 2017). However, the geodetic network in 2011 did not cover the whole source region, lacking information especially on the coseismic rupture in the shallower portion. Thus, it is important to consider seismic and tsunami waveforms, together with geodetic data (e.g., Lay 2018; Uchida and Bürgmann 2021). To visualize (See figure on next page.) Fig. 1 Seafloor motion derived from GNSS-A observations. Average horizontal and vertical velocities for a preseismic period, and three-year cumulative displacements for b 2011-2014, c 2014-2017, and d 2017-2020 are indicated as black arrows and open rectangles, respectively. Selected terrestrial GNSS velocities or cumulative displacements were extracted from the F3 solution of the GEONET sites (Nakagawa et al. 2009). Onshore velocities during 2007-2011 are shown in (a). Purple patches indicate possible afterslip regions, but there is little or no resolution for their spatial spread. Brown contours and colored tiles indicate geodetically derived (Iinuma et al. 2012) and tsunami-derived (Satake et al. 2013) coseismic slip distributions of the 2011 Tohoku-oki earthquake, respectively. Green circles denote repeating earthquakes that occurred in each period (Igarashi 2020). CMT solutions for the use of displacement correction are shown in each panel. Navy and blue lines indicate 2-m and 4-m slip contours of historical earthquakes in the northern area (Nagai et al. 2001). Green rectangles indicate patches with a slip of 20 m for the 1896 tsunami earthquake (Satake et al. 2017)  the difference, Fig. 1 displays both the geodetically and tsunami-derived slip distributions (Iinuma et al. 2012 andSatake et al. 2013, respectively). Tsunami data, which are sensitive to topographical changes of the seafloor, indicate that the tsunamigenic area was extended, especially in the northern (> 39° N) and southern (~ 37° N) areas along the trench ( Fig. 1) (Satake et al. 2013). This feature does not appear in geodetic inversion due to the absence of data. Such difference in the coseismic input would affect the coseismic stress-change field and the postseismic relaxation processes. For instance, Agata et al. (2019) proposed a viscoelastic model using coseismic slip distribution derived from seismic cycle simulation (Nakata et al. 2016). The coseismic input included an additional shallow slip near the trench eastern off-Fukushima (~ 37° N). Their model showed significant seafloor deformation in the off-Fukushima region caused by viscoelastic relaxation, whereas other models based on geodetic coseismic input could not induce enough stress to cause such significant viscoelastic deformation in this region (e.g., Sun et al. 2014;Freed et al. 2017). The difference between the models suggests that we will be able to estimate shallow coseismic slip behavior and following postseismic models by decomposing and discussing the postseismic deformation sources with longer-term seafloor geodetic data, as discussed in Tomita et al. (2020). Some researchers took a different approach to estimate coseismic distribution by constraining their models with postseismic data (e.g., Yamagiwa et al. 2015;Tomita et al. 2020;Fukuda and Johnson 2021).
In this study, decadal seafloor geodetic data that contain the temporal evolutions of surface velocity are used to decompose the deformation sources. Based on the results, we clarify the co-and postseismic slip behaviors in the northern and southern parts of the source region.

Data and methods
To investigate the temporal evolution of seafloor crustal deformation, the Japan Coast Guard (JCG) regularly performs seafloor geodetic observations using GNSSacoustic ranging combined seafloor positioning system (GNSS-A) (Additional file 1: Figure S1), in the Japan Trench region (e.g., . GNSS-A data were obtained at six JCG sites (KAMN, KAMS, MYGI, MYGW, FUKU, and CHOS; Additional file 1: Table S1) from March 2011 to June 2020 using survey vessels, and were processed with GARPOS v1.0.0 (Watanabe et al. 2021d). Note that the present analysis method, which incorporates the estimation process of horizontal gradients in the sound speed structure, has been applied to the previously published data (Watanabe et al. 2014). The JCG also performed GNSS-A observations at five sites (G08, G10, G12, G14, and G17; Additional file 1: Table S1) installed by Tohoku University (TU) (Kido et al. 2015) since 2013, independently of TU. We additionally reprocessed the GNSS-A data before the Tohoku-oki earthquake at five JCG sites (Sato et al. 2013a) using the present analysis method to determine the preseismic seafloor velocities.

Table 1 Summary of previously proposed postseismic deformation models
Models with (1) referencing the seafloor geodetic data, (2) incorporating the afterslip models, (3) modeling the subducting cold slab, and (4)  All GNSS-A data used in this study are available at Zenodo (Watanabe et al. 2021a, b). Note that the GNSS-A data before 2009 were obtained by the old "drifting observation" system, where the ship engine had to be turned off during the observation; thus, the ship track cannot be controlled (Additional file 1: Figure S1a). The data obtained by the "drifting observation" system are less precise than the "sailing observation" system, where the ship moves along a planned track (Sato et al. 2013b;Ishikawa et al. 2020).
To extract the annual-scale velocity changes, we took the following steps: We first subtract the effects of aftershocks that can cause coseismic displacement of more than 1 cm at seafloor sites (sources are shown in Fig. 1) by applying the method of Okada (1992). We used the Centroid Moment Tensor (CMT) solution catalogue provided by the Japan Meteorological Agency, which are available online (https:// www. data. jma. go. jp/ svd/ eqev/ data/ bulle tin/ index_e. html). We then smoothen the time series of postseismic displacements, x(t) , where t denotes the time after the event, by fitting with a function which is modified from the fitting function of Tobita (2016), i.e., where x 0 , v , α 1 and α 2 are estimation parameters. For the time constants τ 1 and τ 2 , we applied the values of Fujiwara et al. (2021), i.e., τ 1 = 2.1176 day and τ 2 = 287.45 day , because the GNSS-A observation frequency is as low as several times per year per site which is insufficient to determine these parameters. Additionally, we approximated the exponential term in their formulation to a linear component, vt , since an extremely large time constant was found by Fujiwara et al. (2021). It should also be noted that we put α 1 = 0 for the TU sites where the observation started in 2013.

Results
Time series of post-and preseismic seafloor displacements with respect to the Okhotsk plate of NNR-MOR-VEL56 model (Argus et al. 2011) within a framework of the International Terrestrial Reference Frame 2014 (Altamimi et al. 2016(Altamimi et al. , 2017 are shown in Figs. 2 and S2 (Additional file 1), respectively. The values of the displacement, without corrections for aftershocks, are available at Zenodo (Watanabe et al. 2021c). To discuss the motion changes over several years, we extracted 3-year cumulative movements from the fitted curves for the periods of Apr. -Apr. 2014, Apr. 2014-Apr. 2017, and Apr. 2017-Apr. 2020 (Table 2, Fig. 1). For the preseismic period, we only consider the average velocity (Table 2, Fig. 1a). (1) The GNSS-A results at the TU sites independently observed by the TU research group (Kido et al. 2015;Honsho et al. 2019) are simultaneously plotted in Fig. 2g-k for comparison. Because the absolute positions estimated in the different GNSS-A systems tend to be biased due to the differences in the GNSS's reference or the assumed sound speed model, there are offsets in the results between two observation systems, i.e., JCG and TU. The offsets were estimated and corrected as follows: we first estimated the linear trend of JCG's results in the period of 2013-2017 as a reference. Offsets were calculated from the average of the residuals of the TU's time series for 2013-2017 relative to the JCG's trend. Text S1 (Additional file 1) describes the validations for data at G17, where the positioning precision is lower than the other sites.

Discussion
Based on the temporal changes of seafloor movement (Fig. 1), we discuss the expressions of viscoelastic relaxation and afterslip in the northern (off-Kamaishi), central (off-Miyagi/main rupture), and southern (off-Fukushima and off-Choshi) parts of the source region, which are conceptually illustrated in Fig. 3.
In previous studies that analyzed the data until 2014 (Watanabe et al. 2014;Sun et al. 2014), landward movements at rates larger than the Pacific plate subduction were detected at the sites located above the main rupture, i.e., KAMS and MYGI (Fig. 1b). This indicated the dominance of the effects of viscoelastic relaxation in the oceanic asthenosphere beneath the Pacific plate, and the interplate backslip, assuming that the interplate coupling was restored (Fig. 3b). Crustal deformation in this area was consistent with the quantitative models incorporating the viscoelastic response to the geodetically constrained coseismic input (Sun et al. 2014;Sun and Wang 2015;Freed et al. 2017;Wang et al. 2018). The large landward movements at these sites continued with a slight decay over the whole period, as well as at G08, G10, G12, and G14 (Fig. 1c). The decay of the landward motion can be explained as the time-dependent viscoelastic deformation (Fig. 3b).
Little temporal change in the present decade was found in the horizontal movement at MYGW near the downdip edge of the main rupture ( Fig. 1b-d). If the interplate coupling in the main rupture had been restored, its landward motion of about 4-5 cm/year should be canceled by a trench-ward motion driven by viscoelastic relaxation in the mantle wedge (e.g., Fig. 4a of Hu et al. 2016) or afterslip to maintain balance for almost 10 years. Although the landward motion cannot be clearly detected, it seems to have been slightly restored after 2017. This might indicate a decrease in the dominance of relaxation processes similar to those in the main rupture area (Fig. 3b).
At the northern edge of the main rupture, little horizontal displacement was observed at KAMN until 2014 (Fig. 1b). Because viscoelastic relaxation is mainly driven by the stress induced in the low-viscosity layer beneath the lithosphere, it tends to cause almost the same movements at KAMN and KAMS which are only 30 km apart. To explain the velocity contrast at KAMN and KAMS, some postseismic models require afterslip to reach the trench on the northern side of the main rupture to cause a relative trench-ward motion Table 2 Displacements with respect to the Okhotsk plate from the fitted curves with variance-covariance information at KAMN with respect to KAMS (Sun and Wang 2015;Freed et al. 2017).

Site Displacement rate (cm/year) Variance-covariance ((cm/year) 2 ) E-ward N-ward U-ward V(E,E) V(E,N) V(N,N) V(U,U)
After 2014, landward motion significantly accelerated at KAMN and had almost the same velocity at KAMS and MYGI (Table 2c and Fig. 1c). With the preseismic velocities (about 3-4 cm/year at both KAMN and KAMS though with large uncertainty; Table 2a and Fig. 1a) taken into account, the consistency in the movements at KAMN and KAMS after 2014 indicates that the two sites have been similarly influenced by long-term viscoelastic relaxation. Because the spatial pattern of the viscoelastic deformation has not significantly changed in the present decade, the viscoelastic relaxation is expected to have caused almost the same displacements at these sites  (Iinuma et al. 2012) and tsunami-derived (Satake et al. 2013) coseismic slip areas of the 2011 Tohoku-oki earthquake, respectively. Blue circles and bars indicate the locations of GNSS-A sites and their projections to the plate interface, respectively. Black rectangles show the locations for cross sections illustrated in subsequent panels. Note that little or no information for the spread of afterslip regions can be added by the present GNSS-A data. b Cross sections for each region with one of possible qualitative explanations for the contributions of each deformation source to the motion at GNSS-A sites. Black, red, purple, and yellow arrows indicate the observed motion and components due to viscoelastic relaxation, afterslip, and interplate coupling, respectively before 2014. This supports that the afterslip in the off-Kamaishi region actually caused the relative trench-ward motion of 10 cm/year in three-year average at KAMN before 2014 to cancel the landward motion driven by the viscoelastic relaxation and interseismic backslip (Fig. 3b). Furthermore, the temporal evolution of KAMN's movement (Fig. 2a) confirms that the afterslip in that region had decayed sufficiently in 3 years. After the decay of the afterslip, the viscoelastic response and interplate coupling became dominant for the crustal deformation around KAMN, similar to the case at neighboring GNSS-A sites (Fig. 3b).
To reproduce the difference in the average displacement rate of about 10 cm/year between KAMN and KAMS before 2014, afterslip in the northern region outside the coseismic rupture with an average displacement rate on the order of meters per year is required (see Additional file 1: Text S2 for detail). However, the slip magnitude depends on the afterslip distribution, which cannot be geodetically constrained. Taking into account the tsunami-derived shallow coseismic slip in the off-Kamaishi area, the afterslip would not be expected to have reached the trench at 39° N (Fig. 1b).
In contrast to the off-Miyagi region, rapid trenchward movements were observed at FUKU and CHOS in the southern region (< 37.5° N) especially in the first 1-2 years after the Tohoku-oki earthquake (Fig. 2e, f ). The trench-ward motion became much smaller after 2013. Almost no significant horizontal movement was found at G17, located on the trench side of FUKU, despite the low positioning accuracy due to instrumental malfunction (< 5-6 cm/year for three-year average displacement rate; see Additional file 1: Text S1 for detail).
For the horizontal movement, it is reasonable to assume that the afterslip caused the trench-ward motion with rapid decay of 1-2 years, as shown in most of the postseismic models (Sun and Wang 2015;Iinuma et al. 2016;Freed et al. 2017;Agata et al. 2019). Actually, it had been indicated that viscoelastic deformation cannot cause significant motion at FUKU in some of the finite element models (Sun et al. 2014;Freed et al. 2017), which used the coseismic input based on the geodetic inversion with a single peak beneath FUKU, such as the model by Iinuma et al. (2012) (Fig. 1b). Therefore, such models require to reproduce both the horizontal and vertical motion at FUKU by only the afterslip, leading to the assumption that strong afterslip in the shallow portion is necessary to explain the rapid subsidence observed at FUKU.
However, the subsidence at FUKU continued at almost a constant rate of about 4 cm/year, even after 2014 when the trench-ward movement had almost ceased. This result conflicts with models which try to explain both the trench-ward motion and the subsidence with shallow afterslip. A single afterslip cannot cause such persistent subsidence without significant trench-ward movement at FUKU and G17, because of the low dip angle of the plate boundary (see Additional file 1: Text S3 for detail). Therefore, we should assume another input for the viscoelastic deformation rather than the geodetic inversion model to reproduce the persistent subsidence at FUKU. Agata et al. (2019) used a coseismic slip model derived from an earthquake cycle simulation (Nakata et al. 2016), which incorporates an additional peak of coseismic rupture near the trench at 37° N. Although the source model has not reproduced the coseismic seafloor uplift observed at FUKU , it can provide a practical example for the viscoelastic deformation pattern associated with the shallow coseismic slip. Their viscoelastic finite element model demonstrated sufficient subsidence at FUKU. According to the comparisons of the numerical examples of viscoelastic deformation to the coseismic slip distribution, as illustrated in Fig. 6 of Sun and Wang (2015), the hinge line of horizontal deformation and the peak of viscoelastic subsidence were located above the downdip side of the major slip. For this reason, the coseismic slip near the trench caused the subsidence at FUKU in the model of Agata et al. (2019). The tsunamiderived coseismic slip distribution in the off-Fukushima region (Satake et al. 2013) has two peaks in the along-dip direction (Fig. 1b). Based on the postseismic model of Agata et al. (2019), the viscoelastic relaxation driven by the shallower coseismic slip can cause long-term subsidence at FUKU (red arrows in Fig. 3b).
The discussion above suggests that the coseismic models that include constraints from tsunami observations (e.g., Yokota et al. 2011;Minson et al. 2014;Romano et al. 2014;Bletery et al. 2014;Melgar and Bock 2015) should be adopted in the viscoelastic relaxation model to reasonably explain the spatiotemporal variation of seafloor deformation. In this case, both the viscoelastic relaxation and possible interseismic backslip are predicted to cause landward motion at G17 (red and yellow arrows in Fig. 3b, respectively). However, it cannot be well detected because of the low accuracy of observations at G17. The data are consistent for both cases where G17 actually moves toward the land and where the landward motion is canceled or weakened by the remaining afterslip. Therefore, the data cannot constrain the degrees of contribution of viscoelastic response and afterslip to the slight trench-ward motion at FUKU in 2014-2017.
In any cases, we can consider that the trench-ward movements at FUKU and CHOS for the first 1-2 years were mainly caused by afterslip in the southern region (purple arrows in Fig. 3b). Although we cannot constrain and estimate the spatial extent of the southern afterslip region because of the low spatial density of geodetic observation sites at that time, annual-scale afterslip was a dominant deformation source in the southern region except for the vertical component at FUKU. Recalling the tsunami-derived rupture distribution at 37° N, for reproducing the trench-ward motion at FUKU, it might be reasonable to put the afterslip on the fault between the two peaks of the coseismic slip (Fig. 1b) rather than putting strong afterslip patches only in far south of FUKU (< 37° N) as considered in the forward slip model shown by Agata et al. (2019).
Considering that afterslip occurs on a fault with aseismic frictional properties, the assumed afterslip regions on the northern and southern sides of the main rupture are considered to have behaved as barriers to the rupture propagation. It is plausible that the shallow tsunamigenic slip in the off-Kamaishi and off-Fukushima areas (Satake et al. 2013) additionally loaded stress on the downdip side, helping to drive the afterslip. The northern afterslip occurred between the rupture zones of major earthquakes, i.e., the 1968 Tokachi-oki earthquake (M w 8.3), the 1994 offshore Sanriku earthquake (M w 7.7) (Nagai et al. 2001), the 1896 Meiji Sanriku tsunami earthquake (M w 8.1) (Satake et al. 2017), and the 2011 Tohoku-oki earthquake. The gap between the major earthquake ruptures is characterized by relatively low seismicity (Mochizuki et al. 2005), including repeating earthquakes (Uchida and Matsuzawa 2013;Igarashi 2020). Slow earthquake activity has been reported in this area (Nishikawa et al. 2019) as well. Although the northern and along-dip extensions of the afterslip zones cannot be resolved due to the absence of geodetic instruments, these features are consistent with the aseismic frictional properties. In the off-Fukushima region, co-and postseismic deformation sequences where the postseismic moment release due to the afterslip exceeded the coseismic moment of the mainshock were observed for the 2008 and 2010 Fukushima-ken-oki earthquakes (M j 6.9 and M j 6.7, respectively) (Suito et al. 2011). Although the spatial resolution for these sequences at that time, estimated only with terrestrial GNSS, is insufficient to compare with the afterslip region in the present study, it is possible that the proposed afterslip in the off-Fukushima region following the Tohoku-oki earthquake occurred in a region near the previous afterslip episodes.
The model proposed by Fukuda and Johnson (2021) seems reproduce the spatial distribution of afterslip inferred in the present study (Fig. 1b). However, there remain deviations in the temporal evolution, especially in the persistent subsidence observed at FUKU. Although it must be difficult to thoroughly reproduce the observations, GNSS-A results suggest that long-term viscoelastic relaxation in the oceanic asthenosphere should have a larger contribution to the movement at FUKU. In addition, because the precisions of preseismic seafloor data are insufficient (as shown in Table 2a and Fig. 1a), effects of the interplate coupling should also be evaluated simultaneously using the postseismic data.

Conclusions
Based on the GNSS-A observations and the above interpretations, the slip behaviors in the northern and southern areas can be summarized as follows.
(1) Afterslip occurred in the northern and southern regions away from the main rupture at depths near the hypocenter, which is consistent with the shallower tsunamigenic slip inducing stress on the downdip aseismic sections of the subduction thrust. It is plausible that these aseismic afterslip regions stopped the north-south rupture propagation (e.g., Nakata et al. 2021). (2) By at least 2014, the afterslip in the off-Kamaishi, off-Fukushima, and off-Choshi regions had almost fully decayed, though there is less evidence for the off-Fukushima region due to the insufficiency of information for determining the viscoelastic deformation patterns. (3) Shallow tsunamigenic slip in the south resulted in subsequent viscoelastic deformation that caused persistent seafloor subsidence, which was captured by postseismic seafloor geodesy. Additionally, the GNSS-A results indicate that the longterm viscoelastic relaxation process is continuing and is dominant in the off-Miyagi and off-Kamaishi regions. It also plays an important role in the off-Fukushima region, although its contribution cannot be well resolved. These long-term behaviors should be investigated by continuing and expanding seafloor geodetic observations.