Structural control and system-level behavior of the seismic cycle at the Nankai Trough

The Nankai Trough in Southwest Japan exhibits a wide spectrum of fault slip, with long-term and short-term slow-slip events, slow and fast earthquakes, all associated with different segments down the plate interface. Frictional and viscous properties vary depending on rock type, temperature, and pressure. However, what controls the down-dip segmentation of the Nankai subduction zone megathrust and how the different domains of the subduction zone interact during the seismic cycle remains unclear. Here, we model a representative cross-section of the Nankai subduction zone offshore Shikoku Island where the frictional behavior is dictated by the structure and composition of the overriding plate. The intersections of the megathrust with the accretionary prism, arc crust, metamorphic belt, and upper mantle down to the asthenosphere constitute important domain boundaries that shape the characteristics of the seismic cycle. The mechanical interactions between neighboring fault segments and the impact from the long-term viscoelastic flow strongly modulate the recurrence pattern of earthquakes and slow-slip events. Afterslip penetrates down-dip and up-dip into slow-slip regions, leading to accelerated slow-slip cycles at depth and long-lasting creep waves in the accretionary prism. The trench-ward migrating locking boundary near the bottom of the seismogenic zone progressively increases the size of long-term slow-slip events during the interseismic period. Fault dynamics is complex and potentially tsunami-genic in the accretionary region due to low friction, off-fault deformation, and coupling with the seismogenic zone.


Introduction
Subduction megathrusts accumulate and release mechanical energy through a variety of slip modes, including earthquakes, low-frequency earthquakes, slow-slip events, and stable creep. These different types of events can be discriminated in seismo-geodetic observations by their distinct peak slip velocities, frequency spectra, durations, and recurrence intervals (e.g., Wallace and Beavan 2006;Ito et al. 2007;Araki 2017). Slow earthquakes involve slip velocities and rupture speeds fast enough to be detected by broadband seismometers, but not fast enough to generate high-frequency seismic waves (Barbot 2019b). Slow-slip events are too slow to be detected by seismometers due to their low frequency content, but may be monitored with geodetic instruments. Hence, in this paper, we refer to slow earthquakes as any of low-frequency earthquakes, very-low-frequency earthquakes and non-volcanic tremors, but exclude slowslip events from this denomination. In contrast, we refer to slow-slip event as either long-term or short-term slowslip events. At subduction zones, each slip mode is commonly associated with distinct down-dip segments of the megathrust (Oleskevich et al. 1999;Rogers and Dragert 2003;Peng and Gomberg 2010;Hirose and Obara 2005;Obara and Kato 2016;Gao and Wang 2017). Structural control also operates along strike in various ways (e.g., Kodaira et al. 2006;Sieh 2008;Contreras-Reyes et al. 2010;Contreras-Reyes and Carrizo 2011;Moreno 2012;Lin 2013;Metois et al. 2016;Bassett et al. 2016;Nocquet 2017;Philibosian 2017;Vaca et al. 2018;Van Rijsingen et al. 2018;Nishikawa et al. 2019), but the down-dip segmentation of slip behavior seems to be an ubiquitous feature of subduction megathrusts. The Nankai Trough offers a striking example of down-dip segmentation, with shallow low-frequency and very-low-frequency earthquakes, large megathrust earthquakes, long-term and short-term slow-slip events, the latter accompanied with deep low-frequency tremors (Fig. 1).
This study represents an attempt to reconcile geological and seismo-geodetic observations at the Nankai Trough, including insights from other subduction zones. We first review the general characteristics of subduction megathrusts, in particular, the seismic potential and source characteristics at each structural segment. The shallowest fault segment, overlain by the accretionary prism, was once thought to undergo stable sliding due to the low seismicity near the trench late in the seismic cycle and the presence of poorly consolidated sediments. However, accretionary prisms are known to host tsunami earthquakes (Kanamori 1972;Pelayo and Wiens 1992;Satake and Tanioka 1999;Bilek and Lay 2002;Geersen 2019), for example, along the Sunda arc (Ammon et al. 2006;Kanamori et al. 2010;Lay et al. 2011;Newman et al. 2011;Bilek et al. 2011;Satake et al. 2013), the Hikurangi subduction zone (Bell et al. 2014), Nicaragua (Kanamori and Kikuchi 1993;Satake 1994), the Aleutian Islands (Johnson and Satake 1997), the Peru-Chile trench (Pelayo and Wiens 1990), and the Kuril-Kamchatka trench (Kanamori 1972). The rupture of giant earthquakes also propagates to the trench, for example, the 1960 Valdivia, Chile, the 1964 Alaska, the 2004 Aceh-Andaman, the 2010 Maule, Chile, and the 2011 Tohoku-Oki, Japan Deep slow-slip events are located down-dip extension of the seismogenic zone. Deep tremors are denoted by dark red dots clustered within a belt-like region down-dip to the deep slow-slip events. The 20-km-interval contours of the depth from 0 to 100 km (Hayes et al. 2018) and trench (thick black toothed lines) depict the geometry of the Philippine Sea slab that is subducting in the direction represented by the black arrow. P-P' represents a reference profile for the two-dimensional model. The Median Tectonic Line (thick solid line) (Kawamura et al. 2003;Ito et al. 2009) and the landward limit of the accretionary prism (solid-dashed line) (Leggett et al. 1985;Kodaira et al. 2000;Takahashi et al. 2003) are the surface expression of the material boundaries shown in Fig. 11. The spatial distribution of deep tremors is from Maeda and Obara (2009) and Obara et al. (2010); great earthquakes, slow-slip events and shallow low-frequency earthquakes from Obara and Kato (2016) earthquakes (Plafker 1965;Ishii et al. 2005;Lorenzo-Martín et al. 2006;Yue et al. 2014;Tomita et al. 2017), with the notable exception of the 2005 Nias, Sumatra earthquake (Hsu et al. 2006;Konca et al. 2007). Several mechanisms have been proposed to explain the failure of the shallow megathrust in giant ruptures. Noda and Lapusta (2010) use thermal pressurization to facilitate dynamic weakening in a creeping segment. Kozdon and Dunham (2013) invoke wave-mediated propagation of rupture into a shallow velocity-strengthening segment. In these models, seismic ruptures nucleate in the seismogenic zone and propagate up-dip to the velocity-strengthening segment. However, discoveries of low-frequency earthquakes and tectonic tremors within the shallow fault segment underlying the accretionary prism, such as in the Nankai Trough (Sugioka 2012;Yamashita 2015;Toh et al. 2018;Nakano et al. 2018), Costa Rica (Jiang et al. 2012;Dixon et al. 2014) and New Zealand (Wallace and Beavan 2006;Wallace 2016;Wallace et al. 2017), indicate that slip instabilities can nucleate at shallow depth. Laboratory experiments on clay minerals indicate a low static friction of accretionary sediments and a complex temperature and velocity dependence of dynamic friction (Saffer and Marone 2003;Ikari 2019) that includes velocity-weakening behavior at tectonic loading rates (Ikari and Kopf 2017). The seismogenic zone is a major section of the megathrust overlain by the continental crust or arc crust, depending on the tectonic setting, that generates large earthquakes (e.g., Pacheco et al. 1993;Heuret et al. 2011;Hayes et al. 2012). The down-dip limit of the seismogenic zone is thought to coincide with the upper-plate Moho or a shallower isotherm. The frictional properties of quartzrich crystalline rocks explain the firmly velocity-weakening properties of this region (Scholz et al. 1972;Tullis and Weeks 1986;Kilgore et al. 1993;Blanpied et al. 1991Blanpied et al. , 1995Mitchell et al. 2013). The low-angle intersection of the megathrust with the upper-plate crust at subduction zones forms the source region for Earth's largest earthquakes. The oblique convergence at the trench is often accommodated by strike-slip faulting in the forearc, for example, the Median Tectonic Line in Southwest Japan, the Sumatran Fault in Indonesia, the Linquine-Ofqui fault system in South Chile, and the Fairweather and Denali faults in Alaska (Schellart and Rawlinson 2013). Earthquake source properties vary among subduction zones (e.g., Ye et al. 2016;Hayes 2017;Meier et al. 2017;Ye et al. 2018), but a remaining puzzling behavior is the full and partial ruptures of the seismogenic zone in supercycles. For example, coral, turbidite, and paleo-tsunami records indicate a wide range of earthquake sizes in the Sunda trench (Sieh 2008;Patton 2015;Philibosian 2017;Rubin et al. 2017). In the Japan trench, only earthquakes of moment magnitude 7 preceded the 2011 Tohoku megaquake in the last century (Yamanaka and Kikuchi 2004;Nakata et al. 2016). Partial rupture of the seismogenic zone also occurs in a continental collision setting, for example, along the Main Himalayan front (Sapkota et al. 2013;Bollinger et al. 2013;Bilham 2015;Mencin 2016;Hubbard et al. 2016;Feng 2017). Herrendörfer et al. (2015), Michel et al. (2017), and Barbot (2019b) argue that the super-cycle is a self-emergent behavior from the nonlinear dynamics of fault slip favored by wide seismogenic zones. Qiu et al. (2016) and Hubbard et al. (2016) suggest the important role of morphological gradients, i.e., local changes in the geometry of the plate interface. How the earthquake recurrence and super-cycle patterns are affected by surrounding activity, notably slow-slip events and asthenospheric viscoelastic flow is another important unknown.
At most subduction zones, the seismogenic zone is followed down-dip by a stability transition where the megathrust slips in slow, dominantly aseismic episodes (e.g., Rogers and Dragert 2003;Schwartz and Rokosky 2007;Szeliga et al. 2008;Gomberg and Group 2010;Beroza and Ide 2011;Kobayashi 2014;Obara and Kato 2016;Gao and Wang 2017;Bürgmann 2018). Slow-slip events are thought to occur in a narrow metamorphic belt of serpentinized rocks where the development of frictional instabilities is inhibited by high temperature and high pore-fluid pressure (Kodaira et al. 2004;Katayama et al. 2013;Audet and Bürgmann 2014;Frank et al. 2015;Goswami and Barbot 2018;Audet and Schaeffer 2018;Nakajima and Uchida 2018;Cruz-Atienza et al. 2018). The long-term slow-slip events, characterized by durations of months or years and recurrence time of several years (Kobayashi 2014), take place immediately underneath the seismogenic zone, in the mantle wedge corner, representing either slip on the megathrust or distributed flow in a narrow metamorphic shear zone (Platt et al. 2018). In contrast, the short-term slow-slip events, with durations of days or weeks and recurrence times of several months Obara 2005, 2006;Obara et al. 2004), take place in the upper mantle and are accompanied by substantial tremors and low-frequency earthquakes (Obara 2002;Obara et al. 2004;Shelly et al. 2006;Ito et al. 2007), a phenomenon referred to as episodic tremor and slip. Examples include the Cascadia (Rogers and Dragert 2003) and Nankai (Obara and Kato 2016) subduction zones (Fig. 1). The slow-slip phenomenon can be reproduced in a laboratory setting within the stable-weakening conditions promoted by low normal stress (Leeman et al. 2016;Scuderi et al. 2017;Mclaskey and Yamashita 2017;Ikari 2019). Conditional stability, a situation that occurs when the size of nucleation is commensurate with the size of the fault segment, can be invoked to simulate slow-slip events (e.g., Liu and Rice 2005;Rubin 2008;Liu 2016, 2017). Conditional stability occurs over a wide range of parameters, especially considering near velocity-neutral conditions, where most of the available energy is consumed by fracture (Wu and Chen 2014;Barbot 2019b). Furthermore, stability can be promoted by activation of viscoelastic flow in a semi-brittle matrix (Goswami and Barbot 2018), velocitydependent frictional properties (Shibazaki and Shimamoto 2007;Matsuzawa et al. 2010;Shibazaki et al. 2011) or dilatancy effects (Liu and Rubin 2010;Segall et al. 2010;Segall and Bradley 2012). The down-dip limit of slowslip is an isotherm corresponding to the stable creep of olivine and pyroxene. Below this depth, the plate interface is stable, creeping at the plate rate, only modulated by stress perturbations from earthquakes and slow-slip events (Barbot 2018b). As the recurrence times of slowslip events and large earthquakes are so widely different, the mechanical coupling between the seismogenic zone and the slow-slip region is poorly known. However, slowslip events can be triggered by distant earthquakes (Araki 2017; Wallace et al. 2017;Wei et al. 2018), act like structural barriers (Vaca et al. 2018;Rolandone 2018), or trigger other slow-slip events (Payero et al. 2008;Kano et al. 2019). Along-strike segmentation of slow-slip events in the Nankai Trough follows the spatial pattern of coupling in the seismogenic zone (Takagi et al. 2019). Numerical simulations indicate that the recurrence pattern of slow slip is modulated by the earthquake cycle .
The subduction of an oceanic slab below another plate is resisted by frictional drag along the megathrust that results in extremely localized brittle deformation. Consequently, the two plates are separated by the megathrust at geological time scales from the trench down to a depth of about 80 km (Wada and Wang 2009), which also marks the down-dip limit of stable creep, the brittle-ductile transition, and the termination of the megathrust. Below this depth, the subducting slab is permanently coupled with the upper mantle rocks, entraining a convection cell that drives local circulation of volatiles, culminating in the formation of a volcanic arc at most subduction zones (e.g., Spiegelman and McKenzie 1987;Kawakatsu and Watada 2007;van Keken et al. 2011). The termination of frictional heating associated with megathrust slip seems required to explain the pattern of surface heat flow at most subduction zones (Wada and Wang 2009). As a result, 50% of all frontal arc volcanoes lie between 85 and 119 km above the subducting slab (Syracuse and Abers 2006). The volatiles originate from the hydrous metamorphism of upper mantle rocks at the spreading center and from bend-related faulting at the outer rise (e.g., Ranero et al. 2003;Rüpke et al. 2004;Naif et al. 2013Naif et al. , 2015. The Nankai Trough offshore Honshu is a notable exception, not being associated with a volcanic arc, presumably due to the shallow dip angle of the Philippine Sea slab that inhibits efficient corner flow (Nakajima and Hasegawa 2007;Hirose et al. 2008;Hayes et al. 2018) and the influence of the Shikoku back-arc basin behind the Izu arc. The deformation of the mantle wedge above the slab is accommodated by the viscoelastic flow of peridotite, facilitated by the weakness of olivine at this temperature and pressure (Karato et al. 1986;Mei and Kohlstedt 2000;Van der Wal et al. 1993;Chopra 1997;Karato and Jung 2003;Masuti et al. 2019). Thermally activated viscoelastic flow in the asthenosphere facilitates plate tectonics and the Wilson cycle (e.g., Burke 2011). Asthenospheric flow is modulated by the seismic cycle (Barbot 2018b), creating large-scale deformation during the postseismic period (Savage and Prescott 1978;Hirahara 2002;Hu et al. 2004;Wang 2007;Pollitz et al. 2008;Wang et al. 2012;Sun 2014;Masuti et al. 2016;Klein et al. 2016;Li et al. 2018;Qiu et al. 2018;Weiss et al. 2019). Viscoelastic flow and fault slip are mechanically coupled during postseismic deformation (Agata et al. 2019;Muto et al. 2011). However, the impact of viscoelastic flow on the recurrence pattern of earthquakes and slow-slip events is poorly known.
Schematics of heat, mass, and fluid transport at subduction zones abound (e.g., Hyndman et al. 1997;Stern 2002;Agard et al. 2009, and references therein). However, how lithology, temperature, and structure affects the seismic cycle is not firmly established. Scholz (1998) invokes the temperature dependence of the friction properties of granite to define the boundary of the seismogenic zone. Lay et al. (2012) describe how other segments can be associated with tsunami earthquakes and short-period seismic radiations. Gao and Wang (2017) highlight the rheological separation of the megathrust into regions of different rupture styles. But the first-order, structural control on the dynamics of fault slip has not been properly asserted. Here, we describe a model of the seismic cycle at the Nankai subduction zone where the style of fault dynamics is controlled by the structure and composition of the overriding plate. The intersection of the megathrust with thoroughly different terranes creates structural domain boundaries that profoundly affect the style of rupture. Using the Nankai Trough as a backdrop because of the manifestation of many rupture styles, we describe how each region of the megathrust-the accretionary prism, arc crust, metamorphic belt, mantle lid, and upper mantle-can be associated with characteristic physical properties and dimensions that have a long-lasting influence on the unfolding of the seismic cycle. We highlight how the nonlinear dynamics of fault slip and the interplay between each domain introduces additional complexity, such as full and partial earthquake ruptures, modulation of slow-slip during the seismic cycle, and aperiodic super-cycles, that overprint the structural control. We also discuss the mechanical coupling between the brittle domain and the viscoelastic substrate.
The manuscript is organized as follows. In the following section, we propose a structural model of the megathrust at the Nankai Trough along a representative cross-section. We describe a numerical model that reproduces the salient features over the following sections. In the third section, we introduce the constitutive laws for fault slip and viscoelastic flow and the relevant governing equations. In the fourth section, we discuss the physical and structural properties of the model. In the fifth section, we present numerical simulations of the systemlevel dynamics. In the following sections, we focus on the mechanical coupling between earthquakes and the shortterm and long-term slow-slip events. We also discuss the impact of viscoelastic flow in the oceanic asthenosphere and the mantle wedge on the recurrence pattern of earthquakes. This is followed by a discussion and conclusion.

Structural and lithological control on megathrust dynamics
We describe a representative cross-section of the Nankai Trough where the major rupture styles are associated with structural boundaries of the overriding plate. From the trench to the asthenosphere, we describe five megathrust segments delineated by structural boundaries or an isotherm (Fig. 2). Each of these segments is associated with physical characteristics and a length scale that control fault dynamics at first order. We do not resolve the interactions with the Median Tectonic Line because pure strike-slip and thrust faults are mechanically decoupled in two-dimensional approximations. Segment A. The top segment represents the intersection of the megathrust with the accretionary prism, a sedimentary sink characterized by intense off-fault deformation (Moore et al. 1990;Hubbard et al. 2015;Sathiakumar et al. 2020), low rigidity (Sallarès and Ranero 2019), and particularly low static friction (e.g., Byrne and Fisher 1990;Cubas et al. 2013). The shortening taken up by secondary faults and folds in the upper plate builds topography and reduces the long-term loading rate on the megathrust, resulting in a low rate of seismicity on the plate interface. The low static friction coefficient of clay minerals and derived sedimentary rocks greatly affects the rupture style, resulting in anomalously long, complex ruptures (Barbot 2019b). The many high-angle ramps and splay faults of the frontal prism are responsible for the tsunami-genic potential. The dominantly velocity-weakening properties of the plate interface at these depths (Saffer and Marone 2003;Ikari and Kopf 2017;Ikari 2019) allows for the development of slowslip events, slow earthquakes, and tsunami earthquakes. Deep ruptures may also propagate through the accretionary prism and arrest near the trench. Segment B. The seismogenic zone corresponds to the cold section of the low-angle intersection between the megathrust and the arc crust. The rupture of the seismogenic zone is controlled by the frictional properties of quartz-rich crystalline rocks, which are considered velocity-weakening at temperatures up to at least 350 °C (Blanpied et al. 1991(Blanpied et al. , 1995Mitchell et al. 2013Mitchell et al. , 2016. The large intersection of the megathrust down to this isotherm allows great and giant interplate ruptures. The width of the seismogenic zone provides a characteristic length scale for rupture size, but full and partial ruptures of the seismogenic zone take place spontaneously because of the nonlinear dynamics of the system (Kato 2012;Wu and Chen 2014;Michel et al. 2017;Cattania and Segall 2018;Barbot 2019b;Cattania 2019). Giant earthquakes may also occur if ruptures propagate across segment boundaries. Segment C. We associate the long-term slow-slip events with a metamorphic belt within the mantle wedge corner, a narrow serpentinized shear zone below the arc crust. The slow-slip events develop in antigorite-rich serpentinite, which is velocitystrengthening at the low temperatures and pressures of the seismogenic zone, but velocity-weakening at and above 450 °C Katayama et al. 2013;Okazaki and Katayama 2015). The low activation energy of the viscoelastic flow of serpentinite (Hilairet et al. 2007) may hinder the development of frictional instabilities in this region (Goswami and Barbot 2018). The depth of the metamorphic belt coincides roughly with the arc lower crust farther inland, but the shear zone of the mantle wedge corner and the lower crust differ in composition. Segment D. The episodic tremor and slip phenomenon coincides with the mantle lid, i.e., the brittle layer of the upper mantle, where olivine-rich peridotite exhibits velocity-weakening frictional properties at temperatures up to 600-800 °C (Boettcher et al. 207;King and Marone 2012). The proximity to the stability transition induced by high temperature creates the conditions for chaotic short-term slow slip where a macroscopic slow rupture is accompanied by intermittent slow earthquakes (Viesca 2016a, b;Barbot 2019b). The mix-mode failure with coincident slow and fast rupture during the same event may provide a first-order explanation for the intense tremor activity of this section. Segment E. The deep velocity-strengthening segment of the megathrust corresponds to temperatures between that of the stability transition of peridotite (600-800 °C) and that of the brittle-ductile transition (1000 °C or 80 km depth). Beyond this depth, the downgoing slab is bound to the mantle wedge rocks and fault slip is prevented or greatly limited. This segment plays an important role during postseismic relaxation, releasing the coseismic stress change rapidly with afterslip, and during the interseismic period towards loading of the seismogenic zone. a Synoptic look at the relationship between brittle and ductile processes with a subduction zone geometry. The megathrust is segmented into the accretionary prism (segment A), seismogenic zone (segment B), metamorphic belt (segment C), mantle lid (segment D), and mantle wedge (segment E). The friction properties of these regions are controlled by the thermo-mechanical properties at the plate interface. The brittle-ductile transition occurs seaward of the volcanic arc. Corner flow in the mantle wedge starts landward of this boundary, in the mantle wedge asthenosphere. b Friction parameters a − b and L of segments A to E, and c corresponding R u and R b numbers. d Static friction coefficient and loading rate along the megathrust. e Rheological parameters of the oceanic asthenosphere at steady state. f Rheological parameters of the mantle wedge at steady state. g Thermal gradient in the upper mantle In the following sections, we describe a two-dimensional model that implements this idealized structural segmentation of the megathrust. Compatible with previous findings, we show that megathrust dynamics is greatly controlled by the structure and composition of the upper plate, but that additional complexity and non-periodic recurrence patterns result from the nonlinear dynamics of the system. This results in a complex history of earthquakes and slow-slip events with large-scale and long-period characteristics that are predetermined by the structural layout.

Constitutive laws and governing equations
In this section, we describe the assumed constitutive relationships for fault friction and viscoelastic flow. We describe a modeling framework that allows us to simulate the dynamics of fault slip in a viscoelastic medium that includes the mechanical coupling between brittle and ductile deformation during the seismic cycle. We then discuss important take-aways from a dimensional analysis of the governing equations that allows us to discuss the emergence of a wide spectrum of rupture styles. These considerations provide the building blocks for simulating the short-term dynamics of subduction zones during the seismic cycle.

Constitutive framework for fault slip
The evolution of fault friction on a fault during all stages of the seismic cycle can be captured by the constitutive framework of rate-and-state friction introduced by Dieterich (1979), Ruina (1983), and Rice and Ruina (1983). The constitutive law is motivated by a wealth of laboratory experiments (Dieterich 1978;Marone et al. 1990;Marone and Kilgore 1993;Nakatani 2001;Yamashita et al. 2014;Lyu et al. 2019, and references therein) and is supported by various micro-physical models for the evolution of strength in fault gouge or bare contacts (Chester 1989; Sleep and Blanpied 1992;Sleep 1995Sleep , 2006Barbot 2019a). The framework of rate-and-state friction has been useful to model a wide spectrum of rupture styles, including slow and fast ruptures Lapusta and Rice 2003;Hori et al. 2004;Liu and Rice 2005;Rubin 2008;Chen and Lapusta 2009;Kaneko et al. 2010;Barbot et al. 2012;Wei et al. 2015;Veedu and Barbot 2016;Lambert and Barbot 2016;Lui and Lapusta 2016;Salman et al. 2017;Yu et al. 2018;Ong et al. 2019) and the dynamics of quasi-static deformation in dominantly aseismic periods Avouac 2004, 1992;Barbot et al. 2004;Bruhat et al. 2011;Rousset et al. 2012;Rollins et al. 2015).
Several forms of constitutive equations that fall into the rate-and-state framework have been proposed (Ruina 1983;Linker and Dieterich 1992;Rice and Ben-Zion 1996;Kato and Tullis 2001;Bizzarri 2011;Nagata et al. 2012), in particular due to uncertainties about the evolution law for the state variable and issues regarding regularization for vanishing velocities. In this study, we adopt the formulation proposed by Barbot (2019a), where the velocity across the fault zone is controlled by the area of the contact junctions that support the shear and normal loads. The model assumes the following constitutive relationship where V is the sliding velocity, τ is the norm of the shear traction resolved on the fault plane, T is the local absolute temperature, µ 0 is the static coefficient of friction at velocity V 0 and temperature T 0 , A is the real area of contact divided by the nominal surface area, χ is the indentation hardness of the fault material, and a ≪ 1 is a power exponent. Fault slip is thermally activated with the activation energy Q; R is the universal gas constant. The real area of contact depends primarily on the applied normal stress, but it evolves during the seismic cycle, modulated by changes in the morphology of the micro-asperities. Flattening of the micro-asperities increases the area of contact junctions and strengthens the fault. Comminution during fault slip reduces the average contact area and weakens the fault. This behavior is captured by a dependence on effective normal stress and a state variable where σ is the effective normal stress, L is the characteristic weakening distance, b ≪ 1 is a power exponent, and θ is a state variable that represents the age of the contact, or equivalently, the duration of contact flattening (Barbot 2019a). Combining (1) and (2), we obtain the constitutive relationship To capture the evolution of the topography of the fault surface at the micro-scale, the state variable follows the evolution law (Barbot 2019a) where H is the activation enthalpy for contact healing (Chester and Higgs 1992;Chester 1994). In this study, we assume isothermal conditions, and, in the limit of this approximation, Eqs. (3) and (4) represent a complete constitutive framework for the dynamics of fault slip. Choosing T = T 0 , Eq. (4) simplifies to the aging law proposed by Ruina (1983). The parameters a, b, and L have the same function as in the classical formulation of Ruina (1983), despite the multiplicative form of Eq.
(3). The spatial distribution of the constitutive parameters µ 0 , a, b, and L and the effective normal stress σ along the megathrust as a function of lithology and temperature will shape the dynamics of the seismic cycle. The constitutive law of Eq.
(3) takes a multiplicative form, differing in that regard from the additive, logarithmic form used elsewhere (e.g., Ruina 1983, and references therein). A truncated power series expansion of the multiplicative form, which remains valid for a wide range of velocities around V 0 = 1 μm/s, reduces to the additive form (Barbot 2019a). However, the logarithmic form is ill-posed for vanishing velocity and requires regularization, such as the one with a hyperbolic sine motivated by thermo-dynamic considerations proposed by Rice and Ben-Zion (1996). In contrast, the multiplicative form is valid for vanishing velocities and does not require additional regularization. The multiplicative form of rate-andstate friction described above produces a different style of rupture dynamics in weak faults ( µ 0 ≤ 0.1 ), resulting in exceedingly long, multi-pulsed ruptures resembling that of tsunami earthquakes (Barbot 2019a, b).

Constitutive framework for viscoelastic flow
The deformation of crust and mantle rocks, inferred from laboratory experiments on aggregates or single crystals, is accommodated by several micro-mechanisms that include diffusion creep, dislocation creep, and dislocation glide (Karato et al. 1986;Hirth and Tullis 1992;Gleason and Tullis 1995;Chopra 1997;Mei and Kohlstedt 2000;Rybacki and Dresen 2000;Karato and Jung 2003;Hirth and Kohlstedt 2003;Dimanov and Dresen 2005). During steady-state, i.e., at constant stress or strain-rate, a constitutive relationship of the form where A is a constant pre-factor, σ is the norm of the deviatoric stress, and m, n, and r are power exponents, predicts the deformation-rate ǫ based on in situ physical conditions, including water content C OH , grain-size d, pressure p, activation volume , activation energy E, and temperature T. For diffusion creep and pressuresolution creep, the power exponent is n = 1 , corresponding to linear viscoelasticity. For dislocation creep, n > 1 and m = 0 , corresponding to a power-law rheology.
Dislocation glide has a nonlinear response with n ∼ 2 , but is also sensitive to grain-size ( m = 0) (Karato 2008). During periods of rapid change of stress or strain-rate, other deformation processes take place within the rock or mineral structure, leading to a short phase of hardening that transitions to the steady-state response, usually referred to as transient creep (Chopra 1997;Wu et al. 2016). In contrast to steady-state, the micro-mechanics of transient creep are poorly understood and different constitutive relationships have been proposed (Sherburn et al. 2011;Thieme et al. 2018;Holtzman et al. 2018), but plastic anisotropy within single crystals may be an important factor (Masuti et al. 2019). Because of the rapid change imposed by the seismic cycle, transient creep may be operating during postseismic relaxation (Pollitz et al. 2008;Freed et al. 2010;Hoechner et al. 2011;Tang et al. 2019;Muto et al. 2011). The strain-rate of transient creep may be represented by an additional strain component (Masuti et al. 2016) where A K is a pre-exponential factor, q = σ − 2G K ǫ K is the effective stress, and G K is strain-hardening coefficient. This rheological framework has been useful to model postseismic deformation (e.g., Nur and Mavko 1974;Savage 2000;Hirahara 2002;Freed and Bürgmann 2004;Hu et al. 2004;Hetland and Hager 2005;Wang 2007;Pollitz et al. 2008;Johnson et al. 2009;Suito and Freymueller 2009;Wang et al. 2012;Sun 2014;Hu et al. 2016;Klein et al. 2016;Muto et al. 2011, and references therein). In this study, we approximate the rheological response by the steady-state, as the details of sub-annual evolution of viscoelastic relaxation is not of central importance to our discussion.

Governing equations
We consider the dynamics of fault slip and viscoelastic flow during the seismic cycle within the quasi-dynamic approximation Rice and Tse 1986;Rice et al. 2001), whereby the wave-mediated stress transfer is ignored. In a viscoelastic medium, the traction along the fault depends on both the velocity of the surrounding fault elements and the anelastic strain-rate in the surrounding rocks. The rate of change of shear traction due to fault slip (e.g., Maruyama 1964;Rice and Tse 1986;Shibazaki and Matsu'ura 1992;Kato 2003;Ziv and Cochard 2006;Fukuyama 2007) is extended to include the effect of viscoelastic flow as follows where K is the Green's function tensor for self-interactions on the fault surface ∂Ŵ (e.g., Chinnery 1963;Okada 1992;Segall 2010;Gimbutas et al. 2012;Nikkhoo and Walter 2015), v is the instantaneous velocity vector, v L is the local loading rate on the fault, ǫ represents the anelastic strain-rate tensor in the ductile domain , and J is the Green's function tensor describing how distributed deformation affects the shear stress on the fault (Faivre 1969;Chiu 1978;Nozaki and Taya 2001;Kuvshinov 2008;Barbot et al. 2017;Barbot 2018a). The last term on the right-hand side corresponds to the radiation-damping term with the rigidity G, shear wave velocity V S and slip acceleration V (Rice and Tse 1986; Tse and Rice 1986; Rice et al. 2001). Viscoelastic deformation is controlled by the amplitude and orientation of stress. Following a decomposition of the total strain into the elastic and anelastic components (Andrews 1976; Barbot and Fialko 2010), the evolution of the deviatoric stress tensor in the viscoelastic substrate can be written where M is the Green's function for self stress in the ductile substrate and L represents the stress impacted by fault slip onto the ductile rocks. Therefore, the Green's functions J and L represent the mechanical coupling between brittle and ductile deformation. The terms v L and ǫ L represent the long-term slip-rate on the fault and the long-term strain-rate in the viscoelastic substrate, respectively, and are considered constants imposed by long-term geodynamics. Combining the constitutive laws (3), (4), and (5) with the governing equations (7) and (8) allows us to simulate the dynamics of fault slip coupled with viscoelastic flow.
Several numerical approaches have been proposed to incorporate nonlinear viscoelastic flow in seismic cycle simulations, including the finite-difference and finite-elements methods (Herrendörfer et al. 2015;Erickson et al. 2017;Sobolev and Muldashev 2017;Biemiller and Lavier 2017;Allison and Dunham 2018;Tong and Lavier 2018). Here, we use the integral method proposed by Lambert and Barbot (2016) and Barbot (2018b), whereby the fault surface is discretized into line elements and the ductile domain is meshed with surface elements in crosssection Barbot 2018a). The integrals in Eqs. (7) and (8) can then be evaluated numerically with algebraic relationships. The approach allows us to include non-planar fault geometry (Ong et al. 2019) and to accommodate nonlinear rheological laws. The dynamics of fault slip coupled with viscoelastic flow is then obtained based on the Runge-Kutta method with adaptive time steps (Press et al. 1992). Some aspects of the numerical method have been validated against other tools (Erickson et al. 2020).

Take-away from dimensional analysis
A dimensional analysis of the governing equation for fault slip can help identify four dimensionless variables that control the rupture style and other characteristics of the seismic cycle in isothermal conditions (Barbot 2019b). Non-dimensional parameters are important to consider because they reduce the parameter space considerably: various combinations of physical parameters producing the same set of non-dimensional parameters are expected to produce similar rupture styles, only with different recurrence times or overall size. The static friction parameter µ 0 controls the overall fault strength and affects the rupture style. In particular, weak faults can rupture in anomalously long and complex ruptures. Another parameter V th = 2V S (b − a)σ /G , linked to seismic wave radiation, controls the peak velocity and the rupture velocity during fast ruptures (Barbot 2019b). This parameter affects the accretionary prism in subduction zones particularly because of the low shear wave speeds of water-saturated sedimentary muds. The remaining two non-dimensional parameters R u and R b represent important degrees of freedom because they vary wildly in nature and play important roles in controlling fault dynamics. The Dieterich-Ruina-Rice number where G is the rock rigidity, W is a representative length scale of the fault segment, and a and b are defined in Eqs. (1) and (2), controls the importance of non-local stress transfer. The R u number can be shown to represent a ratio of segment size to nucleation size (Barbot 2019b). The R b number controls the relative importance of the transient evolutionary effects, i.e., the weakening and strengthening with contact rejuvenation and contact aging, respectively. The R b number also controls the ratio of static and dynamic stress drops during rupture, and hence the amount of fracture energy consumed. The R u and R b numbers form a two-dimensional phase space, each coordinate being associated with a potentially different rupture style. Large R u numbers promote complexity of the seismic cycle, with a wide range of rupture sizes (Wu and Chen 2014;Cattania and Segall 2018;Cattania 2019). Low R u numbers correspond to increasing stability of the fault patch, i.e., a propensity for cycles of slow earthquakes, slow-slip events, or creep, depending on the actual value. Low R b numbers also promote macroscopic stability, but a combination of large R u and low R b leads to chaotic cycles of slow-slip events encompassing slow earthquakes (Barbot 2019b). During model design, we choose the physical parameters considering the predicted dynamics associated with the corresponding set of non-dimensional parameters and their coordinates in phase space.

Model description and physical parameters
We design a numerical simulation of short-term subduction dynamics that incarnates the structural model presented in the previous Section. The model allows us to discuss the emergence of a wide range of rupture styles, each associated with a specific megathrust segment, and to illuminate the nonlinear coupling among the various domains of the megathrust and surrounding viscous substrate. For numerical convenience, we consider a two-dimensional approximation in plane-strain conditions along a trench-perpendicular cross-section of the Nankai Trough (Fig. 2). The model is complex enough to discuss general aspects of fault dynamics in this region, but the reader is referred to even more sophisticated models of the Nankai Trough (Hori et al. 2004;Hori 2006;Matsuzawa et al. 2010;Shibazaki et al. 2012;Matsuzawa et al. 2013;Hyodo and Hori 2013). We follow the megathrust geometry from the Slab1.0 model (Hayes et al. 2012) along the profile P-P' in Fig. 1. We extend the slab geometry down to more than 200 km depth, even though the deepest section is not observed seismologically (Nakajima and Hasegawa 2007;Hirose et al. 2008;Hayes et al. 2018). We divide the megathrust into five segments with piecewise constant dip angles. The physical parameters at each segment are determined considering the lithology and temperature at the plate interface as well as the expected dynamics based on the resulting non-dimensional number in the governing equations. Because of numerical limitations, we adopt a uniform rigidity G = 30 GPa , shear wave speed V S = 3000 m/s , and Poisson's ratio ν = 0.25 . This prevents us from using more realistic variations, in particular, in the accretionary prism, where low rigidity plays an important role (Sallarès and Ranero 2019). The spatial distribution of the relevant physical parameters is compiled in Table 1.
For segment A, we use a low static friction coefficient µ 0 = 0.1 , compatible with the presence of poorly consolidated material in the overlying sediment pile-up. Clay minerals also exhibit velocity-weakening behavior under low normal stress and low loading rate. For simplicity, we conflate the accretionary prism and the frontal prism into one single region of uniform properties. We also ignore the change of frictional properties associated with the smectite to illite transition (Saffer and Marone 2003). The local loading rate on the megathrust V L = 1.2 × 10 −9 m/s , is smaller than the plate convergence of V pl = 1.8 × 10 −9 m/s to account for the partitioning of deformation among the megathrust, thrusts, and folds of the sedimentary wedge (Sathiakumar et al. 2020). Using the down-dip width of the megathrust that intersects with the accretionary wedge, i.e., W = 41.75 km , and the frictional properties shown in Table 1, the resulting non-dimensional parameters R u = 1.82 and R b = 0.025 are associated with a complex long-term slow-slip regime (Barbot 2019b). All other seg- In the seismogenic zone (segment B), we select a set of physical parameters that produces large ruptures within the recurrence period of the Nankai earthquakes, which vary within 140 ± 42 years (Ando 1975;Rikitake 1976;Thatcher 1984;Kumagai 1996). In addition, we choose the frictional properties to obtain complex seismic sequences with full and partial ruptures of the seismogenic zone, which occurs for R u ≥ 16 for a single asperity embedded in a planar fault (Barbot 2019b). Using W = 148.75 km as the characteristic length scale of the seismogenic zone and with the combination of effective normal stress, characteristic weakening distance, and kinematic friction coefficients, we obtain R u = 29.7 , a value high enough to generate complex seismic cycles with a variety of rupture sizes.
In the metamorphic belt (segment C), we choose a set of frictional properties that reflects stable weakening and the emergence of long-term slow-slip events with recurrence times of about 2 years and rupture durations of several months. As discussed in the introduction, this can be obtained with a variety of mechanisms that permit nucleation but inhibit fast slip. For simplicity, we adopt a simple frictional behavior with physical properties such that R u = 1.75 and R b = 0.18 , a near-stable regime of periodic slow slip (Barbot 2019b). The low R b number associated with near velocity-neutral conditions is justified by the gradual transition from velocity-weakening to velocity-strengthening properties as temperature increases with depth. After we added other velocityweakening segments to the model, the recurrence time of the long-term slow-slip events can be extended to about 2 to 6 years, consistent with the observations in different regions along the strike of the Nankai Trough (Hirose et al. 1999;Hirose and Obara 2005;Yarai and Ozawa 2013;Kobayashi 2014).
In the mantle lid (segment D), we consider frictional properties that produce a chaotic form of slow-slip events, whereby macroscopic slow slip is accompanied by slow earthquakes. This phenomenon emerges from the nonlinear dynamics of fault slip and was identified 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 in studies of rupture nucleation (Viesca 2016a, b) and seismic cycles (Barbot 2019b). Chaotic slow-slip cycles occurs in near velocity-neutral conditions, i.e., close to the stability transition forced by high ambient temperatures. We also consider kinematic friction properties and effective normal stress to produce slow-slip events every about 1 year with ruptures that last several weeks. Using the width of the velocity-weakening intersection of the megathrust with the upper mantle, W = 32 km , we arrive at R u = 6.1 and R b = 0.09 . The lower R b number than in the up-dip segment is justified by a greater proximity to the velocity-weakening to velocity-strengthening transition, which is likely an isotherm. This combination of the R u and R b numbers produces complex slow-slip events encompassing slow earthquakes (Barbot 2019b).
In section E, further down-dip, the slab dives deeper into the mantle and the increasing temperature results in a long ( W = 53.3 km ) velocity-strengthening segment. The frictional contact ends at the brittle-ductile transition at a depth of 71 km (the difference in depth from the volcanic arc originates from the absence of topography in the numerical model). For simplicity, we assume the same viscoelastic properties for the mantle wedge and the oceanic mantle, ignoring a more realistic contrast in water content and temperature. The temperature profile follows the cooling half-space model with a plate age of 60 Myr, and a basal mantle temperature of 1673 K. We use the viscoelastic properties for dislocation creep of wet olivine with 1000 ppm H/Si (Hirth and Kohlstedt 2003). We ignore viscoelastic flow in the arc lower crust. We drive viscoelastic flow with the non-zero background shortening rate ǫ 22 = −10 −15 /s , where the subscripts 2 and 3 represent the x 2 and x 3 coordinates in the horizontal and depth direction, respectively.
We discretize the fault with a mesh fine enough to resolve the smallest characteristic lengths of the problem, in particular nucleation size and cohesion size. We choose the final, uniform fault sampling size of 125 m with a numerical convergence test. In the ductile substrate, we refine the mesh to resolve the spatial distribution of the stress change caused by megathrust earthquakes. We use 10-km-sized cross-section volume elements. The ductile substrate is represented by an unstructured mesh of triangular elements (Barbot 2018a, b).

System-level dynamics of a subduction megathrust
The two-dimensional subduction model simulates the dynamics of the megathrust and the surrounding asthenosphere during seismic cycles, including full and partial ruptures of the seismogenic zone, foreshocks and aftershocks, shallow slow-slip events, deep long-term and short-term slow-slip events, creep and afterslip, and viscoelastic relaxation in the ductile substrate (Figs. 3, 4, 5, 6).
Megathrust dynamics in the fault area within or adjoining the accretionary prism is perhaps the most complex of the model. The slip deficit is released by a wide range of slip styles that include long-term creep waves, partial coupling, slow and fast creep, afterslip, slow-slip events, and slow earthquakes (Figs. 3 and 4). The 450-year window into the seismic cycle shows no cyclic pattern whatsoever. Every full or shallow rupture of the seismogenic zone is followed by a long-period afterslip in segment A that endures across multiple earthquake cycles, including across super-cycles that is defined as the period between two full ruptures. The long afterslip response is better characterized as a long-term creep wave. In the interseismic period, creep propagates into the seismogenic zone from the top, resulting in the occurrence of slowslip events of increasing sizes. Some events culminate to a partial shallow rupture of the seismogenic zone. During the few years that precede full or shallow earthquakes, the shallow megathrust appears to creep near or at a fraction of plate rate. In general, the creep velocity is not uniform in the accretionary region, exhibiting a complex spatial pattern of partial locking.
Rupture dynamics in the seismogenic zone is also complex, inherently due to the high R u number, but also because of interactions with the quasi-stable accretionary prism and the slow-slip region. The megathrust generates earthquake super-cycles, defined here as the period separating two full ruptures of the seismogenic zone, with repeating periods varying from 30 to 180 years, with most super-cycles taking 120 years. Within a super-cycle, the recurrence time of full and partial ruptures-aftershocks excluded-is irregular, ranging from 10 to 75 years. Most earthquakes nucleate near the top and bottom boundaries of the seismogenic zone due to concentrated shear stress caused by creep or slow slip in the neighboring segments. A shallow foreshock sequence sometimes precede a full rupture, for example, the shallow events leading up to full rupture 54 (Fig. 6). All earthquakes are followed by afterslip on the deep megathrust and viscoelastic relaxation in the oceanic asthenosphere and the mantle wedge. Event 46 represents an aftershock of a full rupture (event 45).
Long-term and short-term slow-slip events recur every 4 years and 11 months, respectively. But the overall dynamics of slow-slip events is complex due to the systematic interaction of the short-term and long-term slow-slip regions, the trench-ward migration of the coupling boundary, and the long-term modulation by seismogenic zone ruptures. The model illuminates a range of interactions between the seismogenic zone and the slowslip region. Most of the great earthquakes that rupture the entire seismogenic zone penetrate through the down-dip slow-slip region. However, the full ruptures that initiate at  (Fig. 6). During a supercycle, the shallow partial ruptures arrest within the seismogenic zone and all deep partial ruptures penetrate into the slow-slip region. All ruptures that propagate in the slow-slip region greatly perturb the recurrence pattern of slow-slip events, sometimes for decades. The long-term slow-slip region is the most frequently and most durably affected due to its proximity with the seismogenic zone. The recurrence pattern of slow-slip events is perturbed by interactions between the short-term and long-term slow-slip events and interseismic migration of the locking depth. As the locking depth of the seismogenic zone slowly migrates towards the trench, the width of the creeping section is expanded up-dip of segment C, leading to expansion of the long-term slow-slip ruptures (Figs. 3e,  4c). The short-term slow-slip events have stable peak slip rates and recurrence times in the early stage. However, the slowly expanding long-term slow-slip events eventually affects the short-term slow-slip events, which exhibit greater variability in slip rates and recurrence intervals in the later stage of the interseismic period (Figs. 3e, 4d). As a result, the recurrence pattern of long-term and shortterm slow-slip events are opposite: the long-term slowslip events, closer to the seismogenic zone, are suppressed in the early stage of the seismic cycle, and correspondingly, the short-term slow-slip events are quasi-periodic. Later in the seismic cycle, long-term slow-slip events are well developed, and the short-term slow-slip events feature a complex, aperiodic cycle.
Due in part to the chaotic dynamics of slow slip in segment D induced by the coordinates of the physical parameters in phase space, the short-term slow-slip sequences exhibit tremendous complexity. This is exemplified in two slow-slip sequences that start at about 930 and 1093 years after the beginning of the simulation (Fig. 5). As the locking depth of the seismogenic zone migrates up-dip, the long-term slow-slip events in segment C are facilitated, triggering larger stress perturbation in segment D when it is expanding up-dip and down-dip (Fig. 5a, e). The slow slip reaching the segment boundary breaks into segment D, generating a relatively rapid slow-slip event that propagates faster and bilaterally, lasting half an hour for the first sequence (Fig. 5d) or one day for the second sequence (Fig. 5h). A faster rupture nucleates immediately thereafter and breaks the entire locked section of segment D, eventually triggering another slow-slip event (Fig. 5b, f ) in segment C. The evolution of these two sequences of slow slip and slow earthquakes demonstrates that the migration of deep slow slip can facilitate the nucleation of adjacent slow earthquakes. Conversely, the propagation of slow earthquakes affects the pattern of long-term and short-term slow-slip events. These sequences are reminiscent of the slow slip transient that was excited during an episodic tremor and slip event at the down-dip extension of the locked zone beneath western Shikoku Island in Southwest Japan (Kano et al. 2019).
At greater depth, fault dynamics is much simpler, dominated by continuous creep at plate rate during most of the seismic cycle, but occasionally interrupted by afterslip following full and partial ruptures (Figs. 3a, 4f ). The afterslip period consists of a short pulse of rapidly decaying creep that lasts a few years, as observed after all suitably instrumented large subduction earthquakes (e.g., Hsu et al. 2006;Chlieh 2007;Feng et al. 2016;Tsang et al. 2016;Bedford et al. 2016;Klein et al. 2016;Hu et al. 2016;Qiu et al. 2018;Tang et al. 2019;Muto et al. 2011).
In the upper mantle, a transient flow acceleration is triggered by the sudden stress perturbation induced by a full rupture or a deep partial rupture (Fig. 3b, c). Due to variable earthquake sizes and down-dip moment distribution, the peak postseismic strain-rate varies substantially. The full postseismic relaxation takes 50-150 years and during some short interseismic periods the induced stress perturbation is not fully relaxed. During the postseismic transient, the effective viscosity decreases by one order of magnitude from its steady-state value (Figs. 2a,7a), consistent with the stress dependence of effective viscosity. The postseismic strain-rate components ε 22 , ε 23 and ε 33 increase to 50 times larger than the steady-state values, while the locations of the peak strain rates and their directions remain the same except for the region at a depth of 100 to 200 km below the deeper portion of the seismogenic zone and the slow-slip segments (Fig. 7).
(See figure on next page.) Fig. 4 a History of slip velocities at the center of each megathrust segment. b Afterslip, slow earthquakes and slow-slip events in the accretionary prism. c Full and partial ruptures break the locked seismogenic zone with varying recurrence intervals. d Long-term slow-slip events in the metamorphic belt at down-dip of the seismogenic zone. e Short-term slow-slip events in the mantle lid. f Creeping in the mantle wedge corner interrupted by afterslip following full and partial ruptures. g Overall peak slip velocity on the megathrust, which captures events that do not rupture the segment centers. h The strain rate evolution of the volume element shown as the gray triangle in the mantle wedge in Fig. 3. i The strain rate evolution of the volume element shown as the gray triangle in the oceanic mantle in Fig. 3  This feature is consistent with the deep SSEs frequently loading the viscoelastic substrate. In the particularly long earthquake cycles, the relaxation is followed by a slower flow with strain-rates lower than the background loading rate, which preserves the long-term imposed strain-rate on average in the asthenosphere.

Surface dynamics
Recent studies have illuminated the different contributions of megathrust slip and mantle flow in the dynamics of surface deformation (e.g., Barbot 2018b;Hu et al. 2016;Noda et al. 2018;Suito 2017;Sun 2014;Qiu et al. 2018;Agata et al. 2019;Weiss et al. 2019). In particular, the viscoelastic relaxation of the oceanic asthenosphere in these models generates retrograde horizontal displacements in the forearc during viscoelastic relaxation. We describe the surface dynamics by tracking the surface deformation deficit, i.e., the difference between the cumulative displacement and the long-term average throughout several earthquake super-cycles (Fig. 8). We compute the space-time distribution of surface displacement produced by fault slip and viscoelastic flow separately by multiplying the slip or strain in the regions of interest by the respective Green's function in a post-processing step. Overall, the fault slip contribution dominates the surface displacements above the megathrust, but its effect is concentrated in the forearc region and is reduced in the interseismic period because of re-locking of the seismogenic zone. Partial ruptures also produce sudden surface displacements. The complex behavior of the shallow segment generates long-term displacements in the accretionary prism and in the outer rise. Viscoelastic surface deformation is most rapid in the several years following large earthquakes, after which viscoelastic flow gradually returns to background rates in the following 50 to 80 years. The surface displacements produced by viscoelastic flow are more distributed than those produced by fault slip, spreading across the outer-rise, forearc, and backarc regions. All earthquakes, afterslip, slow-slip events and creep produce seaward surface displacement. The viscoelastic relaxation produces retrograde motion in the forearc, most notably above the seismogenic zone. Note that the bottom of the seismogenic zone often coincides roughly with the continental slope, so this phenomenon may only be measured with seafloor geodesy. The vertical component of the surface deformation shows spatial diversity that highly depends on the location of the source. The fault slip of a great earthquake produces substantial uplift above the major rupture area and subsidence in the region landward of the uplifted region (Fig. 8b, indicated by colors). The boundary of uplift and subsidence regions is about 150 km landward from trench. Partial ruptures produce a similar uplift and subsidence pattern, but its location varies with the distribution of coseismic slip. For example, the surface uplift and subsidence pattern created by event 47, which occurs in the down-dip portion of the seismogenic zone, are farther from the trench than the patterns created by event 42 that occurs in the up-dip portion and the full rupture 43 (Fig. 8b). The deformation caused by viscoelastic flow is significant during the postseismic period, producing uplift within the region between 0 to 150 km landward from the trench and the volcanic arc region, and subsidence in the outer rise and within the region between 150 to 250 km from the trench (Fig. 8c). Some important details concerning the dynamics of surface deformation may vary among subduction zones depending on megathrust geometry, the age of the subducting plate, and the rheology of the mantle wedge, all of which affecting the mechanical response of the brittle-ductile system.

System-level dynamics and impact of mechanical coupling
In this section, we compare our coupled subduction model with simpler models where only one or a pair of velocity-weakening segments are included, assuming that the surrounding fault region is velocity strengthening. This allows us to identify the effect of the mechanical coupling between various nearby segments. First, we consider the case of a shallow velocity-weakening region that shares the characteristics of segment A, the accretionary prism. In this case, the shallow megathrust produces creep waves and long nucleations culminating in a slow earthquake (Fig. 9a). This behavior is compatible with the chaotic slow-slip events that emerge at low R b numbers, producing macroscopic slow-slip interspersed by short-duration slow earthquakes (Barbot 2019b). In contrast, the coupled model produces a wider range of behaviors at these depths, including slow-slip with slow  . 7 Effective viscosity and plastic strain-rate during the earthquake cycle. The postseismic relaxation following the full rupture event 43 takes longer time than that following the partial rupture event 47. a Effective viscosity of the asthenosphere after the full rupture. b Plastic strain rate corresponding to four stages in a. c Effective viscosity of the asthenosphere after the partial rupture. d Plastic strain rate corresponding to four stages in c earthquakes, but also regular slow-slip events, long-term creep waves, and afterslip. This comparison illustrates the disruptive effect of the nearby seismogenic zone, which forces the relaxation of slip deficit in various modes during the interseismic period. Next, we compare the coupled model with a simpler case where only the seismogenic zone is included. This is done by changing the friction property of segments A, C, and D from velocity-weakening to velocity-strengthening, keeping the other friction parameters the same. The uncoupled model produces super-cycles of great regularity, where every event is followed by another in a complex, but quasi-periodic repeating sequence (Fig. 9b). The recurrence intervals of great earthquakes in the uncoupled model are about 95 years, while those in the coupled model vary from 70 to 115 in the period considered (Fig. 9d). The uncoupled model produces recurring partial ruptures during a super-cycle, while partial ruptures are not systematic in the coupled model. For example, there is no partial rupture during the second super-cycle of the period shown in Fig. 9d. The uncoupled model illustrates the complexity that spontaneously emerges from the nonlinear dynamics of fault slip for high R u numbers (e.g., Barbot 2019b; Cattania 2019).   Fault slip produces seaward displacement above the megathrust and small landward displacement in the outer rise (denoted by arrows). Fault slip also produces significant uplift in the forearc region above the major rupture area and subsidence outside this region. Small coseismic subsidence close to the landward side the trench is caused by the locking in the shallowest portion of the megathrust (background color). c Postseismic viscoelastic flow in the asthenosphere following full ruptures produces landward displacement (white arrows) on the seaside of the mantle wedge corner and seaward displacement on the other side. Uplift is found above the rupture area and most of the mantle wedge. Viscoelastic subsidence is found mostly above the metamorphic belt and in the outer rise Time (  Time (  Additional complexity, in particular, aperiodic cycles, is obtained by mechanical coupling with the nearby velocity-weakening segments. Finally, we examine the dynamics of another simple model where only the slow-slip events of segments C and D are included. The uncoupled model converges to a repeating period-two cycle of short-term and longterm slow-slip events (Fig. 9c). In comparison, the coupled model shows a reduced amplitude of the slow-slip events, possibly due to the locking of the seismogenic zone (Fig. 9e, f ). Following great earthquakes or deep partial ruptures, afterslip propagates down-dip into the slow-slip region, weakening the fault segment and releasing the strain energy. The long-term slow-slip events are suppressed for 10 to 50 years, depending on the size of the earthquake (Fig. 9g). The short-term slow-slip events are accelerated for several years, after which they resume their natural recurrence pattern (Fig. 9h). The acceleration of short-term slow-slip events following a great earthquake is reminiscent of the increased tremors and slow-slip activity after the 2011 Tohoku-Oki earthquake (Nishikawa et al. 2019). As the long-term slow-slip events become more significant in the later stage of an earthquake cycle, they produce larger stress perturbation in the neighboring segments, which in turn leads to the variability of the recurrence time and amplitude of short-term slow-slip events. That slow-slip events in the metamorphic belt are inhibited in the decades following large earthquakes may explain the absence of recurring slow-slip observations in the Sunda trench. The series of recent Sumatra earthquakes, the 2004 Aceh-Andaman, the 2005 Nias, the 2007 Bengkulu earthquakes, may indeed have suppressed the slow-slip phenomenon in the region for several decades. If the model applies to the Sunda trench, slow-slip events should resume in the next decades.

Seismic cycles modulated by viscoelastic flow
With the model assumption of the viscoelastic flow presented in Table 1, the viscoelastic flow in upper mantle gradually modulates the earthquake super-cycle and the evolution of slow-slip events. This effect is illustrated by comparing the dynamics of the coupled model with another one where the upper mantle flow is preempted. In the latter, the fault slip is the only deformation mechanism included. We build earthquake catalogs for the two models and both of them show high variability in the size and recurrence time of earthquakes (Fig. 10a, b).
To represent the evolution of the super-cycle, the slip deficit at the center of the seismogenic zone, which is the difference between the cumulative slip and the long-term loading, are plotted for both models in Fig. 10c. The effect of viscoelastic flow accumulates for 5 super-cycles (560 years), after which it starts to modulate the short-term slow-slip events that are located in the deepest velocityweakening segment, also the closest to upper-mantle flow (Fig. 10f ). The modulated short-term slow-slip behavior in turn decelerates long-term slow-slip events and influences earthquake cycles in the up-dip segments, amplifying the impact originated from viscoelastic flow (Fig. 10d,  e).
The statistics of recurrence intervals of earthquakes and lengths of super-cycles are shown in Fig. 10g, h, respectively. With the cumulative effect from viscoelastic flow, the earthquake recurrence time shows a more random pattern, even though the most frequent recurrence time of ∼ 55 years is the same in both cases. Surprisingly, the super-cycles modulated by asthenosphere flow are more concentrated, around 75 and 120 years. This could be related to the interactions between deep partial ruptures and slow-slip events. When long-term and short-term slow-slip events are decelerated by viscoelastic flow, the buildup of shear strain neighboring the slow-slip region may be slackened, which might delay a partial rupture that nucleates near the bottom of the seismogenic zone to the next super-cycle. For example, when viscoelastic flow in the asthenosphere is included, the full rupture events 32 and 33 occur consecutively (Fig. 10a) and the supercycle lasts a shorter period. When viscoelastic flow is not activated, partial ruptures always occur between two full ruptures (Fig. 10b), and the super-cycles span a longer period.

Discussion
Fault dynamics at subduction megathrusts is impacted by the structure, composition, and temperature of the rock assembly at the plate interface. A structural model such as the one described in the second Section provides the constitutive and boundary conditions for the seismic cycle, but tremendous complexity emerges from the nonlinear dynamics of frictional sliding, the interactions among fault segments, and the coupling between brittle and ductile deformation. The simulations of megathrust dynamics illustrate the system-level behavior of a complex mechanical system. Analyses of simpler architectures provide useful insights, for example, how ruptures styles can be associated with certain combinations of physical and geometrical parameters. However, some unique behaviors only unfold with added structural complexity. The presence of fault bends, the triggering of neighboring segments, and the long-range interactions between the brittle lithosphere and the ductile asthenosphere play an important role.
Despite the emerging complexity, the structural boundaries provide a virtually permanent setup that bounds the overall mechanical behavior. The slow-slip events are  Slip rate (m/s) Fig. 10 Impact of asthenospheric flow on the seismic and slow-slip cycles. a The earthquake catalog from the simulation with coupled brittle and ductile deformation. Great earthquakes (red stars with event numbers) corresponding to full rupture of the seismogenic zone are connected by black lines with partial ruptures (orange stars) and slow earthquakes (yellow stars) in order of occurrence. b The earthquake catalog from the simulation with only brittle deformation along the megathrust segments. Events are connected by a blue line. c Slip deficit (black and blue for models with and without viscoelastic flow, respectively) at the center of the seismogenic zone (segment B). Earthquake recurrence time starts to be modulated by viscoelastic flow after several super-cycles (about 560 years). a-c share the same x-axis. d The enlarged view of earthquake cycles of the time period between 500 and 600 years in c. e and f share the same x-axis as d. e Effect of viscoelastic flow on long-term slow-slip events. f Effect of asthenosphere flow on short-term slow-slip events. g Distribution of recurrence times for earthquakes with force larger than 10 8 N (all ruptures) with (gray) our without (blue) coupling with viscoelastic flow. h Distribution of recurrence time of great earthquakes with force larger than 10 10 N (only full ruptures) modulated by nearby ruptures, but the slow-slip cycle endures over multiple seismic cycles or super-cycles. Different types of ruptures, slow and fast, occur near the boundaries of the seismogenic zone, but full and partial ruptures still represent the dominant style of rupture in the forearc. The cycles are aperiodic, but the types of ruptures are pre-determined. A first-order understanding of the down-dip segmentation of the seismic cycle at subduction zone benefits from integrating geological and theoretical insights. The variability of earthquake sizes results from the nonlinear dynamics of fault slip, favored by a short weakening distance and a large seismogenic zone. Megathrusts intersect the cold, brittle crust at low angle, providing a large area for great and giant ruptures, but also the conditions for complex earthquake cycles. This has important implications for understanding seismic hazards, because no temporal window of observation may capture a representative seismic cycle or super-cycle. The interactions between the seismogenic zone and the neighboring segments introduce dramatic changes in the seismic cycle and the mechanical system never converges to a limit cycle. In addition, viscoelastic flow in the asthenosphere produces a subtle modulation of the recurrence times and, in some case, rupture styles, on the megathrust.
The proposed mechanical model of the Nankai Trough highlights some important differences with other subduction zones. For example, the Cascadia subduction zone does not feature the type of long-term slow-slip events found in the Nankai. The structural layout of the Middle American trench differs crucially, with longterm and short-term slow-slip events taking place along a large-scale ramp-flat section of the megathrust. In the Japan trench, giant ruptures break all the way to the trench, a type of event that does not appear in the current simulations. The development of long-term slow-slip events in our model requires the presence of a wide metamorphic belt in the mantle wedge corner. Perhaps due to its temperature profile and degree of metamorphism, this region may be under-developed in Cascadia or along the Chilean trench.
The model does not explain the formation of tsunami earthquakes in some accretionary regions. Tsunami earthquakes may occur in a more sophisticated model that better resolves the internal stratification of the accretionary wedge (Leggett et al. 1985;Moore et al. 1990;Park et al. 2000;Kodaira et al. 2000;Takahashi et al. 2003;Ito et al. 2009), with different properties in the frontal prism, e.g., high-angle ramps, smectite and illite layers (Saffer and Marone 2003;den Hartog et al. 2012). The model does captures the unstable nature of the accretionary wedge, but many different rupture styles and types of seismic cycles can be obtained in an unstable domain, depending on the details of the parameters and where they fall in phase space. In addition, even with the same frictional parameters, the dynamics of two accretionary prisms may be grossly dissimilar, simply based on their geometry, i.e., width on the megathrust, dip angle of basal décollement, distance to the free surface (sediment thickness), and degree of internal deformation. Therefore, we expect different behaviors in a sediment-thick Lesser Antilles subduction zone and a sediment-starved northern Chile margin, for example. The fabric of the upper plate also varies greatly within the Nankai Trough, with a mega-splay fault bounding the accretionary prism from the forearc offshore the Kii Peninsula (Park et al. 2002), with many reverse faults occupying the entire outer wedge (Park et al. 2010;Kamei et al. 2012;Tsuji et al. 2014Tsuji et al. , 2015Tsuji et al. , 2017Shiraishi et al. 2019). We therefore anticipate important along-strike variations with different rupture styles in the accretionary region even within the Nankai Trough.

Conclusions
The down-dip segmentation of rupture styles at subduction zones indicates a strong control from the thermomechanical structure and composition of the upper plate. At the Nankai Trough, the development of slow earthquakes, large earthquakes, and long-term and shortterm slow-slip events can be reconciled with a structural model that considers the intersection of the megathrust with the accretionary prism, the arc crust, a metamorphic belt, the mantle lid, and some important isotherms coinciding with the stability transition of granitic rocks and peridotite (Fig. 11). The geological structure provides the boundary and physical conditions for the unfolding of the seismic cycle. Using the integral method within a two-dimensional approximation, we simulate the seismic cycle in the viscoelastic medium of the lithosphereasthenosphere system. With a choice of parameters informed by the geological setting, laboratory experiments, and theoretical insights, we reproduce several key rupture styles with their structural setting, providing new insights into the coupled dynamics of slow and fast ruptures and distributed deformation at subduction zones.
Our model illustrates how the geological structure of the overriding plate controls the down-dip segmentation of ruptures styles, but also how the nonlinear behavior of fault dynamics and the interactions of various parts of the mechanical system induce additional complexity. For example, the recurrence time and magnitudes of earthquakes and slow-slip events exhibit great variability. The short-term slow-slip events can be accelerated for several years following a great earthquake, and the long-term slow slip can be suppressed for decades in the early stage of a super-cycle. This suggests the possible sudden reactivation of slow slip after a long period of creep following a great earthquake. The propagation of creep waves in the shallowest segment reflects the instability of shallow slip and a potential for tsunami-genic ruptures.
Our modeling approach is well suited to simulate fault dynamics within a complex geological framework. Future work will expand the structural complexity. For example, coupling in-plane and anti-plane fault dynamics may be important to represent the oblique convergence at a plate boundary, which is partitioned into the megathrust and an along-arc strike-slip fault in many regions. Addition of splay faults, backthrusts, and normal faults may better represent the seismicity at subduction zones. Eventually, three-dimensional models will be required to capture along-strike segmentation and to produce more realistic simulations of seismological and geodetic data.

Funding
This study was funded by the National Science Foundation under grant no. EAR-1848192.

Availability of data and materials
The code and input files used in the study are available upon request.

Competing interests
The authors declare that they have no competing interests.  The full spectrum of simulated slip modes including full and partial ruptures, slow earthquakes and slow-slip events are distributed within down-dip segments along the megathrust fault. The interaction between neighboring segments leads to additional complexity, including triggered slow slip and slow earthquakes, and aperiodic recurrence patterns of earthquakes and slow-slip events