Frictional and structural controls of seismic super-cycles at the Japan trench

The diverse mechanical behaviors of subduction zones during the seismic cycle emerge from the nonlinear dynamics of a complex mechanical system with interacting brittle and ductile deformation. The 2011 Tohoku mega-quake represents the culmination of a super-cycle of partial and full ruptures of the plate interface, but the physical controls on the down-dip segmentation of the megathrust remain unclear. Here, we propose a two-dimensional rheological model of the Japan trench to explain the variability of earthquake size at the Miyagi section, in northern Honshu, during the last century. We simulate seismicity in a continuum with a physics-based rate- and state-dependent constitutive law for fault slip, producing aperiodic earthquake sequences with a power-law distribution of rupture sizes. Although some partial ruptures of the megathrust are the result of self-emergent behavior, others are structurally controlled. The 1978 and 2005 Mj∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}7 interplate earthquakes took place in a metamorphic belt in the mantle wedge corner bounded up-dip by the arc Moho that constitutes a permanent structural boundary. We explain the succession of the great 1981 Mw = 7.1 and 2003 Mw = 6.9 earthquakes within the forearc and the 2011 giant earthquake as the natural response of a large continuously velocity-weakening fault with a small nucleation size. The shallow segment below the frontal prism only slips in giant through-going ruptures that unzip the whole velocity-weakening interface. The model consistently explains the size, recurrence time, and hypocenter location of historical earthquakes in the Miyagi segment, the slow-slip and foreshock preparatory phase of the 2011 Tohoku earthquake, the large slip near the trench during the giant rupture, and important features of its postseismic deformation. The complex patterns of seismicity at the Japan trench can be better understood by assimilating geological and geophysical observations at various periods of the seismic cycle within an explicative physical framework.


Introduction
Natural earthquakes at plate boundaries result from frictional instabilities on faults that rupture to accommodate the long-term relative movement of tectonic plates. Earth's largest earthquakes occurred at subduction zones where the low-angle megathrust offers expansive area for rupture propagation. The Japanese archipelago is an island arc formed by intra-oceanic subduction zones between the Philippines Sea and Eurasian plates in the Nankai trough, and between the Pacific and the Okhotsk plates at the Japan trench in northern Honshu. The seismic cycle along the Nankai trough is characterized by a down-dip segmentation of rupture styles with slow-slip events, low-frequency earthquakes, and large ruptures that is typically persistent along strike (Hirose et al. 2010;Obara and Kato 2016;Gao and Wang 2017). The dynamics of the seismic cycle at the Nankai trough is discussed in a companion paper in the same volume (Shi et al. 2020). The Japan trench was the epicenter of the 2011 moment-magnitude (Mw) 9.1 Tohoku mega-quake that ruptured the megathrust up to the surface, generating a giant tsunami (Tajima et al. 2013;Satriano et al. 2014;Lay 2018). Previous earthquakes in the Miyagi region were confined to smaller offshore regions (Fig. 1). Barbot Earth, Planets and Space (2020) 72:63 The 1978 and 2005 Mj∼ 7 earthquakes ruptured overlapping regions of the megathrust close to the coast (Okada et al. 2005;Yaginuma et al. 2006;Umino et al. 2006). The 1981 Mw = 7.1 and 2003 Mw = 6.9 earthquakes ruptured adjoining sections of the megathrust farther offshore (Yamanaka and Kikuchi 2004;Nakata et al. 2016).
The large variability of earthquake size is understood as the hallmark of a seismic super-cycle with partial and full ruptures of the megathrust that has operated for at least 4000 years, culminating in the 869 Jogan, 1454 Kyotoku, 1896 Sanriku, and 2011 Tohoku giant earthquakes (Usami et al. 2018). Seismic super-cycles are a relatively recent discovery. For example, coral, turbidite, and paleo-tsunami records indicate a broad range of earthquake sizes in the Sunda trench (Sieh 2008;Patton 2015;Philibosian 2017;Rubin et al. 2017). Along the Cascadia subduction zone, turbidite paleoseismology indicates 19 full-length ruptures and about as many partial ruptures (Goldfinger et al. 2003a), some of them linked to tsunami deposits in Japan (Goldfinger et al. 2003b). Partial rupture of the seismogenic zone also occurs at continental collisions, for example, during the recent 2015 Mw = 7.8 Gorkha, Nepal earthquake (Mencin 2016;Hubbard et al. 2016;Qiu 2016;Feng 2017). Understanding the origin of super-cycles is of fundamental importance for seismic hazards mitigation in Japan, as only Mw~7.0 earthquakes took place in the Miyagi section in the century preceding the 2011 Tohoku megaquake (Yamanaka and Kikuchi 2004;Nakata et al. 2016), providing little insights at the time into the imminent devastation. There are several flavors of super-cycles, characterized by either along-strike or down-dip segmentation of locked asperities. In this study, we focus on the emergence of downdip partial ruptures, perhaps best exemplified outside of Japan in a subduction setting by the sequence of the 2007 Mw = 8.4 Bengkulu earthquake shortly followed by the 2010 Mw = 7.8 Mentawai tsunami earthquake in Sumatra, Indonesia (Konca 2008;Tsang et al. 2016;Feng et al. (Nakata et al. 2016), and the 2011 Mw = 9.1 Tohoku earthquake (Bletery et al. 2014), are shown in colored contours. The continental shelf (dashed profile) provides an indication on the extent of the Japanese island arc (volcanos as red triangles). The triangles represent GEONET GPS stations on land (Sagiya et al. 2000;Ozawa et al. 2011), and seafloor geodesy monuments offshore (Watanabe et al. 2014;Tomita et al. 2015). Only the geodetic data within 50 km of the two-dimensional cross-section A-A' are considered (black triangles) Barbot Earth, Planets and Space (2020) 72:63 2016). Why some large buried ruptures would remain confined at great depth while other ruptures break all the way to the surface, generating a giant earthquake? Some investigators argue that the succession of partial and full ruptures is a self-emergent behavior from the nonlinear dynamics of fault slip favored by wide seismogenic zones (Herrendörfer et al. 2015;Michel et al. 2017;Barbot 2019b). Others suggest the important role of local changes in the fault geometry (Qiu 2016;Hubbard et al. 2016;Sathiakumar et al. 2020). At subduction zones, where the megathrust traverses various geological formations, we expect a combination of these two factors to operate a strong control on the seismic cycle. What mechanism controls the down-dip segmentation of the megathrust at the Japan trench remains unclear, in particular, because the frictional behavior and seismogenic potential of the shallow megathrust underneath the accretionary prism are poorly understood. Assuming that the shallow megathrust is not amenable to rupture nucleation, thermal pressurization of the shallow interface or wave-mediated stress transfer in weak sediments can be invoked to explain trench-breaking ruptures (Mitsui et al. 2012;Kozdon and Dunham 2013;Noda and Lapusta 2013;Cubas et al. 2015;Shibazaki et al. 2019). However, the shallow megathrust section below accretionary prisms may host tsunami earthquakes (Kanamori 1972;Satake and Tanioka 1999;Bilek and Lay 2002;Geersen 2019), low-frequency earthquakes (Toh et al. 2018), and tectonic tremors (Sugioka 2012;Yamashita 2015;Nakano et al. 2018). Laboratory friction experiments also unravel complex frictional properties of accretionary sediments that is not incompatible with an unstable fault behavior (Saffer and Marone 2003;Ikari 2019). In particular, Ikari and Kopf (2017) highlight the velocity-weakening behavior of sedimentary muds at plate tectonic slip rates, which is a necessary condition for rupture self-nucleation. If the shallow layer is assumed frictionally unstable, other mechanisms may be invoked to explain segmentation, for example, the presence of rheological barriers (Kaneko et al. 2010) or morphological gradients (Qiu 2016;Ong et al. 2019;Sathiakumar et al. 2020). Within this view, giant surface-breaking ruptures may represent the synchronous failure in a single event of otherwise isolated asperities (Yamanaka and Kikuchi 2004). Considering the megathrust as a deterministic nonlinear system, segmentation may also happen spontaneously due to residual stress left by previous ruptures and non-uniform loading of the fault, even in a homogeneous velocity-weakening domain (Lapusta and Rice 2003;Wu and Chen 2014;Cattania and Segall 2019;Barbot 2019b;Cattania 2019). The remaining challenge is to determine the respective roles of the frictional and structural controls on the seismic super-cycles at the Japan trench. Previous efforts have already focused on understanding how the geological setting may influence seismicity (Shimizu 2014;Keren and Kirkpatrick 2016;Murphy 2018). Here, we extend the work to seismic cycle simulations and modeling of geodetic data throughout the seismic cycle.
We will argue that the occurrence of deep Mw~7 interplate earthquakes below the arc Moho can be attributed to a permanent structural boundary traversed by the megathrust, but that the super-cycles of forearc Mw~7 and Mw~9 ruptures can take place in an uninterrupted velocity-weakening domain stretching from the arc Moho to the trench, without invoking the presence of structural barriers or complex rheological behavior. To support these assertions, we present a two-dimensional structural and rheological model of the Japan trench at the Miyagi section. In the next section, we identify the important structural boundaries that shape the down-dip variation of the frictional and other physical properties that control fault dynamics on the megathrust. Then, we describe the governing equations and physical parameters considered in simulations of earthquake, aseismic slip, and viscoelastic flow that implement the structural model. We then present a numerical simulation that reproduces many important aspects of the seismic cycle in the Miyagi segment during the last century. This is followed by concluding remarks.

First-order structural and rheological model of the Japan trench
Northern Honshu represents an old intra-oceanic subduction zone separating the Japan trench from the Sea of Japan in the backarc (Wakita 2013). As the Pacific plate dives underneath the Okhotsk plate, the megathrust intersects different terranes, which, along with systematic changes in temperature and pressure, affect the frictional and rheological properties of the plate interface. To describe a first-order model of the structural boundaries, we consider a representative two-dimensional trenchperpendicular cross-section of the Japan trench along the profile A-A' in Fig. 1. The geological structure of the Japan trench is constrained by field observations (Taira et al. 2016), complemented with seafloor bathymetry and gravity (Bassett et al. 2016), active source seismic profiling , seismic tomography (Nakajima et al. 2001;Hasegawa et al. 2005;Zhao 2015;Zhao et al. 2015;Wang et al. 2017), and seismic source imaging Nishikawa et al. 2019), among others. Each dataset provides insights into a particular domain of the plate interface with limited resolution, especially at great depths. Despite this limitation, the collective set of observations at the Japan trench combined with insights from geodynamics at other subduction zones allows us to build a tentative interpretation of how the geological structure may control the seismic cycle.
The megathrust extends from the trench to a depth of about 80 km, where the down-going slab plunges into the asthenosphere, marking the brittle-ductile transition (Muto et al. 2013). Below this depth, the down-going slab entrains the rocks of the hanging wall, the locking between the two plates is permanent, and anelastic deformation in the upper plate is accommodated by the plastic flow of peridotite. This major rheological boundary is ostensibly evidenced by not only arc volcanos further inland, scattered between 80 and 100 km above the slab (Syracuse and Abers 2006;Hayes et al. 2018), but also the change to more dominantly intra-slab seismicity , and the pattern of surface heat flow (Wada and Wang 2009).
At shallower depths, the deformation is localized along the megathrust. Consideration of complementary geophysical datasets offers a cohesive picture of structural variations along the megathrust in relation to seismicity (Liu et al. 2014;Liu and Zhao 2018). The frontal prism (Kodaira et al. 2012;Nakamura et al. 2013) represents a wedge-shaped block of the accretionary prism below the middle and lower slopes with intense internal deformation bounded by a mega-splay that forms a backstop Tsuji et al. 2011;Boston et al. 2017). The frontal prism is packed with fluid-filled sediments characterized by low rigidity and low frictional strength. Indeed, experimental data on clay-rich minerals indicate a consistently weak strength, with steadystate coefficients of sliding friction lower than 0.35 (Ikari et al. 2009;Tesei et al. 2012;Ikari et al. 2013;Tesei et al. 2015). Faulting and folding within the wedge to accommodate shortening result in reduced long-term slip velocity along the megathrust (Sathiakumar et al. 2020). As the long-term loading rate exerts a direct, first-order control on the recurrence time of frictional instabilities, the reduced loading rate at the base of the frontal prism may represent an important factor for down-dip segmentation. Seismic reflection and differential bathymetry data show clear evidence of fault rupture at the trench axis during the 2011 Tohoku earthquake Kodaira et al. 2012). However, this section did not participate in other great earthquake ruptures in the past century.
The arc crust constitutes an important structural domain that delineates down-dip segment boundaries on the megathrust. Wide-angle reflection and refraction imaging Ito et al. 2005) and seismic tomography (Yamamoto et al. 2006;Liu et al. 2014) provide a synoptic view of hanging wall fabric at crustal depths, revealing that the arc Moho extends offshore beyond the shelf to the epicentral region of the Tohoku earthquake, between the rupture areas of the 2003 and 2005 earthquakes, about 120 km from the trench axis. This indicates that the shallow 1981 Mw = 7.1 and 2003 Mw = 6.9 earthquakes ruptured a segment of the megathrust confined to the forearc region underneath thick sedimentary rock layers. The up-dip boundary of this seismogenic zone may coincide with a forearc section of the accretionary prism, the forearc upper crust, that is less actively deforming than the frontal prism, but also pilling poorly consolidated sediments. In contrast, the 1978 and 2005 Mj~7 earthquakes that took place some 150 km from the trench, below the shelf of Japan, close to the coast, presumably originated from a metamorphic belt below the arc Moho, where the megathrust intersects the mantle wedge corner Ito et al. 2005). The rupture areas of the 2003 and 2005 Miyagi interplate earthquakes may therefore be separated down-dip by a permanent structural boundary: the arc Moho. The different geological settings of the 2003 and 2005 earthquakes offshore Miyagi may be generalized to Northern Honshu, where the historical seismicity can be divided into a shallow thrust zone hosting great earthquakes and a deep thrust zone where large earthquakes occur (Kawakatsu and Seno 1983). The downdip limit of the metamorphic belt marks the transition to stable creep on the megathrust, associated with the velocity-strengthening friction of olivine aggregates that may persist up to about 600 • C or more (Boettcher et al. 2007;King and Marone 2012). Here, the transition from velocity-weakening to dominantly velocity-strengthening behavior is evidenced by the down-dip limit of the deep Mj∼ 7 interplate earthquakes (Okada et al. 2005;Yaginuma et al. 2006;Umino et al. 2006). Below these depths, the megathrust cuts through upper mantle rocks and fault slip is largely aseismic, only modulated by the stress transfer from up-dip seismic events, leading to creep during interseismic periods and afterslip during postseismic periods. The transition from velocity-weakening to velocity-strengthening friction may involve small unstable patches of 0.1 to 1 km radius in the nominally velocity-strengthening region, explaining the presence of small magnitude continuously repeating earthquakes down to 60 km depth (Igarashi et al. 2003) and deep aftershocks after the 1978 Mj∼ 7 Miyagi-Oki earthquake (Seno et al. 1980).
A first-order representation of the structural layout of the Japan trench along the A-A' profile is shown in Fig. 2. The rate of convergence between the Pacific and Okhotsk plates varies depending on estimates, for example, from 73.50 mm/year for REVEL 2000 (Sella et al. 2002) to 92.45 mm/year for MORVEL 2010 (DeMets et al. 2010) with intermediate values for GEODVEL, GSRM, ITRF, and NUVEL 1, leading to a range of 79.9 Barbot Earth, Planets and Space (2020) 72:63 ± 12.3 mm/year. Throughout this work, we will assume a trench-perpendicular convergence rate of 84mm/year, consistent with previous crustal dynamics studies (Kato et al. 2011;Usami et al. 2018). The model is limited to brittle deformation localized on the megathrust and viscoelastic deformation in the asthenosphere (Barbot 2018b) and ignores internal deformation in the hanging wall or small-scale viscoelastic deformation at crustal depths below the volcanic arc (Hu et al. 2014;Muto et al. 2016Muto et al. , 2019. From the trench to the brittle-ductile transition, the megathrust is divided into characteristic segments delineated by permanent structural boundaries associated with different provinces of the hanging wall. The top megathrust segment underlies the frontal prism near the trench and is characterized by a low steady-state friction coefficient and a low long-term slip rate. The shortening rate accommodated by splay faults and folding in the frontal prism is poorly known, but perhaps amounts to 40% of the total convergence, due to the presence of many reverse faults and offset terrasses . We therefore assume a local loading rate of V L = 52mm/year in the frontal prism megathrust segment, representing the remaining 60% of the plate convergence, and a basal coefficient of sliding friction of µ 0 = 0.3 . To represent the contractional nature of the frontal prism in the model, we introduce a single mega-splay fault. Farther down, the sub-horizontal megathrust intersects the accretionary sediments of the forearc upper crust and the interface is also characterized with low friction. However, the sediment pileup landward of the backstop is not actively deforming. Hence, the long-term slip rate at this segment and at those following down-dip equals the convergence rate between the Pacific and the Okhotsk plates, i.e., V L = 84mm/year. The next segment down-dip is the forearc basement, where the megathrust intersects the sedimentary rocks of the forearc at a higher dip angle. Due to the increased depth of burial, rocks are more cohesive and their friction properties change accordingly. The forearc basement exhibits unstable weakening properties and hosts recurring Mw∼ 7 earthquakes. The next megathrust segment is the intersect with the mantle wedge corner, located offshore below the arc Moho. The composition, frictional properties, and  Fig. 1. The velocity-weakening and velocity-strengthening segments of the megathrust are shown in yellow and gray color, respectively. Coupling between the two plates is permanent below the brittle-ductile transition. The ambient temperature in the oceanic asthenosphere and in the mantle wedge asthenosphere is shown with the background color. The triangle mesh is used for sampling of the viscoelastic deformation. Contours of steady-state viscosity for 10 20 Pa s, 10 21 Pa s, and 10 22 Pa s are shown with the yellow, green, and blue profiles, respectively. The small numbers at the corner of fault segments indicate their starting depth in the model ambient temperature of the metamorphic rocks at these great depths are unclear. However, the recurring Mj∼ 7 interplate earthquakes at this location clearly indicate unstable weakening. Velocity-weakening friction at high temperature is compatible with the frictional properties of antigorite serpentinite, which are velocity-strengthening at low temperature (Reinen et al. 1991), but velocity-weakening between 400 and 650 • C in wet conditions Sone et al. 2012;Okazaki and Katayama 2015;Brantut et al. 2016). We therefore envision that the megathrust crosses a serpentinized metamorphic belt in the mantle wedge corner, characterized with velocity-weakening properties at these particular hydro-thermal conditions (Goswami and . The arc Moho separates the forarc basement from the mantle wedge corner, providing a permanent structural boundary. We materialize this boundary in the model by a narrow stretch of velocity-strengthening rocks between the two domains, in the forearc lower crust. This barrier is strong enough to allow separate Mw∼ 7 ruptures above and below the arc Moho, but weak enough to allow those rare giant ruptures to break through (Kaneko et al. 2010). We model the mantle wedge corner segment as a homogeneous velocity-weakening region, ignoring internal variations that may explain the complex rupture sequence of the 1978 Mj∼7.0 Miyagi-Oki earthquake (Seno et al. 1980). Below the mantle wedge corner, the megathrust cuts the cold nose of the mantle wedge and the upper mantle region. These two regions are characterized by velocity-strengthening properties, but gradients of frictional properties may operate due to continuously varying temperature, degree of serpentinization, or fluid pressure. We therefore allow distinct properties in these two regions. We include the mega-splay fault in the model, but there are no definitive constraints about its frictional behavior or long-term slip rate. The long-term slip accumulation on the mega-splay fault is presumably a fraction of plate convergence, which we take as 40%. We treat the mega-splay fault as potentially seismogenic, but our numerical simulations indicate that little slip accrues on this fault over the period considered and that ruptures propagating on the megathrust rarely branch onto the mega-splay fault. As a result, the current study does not offer new insights into the mechanics of this structure.
These considerations provide important guidelines for the spatial distribution of frictional properties on the megathrust. More quantitative constraints on the many parameters involved come from the size and recurrence patterns of seismic events and require considering the dynamics of the megathrust as a complex mechanical system. The constitutive properties for the viscoelastic flow of olivine are better constrained from laboratory experiments (Karato and Wu 1993;Karato and Jung 2003;Hirth and Kohlstedt 2003) and from the few studies of postseismic relaxation that implement a compatible rheological framework (Masuti et al. 2016;Sobolev and Muldashev 2017;Barbot 2018b;Muto et al. 2019;Agata et al. 2019).

Frictional and rheological properties
We describe the building blocks of a seismic cycle simulation that incorporates the structural model described in the previous section. The model incorporates the constitutive relationships governing fault slip along the megathrust and viscoelastic flow in the oceanic asthenosphere and mantle wedge. For the frictional constitutive behavior, we assume a micro-physical model of rate-and state-dependent friction in isothermal conditions where the frictional strength is controlled by the real area of contact and the velocity of sliding (Barbot 2019a), following the constitutive relationship where τ is the norm of the shear traction resolved on the fault plane, V the instantaneous sliding velocity, A the real area of contact divided by the nominal surface area, χ the indentation hardness of the fault material, µ 0 a reference coefficient of friction at sliding velocity V 0 , and a ≪ 1 a power exponent. The real area of contact depends on the effective normal stress and is modulated by changes in the shape of the contact junctions Kilgore 1994, 1996;Sleep 2006) and the quality of contact (Li et al. 2011;Liu and Szlufarska 2012;Renard et al. 2012). Flattening of the micro-asperities increases the real area of contact junctions and increases the strength of the fault. The average contact area is reduced during rapid sliding, thereby weakening the fault. These effects may be captured by a dependence of the effective real area of contact on effective normal stress and a state variable that represents the age of contact where σ is the effective normal stress accounting for the effect of pore pressure, L the characteristic weakening distance, b ≪ 1 a power exponent, and θ the state variable. Combining Eqs. (1) and (2) leads to the multiplicative form of rate-and state-dependent friction in isothermal conditions (Barbot 2019a) ( where the parameters a, b, and L play the same role as in the alternative formulation of Ruina (1983). The parameters µ 0 and V 0 represent reference values such that the frictional resistance at steady-state sliding velocity V 0 equals µ 0σ . The state variable follows an evolution law, also assumed isothermal, that allows healing at stationary contact The spatial distribution of the physical parameters µ 0 , a, b, L, and σ controls the dynamics of fault slip. We assume drained conditions around the fault, such that the spatial distribution of these properties is not affected by heat and fluid transport. With the assumptions of the model, frictional instabilities may only develop when the frictional resistance is velocity-weakening at steady state, which occurs for a − b < 0.
The exact nature of the plate interface at great depths, including the rock type, fault fabric, and the surrounding hydrothermal conditions, is poorly known. However, theoretical and numerical studies of the seismic cycle in a single-asperity fault show that a wide range of styles of seismic cycle can be attributed to a coordinate in a two-dimensional phase space associated with two non-dimensional parameters (Barbot 2019b). The Dieterich-Ruina-Rice number, defined in plane strain condition as where G is the shear modulus, ν Poisson's ratio, and W the width of the velocity-weakening patch, controls the relative importance of elastic interactions, with 0 < R u < ∞ in a velocity-weakening patch. As the R u number represents a ratio of the fault dimension to a characteristic nucleation size (Ruina 1983;Rice and Ruina 1983;Rice 1983;Dieterich 1992Dieterich , 1994, fault slip is stable for R u ≪ 1 , quasi-stable for R u ∼ 1 , and firmly unstable for R u ≫ 1 . Drawing parallels with other fields, the role of the R u number is akin to the Rayleigh number in a Rayleigh-Bénard convective system, due to its association with a mechanical instability (Howard 1966), or with the Reynolds number in fluid dynamics (Horowitz and Ruina 1989), for its association with emerging deterministic chaos at high numbers (Kolmogorov 1962). The second non-dimensional number represents a ratio of static to dynamic stress drop and controls the importance of the local evolutionary effects, with 0 < R b < 1 for a weakening patch. A low positive R b number is associated with conditions approaching velocity-neutral, which are relevant near any velocity-weakening to velocity-strengthening transition associated with increasing temperature. A large R b number approaching 1 from below corresponds to strong weakening, unless a is abnormally large. Values of R b exceeding 1 are possible in case of work-hardening behavior where b is negative (e.g., Blanpied et al. 1998;Sone et al. 2012;Carpenter et al. 2015). An important take-away is that the overall style of the seismic cycle is largely controlled by the pair of nondimensional numbers R u and R b , regardless of how the parameters are combined in Eqs. (5) and (6). Different combinations of frictional parameters with the same R u and R b numbers may lead to different recurrence times or amplitude of slip, but the style of rupture, the fault behavior in the interseismic period, and the recurrence patterns would be similar. A preliminary attempt to map the seismic behavior as a function of coordinates in the two-dimensional phase space formed by R u and R b on semi-infinite and finite faults was performed by Barbot (2019b), ignoring the wave-mediated stress transfers. For intermediate R b numbers, R u ∼ 1 leads to cycles of periodic slow-slip events, R u ∼ 7 , to periodic bilateral ruptures. For R u > 17 , the cycles exhibit sequences of partial and full ruptures that propagate unilaterally. For R u > 50 , mainshock/aftershocks sequences spontaneously emerge. Larger R u numbers lead to a type of deterministic chaos with aperiodic sequences of fast ruptures of varying size that include aftershocks. Cattania (2019) even suggests that Omori's law for aftershocks and Gutenberg-Richter statistics emerge in single-asperity faults with R u > 400 . In contrast, low R b numbers promote chaotic seismic cycles in the slow-slip regime. These non-dimensional parameters are clearly defined for a homogeneous velocity-weakening patch on a planar fault, outside of external perturbations. However, we anticipate that they will also control the fault behavior in more complex settings.
The seismic cycle at the Miyagi section of the Japan trench exhibits some of the signature behavior of a velocity-weakening patch associated with an intermediate R b number and a large R u number (Barbot 2019b). We therefore assign velocity-weakening properties to the megathrust segments of the forearc and the frontal prism. We assume that the down-dip boundary of the forearc basement is associated with a thermal boundary, even though temperature is not explicitly described in the brittle domain. The dimension of the system is well constrained by structural data. We may then obtain cycles of partial and full ruptures of the megathrust for a sufficiently low nucleation size, which is controlled by a combination of the characteristic weakening distance, the effective confining pressure that accounts for the effect of pore pressure, and the parameters controlling the velocity dependence of friction (e.g., Ruina 1983;Rubin and Ampuero 2005). Using the typical parameters a = 1.0 × 10 −2 and b = 1.4 × 10 −2 for the steady-state kinematic parameters in a seismogenic zone (Blanpied et al. 1991;Scholz 1998), a tradeoff remains between the effective confining pressure and the characteristic weakening distance. We select the effective confining pressure σ = 110MPa in the forearc segments, corresponding to near-lithostatic pore pressure, to obtain a recurrence time compatible with historical earthquakes at Miyagi segment. Finally, we choose the characteristic weakening distance L = 1 cm in the forearc basement and L = 5 mm in the forearc upper crust segment to obtain complex sequences of partial and full ruptures. If the forearc basement or the forearc upper crust segments were isolated asperities, the R u numbers associated with these parameters and fault dimensions would be 85 and 70, respectively, with R b = 0.28 . Because these segments are connected, we can consider their average frictional properties and cumulative dimension, leading to R u = 240 . These coordinates in phase space are associated with chaotic sequences of aperiodic partial and full ruptures (Barbot 2019b). In the megathrust segment crossing the mantle wedge corner, we adjust the frictional properties following a similar approach to obtain Mw∼ 7 earthquakes with recurrence patterns similar to the history of deep interplate earthquakes at the Miyagi section, leading to σ = 20MPa, b − a = 2.0 × 10 −3 , L = 5mm. These parameters are associated with R u = 13 , corresponding to cycles of unilateral ruptures. The resulting dynamics of the seismic cycle on the megathrust is expected to be complex based on insights from simplified systems, but we expect unique behaviors due to interacting segments with different frictional properties and the step-wise variations of fault dip from the trench to the brittle-ductile transition.
For the viscoelastic substrate, we adopt a rheological framework that captures the transient creep and the steady-state creep of olivine aggregates. The bulk deformation of any rocks parcel is modeled with a Burgers assembly where the dashpots obey a nonlinear stress/ strain rate relationship (Masuti et al. 2016). The total anelastic strain rate is given by ǫ i =ǫ K +ǫ M , where ǫ K and ǫ M are the instantaneous strain rates in the Kelvin element and the Maxwell element, respectively. The stress/ strain rate relationship in the Maxwell dashpot is given by (Hirth and Kohlstedt 2003) where A is the pre-exponential factor, C OH the water concentration, r the water content exponent, σ the norm of the deviatoric stress, and n the stress exponent. Plastic flow is thermally activated with the activation energy E, the activation volume , the confining pressure p, and the temperature T, with the universal gas constant R. The strain rate in the Kelvin dashpot is a function of the effective stress where ε K is the cumulative strain in the Kelvin element, and G K is a work-hardening coefficient. The cumulative strain in the Kelvin dashpot can be viewed as a state variable for viscoelastic flow. The product 2G K ε K represents the internal stress associated with strain partitioning among several dislocation planes within a representative volume element. The strain rate in the Kelvin dashpot is then given by In general, the constitutive parameters for the Maxwell and Kelvin elements may differ, but because they are poorly known, we use the same values. The rheological properties of dry or hydrous olivine at steady state, particularly the activation energy, depend on the water content (Karato and Jung 2003). The effective viscosity is thermally activated and depends on the ambient temperature in the lithosphere-asthenosphere system.
We build a simplified two-dimensional thermal model for the Japan trench. In the Pacific plate, heat transport is accommodated not only by conduction and horizontal advection below the ocean basin, but also by downward advection past the trench. We base the thermal model on a cooling half-space with a basal temperature of 1380 • C and a thermal age of 120 million years (Myr) where the top of the oceanic crust is used as the reference depth. The young thermal age is to compensate for the misfit of the cooling half-space model with heat flow and bathymetry data for ocean basins older than 100 Myr (Stein and Stein 1992). For the mantle wedge, we base the thermal model on a cooling halfspace with a thermal age of 53 Myr and a basal temperature of 1380 • C. However, we reduce the temperature exponentially by up to 900 • C as a function of proximity to the plate interface. The resulting thermal model is shown in Fig. 2. The different temperature profiles across the plate boundary represent the dominant rheological contrast in the model. In addition, we assume a gradient of fluid content in the mantle wedge, with a large water concentration of C OH = 4000 ppm H/ Si in the mantle wedge corner below the volcanic arc and C OH = 1000 ppm H/Si farther away. The activation of viscoelastic flow depends on the water content, as water fugacity is also thermally activated (Karato and Jung 2003). This implies an associated activation energy of E = 490 kJ/mol in the mantle wedge corner and E = 510 kJ/mol below the Sea of Japan. We use a water content of C OH = 1000 ppm H/Si in the Pacific plate asthenosphere, a value also inferred in the Indian Ocean asthenosphere based on analysis of postseismic data (Masuti et al. 2016). The viscoelastic flow is driven by a background anelastic strain rate of horizontal uniaxial shortening of ǫ 22 = −10 −14 /s, where x 2 is the trench-perpendicular coordinate, compatible with plate convergence distributed over 260 km, encompassing the lithosphere-asthenosphere system.

Numerical solution with the integral method
We simulate the dynamics of fault slip and viscoelastic flow during the seismic cycle with the integral method (Lambert and Barbot 2016;Barbot 2018b;Shi et al. 2020). The approach extends the boundary integral method to include both surface and volume elements. As anelastic strain accrues in the viscoelastic domain, the change of fault traction depends on the surrounding fault slip and distributed strain. Accordingly, the rate of change of shear traction along the fault can be written as where K is the stress kernel for slip on the megathrust surface ∂Ŵ , v is the instantaneous slip velocity vector, v L is the local loading rate on the fault, ǫ represents the anelastic strain rate tensor (Andrews 1976; Barbot and Fialko 2010) in the ductile domain , and J is another stress kernel connecting distributed deformation to fault traction. The last term on the right-hand side of Eq. (10) corresponds to radiation-damping with the rigidity G, shear wave velocity V S , and slip acceleration V (Rice 1993). Similarly, the evolution of stress in the ductile region can be written as where M is stress kernel associated with viscoelastic flow and L another stress kernel associated with fault slip. The dynamics of deformation includes four directions of (10) coupling captured by the stress kernels M for fault slip self-interactions, J for strain reloading of the fault, M for viscoelastic relaxation, and L for post-seismic deformation. The formulation assumes a background slip rate v L and strain rate ǫ L , both possibly non-uniform, at which no stress is accumulating. Combining the constitutive laws for frictional resistance and viscoelastic flow with the above equations forms a closed system that can be evaluated numerically. The integral equations (10) and (11) are approximated with finite sums over surface elements and volume elements. The traction and stress in a half-space caused by fault slip and distributed strain are obtained with the solutions of Okada (1992) and Barbot (2018a), respectively. The time stepping is then evaluated numerically using the fourth/fifth-order Runge-Kutta method that also provides adaptive time steps (Press et al. 1992). The calculations are based on the Unicycle software (Barbot 2018b(Barbot , 2019b, which has been partially validated against other numerical methods using two-dimensional benchmarks in anti-plane strain conditions (Erickson 2020), part of a community effort to develop and share reliable software for crustal dynamics.

Numerical simulations of seismic super-cycles at the Japan trench
We simulate the seismic cycle in viscoelastic media using the integral method, combining surface elements and volume elements to capture the coupling between fault dynamics and mantle flow. We use a two-dimensional cross-section and model the deformation within the in-plane strain approximation (Barbot 2018b). Threedimensional models would be more appropriate, allowing us to capture the trench-parallel structure. However, modeling fault dynamics at high R u number is a challenging task with current numerical methods and including viscoelastic flow increases the computational burden further more. We therefore proceed with a two-dimensional approximation to capture the first-order behavior of the system. The computational framework is numerically efficient because the stress interactions through elastic coupling are based on analytic solutions Barbot 2018a). The model makes a number of additional approximations for numerical convenience, such as a planar Earth's surface, uniform elastic moduli, and isothermal conditions. Yet, evaluating 1000 years of seismic cycles, which encompasses more than one hundred seismic events, still requires a minimum of 15 h of computational time using a single-node 16-core computer, and every simulation is followed by additional computationally intensive data analysis.
We approximate the megathrust geometry from Slab2 (Hayes et al. 2018) by piecewise planar segments and we discretize the mega-splay geometry based on seismic imaging . We discretize the megathrust and the mega-splay fault with planar patches of 200 m width. For the viscoelastic domain, we use triangular cross-sectional elements with a circumsphere radius of about 10 km forming an unstructured mesh (Barbot 2018a, b). We consider a Poisson solid with a rigidity of 30 GPa. Building up on the considerations described in the previous sections, we further optimize the frictional and viscoelastic properties by trial and error within existing constraints from laboratory friction experiments and mineral physics (e.g., Bürgmann and Dresen 2008, and references therein) in dozens of numerical simulations to explain the recurrence pattern of interplate earthquakes at the Miyagi section and the geodetic data associated with the 2011 Tohoku earthquake. Simulations at high R u number in a different tectonic context which also showcase sequences of partial and full ruptures indicate that the characteristics of the first few events, part of the first super-cycle, differ from the subsequent sequences, due to bias from the initial conditions (Erickson 2020). We therefore exclude from the analysis the entire first super-cycle, which in our case represents the first 100 or so ruptures, depending on the model parameters considered. We compare the model predictions for a 3000-year period, as considering longer simulations within an optimization scheme would be impractical. We do observe slightly longer recurrence times at longer time scales, as the system behavior evolves slowly over many supercycles, but we have not observed qualitatively different rupture behaviors in the second and following supercycles. In the following, we only present the result of the preferred model, based on the distribution of physical parameters summarized in Table 1 for the brittle section and in Table 2 for the viscoelastic substrate.
We first present the dynamics of fault slip and viscoelastic flow for a period of 30,000 years, representing close to 2500 earthquakes of varying size (Fig. 3). Seismic events are defined as episodes during which the peak velocity has exceeded an initial threshold of 1 cm/s and then maintained a peak velocity above 0.1 mm/s. The peak velocity represents the maximum slip speed attained anywhere in the model at a given time. The simulation features seismic super-cycles, defined as the period separating two giant ruptures that bookend small-size ruptures. During the last 10,000 years of the simulation, the super-cycles span from 1058 to 1344 years, encompassing from 5 to 7 smaller single-or multiple-segment breaking ruptures and overall 100 seismic events of varying size, location, and source properties. Over the period considered, the model never converges to a limit cycle, compatible with the emergence of deterministic chaos with sequences of mainshocks and aftershocks of varying size at high R u number (Barbot 2019b;Cattania 2019). Although our assumptions of continuously velocity-weakening properties from the arc Moho to the trench make the through-going ruptures not surprising, it may be less intuitive that the same region occasionally breaks in smaller ruptures. However, this

Table 1 Summary of the physical parameters controlling fault dynamics
The model assumes uniform elastic properties with a shear modulus of G = 30 GPa, Poisson's ratio ν = 0.25 , and shear wave speed V S = 3 × 10 3 m/s. The reference slip velocity is V 0 = 10 −6 m/s behavior is typical of where the physical parameters fall in the two-dimensional phase space formed by R u and R b , despite the additional complexity brought by varying frictional properties, fault geometry, and coupling with viscoelastic flow. Some other characteristic features include long-term creep episodes that develop in the velocity-weakening region after partial ruptures of the forearc, near the tip of a previous rupture or across a fault bend at a segment boundary. The hypocenters of the through-going ruptures occur sometimes at the bottom of the forearc or somewhere in the middle of the velocity-weakening region. The megathrust below the lower arc is systematically creeping during the interseismic period. Because they are structurally bounded, i.e., surrounded by permanent velocity-strengthening regions, the deep ruptures below the mantle wedge corner exhibit more regularity (Fig. 4a). When the deep earthquakes are excluded from the 30,000-year catalog, the size of the remaining events approximately follows a Pareto distribution, a type of power-law probability distribution with the probability F of the moment per unit length of an event to exceed M given by F = (M l /M) α if M ≥ M l and F = 1 otherwise. The Pareto distribution is related to the Gutenberg-Richter distribution of earthquake size, but the two-dimensional approximation prevents us from estimating moment magnitudes and their associated b value. Comparing the Pareto distribution to a cumulative normalized histogram of event sizes, we find M l = 5 × 10 7 N and α = 0.4 by manual adjustment (Fig. 4b). The large variability of rupture size in the model can be attributed to the large effective R u number, the geometrical complexity, and the rheological gradients. The source properties of simulated events also feature a great variety of source properties, including slow and fast ruptures. The partial and full ruptures of the forearc region seem to follow the expected square-root scaling between moment per unit length and duration, given the two-dimensional approximation (Fig. 4b). In contrast,

Table 2 Summary of the physical parameters for viscoelastic flow
The thermal model is based on a cooling half-space with a thermal age a and a basal temperature T 0 , where the top of the slab is the reference depth, to account for advection. The mantle wedge is cooled as a function of proximity to the down-going slab, following the perturbation �T e −z/z0 , where T = −900 • C, z is the distance to the slab, and z 0 = 25km  . 3 Numerical simulations of 8000 years of seismic super-cycles at the Japan trench (Miyagi section) that was preceded by a 22,000-year period. a Seismic cycle from 22,000 to 24,000 years. b Period from 24,000 to 26,000 years. c Period 26,000 to 28,000 years. d Period 28,000 to 30,000 years. Super-cycles (thin dashed lines) are defined as the period separating two full ruptures of the forearc. Long-term creep in the velocity-weakening domain may take place near fault bends or near the rupture tip of previous partial ruptures. The domain boundaries are shown by the white dashed lines. The stars indicate the hypocenter and magnitude (moment per unit length) of the simulated ruptures. The event number indicates simulated seismic events during which the peak velocity has exceeded an initial threshold of 1 cm/s and then maintained a peak velocity above 0.1 mm/s slow earthquakes, and deep, slow aftershocks, both with a peak velocity lower than 10 cm/s, exhibit a different scaling. Here, slow aftershocks are defined as slow earthquakes that occur in a postseismic sequence. Within the period simulated, the recurrence pattern of seismic ruptures is aperiodic and the spectrum of rupture size, in terms of moment per unit length, covers five orders of magnitude.
We then describe the dynamics of the seismic cycle for the same physical parameters considering a smaller simulation period (Fig. 5), which allows us to discuss key similarities with geophysical data at the Japan trench. We exclude from the analysis the entire first super-cycle, which in this case represents the first 87 ruptures. Over the remaining period considered, from 560 to 3000 years, the super-cycles differ in the details, recurring aperiodically between 520 and 800 years, i.e., slightly shorter than after more super-cycles have taken place (Fig. 3). The super-cycles share similar broad characteristics, as they all include quasi-periodic ruptures of the mantle wedge corner every about 20 years, partial ruptures of the forearc separated by about 100 years of quiescence at the beginning of the super-cycle, and full ruptures of the megathrust from the trench to the mantle wedge corner. The moments per unit length of simulated events cover more than five orders of magnitude, from 8.7 × 10 5 N to 2.1 × 10 11 N. Every large rupture in the forearc, partial or full, is followed by aftershocks nucleating down-dip, and by a long-lived postseismic transient that includes afterslip on the megathrust and viscoelastic flow in the Fig. 4 Scaling relationships and statistical properties of simulated ruptures within seismic super-cycles at the Japan Trench. a Histogram of event moment per unit length for a simulation duration of 30,000 years that includes more than 2500 events. The recurring interplate ruptures below the mantle wedge (labeled as deep events) repeat with more regularity than other events. b Normalized cumulative histogram for all events, and for a catalog that excludes the deep earthquakes showing the probability of events to exceed the moment per unit length of the absissa. The latter approximately follows a Pareto distribution (dashed line) with a power exponent α = 0.4 for moment per unit length above 5 × 10 7 N. c Moment per unit length versus duration of simulated ruptures. The peak velocity corresponds to the maximum slip velocity reached at any point and time during a rupture. Slow earthquakes and deep aftershock follow a different scaling from the partial and full ruptures of the forearc. Linear and square-root scaling are shown for reference Fig. 5 Dynamics of the Japan Trench during a seismic super-cycle. a, b Evolution of slip velocity down-dip of the mega-splay fault and the megathrust, respectively, for 1000 years. The horizontal dashed profiles indicate segment boundaries. The stars indicate the hypocenters of earthquake nucleation, colored by the peak slip velocity attained during the event. Partial ruptures are confined to the forearc basement. During giant ruptures, slip propagates to the trench and to the mega-splay fault. Afterslip follows all ruptures. The cluster of events before and after giant ruptures corresponds to foreshocks and aftershocks, respectively. c Evolution of slip velocity at the center of each velocity-weakening segment and peak velocity anywhere along the megathrust/mega-splay system. The recurrence time of ruptures in the mantle-wedge corner increases during a super-cycle. d Peak strain rate in the oceanic asthenosphere and in the mantle wedge. e Sequence of full and partial ruptures, foreshocks and aftershocks, and style of rupture. Each event is indicated by a star colored by the peak slip velocity attained during rupture propagation (See figure on next page.) Barbot Earth, Planets and Space (2020) 72:63 oceanic asthenosphere and in the mantle wedge ( Fig. 5bd). The partial ruptures in the forearc and in the mantlewedge corner are also followed by afterslip, compatible with observations of afterslip following the 2005 Miyagi-Oki earthquake (Miura et al. 2006). The recurrence time of the deep ruptures is modulated by the giant ruptures and their postseismic relaxation, with measurably shorter recurrence intervals early into a super-cycle, and longer ones later on. The short recurrence time of Mw∼ 7 earthquake after a giant rupture of the Japan trench is compatible with the more sophisticated three-dimensional simulations of Nakata et al. (2016), illustrating that first-order features may be captured in a two-dimensional model. While the ruptures of the mantle-wedge corner all spread over similar distances, the slip distributions of the forearc ruptures show larger variability (Fig. 6b). While some forearc ruptures may propagate through the upper crust, only the giant through-going ruptures break through the frontal prism and reach the trench. The mega-splay fault only breaks during these full ruptures, which occurs rarely, and the mega-splay ruptures do not reach the Earth's surface, compatible with findings from seismic reflection imaging of the frontal prism (Kodaira et al. 2012). During the interseismic  (Iinuma 2012). c Surface horizontal displacements induced by two simulated giant ruptures 2137.2 (black profile) and 2656.9 (red profile) years after the start of the simulation and geodetic coseismic displacements within 50 km of profile A-A' in Miyagi for the Tohoku mega-quake. The on-land measurements originate from GEONET (Sagiya et al. 2000;Ozawa et al. 2011). The offshore measurements include GPS-acoustic and ocean bottom pressure gauges data summarized by Iinuma (2012). d Same as above panel, but for the vertical displacement period, the megathrust underneath the frontal prism is always locked. We focus our attention on the two giant ruptures that occur 2137 and 2656 years after the beginning of the simulation, representing the 257th and 313rd events, respectively. Before 2137 years, 3 through-going ruptures separated by two super-cycles already took place, in addition to 7 other ruptures that spread over the entire forearc region, and another 7 that unzipped the entire forearc basement. The remaining 239 ruptures represent aftershocks and small isolated events. Even though we do not expect any two events to be exactly the same in Fig. 7 Fault dynamics before, during, and after the giant rupture that occurs after 2656.9 simulated years. The few days before the mainshock exhibit acceleration of slow slip above the background convergence rate. A few slow-earthquake precursors also precede the mainshock. After the giant rupture, the forearc segments relock, aftershocks take place immediately up-dip of the lower arc and in the mantle wedge corner, and afterslip propagates in the lower arc, cold nose, and upper mantle segments. The peak velocity is the maximum velocity reached anywhere on the fault at a given moment the simulation, we assume that these two giant ruptures form a representative sample of the model outcome. The slip distributions of these events showcase up to 65 m and 45 m of slip near the trench, respectively (Fig. 6), and show some similitudes with geodetically inferred coseismic slip distributions for the 2011 Mw = 9.1 Tohoku earthquake (e.g., Iinuma 2012), even though the latter suffers large uncertainties (e.g., Sun et al. 2017). The amplitude of shallow slip for these events falls within the range of various estimates based on differential bathymetry and seafloor geodesy Sato et al. 2011;Ito et al. 2011;Sun et al. 2017;Fujiwara et al. 2017). The surface displacements associated with the simulated coseismic slip explain the geodetic observations for the Tohoku earthquake within 50 km of profile A-A' , both on land and offshore, in the horizontal and vertical directions (Fig. 6c, d).
The through-going rupture of the simulated year 2656 (event 313) shows some other striking similarities with the Tohoku mega-quake. For example, this event is preceded by two simultaneous long-term creep episodes of increasing size, one close to the epicenter and another at the boundary between the forearc upper crust and forearc basement (Fig. 5b), and by two large earthquakes in the previous decade, including one in the mantle wedge corner (Fig. 7a). The thin strip of interseismic creep aligned with the tip of a previous rupture eventually becomes the nucleation point of the giant rupture. The nucleation was also preceded by a foreshock sequence close to the epicentral region in the preceding hours and a few aftershocks in the following days (Fig. 7b). We note that the hypocenter location is deeper than for the Tohoku earthquake, but the simulation features other giant ruptures with a shallow nucleation point, for example, the full rupture of event 257 (Fig. 5b). Therefore, a shallow nucleation of giant ruptures is allowed in the model. The rise of interseismic creep far into a velocity-weakening region can occur due to the residual stress left by previous ruptures or due to spatial variations in fault geometry that affect the loading rate. Both effects have been shown in seismic cycle simulations on the Ventura thrust system, California (Ong et al. 2019) and in numerical models of earthquake cycles in fault-bend folds (Sathiakumar et al. 2020). However, interseismic creep in a velocity-weakening region is more prominent at fault bends, possibly occurring during the entire interseismic period of a super-cycle, depending on the change of dip at the bend (Sathiakumar et al. 2020).
The inclusion of the mechanical coupling between brittle and ductile deformation in the model allows us to investigate the postseismic relaxation following a giant rupture (Fig. 8). We compare the displacements of the giant event of the simulated year 2656 (event 313) with the deformation that followed the Tohoku earthquake between April 2011 and December 2011, i.e., between 21 and 265 days after the mainshock. The model explains the wide range of forward and retrograde motion and the spatial distribution of postseismic horizontal displacements, including on land and offshore, except for short-wavelength displacements centered around the volcanic arc, as the model ignores viscoelastic flow in the lower arc. The postseismic deformation is explained by a large amplitude afterslip contribution with a dynamic range of about 60 cm and a long-wavelength viscoelastic deformation with a dynamic range of 10 cm. Both mechanisms of deformation therefore contribute to retrograde motion in the forearc, compatible with seafloor geodetic data (Watanabe et al. 2014;Sun 2014). The role of viscoelastic relaxation in the oceanic lithosphere on the postseismic deformation in the forearc is well understood (Suito 2017;Noda et al. 2018;Barbot 2018b). That afterslip on the megathrust would also create retrograde motion at the surface is less intuitive. Indeed, deep afterslip on a planar megathrust would only create forward motion in the hanging wall. However, as the megathrust incrementally bends from shallow to deep segments, the deep afterslip occurs on a plane that projects halfway between the trench and the shelf boundary at the Earth's surface. The direction of motion caused by deep afterslip below this imaginary plane is compatible with displacements in the footwall, thereby explaining these observations. A vertical cross-section of the displacement field across the plate boundary due to deep afterslip that explains these results is also shown by Agata et al. (2019).
The simulated postseismic displacement is strongly affected by the aftershock sequence that rapidly follows the mainshock (Fig. 5b). Identifying the contribution of these aftershocks on the postseismic relaxation (See figure on next page.) Fig. 8 Postseismic surface displacements in the year following the mainshock. a Geodetic horizontal displacements that accumulated between April 2011 and December 2011, i.e., between 21 and 265 days after the Tohoku earthquake (Muto et al. 2019). Subsidence of the forearc is more accentuated in subsequent years (Watanabe et al. 2014). The black profile shows the modeled displacement for the event that occurs after 2656.9 simulated years. The red and blue profiles show the relative contributions of fault slip and viscoelastic deformation. The simulated displacements include the effects of 4 aftershocks that occurred during that time. b Simulated horizontal displacement (solid profile), surface displacement when the coseismic surface displacements of aftershocks have been removed (dashed profile), and surface displacement when the coseismic and immediate afterslip displacement of aftershocks have been removed (dot profile). c Same as above panel but for the vertical displacements is challenging within the two-dimensional approximation. Indeed, all ruptures assume an infinite length in the trench-parallel direction, giving small events a disproportionate impact on surface displacements. Since aftershocks are impulsive, they can be excluded from the postseismic displacements, which does not greatly impact the quality of fit to the geodetic data (Fig. 8b). However, aftershocks also trigger a postseismic relaxation of their own that contributes significant postseismic displacement, i.e., more than they would in a three-dimensional model, particularly in the vertical direction (Fig. 8c). We conclude that two-dimensional models of postseismic deformation within seismic cycles that include aftershocks are not ideally suited to constrain the viscoelastic properties of the asthenosphere. Other approaches based on geodetic imaging (Tsang et al. 2016;Moore 2017;Qiu et al. 2018;Tang et al. 1999;Weiss 2019), or forward modeling of the postseismic phase alone (e.g., Rousset et al. 2012;Sun 2014;Hu et al. 2016;Freed et al. 2017;Agata et al. 2019;Muto et al. 2019) may be more befitting.

Conclusions
The structural layout of the Japan trench and the historical seismicity offshore Miyagi in the last centuries can be reconciled in a consistent rheological model of the subduction zone (Fig. 9). With the exception of a thin segment in the arc lower crust, the megathrust may be characterized by unstable weakening friction in a wide region spanning from the trench to the coast, allowing seismic super-cycles. Down-dip variations in frictional properties on the megathrust are due to the varying fabric of the hanging wall above the plate interface and systematic changes in temperature and pressure that affect metamorphism and the dominant deformation mechanism.
The mantle wedge corner delineates a segment of the megathrust that breaks regularly in Mw∼ 7 earthquakes Fig. 9 Conceptual rheological model of the Japan trench at Miyagi. a The down-going Pacific slab intersects different terranes of the Japanese arc at various dip angles, introducing rheological and morphological gradients that shape the seismic super-cycle. The megathrust occupies the brittle regime of the plate interface, stopping down-dip at the brittle-ductile transition. In the mantle wedge, the coupling between the two plates is permanent, entraining corner flow below the volcanic arc. Fluids enter the mantle wedge underneath Honshu from the Pacific lithosphere and the Sea of Japan asthenosphere. b The megathrust is segmented into the weak frontal prism and the forearc upper crust, the forearc basement, the mantle-wedge corner, cold nose, and upper mantle. The segments are named after the hanging wall block that they underlie. The mega-splay fault delineates the active frontal prism. Large Mw∼ 7 earthquakes are mostly confined in a single segment, with rupture nucleation close to structural boundaries. Giant Mw∼ 9 ruptures break the entire offshore section, reaching the trench. Megathrust segments in the upper mantle exhibit stable creep. The cross-section area in b) is shown as a dashed box in the top panel with a recurrence period that increases during a seismic super-cycle. These isolated ruptures take place because of a permanent barrier in the arc lower crust above the arc Moho. In the forearc basement and the shallower regions, the rupture sequence is more complex because of the intrinsically nonlinear dynamics of fault slip when the unstable area is much larger than the nucleation size. In these conditions, a hallmark of the seismic cycle is a sequence of partial and full ruptures of the velocity-weakening region, even within a single velocity-weakening asperity. This regime of frictional parameters also leads to a wide range of earthquake sizes, making the sequence aperiodic. The full ruptures correspond to trench-breaking events. The presence of active faults and folds in the frontal prism facilitates efficient uplift of the seafloor, such that shallow ruptures are also tsunamigenic (e.g., Ma 2012;Lotto et al. 2017). The forearc partial ruptures vary in size from one super-cycle to another and even within the same super-cycle because they are not structurally controlled. Instead, their spatial extent is dictated by concentrated stress residuals from previous ruptures, which is a dynamic feature. The bending of the subducting slab below the forearc and the stress concentrations at previous rupture tips promote long phases of aseismic creep in the middle of locked patches, sometimes next to the epicenter, as was observed before the Tohoku earthquake (Ito 2013). The dynamics of fault slip at the Miyagi segment of the Japan trench is therefore both structurally and frictionally controlled. The down-dip segmentation is due, first, to a large velocity-weakening area stretching from the arc Moho to the trench combined with a small nucleation size that makes the mechanical system fall into a chaotic regime of partial and full ruptures and, second, to the presence of a permanent structural boundary around the arc Moho, the arc lower crust, that allows Mw∼ 7 interplate earthquakes to frequently repeat below the mantle wedge corner.
This study demonstrates how to reconcile structural and seismological data in northeast Japan at the Miyagi section without invoking thermal pressurization of the shallow sediments during rapid slip (e.g., Mitsui et al. 2012;Noda and Lapusta 2013;Cubas et al. 2015), evolving frictional properties with sliding velocity , nor tying all historical ruptures with permanent asperities (e.g., Hori and Miyazaki 2011;Kato et al. 2011). These possible end-member models predict different interseismic behavior near the trench. Our model implies a systematically locked frontal prism during the interseismic period, but this cannot be verified with current seafloor geodesy data (Yokota et al. 2016;Tomita et al. 2017), because of insufficient resolution near the trench (Sathiakumar et al. 2017). We concede that a shallow velocity-strengthening segment below the frontal prism is not entirely excluded, but it should be weak enough to allow trench-breaking ruptures with substantial coseismic slip. Furthermore, it is possible that the frictional behavior near the trench may not be well described by considering solely the velocity dependence of fault strength if the material is temperature weakening (Chester 1995;Blanpied et al. 1995;Chen et al. 2015). More comprehensive models that include a possible positive feedback between shear heating and a temperaturedependent frictional strength should be explored in future work.
Our two-dimensional model provides a first-order explanation for the size, recurrence time, and hypocenter location of earthquakes in the Miyagi segment during the last century, for the precursory creep and foreshocks activity of the 2011 Tohoku earthquake, and for several aspects of its coseismic and postseismic deformation. However, the Japan trench exhibits large trenchparallel changes in seafloor morphology and hanging wall fabric (e.g., Bassett et al. 2016;Nakajima et al. 2001;Hasegawa et al. 2005;Zhao et al. 2015;Wang et al. 2017;Uchida et al. 2016;Nishikawa et al. 2019), implying additional segmentation of the megathrust along the trench axis. Incorporating both along-strike and down-dip segmentation in three-dimensional numerical simulations of the seismic cycle may open the door to more comprehensive models of the seismic cycle at the Japan trench. However, these simulations need to resolve large unstable fault regions and small nucleation sizes to capture the complex dynamics that emerges in this numerically challenging parameter regime. Major improvements in methodology are needed to make these simulations more numerically efficient and widely available across the scientific community.