Stress evolution during the megathrust earthquake cycle and its role in triggering extensional deformation in subduction zones

Abstract Great (moment magnitude Mw ∼8.0 and larger) subduction megathrust earthquakes are commonly followed by increased rates of normal faulting seismicity. Extensional activity within the subducting slab is amplified when megathrust slip propagates close to the trench, and forearc extension is triggered by the largest magnitude (Mw 8.5 and larger) events. To better understand these observations, we develop an earthquake cycle model with a realistic slab geometry and stresses that are in balance with plate interface slip and bulk viscous relaxation. The modeled stresses represent perturbations to the long-term background tectonic stress field. The steady-state inter-seismic earthquake cycle stresses are compressive and their magnitudes depend on the interface locking configuration, from 25 MPa for a fully locked seismogenic zone to 5 MPa for discrete asperities on an otherwise unlocked plate interface. The co-seismic slip and corresponding co-seismic stress changes are similar in these models, independent of the locking configuration. The co-seismic stress change magnitudes are up to 5 MPa. This implies that the earthquake cycle stresses after the mainshock would still be compressive in the case of a continuous seismogenic zone, but would be widely reset to zero in the case of discrete asperities. Models with co-seismic slip confined to shallow depths produce tensional stress changes in the slab and the forearc near the trench, whereas events rupturing the base of the seismogenic zone produce tensional stress changes limited to region immediately surrounding the rupture. The locations of normal faulting aftershocks generally correlate well with the tensional co-seismic stress changes greater than 1 MPa. Following the mainshock, rapid afterslip occurring down-dip of the megathrust rupture expands the region of relative extension, but reduces its magnitude and does not promote additional normal faulting events. Bulk viscous relaxation does little to the state of stress in the elastic parts of the model. Continued plate convergence across a re-locked interface returns the system to the pre-earthquake state of compression.


Introduction
The long-term stress field (and corresponding deformation) in subduction zones is largely controlled by stationary or relatively slow-acting tectonic processes (Forsyth and Uyeda, 1975).This background tectonic stress varies systematically across the subduction system: bending dominates the subducting plate (Walcott, 1970); the accretionary wedge and forearc are often compressional (Dahlen, 1990); and global backarc stress conditions range from tensile to neutral to compressive (Heuret et al., 2011).Superimposed onto this background stress state is a time-varying stress field associated with the subduction plate interface (megathrust) earthquake cycle.In general, plate convergence across a locked megathrust results in the accumulation of compressive stresses in between earthquakes, which are subsequently released to drive co-seismic slip when the megathrust unlocks and an earthquake occurs.The rate of extensional deformation throughout a subduction system can therefore abruptly increase after a sufficiently large megathrust earthquake.The upper plate forearc may even transform (at least locally) from a state of compression to extension (Hardebeck, 2012).In this study, we constrain these timedependent stress variations in the subduction system and the corresponding deformation behavior caused by the megathrust earthquake cycle.
A common interpretation for increased extensional deformation occurring after subduction megathrust earthquakes is triggering by co-seismic Coulomb stress changes ( CS; Reasenberg and Simpson, 1992).For example, after the March 2011 moment magnitude (Mw) 9.0 Tohoku earthquake, CS at the locations of extensional aftershocks occurred was up to +5 MPa in the upper plate above the rupture zone and up to +1.5 MPa in the slab (Hiratsuka and Sato, 2011;Toda et al., 2011).Co-seismic extension in the Chile https://doi.org/10.1016/j.epsl.2020.1163790012-821X/© 2020 The Author(s).Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).and Peru forearcs also occurred in regions with CS up to +5 MPa (Loveless and Pritchard, 2008;Aron et al., 2013;Cortés-Aranda et al., 2015).Enhanced extensional deformation in the slab is associated with the positive CS from megathrust ruptures that propagate near to the trench (Lay et al., 2009;Herman et al., 2017;Sladen and Trevisan, 2018).Implicit in these CS studies is that the pre-mainshock (inter-seismic) stresses are unknown or even irrelevant to the triggering process.However, a small tensile stress change, if superimposed onto a much larger compressive stress, is unlikely to result in extensional faulting.Therefore, it is also important to constrain the inter-seismic state of stress prior to the megathrust event.
One way to estimate the magnitude of differential stresses before a megathrust earthquake is by combining the co-seismic stress change calculation with the change in focal mechanisms of smaller events.The rotation in principal stresses inferred from these focal mechanisms (or other stress indicators such as borehole breakouts; e.g., Lin et al., 2011) can be up to 25 • in places where thrust faulting gives way to normal faulting (e.g., Hardebeck, 2012;Xu et al., 2016).In the upper plate offshore Japan, this approach suggests that the pre-2011-earthquake differential stresses were 10-20 MPa (Hasegawa et al., 2012).Onshore Japan, there were fewer premainshock earthquakes and the co-seismic stress changes had smaller magnitudes (≤1 MPa).Normal faulting aftershocks in this region either represent major rotations of low pre-mainshock differential stresses (Kato et al., 2011;Yoshida et al., 2012) or a pre-mainshock stress field that already favored extension and was further induced by the 2011 Tohoku earthquake (Imanishi et al., 2012;Yoshida et al., 2018).Recent studies of the Japan subduction zone also include the effects of post-seismic processes on these stresses, for example concluding that afterslip enhanced the coseismic stress changes over the year following the earthquake, and subsequent bulk viscous relaxation redistributed stresses within the shallower elastic regions (Becker et al., 2018).These stress rotation analyses provide an estimate of the ambient stress (given the underlying assumption that the seismicity consistently reflects this stress field), but do not quantitatively evaluate their source(s) in the time period before the subduction megathrust earthquake.In particular, the contribution to the pre-earthquake stress field from earthquake cycle processes is still poorly constrained.
Some explanations for local co-seismic extension have been proposed, especially in the upper plate forearc close to the trench.McKenzie and Jackson (2012) argue that co-seismic gravitational collapse played a major role in extension of the Japan forearc near the trench.Wang et al. (2019) present a Coulomb wedge model of pre-and post-seismic stresses in the Japan subduction zone to explain the apparent stress rotations in the forearc wedge and from this infer a frictionally weak megathrust corresponding to low differential stresses (<10 MPa) in the wedge.These analyses help to understand the deformation processes acting in the forearc, but do not provide a comprehensive explanation of the temporal evolution of extension throughout the subduction system.
Earthquake cycle models with consistent pre-, co-, and postseismic stresses can help explain some of the observed patterns of extensional deformation.These models show how stresses are accumulated and released in and around stick-slip asperities (Dmowska et al., 1996), and reveal feedbacks between megathrust slip and intraplate deformation (van Dinther et al., 2014).In this study, we use observations of extension occurring near great (Mw 8.0 and larger) megathrust earthquakes to motivate the development of a generic subduction earthquake cycle model.We focus particularly on the impact that plate interface slip behavior over time has on the stress field.This allows us to address the question: what characteristics of the earthquake cycle are compatible with the timing and distribution of extensional deformation in subduction zones?

Observations
We first examine the spatial and temporal patterns of seismicity relative to 13 Mw 8.0 and larger subduction megathrust earthquakes from 2001 to 2015 (Fig. 1; Fig. S1 and Supplement).The magnitudes and slip distributions of the great earthquakes are from the United States Geological Survey (Hayes, 2017).We also include slip models from the SRCMOD Finite-Source Rupture Model Database (http://equake -rc .info/SRCMOD).The locations, magnitudes, and moment tensors of the seismicity comes from the Global Centroid Moment Tensor (GCMT) catalog (Ekström et al., 2012).We take the events from the GCMT catalog in a time window of 10 yr before to 5 yr after each mainshock.Although we explored using local catalogs with lower magnitudes of completeness, the results did not differ significantly from what we present here because the GCMT catalog contains >95% of the moment release.Any event with a normal faulting fraction ≥0.70 is defined as a normal faulting earthquake and events with thrust faulting fraction ≥0.60 are thrust faulting events (Frohlich, 1992).
We determine the plate in which each normal faulting aftershock occurred by comparing its centroid depth (±5 km) to the depth of the subduction plate interface at that location (Hayes et al., 2018;Fig. S1).
Prior to the mainshocks, the subduction margins are generally dominated by normal faulting in the outer rise, thrust faulting in the forearc, and mixed behavior in the arc to backarc, although the detailed distribution differs in each setting (Fig. S1).Out of the 13 earthquake sequences examined, nine have increased normal faulting aftershock activity in the vicinity of the rupture zone within days to weeks, and show elevated extensional moment release (relative to pre-earthquake rates) over the following months to years (Fig. 1A; Fig. S1).The decade before these nine events contains 0-60% of the total normal faulting moment release (from the 15-yr window around the mainshock).In all but two sequences (2007 Sumatra and 2015 Illapel), the pre-mainshock extensional moment release is less than 30% of the total.The moment release of normal faulting aftershocks within 6 months of the mainshock accounts for no less than 30% of the total.In the five years following these mainshocks, the normal faulting moment release rate slows, but stays at a comparable or higher rate than in the decade before the earthquake.
The cumulative seismic moment released by normal faulting earthquakes in the five years after the mainshock is typically 1-4 orders of magnitude smaller than the mainshock moment release, but extension can account for a significant fraction of the total aftershock deformation, especially after the largest megathrust earthquakes.In the aftershock sequences of the three largest events (2004Sumatra, 2010Maule, and 2011 Tohoku), the cumulative extensional seismic moment is 30-40% of the total aftershock moment release.In the other aftershock sequences, normal faulting seismicity accounts for 10% or less of the aftershock moment release.The remaining aftershock seismic moment is dominated by thrust faulting on the megathrust.The main exception is in the 2006 Mw 8.3 Kuril Islands sequence, in which the largest normal faulting aftershock was Mw 8.1.Other than this event, normal faulting accounts for less than half of the aftershock moment budget.
The locations of elevated normal faulting aftershock activity are correlated with the moment magnitude and centroid depth of the mainshock (Fig. 1B; Fig. S1).The centroid depth is related to the down-dip position of slip in the seismogenic zone, and the moment magnitude is proportional to both rupture dimensions and slip magnitude.The four largest events (Mw 8.7-9.1),which had the greatest slip and likely ruptured the full down-dip extent of the locked megathrust, are followed by normal faulting aftershocks in both the upper plate and the slab.The smaller events have nor- Many sequences show an immediate spike in moment release after the mainshock (dashed line), followed by a period of increased extensional seismic activity relative to the pre-mainshock period.Events match Panel (B) and are labeled by the year of the mainshock.(B) Plate of origin of normal faulting aftershocks and extensional surface structures following a great megathrust earthquake as a function of mainshock moment magnitude and centroid depth.Colors indicate which plates the normal aftershocks occur in.Dotted circles indicate events which had active surface structures observed, including faulting and tension cracks.The datasets used to create this figure are shown in the Supplement.(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)mal faulting aftershocks only if their centroid depths are less than ∼20 km, and then these aftershocks only occur within the slab (Sladen and Trevisan, 2018).Megathrust earthquakes smaller than Mw ∼8.5 that rupture only the deeper part of the plate interface do not have normal faulting aftershocks in the GCMT catalog.The exception is the 2013 Mw 8.0 Santa Cruz Islands earthquake; however, this event was also in a complex tectonic setting near a slab edge that hosts a higher than normal rate of background extensional faulting activity (Lay et al., 2013).
Extensional fault slip and surface cracking of the forearc are sometimes observed after some megathrust earthquakes, including some cases where there is no record of normal faulting aftershocks (Fig. 1B).Normal faulting was observed in the Atacama (Chile) fault system immediately following the 1995 Mw 8.0 Antofagasta event (Delouis et al., 1997) and above the down-dip limit of the 2010 Maule earthquake (Aron et al., 2013).An Mw 6.6 normal faulting aftershock of the 2011 Tohoku earthquake nearly 200 km from the trench produced 14 km and 16 km of surface rupture on two segments (Mizoguchi et al., 2012).The finite offsets on these faults are interpreted to reflect repeated activation over thousands of earthquake cycles.Tension cracks exposed throughout southern Peru and northern Chile have broadly uniform orientations parallel to subduction strike (extension perpendicular to strike), and cluster around mapped faults; this is interpreted to reflect an earthquake cycle-related origin (Loveless et al., 2005).The tension cracks were seen to be active immediately after the 2001 Arequipa and 2014 Iquique earthquakes (Scott et al., 2016), despite these events having no normal faulting aftershocks in the GCMT catalog.
The relationship between megathrust earthquakes and extensional forearc structures in the absence of recent activity is more speculative.Several other normal faults in eastern Japan may accumulate their offsets during megathrust events over multiple earthquake cycles (Regalla et al., 2017).Similarly, mapped normal faults in southwestern Mexico (Gaidzik et al., 2016) and southern Iran (Normand et al., 2019) occur in areas otherwise characterized by shortening.Normal faults are also identified in marine seismic surveys offshore New Zealand 40-50 km from the Hikurangi trench (Böttner et al., 2018).Dielforder et al. (2015) even suggest there is preserved evidence of extension caused by megathrust earthquakes in an exhumed accretionary complex from the Alps.However, it remains unclear whether forearc extension is limited to the time period shortly after a large megathrust earthquake or can occur throughout the inter-seismic period (Loveless et al., 2010).
These observations show how the rate of normal faulting seismicity in subduction zones increases significantly following large megathrust earthquakes.The magnitude of the mainshock and the depth extent of its rupture appear to be the primary factors controlling the distribution and the moment release of these extensional earthquakes.The largest events trigger normal faulting within the slab and the upper plate, whereas smaller events trigger normal faulting only in the slab and if the rupture propagates close to the trench.For larger mainshocks, normal faulting makes up a larger percentage of the total aftershock moment release.After the initial spike in extensional deformation rate, normal faulting seismicity continues at an elevated level for years to decades after the mainshock.When the subduction system is in the inter-seismic stage, normal faulting earthquakes can still occur regularly in some locations, such as the outer rise, but do so at a lower rate than in the aftermath of a large megathrust event.

Model setup
The systematic, global occurrence of extension associated with and triggered by megathrust earthquakes suggests that it is partially controlled by common physical attributes of the earthquake cycle, not by region-specific geometric or rheological characteristics.We therefore develop a generic subduction zone earthquake cycle finite element model (FEM) to investigate these relationships, building off of the models of Govers et al. (2018).The mechanical equilibrium equations are solved using the GTECTON FEM platform (Govers et al., 2018), and our results are not sensitive to further mesh and time step size refinements.
The model slab has the shape of the subducting Pacific plate taken along a transect perpendicular to the trench through the centroid of the 2011 Tohoku mainshock (143.05• E, 37.52 • N) (Hayes et al., 2018;Fig. 2).This geometry is generally representative of subducting slabs at the depths where megathrust earthquakes occur, with a low dip near the trench increasing smoothly to the maximum dip (∼30 • ) at depth.The seismogenic zone is the area of the plate interface from 15 to 40 km depth, where plate interface thrust earthquakes usually occur (Hayes, 2017).We extrude after the 2011 Tohoku earthquake.Primary afterslip broadens the zone of 1 MPa extension (relative to the late inter-seismic stage), but reduces its maximum magnitude.The mantle wedge is now highly stressed.(Bottom Row) Bulk viscous relaxation of the mantle wedge strongly reduces stress magnitudes in the wedge, but has a limited effect on the stresses in the elastic upper plate.Re-locking of the plate interface gradually returns the system to horizontal compression.
this cross-section 1000 km along strike (the y-axis) so we can test the effects of lateral variations.
Deformation in the model is driven by kinematic boundary conditions.The front and back of the subducting plate move at 90 mm/yr parallel to the slab dip, while the back of the upper plate at x = −600 km is fixed.The lateral sides of the model are constrained to move perpendicular to strike only.An isostatic restoring force is applied to the base of the slab.In our model, the earthquake cycle is strictly periodic, and each inter-seismic stage lasts for 350 yr, corresponding to 31.5 m of relative motion between the plates.This repeating condition approximates a yield strength criterion.
The elastic properties of the model are taken to be uniform (shear modulus = 40 GPa and Poisson's ratio = 0.25).The 40-kmthick upper plate and the subducting slab are cooler, so we treat these as purely elastic.The warmer mantle wedge beneath the upper plate is visco-elastic, with a linear Maxwell rheology (viscosity = 10 19 Pa s, corresponding to a relaxation time of 8 yr).We do not model the effects of the asthenosphere beneath the slab, because this was shown to have a small influence in Govers et al. (2018).
The plate interface is represented by a narrow fault/shear zone using slippery nodes (Melosh and Williams, 1989).The interface shallower than 15 km, adjacent to the accretionary wedge, is treated as mechanically weak (unlocked) throughout the entire earthquake cycle.The seismogenic zone is locked in the time period between earthquakes to build slip deficit and stresses, and then unlocked at selected times to allow co-seismic slip and stress release.The modeled earthquake has a slip distribution and corresponding deformation field compatible with global observations (Govers et al., 2018).The deeper plate interface (depths greater than 40 km) is modeled as a viscous shear zone (Kneller et al., 2005).At the high strain rates of the co-seismic rupture, the downdip shear zone behaves elastically (locked).After the earthquake, we re-lock the seismogenic zone and simulate relatively rapid post-seismic slip on the down-dip shear zone (Muto et al., 2019) by unlocking this region.We refer to this stage as "primary afterslip."Although we model this process as occurring instantaneously, post-seismic slip on the deeper plate interface (or shearing within a narrow, low viscosity shear zone) likely takes several years after the mainshock.The down-dip shear zone remains unlocked as post-seismic viscous relaxation occurs in the mantle wedge and the model returns to the inter-seismic stage.
The model results represent perturbations to the long-term background stress field.We build up steady state stresses ("earthquake cycle stresses," σ E Q C ) in the model over 15-20 earthquake cycles.Each earthquake in this spin-up stage ruptures the entire length of the seismogenic zone.Once in steady state, each cycle is identical to the previous.The motivation for spinning up the model is that in steady state the stress field perturbations are consistent with the plate interface slip and bulk deformation at all times.After the spin-up is completed, we simulate an Mw ∼9 earthquake (∼400 km long).We then model the primary afterslip following this event and 20 yr (2.5 relaxation times) of early post-seismic relaxation.
The finite element models in this study are based on solving the equations for static mechanical equilibrium with a (visco-)elastic rheology.The elastic response to the imposed boundary conditions inherently includes the effects of flexure, in addition to the deformation produced by megathrust slip and viscous relaxation in the mantle wedge.This is particularly relevant in the subducting slab as it bends into the trench, where flexural stresses are most strongly expressed in the observed seismicity (Craig et al., 2014) and the geological record (Walcott, 1970).We can therefore use our model results to determine the importance of bending stress changes within the slab over the earthquake cycle.
The model as previously described is the reference model.We investigate the sensitivity to varying the plate interface rheology as represented by the pattern of locking and unlocking.The plate interface near the base of the seismogenic zone may be relatively weak across a range of strain rates (by deforming through both brittle and ductile mechanisms or having conditionally stable frictional properties; Peng and Gomberg, 2010;Lay et al., 2012;Scholz and Campos, 2012), so we run a model in which the plate boundary from 40 to 45 km depth remains unlocked throughout the earthquake cycle.Lateral areas on the interface may also exhibit conditionally stable sliding behavior, so we run a model in which the center 50 km of the seismogenic zone is always unlocked.The asperities in the seismogenic zone that lock and unlock may be discrete and isolated rather than broadly continuous, so we run a model with 12 40-km square asperities distributed evenly along strike (i.e., each separated by 40 km).Each asperity corresponds to the size of an Mw ∼7 earthquake, and their size and spacing are consistent with inter-seismic coupling models (Herman et al., 2018).In these models, the permanently unlocked regions near locked sections still accumulate (partial) inter-seismic slip deficit because they are adjacent to the locked patches, and can also slip during the earthquake (Herman et al., 2018).Finally, we test the effects of the mainshock rupturing only part of the seismogenic width with a shallow co-seismic slip model (0-30 km), and a deep co-seismic slip model (30-40 km).

Results
We show the spatial distribution of inter-, co-, and post-seismic σ E Q C in the model along a vertical cross-section taken perpendicular to the trench (Fig. 2) and a horizontal cross-section at 20 km depth (Fig. 3).The potential for seismicity at each point in the model is quantified with the maximum shear stress (the second invariant of the deviatoric stress tensor): which is a measure of the tractions driving fault slip.The type of faulting at each point is determined by treating the modeled principal stress axes as the principal axes of the earthquake focal mechanism and using the ternary classification of Frohlich (1992).

Reference model
In the time step before the earthquake, the σ E Q C mainly encourage thrust faulting (Fig. 2A).From the backstop to the base of the locked seismogenic zone, the upper plate favors thrust faulting with stress magnitudes of 15-25 MPa.The upper plate has a narrow (∼10 km) section immediately above the base of the locked zone where normal faulting is promoted.From the locked zone to the trench, the σ E Q C magnitudes in the upper plate decrease until they are <1 MPa.The mantle wedge at this stage of the cycle is completely relaxed.In the slab shallower than 50 km, the σ E Q C favor thrust faulting with slightly lower magnitudes than in the upper plate (5-20 MPa).Down-dip of the locked zone, the slab is under slab-parallel tension.
The static co-seismic stress changes produced by our model are consistent with previous analyses (Fig. 2B; Fig. 3A; e.g., Hiratsuka and Sato, 2011;Toda et al., 2011).The earthquake promotes normal faulting over a region extending up to 100 km from the rupture zone and spanning the rupture length with stress change magnitudes of up to 5 MPa.This means that the distribution of σ E Q C in the reference model is largely unchanged after the earthquake, because the steady state σ E Q C are generally >10 MPa.Minor visible changes to σ E Q C are caused by the earthquake, namely: the width of the normal faulting region in the upper plate immediately above the seismogenic zone is expanded by ∼15 km; and the ∼20-km long nose of the mantle wedge becomes highly stressed immediately down-dip of the rupture tip.
Primary afterslip on the deeper plate interface broadens the zone of >1 MPa extension in the upper plate by ∼50 km and extends much of the mantle wedge (relative to the pre-earthquake state) (Fig. 2C).At the same time, this deep slip reduces the coseismic tensional stress changes above the base of the seismogenic zone by 1-2 MPa.The σ E Q C in the slab at depths <50 km are unaffected by primary afterslip, but the lobe of extension beneath the base of the seismogenic zone disappears.Again, the background σ E Q C have much greater magnitudes and are visibly unchanged by primary afterslip.After 20 yr (2.5 relaxation times) of bulk viscous relaxation, the stresses in the mantle wedge have decayed to <1 MPa (Fig. 2D).
This viscous relaxation itself produces relatively little effect on the state of stress in the upper plate (despite potentially large trenchward surface velocities, e.g., up to 10 cm/yr in Honshu; Govers et al., 2018;Muto et al., 2019).The upper plate extension relative to the pre-earthquake state does decrease throughout the postseismic stage; this is primarily due to the re-locking of the plate interface and continued plate convergence.Because the relative plate velocity is constant, the return towards the late inter-seismic state of stress is linear in time.

Down-dip unlocked segment
When the plate interface from 40 to 45 km depth always remains unlocked, the resulting σ E Q C magnitudes are 10-15 MPa lower than in the reference model (Fig. S2).This is because much lower stresses build up around the base of the seismogenic zone if this region is persistently weak instead of a sharp rheological contrast.Immediately before the earthquake, the σ E Q C magnitudes in the upper plate are less than 10 MPa, and still mainly promote thrust faulting as in the reference model (Fig. S2A).The narrow zone of extension above the base of the seismogenic zone has given way to compression in this model, while the area near the trench continues to have low σ E Q C magnitudes.The σ E Q C in slab are also reduced to less than 5 MPa, with a similar pattern as in the reference model.The co-seismic stress changes are similar in magnitude and spatial distribution to that of the reference model (Fig. S2B), only without the stress lobes at the base of the rupture zone that are associated with the sharp down-dip slip gradient in the reference model.Because the co-seismic stress changes have magnitudes closer to the σ E Q C in this model, after the earthquake, the σ E Q C As a result, although there is a good correlation between the locations of tensional co-seismic stress changes and normal faulting aftershocks, the steady state stresses favor compression and thrust faulting rather than extension.(B) Results from the discrete asperity model.The area with co-seismic slip magnitudes greater than 5 m is similar to the reference model, hence the similarity in modeled co-seismic extension.In this model, the earthquake reduces the steady state stresses to nearly zero, resulting in an improved correlation between the locations of normal faulting earthquakes and the areas of tensional co-seismic stress changes and low σ EQ C . in the upper plate and much of the shallow slab are reduced to 1-3 MPa in the regions where co-seismic stress changes promote extension.Specifically, a corridor of reduced compressive stresses extends landward from the base of the seismogenic zone to the surface over ∼100 km; a narrow (∼15 km wide) region above the seismogenic plate interface promotes normal faulting; and the earthquake cycle stresses are reduced to zero in much of the slab near the trench.The post-seismic effects are similar to that of the reference model, with lower σ E Q C magnitudes (Fig. S2C-D).

Lateral unlocked segment
In the model with the center 50 km of the interface always unlocked, the σ E Q C magnitudes are ∼5 MPa smaller than those of the reference model, but still 5-10 MPa greater than in the downdip unlocked model (Fig. S2).As in the reference model, the σ E Q C magnitudes are much larger than the stress changes, so they do not change significantly in each stage of the cycle.The style of faulting promoted by the σ E Q C before the earthquake are generally similar to the other models, except in the upper plate above the unlocked section of the seismogenic zone, where they favor strikeslip or oblique mechanisms (Fig. S2E).This is because the cross-section is taken through the central unlocked section, and the edge effects of the adjacent locked areas rotate the stress field.
The style of stress changes (particularly those favoring extension) associated with the central unlocked area does differ from that of the other two models.The co-seismic stress changes that promote extension in the upper plate extend ∼350 km from the trench (Fig. S2F), compared to ∼250 km for the other models.Coseismic extension also occurs in the slab from under the trench to the outer rise, but unlike the other models the region beneath the seismogenic zone does not undergo co-seismic extension.The effects of primary afterslip and post-seismic relaxation are similar to those of the reference model, but the larger spatial footprint of tensional stress changes persists well into the post-seismic period (Fig. S2G-H).

Discrete asperities
The heterogeneous, discrete asperity model combines characteristics of the down-dip and lateral unlocked models and this is reflected in the resulting steady state σ E Q C (Fig. 2E).The magnitude of these stresses is the lowest of any model (<5 MPa) and they still favor dominantly thrust faulting throughout most of the upper plate and much of the shallow slab.The nose of the upper The earthquake stress changes are >1 MPa and favor extensional faulting throughout most of the upper plate, with the exception of the nose near the trench (Fig. 2F; Fig. 3B).The coseismic extension in the slab extends from under the rupture to 100 km beyond the trench.The general pattern of co-seismic stress changes is similar to that of the reference model because of the similarity in slip.This fault slip produces a ∼10 km wide zone of σ E Q C favoring normal faulting above the seismogenic zone, similar to previous models.As in the down-dip unlocked model, the compressive σ E Q C in the corridor from the base of the seismogenic zone to the surface of the upper plate at x = −250 km are reduced, and in this model they are <1 MPa.
After the earthquake, the evolution of stresses in the postseismic period closely mirrors that of the lateral unlocked segment model: primary afterslip broadens the zone of extension (Fig. 2G), then bulk viscous relaxation reduces stresses in the mantle wedge while re-locking recovers inter-seismic compression (Fig. 2H).

Deep and shallow ruptures
In order to avoid the stress concentration that forms at the deep edge of the seismogenic zone with a sharp rheological contrast, we base the deep and shallow rupture models on the down-dip unlocked model.When the earthquake ruptures the plate interface from 30 to 45 km depth, both its slip magnitude and rupture area are reduced compared to the down-dip unlocked model.As a consequence, the co-seismic stress changes promoting normal faulting have a smaller spatial footprint; they only extend to a depth of 20 km in the upper plate (instead of all the way to the surface) and they are confined to small lobes in the slab (Fig. 4A).The outer rise is notably unaffected by co-seismic extension in an earthquake that only ruptures the deep section of the seismogenic zone (Herman et al., 2017;Sladen and Trevisan, 2018).
The shallow earthquake ruptures the plate interface from the trench to 30 km depth.Like the deep earthquake, this shallow earthquake also has reduced slip magnitude and area compared to the earthquake in the down-dip unlocked model.However, because of the position of the rupture, stress changes favoring extension do reach the surface of the upper plate immediately above the rupture zone (Fig. 4B).The upper plate in the wedge within ∼50 km of the trench remains at a low stress level.Tensional stress changes also extend into the shallow (<10 km depth) part of the outer rise, with magnitudes greater than 0.5 MPa up to 100 km oceanward of the trench.

Slab bending
In all of these models, the slab is deflected downward as it enters the trench and moves along the megathrust into the mantle.This deflection results in flexural stresses within the slab, producing tension near the surface and compression at depth (Turcotte and Schubert, 2002).To estimate the accumulation rate of these bending stresses, we run a model for 350 yr (the length of one earthquake cycle) with same velocity boundary conditions, but no resistance on the plate interface (Fig. S3).The bending stresses are largest at the top and bottom surfaces of the slab.The maximum shear stress magnitude associated with this bending is 0.2 MPa.This is smaller than the other earthquake cycle stress changes, which exceed 1 MPa over large areas.These bending stress magnitudes are linearly proportional to the elastic thickness of the oceanic lithosphere, so a thinner slab would result in smaller bending stresses.

Discussion
Because the model is based on the geometry of the Japan subduction zone, we compare the results to the locations of normal faulting aftershocks of the 2011 Tohoku earthquake.The global similarities in extensional aftershock distribution relative to mainshock slip (especially the vigorous activity in the slab for shallow ruptures and the upper plate extension triggered by the largest events; Fig. S1) motivate us to interpret all of these sequences in the same framework.Despite the differences in plate boundary slip behavior, the magnitudes (up to ∼5 MPa) and patterns of co-seismic stress changes are similar in all of our models.There is a good correlation between the locations of normal faulting aftershocks and the patterns of modeled co-seismic stress changes (Fig. 2B and 2F; Fig. 3), as seen in previous studies (e.g., Hiratsuka and Sato, 2011;Toda et al., 2011).Extensional aftershocks dominantly cluster in the regions where the modeled stress changes promote normal faulting and have magnitudes larger than 1 MPa, and are generally absent from regions with lower stress change magnitudes or locations where the stress changes favor other mechanisms.
The correlation between extensional aftershocks and the modeled steady state σ E Q C depends on the pattern of plate boundary locking.The σ E Q C in the reference model have magnitudes up to 25 MPa due to the stress buildup around the sharp rheological contrast at the base of the seismogenic zone, which is required to produce repeated co-seismic slip on the plate interface.The co-seismic stress changes are nearly an order of magnitude too small to overcome these σ E Q C , so the model always favors thrust faulting.Therefore, the normal faulting aftershocks occur in regions where the earthquake cycle has superimposed persistently compressive stresses onto the background stress field (Fig. 2B; Fig. 3A).When we add permanently unlocked areas (representing conditionally stable, mixed brittle-ductile, or persistently weak rheologies) on the plate interface, the steady state σ E Q C magnitudes decrease.Lateral unlocked zones reduce the magnitudes to ≤20 MPa.If the base of the seismogenic zone can slide freely, the σ E Q C magnitudes are ≤10 MPa.Including both down-dip and lateral permanently unlocked areas (i.e., placing discrete asperities on the plate interface) results in the lowest steady state stresses (magnitude: ≤5 MPa).In the discrete asperity model, the σ E Q C are significantly lowered where the co-seismic tensional stress changes are >1 MPa, with the result that normal faulting aftershocks occur where the modeled steady state σ E Q C are ≤1 MPa (Fig. 2F; Fig. 3B).Depending on the precise nature of the background tectonic (deviatoric) stress field, any of our models (which are perturbations superimposed onto these background stresses) could be compatible with the observed spatial and temporal patterns of seismicity.Our models can still provide additional constraints on the range of possible background stress fields.
The deformation in the subducting slab indicates the broad dominance of bending stresses.Normal faults are visible on the seafloor in outer rises globally (Masson, 1991;Ranero et al., 2005), along with shallow extensional seismicity and, in some cases, thrust faulting seismicity just below the normal faulting events (Craig et al., 2014;Fig. S1).It is straightforward to show that the bending-related stresses responsible for this deformation (accumulated over geological timescales of at least 100s of Kyr) are much larger than the earthquake cycle-generated stresses, despite being relatively insignificant over a single cycle (Fig. S3; Fig. S5).For example, the radius of curvature of the Japan slab in our model is ∼1000 km, corresponding to maximum fiber strains of 0.02 (for a 40 km thick elastic beam; Turcotte and Schubert, 2002).This implies maximum shear stress magnitudes of ∼1 GPa (for Young's modulus = 100 GPa and Poisson's ratio = 0.25).This is significantly higher than the yield strength of rocks at lithospheric depths, so the slab is likely at or near its tensional yield strength within 10-15 km of the surface (and may also exceed its compressive yield strength at depths of 25-40 km).The precise yield strength depends on the deformation mechanism, material properties, fluid content, depth, and temperature, but is likely 10s to 100s of MPa (e.g., Turcotte et al., 1978;Chapple and Forsyth, 1979;Melosh and Raefsky, 1980).This is significantly larger than the earthquake cycle related compressive stresses that our models produce in this region (a maximum of 1-10 MPa during the interseismic stage, depending on the asperity configuration; Fig. 2).
Although the megathrust earthquake cycle is not the dominant long-term tectonic process driving deformation in the outer rise, the sequence of inter-seismic loading, co-seismic slip, and postseismic processes modulates the state of stress in this region over decades to centuries, as the varying rates of normal faulting seismicity and our model results indicate.Inter-seismic compression may inhibit shallow normal faulting in the outer rise while promoting deeper thrust faulting.Normal faulting aftershocks are easily triggered in the outer rise whenever the stress changes from the event promote extension in the bending slab; this occurs only when the megathrust slip propagates to shallow depths (Fig. 1B, Fig. 4A; Sladen and Trevisan, 2018).Post-seismic re-locking and convergence then result in the re-accumulation of compression.An example of this full cycle of outer rise faulting behavior was observed in the Kuril Islands; the region oceanward of the trench adjacent to the 2006 earthquake hosted thrust faulting events prior to 2006, normal faulting events for the two years after the earthquake, and by 2009 was back to hosting thrust faulting earthquakes (Lay et al., 2009).Our interpretation is that the megathrust earthquake cycle shifts the slab back and forth across its yield threshold, which is controlled primarily by its long-term bending.The focal mechanisms of intraslab normal faulting aftershocks dominantly reflect these background bending stresses, with T-axes nearly perpendicular to the trench (Fig. 3; Fig. 5).
A relevant distinction is that although a megathrust earthquake can cause horizontal extension in the outer rise of the subducting slab that triggers normal faulting, this is not equivalent to saying that it promotes additional bending.In fact, because the earthquake produces down-dip tractions on the top of the subducting slab where co-seismic slip occurs, the associated torque balances on the slab result in unbending at the outer rise (Fig. S4).The stress changes associated with this unbending are much smaller than the dislocation-generated stress changes and hence negligible compared to the full background bending stresses.Nevertheless, the earthquake itself actually slightly unbends the slab.
In contrast to the relatively uniform focal mechanisms of aftershocks within the slab, the normal faulting aftershocks in the upper plate forearc show significant spatial variability, with Taxes ranging from trench-perpendicular to trench-parallel over distances of <10 km in some places (Fig. 3).One possible source for such a spatially heterogeneous background stress field in the forearc is lateral differences in gravitational potential energy due to geological heterogeneity (e.g., Levandowski et al., 2017).These stresses vary on the scale of density variations, which can be small (1s-10s of km) compared to plate boundary-scale (10-100s of km) tectonic stresses.The dominant style of deformation in the forearc also varies over the earthquake cycle, from inter-seismic compression to co-seismic extension; this suggests that the background stresses are comparable in magnitude to the variations in stress throughout the earthquake cycle.Because deviatoric stresses in the upper plate are poorly constrained, these observations could be explained in the context of our reference model (Fig. 2B) by adding a tensional stress in the forearc exactly counteracting the compressive steady state earthquake cycle stresses.However, there is no clear candidate for such a convenient source of tensional stress.We therefore prefer the simpler interpretation provided by the discrete asperity model (Fig. 2F), which can be directly superimposed onto the heterogeneous background stress.When a large megathrust earthquake occurs and the trench-perpendicular horizontal compression in the forearc is removed, extension in various directions Fig. 5. Synopsis of the relationship between plate interface slip and extensional faulting in subduction zones interpreted from our modeling results.The seismogenic part of the plate interface consists of discrete stick-slip asperities embedded in an area that can slip easily in both the inter-seismic and co-seismic stages of the earthquake cycle.In addition, the region immediately down-dip of the asperities (the transition from brittle to ductile rheology) also slips, flows, or otherwise deforms throughout the earthquake cycle.These slip behaviors result in reduced stress magnitudes associated with the earthquake cycle (∼5 MPa) compared to a scenario where the plate interface is more broadly locked (∼25 MPa).As a consequence, the steady state earthquake cycle stresses are similar in magnitude to the co-seismic stress changes associated with an Mw ∼8.5 and larger earthquake.When such a giant earthquake ruptures the interface, the inter-seismic compression is removed, and widespread normal faulting aftershocks occur in both the subducting slab and the upper plate.The orientations of these aftershocks are controlled by non-earthquake cycle stresses.In the slab, bending stresses dominate, producing broadly uniform trench-perpendicular extension.In the upper plate, slab rollback, subduction erosion, or reorganization of the plate kinematics can produce similarly long-wavelength tensional stresses.However, the heterogeneity in normal faulting aftershock focal mechanisms suggests a short-wavelength source of tensional stresses, e.g., from gravitational forces associated with geological heterogeneity.The triggering of extensional faulting by earthquakes smaller than Mw ∼8.5 depends on the mainshock location on the megathrust.Shallow events that rupture near the trench easily trigger intraslab and outer rise seismicity by pushing the slab past its failure point defined by its bending.Deeper events do not typically trigger significant intraslab normal faulting seismicity and are also inefficient at triggering upper plate extension.They may still trigger small normal faulting aftershocks or other extensional deformation in the upper plate, such as tension cracks.can be favored by the remaining heterogeneous background stress field (Fig. 5).
In addition to generating a short-wavelength stress field in the upper plate, rheological heterogeneity can play an important role in how stress is distributed throughout the earthquake cycle.For example, the mantle wedge viscosity defines the temporal evolution of stresses below the upper plate.However, the earthquake and primary afterslip are largely insensitive to variations in the mantle wedge viscosity (Govers et al., 2018), so we do not expect the stress evolution within the overriding plate to be substantially different from the models in this study.Another potentially important aspect of rheological heterogeneity is the mechanical strength contrast between the upper plate and subducting slab.Deformation will be preferentially transferred to the mechanically weaker plate (as compared to the models presented here, which have uniform elastic properties).Recent megathrust earthquake sequences suggest that this mechanical contrast can play an important role in defining the co-seismic slip (Furlong et al., 2009) and the aftershock behavior (Lay et al., 2009).
Finally, in this study we focused on the rheological characteristics of the plate interface.In the preferred discrete asperity model, the stick-slip asperities are embedded in sections of the megathrust that are always unlocked (Fig. 5).These unlocked regions still accumulate inter-seismic slip deficit at 50% or more of the convergence rate and can slip during the earthquake (Herman et al., 2018), resulting in similar co-seismic stress changes as the models with more broadly locked seismogenic zones (Fig. 3).There are several processes that can explain how the areas in between asperities remain weak during both inter-seismic and co-seismic stages.In the inter-seismic stage, slip could occur in these regions by numerous smaller earthquakes rupturing the megathrust in be-tween full-asperity-rupturing events, producing earthquake supercycles (Sieh et al., 2008); episodic tremor and slow slip (Obara and Kato, 2016); or simply continuous stable sliding as we have modeled it here.To also remain unlocked during the time period we model as the co-seismic rupture, this region may be conditionally stable (Scholz and Campos, 2012), weaken at the high strain rates of the earthquake (Di Toro et al., 2004), or begin its afterslip before the seismogenic zone re-locks (the timing of re-locking is likely less than one year, but poorly constrained; Bedford et al., 2016).

Inter-seismic extension
Based on elastic half-space dislocation models, Loveless et al. (2010) argue that plate interface coupling can produce extension in the upper plate during the inter-seismic stage.Their models specifically suggest that the top 1-2 km of the upper plate could be in a state of horizontal tension.Our models disagree with this result, instead suggesting that the earthquake cycle-generated stresses in the upper plate should generate broad compression during the inter-seismic stage.This is mainly a consequence of the elasticover-viscous rheology (compatible with the thermal structure of the subduction system; Govers et al., 2018) combined with the shallowly dipping geometry of the plate interface, which produces horizontal shortening in the upper plate.Although the upper plate also undergoes flexure because of this geometry, this does not overcome the compression generated by shortening.We do note that the surface of the upper plate ∼300 km from the trench has low inter-seismic stress magnitudes in the discrete asperity model, which may be overprinted by stresses in the arc.In general, our models suggest that (a) extensional deformation observed during the inter-seismic stage is due to sources of stress other than the σ E Q C from earthquake cycle processes such as the gravitational collapse of high elevation, slab rollback, or subduction erosion (von Huene et al., 2004), and (b) in the absence of other evidence, extensional structures in the upper plate are most likely to be active immediately after a large megathrust earthquake.

Temporal evolution of extension in subduction zones
Our models include the effects of primary afterslip on the deeper plate interface (shear zone) and bulk viscous relaxation of the underlying asthenosphere wedge.When we allow afterslip on the deeper plate interface, this broadens the zone of extension in the upper plate, but reduces the maximum magnitude of the tensional stresses that were generated by the earthquake (Fig. 2C, Fig. 2G).Afterslip also produces a small increase in the tensional stresses in the subducting slab.When we allow afterslip to propagate laterally along strike from the rupture zone, we observe additional upper plate and intraslab extension propagating with the afterslip.Because of the small effect of primary afterslip on the tensional stresses in the model, we expect it to have a similarly limited impact on the frequency and distribution of extensional deformation.For example, in the Tohoku aftershock sequence, normal faulting events occurred in approximately the same locations for years after the mainshock (Fig. 2D, Fig. 2H).
Bulk viscous relaxation has little effect on the state of stress in the elastic parts of the model, despite a large kinematic effect visible at the surface (Govers et al., 2018;Muto et al., 2019).Because of its small effect on stresses in the elastic regions we do not expect this process to contribute significantly to the advancement or delay of normal faulting seismicity.In contrast, the relocking of the plate interface (and continued relative plate motion) appears to be the key process that re-establishes inter-seismic compression.This compressive stress accumulation occurs linearly in time, so we expect the decay in normal faulting seismicity rate (and return to compression) after a megathrust earthquake to follow a roughly similar trend (Fig. 1A).If this re-locking occurs rapidly after the earthquake, then the normal faulting events will be only a transient occurrence lasting a short duration (on the timescale of the earthquake cycle) after the megathrust rupture.

Conclusions
The triggering of normal faulting aftershocks and extensional deformation by great subduction megathrust earthquakes has systematic patterns globally.The primary locations of normal faulting moment release depend on the mainshock slip distribution; shallow slip triggers extension in the slab and accretionary wedge, whereas deep slip triggers extension in the upper plate.However, deep slip is less efficient at triggering normal faulting seismicity; deep earthquakes smaller than Mw ∼8.5 trigger little to no extensional aftershocks, whereas Mw ∼8.0 shallow earthquakes trigger outer rise or intraslab events.Our numerical models of the subduction earthquake cycle reiterate that these normal faulting aftershocks are triggered by co-seismic stress changes.These stress changes generally correlate very well with the distribution of extensional seismicity.New in this study is that we consistently model the stress field associated with the megathrust earthquake cycle.These steady state σ E Q C that drive co-seismic slip in the models favor thrust faulting and convergent deformation in most locations throughout the earthquake cycle.These σ E Q C have maximum shear stresses up to 25 MPa when the entire plate interface is locked and unlocked during the earthquake cycle.If instead the seismogenic region of the plate interface has persistently weak (unlocked) sections surrounding relatively small, discrete asperities, as well as an easily deformable brittle-ductile transition zone at the base of the seismogenic zone, then the co-seismic stress changes are similar in magnitude to the stresses built up to drive earthquake slip (∼5 MPa).In all models, subsequent afterslip on the deeper plate interface broadens the volume with stress changes favoring normal faulting perpendicular to the trench, but reduces the magnitude of these stresses.Bulk viscous relaxation of the deeper, warmer regions of the subduction system does not significantly change the state of stress in elastic regions.The original, compressive state of stress is recovered dominantly by the re-locking of the interface, implying that the stresses are imposed essentially linearly over time.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1. (A) Normal faulting moment released over time for 13 Mw 8.0 and larger earthquake sequences, as a fraction of the total normal moment release in the 15 yr shown.Many sequences show an immediate spike in moment release after the mainshock (dashed line), followed by a period of increased extensional seismic activity relative to the pre-mainshock period.Events match Panel (B) and are labeled by the year of the mainshock.(B) Plate of origin of normal faulting aftershocks and extensional surface structures following a great megathrust earthquake as a function of mainshock moment magnitude and centroid depth.Colors indicate which plates the normal aftershocks occur in.Dotted circles indicate events which had active surface structures observed, including faulting and tension cracks.The datasets used to create this figure are shown in the Supplement.(For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Cross-sections at y = 500 km showing the stresses throughout the earthquake cycle.The hue indicates the focal mechanism type corresponding to the modeled stress tensor, and the intensity indicates the magnitude of the stress.The locations of the cross-sections are indicated on the perspective views of the model at the top; colors on the plate interface indicate the locking behavior (red: locked during inter-seismic loading, unlocked during the earthquake; orange: unlocked during the inter-seismic stage, locked during the earthquake; blue: always unlocked).The left column shows results from the reference model and the right column shows results from the discrete asperity model.(Top Row) In the late inter-seismic stage, steady state stresses in the reference model are 10-25 MPa and <5 MPa in the discrete asperity model, and generally favor thrust faulting.(Second Row) Co-seismic stress changes greater than 1 MPa are indicated by hatched areas.Circles indicate the first 6 months of normal faulting aftershocks of the 2011 Tohoku earthquake.In the reference model, the σ E Q C after the modeled earthquake are mostly unchanged.In contrast, the similar magnitude co-seismic stress changes in the discrete asperity model reduce the σ E Q C to nearly zero across wide regions.(Third Row) Circles indicate normal faulting aftershocks from 6 months to 5 yr

Fig. 3 .
Fig. 3. Cross-sections at z = −20 km showing co-seismic stress changes relative to the late inter-seismic stresses.The hue indicates the focal mechanism type corresponding to the modeled stress tensor, and the intensity indicates the magnitude of the stress.Black lines indicate the orientation of maximum co-seismic extension in the model, and red lines indicate the orientations of T-axes from normal faulting aftershocks of the 2011 Tohoku earthquake.Most events to the left of the plate interface occurred in the upper plate and most events to the right occurred in the subducting slab.Note the nearly uniform T-axes of intraslab events and the more heterogeneous orientations of upper plate events.Co-seismic slip contours of 1 and 5 m are shown as dashed lines.(A) Results from the reference model.The steady state stresses are much larger than the co-seismic stress changes.As a result, although there is a good correlation between the locations of tensional co-seismic stress changes and normal faulting aftershocks, the steady state stresses favor compression and thrust faulting rather than extension.(B) Results from the discrete asperity model.The area with co-seismic slip magnitudes greater than 5 m is similar to the reference model, hence the similarity in modeled co-seismic extension.In this model, the earthquake reduces the steady state stresses to

Fig. 4 .
Fig. 4. Modeled co-seismic stresses from earthquakes that rupture only part of the seismogenic width.The section of the interface that slips is shown as a bold line.The hue indicates the focal mechanism type corresponding to the modeled stress tensor, and the intensity indicates the magnitude of the stress.Co-seismic stress changes greater than 1 MPa are indicated by hatched areas and the 0.5 MPa stress change contour is indicated by the dotted line.(A) The modeled deep earthquake produces a lobate pattern of tensional stress changes around the rupture zone.The zone of >1 MPa extension does not extend to the surface (but potentially could if the earthquake has a sufficiently high slip magnitude).The shape of the lobes up-dip of the rupture shows that significant extension does not extend far into the outer rise, irrespective of the earthquake magnitude.(B) In contrast, the modeled shallow earthquake produces >1 MPa extension in the forearc and >0.5 MPa extension up to 100 km oceanward of the trench.plate near the trench has σ E Q C magnitudes <2 MPa and promotes strike-slip or oblique mechanisms, because the cross-section is taken through an unlocked section of the plate boundary (as in the lateral unlocked segment model).The slab on the ocean side of the trench at depths less than 40 km has low σ E Q C magnitudes, which promote extension.