Interplate channel layer controls of slip during and after the 2011 Tohoku-Oki earthquake


 There are remarkable differences between the middle and southern segments of the Japan Trench in terms of the seismic and aseismic slips on the plate interface and seismic velocity structures. The large coseismic slip of the 2011 Tohoku-Oki earthquake was limited to the middle segment, yet the observed negative residual gravity anomaly area in the southern segment corresponds to the postseismic slip area of the Tohoku-Oki earthquake. A model can explain the different slip behaviors of the two segments by considering their structural differences. The model indicated that the plate interface in the south was covered with a thick channel layer, as noted by seismic survey imaging, and this layer resulted in the residual gravity anomaly. Numerical simulations that assumed obvious frictional heterogeneity caused by the layer existing only in the south successfully reproduced M9 earthquakes recurring only in the middle, followed by evident postseismic slip in the south. We suggest that, while the layer makes the megathrust less compliant to seismic slip, it promotes aseismic slip following the growth of seismic slip on the fault in an adjacent region.


Introduction
Many studies of spatiotemporal distributions of various seismic and geodetic events and underground structures have demonstrated that there are remarkable differences between the middle and southern segments of the Japan Trench. The shallow area off the coast of Fukushima Prefecture, the southern segment of the trench, is known to have few earthquakes. Interplate coupling in the southern segment is signi cantly weaker than that off the coast of Miyagi Prefecture-the middle segment of the trench 1 .
Seismic and geodetic data have demonstrated that the large coseismic slip in the shallow area that caused the 2011 Tohoku-Oki earthquake was limited to the middle segment 2,3,4 . Additionally, paleoseismological evidence collected along the Japan Trench demonstrated that deep-sea turbidites that originated from large earthquakes potentially similar to the 2011 Tohoku-Oki earthquake are recorded only in the middle segment 5 . The characteristics of spatial heterogeneity along the trench may not have changed over a long period.
In contrast, postseismic slip was predominant in the southern segment 6,7 . Slow earthquakes include lowfrequency earthquakes, tectonic tremors, and very-low-frequency earthquakes (VLFs), which are distributed in a manner complementary to the coseismic slip and overlap the postseismic slip area of the Tohoku-Oki earthquake 8,9,10 (Fig. 1).
The postseismic slip during April-December 2011 was estimated to be less than 1.2 m 6 , and the one-year postseismic slip following the 2011 earthquake was estimated to be approximately 2.5 m 11 at the shallow portion of the southern segment of the Japan Trench. As of December 2017, analysis of sea oor displacement demonstrated that postseismic slip was continuing in the southern segment, and the average slip rate between January 2015 and December 2017 were halved of the average slip rate between March 2011 and December 2014 off the coast of Fukushima Prefecture (at the FUKU station) 12 .
Researchers have been working to determine what controls these complementary spatial distributions of slips along the Japan Trench. Liu and Zhao 13 showed a relationship between the fault slip behavior in the Tohoku-Oki earthquake and the P-wave tomography of the Tohoku forearc. They suggested that structural heterogeneities in both the overriding upper plate and atop the subducting lower plate controlled the processes of nucleation and rupture in the Tohoku-Oki earthquake. Bassett et al. 14 found that the large coseismic slip zone of the 2011 earthquake is delimited by the evident negative residual gravity anomaly zone in the south. Distributions of postseismic slip from the Tohoku-Oki earthquake and slow earthquakes correspond to areas with large negative values of residual gravity anomalies 14 .
In this study, we focused on the distribution of the sedimentary units observed on multichannel seismic (MCS) sections along the trench 15 . The sedimentary units represent a narrow wedge shape in the middle segment but take on a thin and broad distribution along the plate interface, forming a channel layer 15 in the southern segment. The borders of these two different characteristics in sediment distribution correspond very well to the southern border of the extensive coseismic slip zone, as well as to the boundary delimiting the negative residual gravity anomaly area (Fig. 1). It seems plausible that the existence of the channel layer controls the slip behavior along the plate interface because there may be a positive relationship between the thickness of a fault layer and a frictional parameter controlling the slip weakening.
The channel layer can be regarded as a fault gauge layer. Laboratory measurements and observations have suggested that the thickness of the fault gauge layer is positively correlated with the critical slip distance, a type of frictional parameter in the rate-and state-dependent friction law 16,17,18 . In the framework of the friction law, where critical slip distances are large, aseismic slip is more likely to occur 19,20 . Then, we hypothesize that the large value of the critical slip distance spreads on the plate interface in the southern segment. This hypothesis is consistent with the aforementioned observation, as some types of aseismic slips are observed in the southern segment, where a thick channel layer exists.
Based on these previous studies, the presence of the channel layer can be expressed as an increased critical slip distance in terms of frictional mechanics. If this theory is correct, then the massive coseismic slip during the mainshock of the 2011 event did not occur, whereas evident aseismic postseismic slip took place in the area of negative residual gravity anomaly, re ecting the presence of a channel layer.
In this study, we semi-quantitatively con rmed that the presence of the channel layer is consistent with the observed negative residual gravity anomaly by constructing a crustal density model based on seismic velocity models. We also conducted numerical simulations of earthquake generation cycles along the Japan Trench by considering the presence of a channel layer with a characteristic slip distance value. By attempting to reproduce megathrust earthquake cycles along the trench with four frictional models, we revealed how structural heterogeneity controls the occurrence of earthquakes and complementary spatial distributions of slip in the study region.
Structural Model Expressing Observed Differences Between The Middle And Southern Segments Bassett et al. 14 explained that the evident difference in residual gravity anomaly between the middle segment, where large coseismic slip occurred, and the southern segment, characterized by shallow postseismic slip, could be explained by the large contrast in the density of the hanging wall side of the plate interface. They assumed a single-layered structure for the hanging wall and showed that the density contrast of ~180 kg/m 3 in the layer accounts for the contrast of the observed residual gravity anomaly across the segment boundary.
However, the results of wide-angle active seismic experiments do not support structural differences throughout the hanging wall crust. According to the P-wave velocity (Vp) models by Miura et al. 21,22 obtained for the middle and southern segments, in the deeper part of the crustal layer underlain by the plate boundary, the difference in Vp is insigni cant (Supplementary Fig. 1) and too small compared to the large density contrast to explain the residual gravity anomaly. According to the standard relationship between Vp and the density of the crustal materials 23 , a density contrast of 180 kg/m 3 is equivalent to a Vp anomaly of approximately −10%, which is not consistent with the Vp models using active seismic explorations ( Supplementary Fig. 1).
Here, we hypothesized that it is the presence of the channel layer imaged by seismic surveys in the southern segment that accounts for the observed spatial pattern of the residual gravity anomaly. We con rmed this hypothesis by performing a simple gravity calculation with a horizontally strati ed twodimensional (2D) model composed of two segments, which correspond to the middle and southern segments of the Japan Trench forearc, respectively (Method; Supplementary Fig. 2).
We have presented three density models that explain the observed residual gravity anomaly equally well ( Fig. 2 and Supplementary Fig. 2). One is composed of a single crustal layer (Model A1) similar to that presented by Basset et al. 14 . The other two are based on low-Vp zones identi ed by Miura et al. 21,22 .
There are two signi cant low-Vp zones, which are candidates for layers of a negative density anomaly accounting for the negative residual gravity anomaly at depths, that range from 2 to 5.5 km and 10 to 11 km, respectively ( Supplementary Fig. 1). The model that assumes a shallower layer with a thickness of 3.5 km is solely responsible for the negative residual gravity anomaly (Model A2) has an extreme density contrast of 400 kg/m 3 , which is much larger than the expected anomaly (−200 to −100 kg/m 3 ) from the estimated Vp anomaly (−10%). The assumed density contrasts in Model A3, where the channel layer has a thickness of 1 km, are well within the plausible range of density variations expected from the Vp pro les. We have concluded that a structural model that considers the presence of a low-density channel layer is the most likely, as it can explain the observed residual gravity anomaly as well as the Vp model obtained in the region without any signi cant contradictions.

Effects Of The Channel Layer On The Spatiotemporal Slip Distribution
We performed calculations using numerical simulations across an extensive period (approximately 5000 years) and identi ed 7-8 instances of M9 (M = 8.9-9.1) earthquakes, using four models with different parameter settings (Supplementary Tables 1-3 and Supplementary Fig. 3-4). For all the models, the time intervals of the M9 earthquakes were 540-770 years ( Fig. 3a-b, and Supplementary Fig. 5a-d, black line). However, coseismic slip areas larger than 10 m were quite different in spatial extent ( Fig. 4 and Supplementary Fig. 6). M9 earthquakes in Models B1 and B2 ruptured both middle and southern segments, whereas the rupture areas of Models B3 and B4 were limited to the middle segment of the Japan Trench. The rupture patterns produced by Models B3 and B4 are consistent with surveys of turbidites, which showed depositional events occurring at intervals of 500-900 years over the past 4000 years at the middle segment 5,24 .
Postseismic slip was dominant in the shallow part of the southern segment in the results of Models B3 (Fig. 4) and B4 (Supplementary Fig. 6c) The slip continued for more than ve years in both Model B3 (  Table 3, and Supplementary Fig. 4)).
If we assume a gradual variation in L between the middle and southern segments (Model B1; Supplementary Table 2, and Supplementary Fig. 3), the M9 coseismic slip propagated southward, and postseismic slip did not dominate at the southern segment ( Supplementary Fig. 6a). If we adopted the segmentation boundary with the smaller contrast of frictional parameters (Model B2 in Supplementary  Table 1 and Supplementary Fig. 3d, purple line), the resulting slip distribution was similar to that of Model B1 ( Supplementary Fig. 6b), and postseismic slip at the southern segment quickly terminated ( Supplementary Fig. 5e, purple line and 9).
On the other hand, the extremely large contrast in L between the two segments is not consistent with the postseismic slip observations. If we used a greater value of L in the southern segment (Model B4 in Supplementary Table 1, and Supplementary Fig. 3d, green line), M9 coseismic slip was limited to the middle segment, similar to Model B3 (Supplementary Fig. 6c). However, the propagation of postseismic slip to the south would take a longer time than in Model B3 (Supplementary Fig. 8). Consequently, the postseismic slip started in the southern segment approximately 0.5 years after the simulated M9 earthquake ( Supplementary Fig. 5e, green line). The signi cantly delayed onset of postseismic slip is not consistent with observation of the 2011 Tohoku-Oki earthquake estimated from small repeating earthquake activity, which increased shortly after the mainshock 25 , although the exact timing of the onset is not evident from geodetic observations. Therefore, the temporal distribution of postseismic slip could constrain the contrast of L between the middle and southern segments.
When the contrast of L between the two segments was similar to that in Model B3, we could explain the remarkable difference in the coseismic and postseismic slip distributions along the trench and the long duration of the postseismic slip, as observed. The spatiotemporal distributions of coseismic and postseismic slips of Model B3 were roughly consistent with geodetic observations 3,7 . However, in our simple model, the amount of postseismic slip was considerably larger than that observed in the southern segment 6,11 . However, there are large uncertainties in both observations and simulations. Although the present modeling indicated that most of the slippage occurred immediately after the M9 earthquake, it is very di cult to constrain aseismic fault motion outside of the coseismic rupture area in the very early postseismic period from seismic and geodetic observations. Early postseismic behavior might be subject to treatments of dynamic rupture and/or viscoelasticity of the earth. For example, dynamic weakening processes such as thermal pressurization 26 and viscoelastic relaxation 11,27,28 have to be carefully considered to quantify coseismic/postseismic slips in detail. It may be necessary to include these physical and structural heterogeneities, which are smaller than the resolution of the structural surveys in our model.

Conclusions
We showed that the presence of a thick, low-velocity channel layer in the subducting plate at the southern segment of the Japan Trench is consistent with the characteristics of the gravity anomaly distribution, and accounts for the signi cantly different slip behavior of the middle and southern segments near the trench axis. A structural model that considers the presence of the channel layer is the best way to explain the observed residual gravity anomaly and seismic velocity without any signi cant contradictions. By representing the presence of the channel layer as the value of a characteristic slip distance in the rateand state-dependent frictional law, we successfully reproduced the longer duration of the postseismic slip of the 2011 Tohoku-Oki earthquake, as well as the spatial distributions of the complementary coseismic and postseismic slips at the middle and southern segments, maintained for recurrence intervals of ~ 600 years.

Method Of Estimation Of Density Distribution
In the present study, we sought to create 2D density anomaly models explaining the difference in the residual gravity anomaly across the forearc segment boundary (FSB) presented by Bassett et al. 14 in the southern segment of the Japan Trench, where postseismic slip is evident. We calculated gravity anomalies using three different density distribution models (Models A1, A2, and A3, in Supplementary   Fig. 2) using the classical Talwani method 29 and compared them to the observed residual gravity anomaly. The crustal layer is assumed to have a thickness of 10 km, since the depth to the plate boundary is ~ 10 km below the sea oor in the postseismic slip zone. We included the observed gravity anomaly in the region with a plate boundary depth from 9 to 11 km. We created a 2D pro le along the trench axis and set x = 0 at the location of the FSB. A negative density anomaly was given to the southern half of the models (x < 0). Model A1, the simplest model, had a homogeneous density anomaly of 180 kg/m 3 in the crust, and was similar to the model used by Bassett et al. 14 to interpret the observed gravity anomaly. Model A2 had a thinner layer of density anomaly, emulating the lower sedimentary layer imaged as the low-Vp layer found by seismic explorations (Supplementary Fig. 1). Model A3 had two layers of negative density anomaly: the rst was the sedimentary layer corresponding to that in Model A2, but with a smaller anomaly, and the second was below the base of the crustal layer. The calculated gravity anomaly was adjusted by adding 40 mgal in Fig. 2 because we were interested in the amount of difference in the residual gravity anomaly across the FSB, not in the absolute values of the anomaly.

Method Of Simulation Of Earthquake Generation Cycles
We conducted numerical simulations of earthquake generation cycles using the realistic threedimensional geometry of the subducting Paci c Plate along the Japan Trench 30 . We used the same equations, initial conditions, seismic radiation damping term, plate geometry, and plate convergence rate as in our previous study 31 , which roughly approximated several characteristics of M7 earthquakes (such as the 1978 Miyagi-Oki earthquake) and M9 earthquakes (such as the 2011 Tohoku-Oki earthquake) in the middle segment of the Japan Trench. Seismic and aseismic events were modeled to represent the release of the slip de cit or backslip that accumulates during interseismic periods 32 . Such spatiotemporal variations in the slip velocity were assumed to indicate an unstable slip with a frictional interface. We used a rate-and state-dependent frictional law 16 as an approximate mathematical model for large-scale frictional behavior on the plate interface and a fault constitutive law 33 to determine the slip rate for a given stress and strength value. In addition, we used an aging law 16,34 , which can be considered an evolution law for changes in strength that varies depending on prior slip history. In our simulation of earthquake generation cycles, we used the quasi-dynamic (QD) model 32 with a smaller radiation damping term instead of a fully dynamic (FD) model. Usually, the QD model cannot simulate fast slip during seismic events; however, in models with standard rate-and-state friction and relatively uniform fault properties, the FD and QD simulations produce qualitatively similar fault behaviors, with crack-like ruptures and similar earthquake patterns 35 .
During our simulations, we only changed the numerical values and the spatial distributions of the frictional parameters (A, B, and L), focusing on the slip in the shallow, southern area. Nakata et al. 31 approximated the spatial distributions of frictional parameters for seismic sources with circular patches for M7 earthquakes, which repeatedly occurred at intervals of several decades in the past. In the present study, we removed these circular patches for simplicity and used the model without these unstable patches (Model B1, Supplementary Fig. 3, Supplementary Table 2), instead of Nakata et al. 31 , for comparison. The calculation time was then reduced so that we could calculate multiple cycles of the M9 earthquakes. In a method similar to that of Nakata et al. 31 , the M9 coseismic rupture area located in the middle segment of the Japan Trench was approximated by a rectangle (an area 150 km long along the strike and shallower than 22 km in depth). At the boundary between the southern and northern sides of the M9 source area, the frictional characteristics changed gradually, and differences between the deep and shallow areas were not considered.
For Model B3 (Supplementary Fig. 4, Supplementary Table 3), the values and spatial distribution of the frictional parameters (A, B, and L) were determined with reference to observations of the residual gravity anomaly by Bassett et al. 14 . The M9 coseismic rupture area was approximated by a rectangle (170 km along the strike and shallower than 24 km depth). Here, we divided the surrounding area, which we assumed to have a uniform frictional condition as in Model B1, into three segments: northern, southern, and deep segments ( Supplementary Fig. 4). Difference of frictional parameters at the boundaries between the northern, southern, and middle (the M9 coseismic) segments were then assumed to be larger than those of Model B1 (Supplementary Fig. 3c-d), and their locations were based on an observed residual gravity anomaly < 0. Spatial distributions of the coseismic slip, postseismic slip, seismic velocity, and residual gravity anomaly showed that uniform and smooth assumptions of frictional boundaries may not be appropriate. Boundary settings along the trench between the middle and southern segments were the main differences between Model B1 and Model B3. In addition, based on our hypothesis, the frictional parameter controlling slip weakening (L) at the southern segment is assumed to be larger than that of the middle and deep segments. Through iterative processes of trial and error, ne-tuning of the frictional parameters was undertaken between Models B1 and B3.
To assess the effect of the channel layer, we conducted calculations using four models (Supplementary Table 1). Model B1 was almost the same as in our previous study. This model did not have a channel layer. Models B2, B3, and B4 were the modi ed models in this study. Model B2 had a small L value (the same as Model B1), Model B3 had a medium L value, and Model B4 had a large L value in the southern segment. Model B2 had a channel layer in which contrast between the middle and southern segments was weak. Models B3 and B4 represented models with a high-contrast channel layer. The L values in the southern segment were 2.1, 4.6, and 6.4 times larger than in the middle segment for Models B2, B3, and B4, respectively.
Declarations the design of the study. R.H. conceived and participated in the design of the study, helped draft the manuscript, and conducted density distribution modeling. All authors have reviewed the manuscript.