Megathrust Stress Drop as Trigger of Aftershock Seismicity: Insights From the 2011 Tohoku Earthquake, Japan

Numerous normal‐faulting aftershocks in subduction forearcs commonly follow large megathrust earthquakes. Postseismic normal faulting has been explained by stress changes induced by the stress drop along the megathrust. However, details of forearc stress changes and aftershock triggering mechanisms remain poorly understood. Here, we use numerical force‐balance models combined with Coulomb failure analysis to show that the megathrust stress drop supports normal faulting, but that forearc‐wide aftershock triggering is feasible within a narrow range of megathrust stress drop values and preseismic stress states only. We determine this range for the 2011 Tohoku earthquake (Japan) and show that the associated stress changes explain the aftershock seismicity in unprecedented detail and are consistent with the stress released by forearc seismicity before and after the earthquake.

Previous work explained the stress reversal to result from the stress drop of the megathrust earthquake, i.e., the coseismic decrease in megathrust shear stress due to fault weakening processes (Di Toro et al., 2011;Hardebeck, 2012;Hasegawa et al., 2011;Scholz, 1998).The megathrust shear stress loading in the interseismic period causes compression of the forearc (Lamb, 2006;Wang & He, 1999), without which the forearc would experience deviatoric tension due to gravitational stresses resulting from density contrasts and margin topography.If the megathrust shear stress decreases during an earthquake, the loss in compression can cause a reversal from deviatoric compression to deviatoric tension and trigger normal faulting in the forearc (Cubas et al., 2013;Dielforder, 2017;Wang et al., 2019).However, the detailed stress changes caused by the stress drop of megathrust earthquakes like Tohoku remain unresolved.Consequently, it is unknown whether the stress change can explain the broad aftershock seismicity in the forearc.
Aftershock seismicity of megathrust earthquakes has been further investigated by Coulomb failure stress (CFS) models, which test whether the resulting stress changes of an earthquake promote or suppress failure on neighboring faults (Farías et al., 2011;Terakawa et al., 2013;Toda et al., 2011).CFS models for the Tohoku earthquake showed that the stress change promoted some of the aftershocks and likely increased forearc seismicity rates.Aftershock seismicity not promoted by the stress change was interpreted to have been triggered either by  c, d) Cross-sectional view.Gray dots are aftershock hypocenters from the Japan Meteorological Agency (JMA).Beach balls denote JMA focal mechanism solutions and are shown in profile view.Reverse faulting events in the outer marine forearc above the megathrust are likely poorly located interplate events (Nakamura et al., 2016) or the slab model is not accurate enough to classify them as interplate events.(e) Number of forearc earthquakes.Count includes events with magnitude ≥ magnitude of completeness (M c = 3.4,Text S3 in Supporting Information S1).
an increase in fluid pressure (Terakawa et al., 2013) or the presence of small faults of variable orientation (Toda et al., 2011).CFS models commonly determine the stress change in the forearc solely from the earthquake slip distribution and neglect the total stresses in the forearc resulting from gravitational and megathrust shear stresses.The total stresses determine, however, the forearc stress state (i.e., the magnitude and orientation of principal stresses) and we show in Section 2 that the preseismic stress state affects the stress change in the forearc that results from a stress drop on the megathrust.Thus, assessing whether a megathrust stress drop supports or inhibits failure in the forearc requires detailed information on the preseismic or postseismic stress state.
For Japan, estimates of the megathrust stress drop and preseismic and postseismic stress states exist (Brodsky et al., 2020;Brown et al., 2015;Wang et al., 2019;Yang et al., 2013), but the available data do not allow an unambiguous assessment of forearc stress changes.We therefore pursue an innovative modeling approach that determines the preseismic and postseismic stress states and related stress change compatible with broad aftershock triggering due to the megathrust stress drop.Using the aftershock distribution and fault kinematics as modeling constraints allows determining the precise stress conditions required for triggering the bulk forearc seismicity that occurred in the first months after the Tohoku earthquake.The model results are supported by independent stress drop observations.

Modeling Approach
We use a plane-strain finite-element model of force balance following Wang et al. (2019) to assess the stress change caused by the megathrust stress drop.The model is created with the software ABAQUS and comprises a rigid lower plate in frictional contact with an elastic upper plate representing the forearc (Figure 2a, Text S1 and Figure S1 in Supporting Information S1).ABAQUS computes the total stress tensor in the upper plate resulting from all applied boundary conditions, including gravity, isostasy, and friction along the plate contact (megathrust).We compute gravitational stresses for average densities of 2,800 and 3,300 kg m −3 for crustal and mantle parts, respectively, and seawater load using a density of 1,025 kg m −3 .Gravitational acceleration is 9.81 m s −1 .Friction along the megathrust is generated by displacing the lower plate and computed for the effective friction coefficient μʹ assigned to the megathrust, such that the megathrust shear stress τ is given by standard Coulomb friction (τ = μʹσ n , where σ n is normal stress).
To simulate the megathrust stress drop, the finite-element model includes a preseismic and a postseismic model step between which τ changes.The two model steps describe the stress state shortly before and after the megathrust earthquake, but not the stress state averaged over multiple earthquakes cycles as in other force-balance applications (Dielforder & Hampel, 2021;Lamb, 2006).Because τ also increases during earthquakes, e.g., due to fault strengthening processes downdip of the main rupture zone (Brown et al., 2015;Scholz, 1998), we divide the megathrust into fault weakening and fault strengthening segments (Figure 2a).The model implements weakening and strengthening behavior by decreasing and increasing the μʹ values of fault segments between the preseismic and postseismic model steps, respectively.Because μʹ is adjusted manually, the model does not involve a rate-and-state friction law and the results are independent of slip rate.As such, we do not model the coseismic stress evolution as in seismic cycle models (Sobolev & Muldashev, 2017;van Zelst et al., 2019).
To determine the forearc stress change as function of the preseismic and postseismic stress states, we solve the model for different pairs of preseismic and postseismic μʹ values.We then determine the megathrust stress drop and resultant forearc stress change as difference in megathrust shear stress and forearc stress between the postseismic and preseismic model steps, respectively.Finally, we determine whether the stress change brings the forearc closer to or further from Coulomb failure in each model run.To do so, we first calculate the critical friction coefficient μ c , for which faults in the forearc were critically stressed at the given stresses.Note that μ c is calculated from the preseismic and postseismic model solutions using the Coulomb criterium for a cohesionless fault and is not a model input parameter (Text S2 in Supporting Information S1).If μ c increases from the preseismic to the postseismic model step, the postseismic stresses support failure on stronger faults, which is compatible with an increase in seismicity after a megathrust earthquake.Conversely, If μ c decreases, the stress change inhibits failure and seismicity and is incompatible with a seismicity increase.Accordingly, we can define the failure tendency Λ = (μ c,post − μ c,pre )/μ c,pre to describe whether the stress change supports (Λ > 0) or inhibits failure (Λ < 0). Figure 2b illustrates the stress state in the forearc as function of μʹ for an exemplary point in the outer forearc.The calculations were carried out for a generic forearc setup and varying μʹ uniformly along the megathrust (Text S1 in Supporting Information S1).If μʹ approaches zero, gravitational stresses dominate and the forearc experiences deviatoric tension (plunge of maximum compressive stress σ 1 > 45°).Conversely, as μʹ increases, megathrust shear stresses dominate and the forearc experiences deviatoric compression (plunge of σ 1 < 45°).Note that for common μʹ values of ≤0.06 (Gao & Wang, 2014) the deviatoric stress in the forearc and plunge of σ 1 do not vary linearly with μʹ due to the opposing effect of gravitational and megathrust shear stresses on the stress field (Figure 2b).Thus, the stress change in the forearc and final stress state resulting from a change in τ depend on the initial megathrust shear stress (Figure S2 in Supporting Information S1).For example, a decrease in τ of 10 MPa from 35 to 25 MPa decreases the deviatoric stress by ∼11 MPa while the stress state remains compressive.For comparison, a decrease in τ from 15 to 5 MPa increases the deviatoric by ∼5 MPa but the stress state switches from deviatoric compression to tension (Figure 2b).
Figures 2c-2j illustrates for the generic forearc model that the failure tendency varies substantially with the preseismic stress state, even if the postseismic stress state is identical (note that μʹ differs for the strengthening and weakening segments in these model runs).If the megathrust shear stress and forearc stress are low before the earthquake (Figure 2c), Λ increases in most of the forearc (Figure 2i).The increase in Λ results from two effects: first, the stress drop results in a reversal in the stress state from deviatoric compression to deviatoric tension, which supports failure because normal faulting operates at lower stresses than reverse faulting (Sibson, 1998).Second, the decrease in horizontal compression results in a net stress increase due to gravitational effects in areas of steep topography, such as the outer forearc (Figure 2g).For comparison, if the megathrust shear stress and forearc stresses are higher before the earthquake (Figure 2d), then the stress drop results in a net stress decrease, which tends to stabilize the forearc and causes low values of Λ, despite the concomitant stress reversal (Figure 2j).The stabilizing effect of the stress decrease reflects that the forearc must sustain much higher stresses before than after the earthquake, which makes postseismic failure unlikely.The effect of the stress decrease is also not counterbalanced by an increase in τ along the strengthening megathrust segment, which has a small impact on the stress in the inner forearc only.Consequently, only a narrow range of preseismic and postseismic stress states leads to an increase in Λ across the forearc.

Model Application to Japan and Discussion of Results
We apply our approach to Japan using finite-element models that account for margin topography, crustal thickness, water load, slab morphology, and extent of the seismogenic megathrust along the Sendai and Iwaki transects shown in Figure 1.We adjust the preseismic and postseismic μʹ values of the weakening and strengthening megathrust segments until the following conditions are fulfilled: first, the bulk of aftershocks and seismic moment release occurs in areas of increased failure tendency.Second, the postseismic stress state is consistent with the prevailing fault kinematics in the forearc after the Tohoku earthquake (Figure 1d).These conditions are fulfilled for the preseismic and postseismic stress states shown in Figures 3a and 3b (Figure S3 in Supporting Information S1).The modeled stress change (Figure 3c) causes a failure tendency increase over large areas that encompass ∼98% of the aftershocks and seismic moment release, while all normal faulting occurs in areas under deviatoric tension (Figures 3d and 3e).For comparison, slightly different preseismic and postseismic stress states are incompatible with the observed aftershock distribution and fault kinematics (Figures S4 and S5 in Supporting Information S1).
The modeled megathrust stress drop is spatially heterogeneous and differs for the Sendai and Iwaki transects (Figure 3c).The largest decrease in megathrust shear stress occurs at shallow and intermediate depth along the Sendai and Iwaki transects, respectively, resembling spatial differences in fault slip (Figure 1b).The stress drop along the weakening and strengthening fault segments varies between 18 and −12 MPa, respectively.The net stress drop averaged over the Sendai and Iwaki transects is ∼2 MPa.The modeled stress drop is comparable to independent stress-drop estimates from 40 different slip-distribution models that vary between 30 and −20 MPa and yield a rupture-zone average <5 MPa (Brown et al., 2015;Wang et al., 2019).The megathrust stress drop in our models relates to changes in μʹ along the weakening and strengthening segments of <0.02 (Figures 3a  and 3b).The preseismic μʹ values of the seismogenic megathrust vary between 0.02 and 0.025, in good agreement with previous estimates of ∼0.02-0.03derived from heat dissipation and force-balance models (Gao & Wang, 2014;Lamb, 2006;Seno, 2009;Wang et al., 2019).
The stress drop causes changes in deviatoric stress (σ dev ) from −10 to 20 MPa in the forearc (Figure 3c), while absolute values of σ dev vary between 5 and 60 MPa both for the preseismic and postseismic model steps (Figures 3a  and 3b).Megathrust stress drop, σ dev , and change in σ dev have the same order of magnitude, in agreement with previous estimates derived from stress reversals in the marine forearc (Hardebeck, 2012;Hasegawa et al., 2011).Furthermore, we find that the stress drop does not result in a general forearc stress decrease, but also in local stress increases.Interestingly, ∼78% of the aftershock seismicity and ∼92% of the seismic moment release occur in areas in which σ dev increases (Tables S2 and S3 in Supporting Information S1).The areas lie within areas of increased failure tendency and encompass the marine forearc and coastal region near Iwaki (Figures 3b-3d).The megathrust stress drop and gravitational stresses arising from steep forearc topography control the σ dev increase, particularly in the marine forearc.The topographic effect near Iwaki, however, is due to a local topographic high along the coast, the Abukuma plateau-a feature that is missing near Sendai (Figure S1b in Supporting Information S1).Our models, therefore, indicate that differences in forearc topography together with spatial stress drop differences have a discernible impact on forearc stress changes and likely facilitated the differences in coastal aftershock seismicity along strike of the forearc (Figure 1b).We also find a dearth of aftershock activity in forearc areas that experience a decrease in σ dev , including inland Japan and the shelf region along the Sendai transect.Only the outermost marine forearc offshore Sendai deviates from this trend and shows little aftershock seismicity despite σ dev increases.
Spatial variations in fault strength may govern spatial variability in aftershock seismicity (Wang et al., 2019).We therefore evaluate how the critical friction coefficient μ c,post varies across the forearc (Figure 3f).Parameter μ c,post is calculated from the postseismic stress solutions (Section 2) and is a composite parameter that includes the effects of both intrinsic friction and pore fluid pressure in the fault zone.It therefore provides an estimate for the effective strength of faults that were seismically active after the Tohoku earthquake.The effective fault strength in forearc areas showing no aftershock seismicity is likely higher than μ c,post as they seem not to have failed at the given stresses.
The μ c,post in the outermost forearc offshore Sendai is ∼0.075-0.24,i.e., highest in the forearc.The comparatively high strength may indicate that the outer forearc is generally stronger than more internal parts.Alternatively, the apparent strength may represent a postseismic condition and result from the large coseismic seaward displacement of the outer forearc, which exceeded 20 m above the main slip area close to the trench (Kido et al., 2011;Sato et al., 2011).The strong seaward displacement causes forearc dilation that tends to decrease pore pressures and increase fault strength (Manga et al., 2012).This interpretation is consistent with the higher seismicity in the outer forearc offshore Iwaki, where the seaward displacement was much smaller (∼5 m) and should not lead to a similar dilational strengthening effect.
Most aftershocks in both transects occur between ∼120 and ∼200 km from the trench and locate in the mantle wedge and overlying crust.The associated faults are likely very weak, in accordance with μ c,post values of ∼0.01-0.06(Figure 3f).The low fault strength may relate to mantle wedge serpentinization and pressure buildup from fluids liberated from the subducting slab (Hyndman & Peacock, 2003;Tulley et al., 2022).Thus, the intense seismicity in the area may be linked to the position above the dehydrating slab.
Further inland, μ c,post is <0.02 for all but the upper 5-10 km of the forearc, and suggests that inland Japan hosts weak faults and quasi-lithostatic pore fluid pressures.In contrast to the mantle-wedge setting, high pore fluid pressures are likely locally restricted because of the large distance to the dehydrating slab.The local restriction of high fluid pressures is consistent with previous studies showing that the inland seismicity is linked to hydrothermal fluid migration along mature fault systems (Okada et al., 2011;Yoshida et al., 2017).The inland seismicity was also interpreted to have been triggered by local fluid pressure increases (Terakawa et al., 2013).A pressure increase would raise the failure tendency and may explain the seismicity underneath the volcanic arc not captured in our model (Figure 3e).Nevertheless, the low inland seismicity may be conditioned by very low background stresses that are marginally capable of initiating frictional failure, suggesting that these regions deform mainly viscously.Similar conditions will also apply to mantle areas that show very low or no seismic activity, although the failure tendency partially increases (Figures 3d and 3e).

Comparison of Model Stresses With Static Stress Drop of Forearc Earthquakes
Our models constrain the preseismic and postseismic stresses in the Japanese forearc.We now validate the stress solutions independently by comparison with static stress-drop estimates of earthquakes that occurred in the forearc before and after the Tohoku earthquake (Figure 4a).Earthquake stress drop magnitudes are restricted by the stress along the fault that ruptured and provide a lower bound for deviatoric stress derived from the finite-element models.In addition, spatiotemporal stress drop heterogeneities may hint at variations in ambient stresses (Kemna et al., 2021).
We estimate stress drop values using seismic data from the NIED High Sensitivity Seismograph Network (National Research Institute for Earth Science and Disaster Resilience, 2019) and the Japanese Meteorological Agency (JMA).As individual stress-drop estimates may reflect local stress heterogeneities that are not representative of average stresses, we first determine individual earthquake stress drop values using both single-spectra and spectral-ratio fitting methods.We then calculate average stress drop values (  ∆ Forearc ) at 10-km depth intervals 10.1029/2022GL101320 9 of 12 with 50% overlap (Text S3 and S4 and Figures S8 and S9 in Supporting Information S1).The averaged values better reflect bulk stress conditions in the forearc and are better suited for comparison with modeled stresses (Moyer et al., 2020).
The pre-Tohoku and post-Tohoku  ∆ Forearc values are largely similar and increase gradually with depth from ∼2 to 8 MPa (Figure 4b).The stress drop estimation incorporates a depth-dependent velocity model from the JMA, thus the observed increase with depth is not an artifact of a constant shear velocity assumption (Abercrombie et al., 2021;Hauksson, 2015).A selection of similar event pairs within 5 km hypocentral distance further ensures a depth-dependent attenuation correction for stress drop values obtained using the spectral-ratio fitting method (Abercrombie et al., 2021;Figure S10 in Supporting Information S1).Aside from the depth dependency, the  ∆ Forearc values are spatially homogenous and do not vary significantly along forearc strike (Figure 4a).Likewise, stress drop values vary negligibly in time without a discernible trend (Figure S11d in Supporting Information S1).The range in estimated stress drop magnitudes presented here is consistent with previous estimates for the Japanese forearc (Yoshida et al., 2017).
We compare stress-drop estimates with model results by calculating the depth-dependent maximum values of σ dev from the models for the Sendai and Iwaki transects (Figure 4b) and for the subregions that include most of the aftershocks (Figures 4c-4e).The depth-dependent deviatoric stresses are similar for the preseismic and postseismic states, which reflects that the maximum stresses did not change significantly due to the modeled stress drop.Similarity in the preseismic and postseismic σ dev values resemble the constancy of  ∆ Forearc values.Moreover, the deviatoric stresses are generally higher than the stress-drop estimates, consistent with partial stress release during earthquakes (Brune, 1970;Simpson, 2018), and exhibit a similar depth dependency to  ∆ Forearc .The remarkably consistent correlation across the forearc follows the trend of  ∆ Forearc being ∼10% of σ dev (Figures 4c-4e), suggesting that the average stress release scales with ambient stress.Taken together, we find that the deviatoric stresses obtained from the finite-element models agree favorably with the  ∆ Forearc values obtained from seismic data and validates the significance of the modeled forearc stresses.

Summary and Conclusions
We examine aftershock triggering by the stress drop of a megathrust earthquake and show that triggering is only possible for a narrow range of stress drop values and preseismic forearc stresses.We determine the range of values for the aftershock seismicity of the Tohoku earthquake using numerical force-balance models that are constrained by the observed aftershock distribution and fault kinematics.The model results are consistent with major aspects of the earthquake, including the magnitude and spatial heterogeneity of the stress drop, the strength of the Japanese megathrust, and the partial stress reversal in the forearc.Our models further indicate that ∼78% of the aftershock seismicity and ∼92% of the seismic moment release occurred in forearc areas of net stress increase, and that the detailed aftershock distribution was also governed by spatial variability in fault strength.Finally, we find that the modeled forearc stresses are remarkably consistent with independent estimates of earthquake stress drop values in the forearc before and after the Tohoku earthquake.Given the consistency in the data and its implications, we conclude that the bulk forearc seismicity in the first months after the Tohoku earthquake was triggered by the megathrust stress drop and resulting stress change in the forearc.This stress change will not remain unchanged after the megathrust earthquake but will be altered by viscoelastic relaxation processes, and relocking and reloading of the megathrust in the years after the megathrust earthquake (Bedford et al., 2016;Bürgmann et al., 2016;Sobolev & Muldashev, 2017;Sun et al., 2014).Viscoelastic relaxation and reloading tend to diminish the failure tendency and forearc seismicity with time.Whether the observed power-law like decrease in forearc seismicity already reflects these or other processes remains to be resolved.Independently, our approach helps to determine the stress perturbation in the forearc caused by the megathrust stress drop and can be applied to other great megathrust earthquakes to obtain detailed information on forearc stress and strength.events used in this study, and ABAQUS output databases are available from https://doi.org/10.25835/0072357.We used the mtspec Python wrapper for spectral estimations (Prieto et al., 2009).We process seismic data and create maps, swath profiles, and diagrams using the Python packages Obspy (Krischer et al., 2015), Matplotlib (Hunter, 2007), and GMT (Wessel et al., 2019).Color schemes follow Scientific colour maps (Crameri et al., 2020).The finite-element models were calculated, processed, and plotted using the commercial software packages ABAQUS and MATLAB and the Matlab tool Abaqus2Matlab by Papazafeiropoulos et al. (2017).

Introduction
Text S1 to S4 provide information on the methods used in this study.Figures S1 to S7 provide information on the finite-element models and model results not shown in the main manuscript.Figures S8 to S11 provide information on the stress drop estimation and stress drop data not shown in the main manuscript.Table S1 provides details on the frictional contact in the finite-element models for Japan.Tables S2 and S3 compare the aftershock distribution and seismic moment release, respectively, with results of the finite-element models.

Text S1. Finite-Element Model
To calculate the stress in the forearc, we use a plane-strain finite-element model of an elastic upper plate and a rigid lower plate in frictional contact following earlier studies (Dielforder & Hampel, 2021;Wang & He, 1999;Wang et al., 2019).The models are created and calculated using the commercial software ABAQUS (version 6.14).The basic model setup, boundary conditions, and material parameters are shown in Figure S1.All models are meshed using linear triangular elements with an average element edge length of ~2 km.
ABAQUS computes the total stress in the elastic upper plate, which results from all applied boundary conditions, including gravity, isostasy, and friction along the plate interface (Figure S1a).The gravitational force is calculated for a gravitational acceleration g = 9.81 m s -2 and densities of 2,800 and 3,300 kg m -3 for crustal and mantle parts, respectively.Seawater load is modeled using a pressure boundary condition, with values determined from seawater density (1,025 kg m -3 ) and depth.Friction is calculated for standard Coulomb friction by displacing the lower plate, such that the shear stress along the contact t depends on friction coefficient µ' and normal stress sn (t = µ'sn).Using the friction law allows us to describe the shear stress along the megathrust in a consistent manner.Still, it should be noted that after the earthquake the actual shear stress can be somewhat lower than the one predicted by the friction law.
Basal drag at the base of the lithosphere of the upper plate is not considered in our model, because basal drag is thought to have a small impact on forearc stresses only (e.g., Erdös et al., 2021).As the upper plate is elastic, it has no yield criterion and deforms only elastically, but not permanently due to plastic or viscous deformation.The model calculates no pore fluid pressures and no effective stresses (total stress minus pore fluid pressure).
At the beginning of each model run, isostatic equilibrium is established following the procedure described in Hampel et al. (2019).Afterwards, the model calculates the total stress in the forearc for a preseismic and a postseismic model step.The model steps differ by the µ' values assigned to the frictional contact.The preseismic and postseismic µ' values and resultant megathrust shear stresses are representative for the strength of the megathrust just before and after the megathrust earthquake (cf.Wang et al., 2019).Strengthening and reloading of the megathrust in the postseismic and interseismic periods is not modelled.
The frictional contact comprises three sections that represent different parts of the plate interface as illustrated in Figure S1 and outlined in the following.The upper section of the interface represents the seismogenic megathrust and comprises weakening and strengthening segments to allow for a coseismic decrease and increase in shear stress, respectively.The extent of the respective megathrust segments and their preseismic and postseismic µ' values are variable.
The middle section of the interface represents the upper part of the viscous plate interface, directly downdip of the seismogenic megathrust, where the plate-interface rheology transitions into viscous creep in nature.The stress along the viscous plate interface depends on temperature and strain rate.In the late interseismic period, the strain rate diminishes and the shear stress is low (cf.Wang et al., 2019).During the earthquake, the strain rate and the shear stress increase (e.g., Brown et al., 2015).To account for the low preseismic shear stress and the coseismic stress increase, we assign preseismic and postseismic µ' values of 0.001 and 0.005 to the middle section of the interface.
The lowermost section of the interface represents the viscous plate interface at depth >75 km.In nature, the temperature increases significantly along the viscous plate interface at about 70-80 km depth due to the onset of mantle wedge flow and associated heat advection (e.g., Wada & Wang, 2009).The interface shear stress is assumed to be negligible at depth greater than 70-80 km (e.g., Lamb, 2006).We therefore model the lowermost section of the interface frictionless (µ' = 0), such that the shear stress is zero.
We construct different finite-element models for the generic forearc (Figures 2) and for the Sendai and Iwaki transects across the Japanese forearc (Figures 3).The generic forearc models include a curved plate interface with an average dip angle of 18º, a 5 km deep trench, a trench-coast distance of 125 km, a maximum elevation of 1 km in the volcanic arc area at 250 to 270 km distance from the trench, and a 30-km-thick crust that thickens underneath the area of high elevation to isostatically support the topography.The forearc setup is representative for the average of global ocean-continent subduction zones (see Table 1 in Dielforder & Hampel, 2021).Note that in Figure S2e-h, we present alternative generic models for a flat forearc wedge to illustrate the gravitational effect of margin topography.
The generic forearc models are used to illustrate how 1) the deviatoric stress in the forearc and plunge of s1 vary with µ' and t (stress paths in Figures 2b and S2b and f), and 2) the failure tendency L varies for different pairs of pre-and postseismic µ' values (Figure 2c-j).
For the calculation of the stress paths, the upper interface section encompasses the upper 50 km and the model is solved for µ' values between 0 and 0.12 that apply to the entire upper interface section.Strengthening along the middle interface section is not considered and the µ' value of the middle section is always 0.001.For the calculation of the failure tendency L, the upper interface section encompasses the upper 40 km and the µ' values of the upper and middle interface section are varied to model weakening and strengthening behavior as indicated in Figure 2c-f.The models for the Sendai and Iwaki transects account for the site-specific forearc topography, seawater load, interface geometry, and crustal thickness.The forearc topography is approximated by the mean elevation (Figure S1b), that we calculated along 100-km-wide swath profiles from the ETOPO1 global relief model using TopoToolbox for MATLAB (Amante & Eakins, 2009;Schwanghart & Scherler, 2014).The swath profiles run parallel to those shown in Figure 1.The interface geometry has been adopted from the Slab2 model (Hayes et al., 2018) by fitting an arc with constant curvature through the upper 100 km of the slab model (e.g., Cailleau & Oncken, 2008;Wang et al., 2019).The crust is 25 km thick at the Moho-megathrust intersection and thickens to 30 km inland Japan (Iwasaki et al., 2013).The upper interface section extends down to ~58 km depth, in accordance with a downdip limit of the seismogenic megathrust in Japan at about 55-60 km depth (e.g., Hayes et al., 2018;Gao & Wang, 2014).
The models for Japan allow weakening and strengthening along the upper interface section in accordance with stress-drop distribution models (e.g., Brown et al., 2015).The extent of the weakening and strengthening fault segments and the associated µ' values are listed in Table S1.Both the extent of the fault segments and the µ' values have been obtained during the fitting process, in which we determined the forearc stress change that best explains the observed aftershock distribution (as illustrated in Figure 3).Using a finer segmentation (for example, one segment every 5 km), does not improve the model fit.Using less segments (for example, only one or two) diminishes the fit.

Text S2. Calculation of failure tendency
The finite-element models compute the total stress in the forearc for a preseismic and postseismic model step that represent the time before and after the megathrust earthquake, respectively.The coseismic stress change is the difference in total stress between the two model steps.To evaluate whether the stress change supports or inhibits frictional failure in the forearc, we first determine from the model results the critical friction coefficient µc.Parameter µc is indicative for the maximum effective strength of cohesionless faults that can slip at a given stress (Figure S6a).Hence, µc can be expressed analogous to the effective friction coefficient, that is, where µ is the coefficient of friction and l is the pore fluid pressure ratio, that is, the ratio of pore fluid pressure Pf to lithostatic stress sz.While the lithostatic stress is given by the results of the finite-element models, the pore fluid pressure is not given but obtained by rearranging the Coulomb criterion for a cohesionless fault t = µ(sn -Pf), that is, where t is the shear stress and sn is the normal stress on the fault (Sibson, 1998).The shear stress is calculated from the model solutions for planes optimally oriented for failure with respect to the maximum compressive stress s1 as where s3 is the minimum compressive stress and q is the angle between s1 and the failure planes (Sibson, 1998).The principal stresses s1 and s3 are given by the results of the finite element models.Angle q is  = 0.5tan (& (1/).
(4) The normal stress is calculated as (5) Combining equations 1 and 2 gives Based on equations 1-6, the critical friction coefficient is calculated for the stress solutions of the preseismic (µc,pre) and postseismic (µc,post) model steps.Comparing the µc values for the two model steps indicates whether the stress change supports or inhibits failure.If µc,post > µc,pre, then the postseismic stress can cause failure on stronger faults than the preseismic stress, and the stress change supports failure in the forearc (Figure S6b).Conversely, if µc,post < µc,pre, then the preseismic stress was more favorable for failure and the stress change inhibits failure (Figure S6d).Accordingly, the failure tendency can be defined as Based on equations 1-7, we calculate L for each point in the forearc.We use µ = 0.7 for all calculations.Using lower and higher values for µ slightly decreases and increases the value L, respectively, but has no impact on our findings and conclusions (Figure S7).

Text S3. Earthquake catalogue
We obtain hypocentral solutions and P and S phase arrivals from the Japan Meteorological Agency (JMA) (https://www.data.jma.go.jp/svd/eqev/data/bulletin/index_e.html).We consider events occurring at depths shallower than 1 km above the top of the slab (Hayes et al., 2018) to limit the amount of interplate and intraslab events in the final dataset.We further restrict the dataset of forearc events by using only high-quality hypocentral solutions, as described in the following.We consider earthquakes that occurred between the 01/01/2008 and the 12/31/2011, with JMA magnitude of MJMA ≥ 2.5, with vertical errors ≤5 km, and latitude and longitude errors ≤0.025°.The filtered catalogue consists of 21,665 events (Figure 1a) with an average latitude error of 0.005±0.003°(~0.5±0.3 km), an average longitude error of 0.01±0.006°(~1±0.6 km), and an average vertical error of 1.5±0.9km.The filtered catalogue has a magnitude of completeness (Mc) of 3.4 ± 0.3.

Text S4. Static stress drop estimates
We use the catalogue of forearc events described in Text S3 for static stress drop calculations (e.g., Abercrombie, 2021;Kanamori & Anderson, 1975) from single spectra and spectral ratio fitting methods (Figure S8a).We estimate the displacement spectral amplitude using Thomson's multitaper method (Thomson, 1982) from S-wave windows starting 0.2 seconds before the arrival and containing 90%, 80%, and 70% of the energy at a specific station within a hypocentral distance of 25 km, 25-50 km, and ≥50 km, respectively (Pacor et al., 2016) (Figure S8b).We ensure that each spectrum has a signal-to-noise ratio (SNR) ≥ 3 in a magnitude-dependent frequency band (Figure S8c-d).We define the frequency band for a particular event by calculating corner frequencies corresponding to theoretical stress drop values of 0.01 and 1000 MPa, which we use as lower and upper-frequency limits, respectively.We then fit the spectra (single spectra using a trust-region-reflective minimization algorithm following Newville et al. ( 2016) where W0 is the long-period spectral amplitude, f is the spectral frequency, t is the travel time, Q is the quality factor, fc is the corner frequency, and n is the high-frequency falloff rate (Boatwright, 1980;Brune, 1970).We use a spectral shape constant g of 2 (Boatwright-model), which fits most of the spectra.We then calculate the seismic moment (M0) using the fitted W0 values as where r is the density, c is the S-wave velocity at the depth of the hypocenter, R is the station-event hypocentral distance, and Ufq is the mean radiation pattern for S-waves (Madariaga, 1976).We use S-wave velocities at the hypocentral depth taken from the JMA velocity model (https://www.data.jma.go.jp/svd/eqev/data/bulletin/catalog/appendix/trtime/trt_e.html).We then calculate mean values and 95% confidence intervals for a particular event using a deleteone jackknife-mean (Prieto et al., 2006) when at least five M0 and fc S-wave estimates are available.
Non-source related terms such as site and path effects can bias the estimates of fc and final stress drop.To ensure accurate estimates of fc, we also use a spectral ratio approach.The ratio between two co-located event spectra at a specific station cancels possible site and path effects and allows high-quality fc estimates from one or both events in the pair depending on the frequency range of high SNR (Bakun & Bufe, 1975) (Figure S8a-c).To obtain event pairs, we cross-correlate 3-component full waveforms of events within 5 km hypocentral distance and retain event pairs with a cross-correlation coefficient ≥0.7 and a magnitude difference ≥0.5.An 0.5 magnitude difference ensures the selection of event pairs with fc values with a resolvable difference in the spectral ratio fitting.The displacement spectral ratio Wr(f) between the two event spectra can be written as Where fc1 and fc2 are the corner frequencies of the larger magnitude target and the smaller magnitude empirical Green's function events (eGf), respectively.The spectral shape constant g is set to 2. We require at least five S-wave station ratios for individual event pairs, manually review the spectral ratio fits to ensure high quality, and check whether fc1 and fc2, or only fc1 values are resolvable (Figure S8e).
We use M0 values from single spectra and fc values either from single spectra or spectral ratio analysis to calculate stress drop values (Ds) assuming a circular crack model with radius r = kb/fc as (Eshelby, 1957) where b is the shear wave velocity at the hypocentral depth and k is a constant set to 0.26 assuming a symmetrical circular model with a rupture velocity of 0.8 b (Kaneko & Shearer, 2015).
We use event waveforms from the High-Sensitivity Seismograph Network of Japan (Hinet) (Obara et al., 2005) available from the National Research Institute for Earth Science and Disaster Resilience (NIED, 2019).We downloaded waveforms from a total of 259 Hi-net stations (Figure S9a).In case of events with MJMA ≥ 3.5, we deconvolve short-period Hi-net velocity records by their seismometer response to simulate seismic waves recorded with broadband seismometers (Maeda et al., 2011) and estimate earthquake source parameters.
We calculate moment magnitudes (Mw) (Hanks & Kanamori, 1979) using the M0 values from S-wave single spectra and observe a good fitting with respect to the magnitude values reported in the initial catalogues (Figure S9d).We do not consider events with Mw ≤3 as we observe a systematic trend of lower stress drop values for these events which are likely related to observational constraints.The trend of lower stress drops for lower magnitudes (fc ~> 10-15 Hz) has commonly been related to decreasing SNR values with magnitude and frequency and to bandwidth limitations of the 100 Hz seismometers (Ide et al., 2003).We compare the corner frequency values of events for which both single spectra and spectral ratio estimates are available and find good agreement between them (Figure S9b).The good agreement between single and spectral ratio fitting estimates, together with the good spatial sampling of spectral ratio estimates (Figure S9a), ensures a high quality of the single spectra estimates.In case of events for which we obtain both spectral estimates, we associate the stress drop value obtained with the spectral ratio fitting to the event.
We obtain 5,819 stress drop estimates, including 5,458 from single spectra fitting and 361 from spectral ratio fitting (Figure S9c).The average values of Q and n are 1,121 and 2.5, respectively.Here g is gravitational acceleration, r is density, n is Poisson ratio, E is Young's modulus, and µ' is the effective friction coefficient.The model is isostatically supported by an elastic foundation at the base of the upper plate.The left model side may move vertically and is set 500 km from the trench, that is, farther then shown in Figures 2 and 3, to minimize boundary effects.Following previous studies (Dielforder & Hampel, 2021;Wang & He, 1999;Wang et al., 2019), we adopted an almost incompressible material (n = 0.48) for the upper plate.Note that using a lower Poisson ratio (e.g., n = 0.3) has a negligible effect on the stresses only (Dielforder & Hampel, 2021).The lower plate is modeled rigid, because its role is simply generating shear stress along the frictional contact.(b) Mean elevation along 100-km-swath profiles running parallel to the Sendai and Iwaki transects shown in Figure 1.The mean elevation is used to approximate the forearc topography in the finite-element models for Japan.The gravitational effect causes that sdev and the plunge of s1 do not vary linearly with µ' and t as illustrated by the stress paths in (b).Subfigures (e-h) are analogue to (a-d) but show solutions for a flat forearc wedge without margin topography.The absence of topography removes the gravitational effect and the forearc experiences no deviatoric tension.In that case, the pattern of sdev change does not depend on initial stress state (compare g and h) as sdev varies linearly with µ' and t as illustrated by the stress path in (f).S2).The percentages are calculated for aftershocks with a magnitude larger or equal to the magnitude of completeness (Mc = 3.4).Seismic moment release M0 = 10 (1.5Mw+9.1)(Hanks & Kanamori, 1979), where Mw is moment magnitude.Mw is obtained according to the empirical relationships illustrated in Figure S8d.

Figure 1 .
Figure 1.Seismotectonic setting of NE Japan.(a-d) Forearc seismicity before and after the Tohoku earthquake.(a, b) Map view.Black lines indicate location of cross-sections shown in (c) and (d).Dashed rectangles indicate the width of swaths (200 km) projected into the cross-sections.Se = Sendai, Iw = Iwaki.(c, d) Cross-sectional view.Gray dots are aftershock hypocenters from the Japan Meteorological Agency (JMA).Beach balls denote JMA focal mechanism solutions and are shown in profile view.Reverse faulting events in the outer marine forearc above the megathrust are likely poorly located interplate events(Nakamura et al., 2016) or the slab model is not accurate enough to classify them as interplate events.(e) Number of forearc earthquakes.Count includes events with magnitude ≥ magnitude of completeness (M c = 3.4,Text S3 in Supporting Information S1).

Figure 2 .
Figure 2. Model setup and generic models results.(a) Schematic representation of forearc model.Here, ρ is density, g is gravitational acceleration, μʹ is the megathrust effective friction coefficient, τ and σ n are the shear and normal stresses, respectively, and σ 1 is maximum compressive stress.(b) Deviatoric stress σ dev (black) and plunge of σ 1 (orange) as function of μʹ and average megathrust shear stress.Solutions are for site P1 in (a).DT and DC denote deviatoric tension and deviatoric compression, respectively.(c-j) Model results.(c-f) Deviatoric stress and plunge of σ 1 for preseismic and postseismic model steps.(g, h) Change in σ dev due to megathrust stress drop Δσ.Positive and negative values of Δσ indicate decrease and increase in megathrust shear stress, respectively.(i, j) Failure tendency Λ.

Figure 3 .
Figure 3. Model results for Japan.(a, b) Deviatoric stress and plunge of σ 1 for preseismic and postseismic model steps.(c) Change in σ dev in the forearc due to megathrust stress drop Δσ.Positive and negative values of Δσ indicate decrease and increase in megathrust shear stress, respectively.(d, e) Failure tendency Λ. Gray dots and red, blue, and black solid circles in (e) are hypocenters of forearc seismicity and reverse, normal, and strike-slip faulting events after the Tohoku earthquake, respectively (see Figure 1d).(f) Critical friction coefficient μ c,post , which estimates the effective strength of faults in areas showing aftershock seismicity.Hatched areas indicate areas showing no seismicity, for which fault strength cannot be estimated.

Figure 4 .
Figure 4. Earthquake stress drop values and comparison to model stresses (solid and dashed lines in c-e).(a, b) Earthquake hypocenters with estimated static stress drop.(c-e) Gray dots are individual stress-drop estimates (n = total number of estimates), squares are average stress drop values (  ∆ Forearc ).Error bars represent one standard deviation.Basal events in (e) are events from 1 to 5 km above the slab interface.

Figure S1 .
Figure S1.Model setup and forearc topography.(a) Sketch illustrating the basic setup of the finite-element model with an elastic upper plate and a rigid lower plate in frictional contact.

Figure S2 .
Figure S2.Dependency of forearc stress change on initial stress state.Subfigures (a) and (b) are similar to Figure 2a, b in the main text.Here r is density, g is gravitational acceleration, µ' is the megathrust effective friction coefficient, t and sn are the shear and normal stresses, respectively, and s1 is maximum compressive stress.DT and DC denote deviatoric tension and deviatoric compression, respectively.(c, d) Change in deviatoric stress (sdev) in the forearc due to a decrease in average megathrust shear stress of 10 MPa from 15 to 5 MPa (c) and from 35 to 25 MPa (d).The pattern of sdev change is different for the two cases (compare c and d), because of the deviatoric tension resulting from gravity in the presence of margin topography.

Figure S3 . 10 Figure S4 .
Figure S3.Model results for the Sendai and Iwaki transects, Japan.(a, c) Maximum compressive stress s1 (thick black bars) and minimum compressive stress s3 (thin red bars) for the preseismic (a) and postseismic (c) model steps.Length of bars scaled.Note that the principal stresses are plotted every 5 km, while in Figure 3 s1 is plotted only every 10 km for reasons of illustration.(b, d) Deviatoric stress s dev = 0.5(s 1 -s 3) for preseismic (a) and postseismic (b) model steps.Note that the deviatoric stresses are also shown in Figure 3a and b.

Figure S5 .
Figure S5.Set of model results illustrating the effect of the stress drop magnitude (postseismic µ' value) on the stress state (a, c) and failure tendency L (b, d).The model results are obtained from finite-element models that are identical to the models discussed in the main text but use a different µ' value for the postseismic model step.(a, b) Results for a low stress drop, which is modelled by decreasing µ' along the weakening segments of the megathrust by 0.005 (see TableS1for the segmentation of the megathrust).(c, d) Results for a total stress drop, which is modelled by decreasing µ' along the weakening segments of the megathrust to zero.Both for a low (b) and total stress drop (d), most aftershocks of the Tohoku earthquake occur within areas of increased failure tendency (number in brackets in (c, d), see Figures3 and S4for comparison).However, the modelled stress states do not fit the observed fault kinematics in the forearc after the Tohoku earthquake.For a low stress drop (a), some forearc areas are under deviatoric compression in the model, but experienced normal faulting in nature.For a total stress drop (c), some forearc areas are under deviatoric tension in the model, but experienced reverse and strike slip faulting in nature.

Figure S6 .
Figure S6.(a) Mohr diagram illustrating the derivation of the critical friction coefficient µc.See Text S2 for details.(b-d) Mohr diagrams illustrating changes in failure tendency L due to the stress change between the preseismic (blue) and postseismic (orange) model steps.Shown solutions are exemplary.The dashed failure envelope is shown for guidance only, but does not imply critical or unstable stress states.(b) Increase in failure tendency (L>0).(c) No change in failure tendency (L=0).(c) Decrease in failure tendency (L<0).

Figure S7 .
Figure S7.(a, b) Effect of friction coefficient µ (see equations 1,2, and 6 in Text S2) on failure tendency L. (a) Results for µ = 0.5.(b) Results for µ = 0.9.Note that L decreases and increases for lower and higher values of µ, respectively.See Figure 3d, e for comparison.

Figure S8 .
Figure S8.Single spectra and spectral ratio fitting methods.(a) Schematic representation of single spectra and spectral ratio fitting methods.Gray circles are earthquake hypocenters (scaled by magnitude), white triangles are seismic stations, and dashed lines are source-station ray paths.The blue-and yellow-colored dashed circles denote the respective crustal source volumes considered in a single iteration of single-spectrum and spectral-ratio fitting, respectively.(b) Waveform recordings of a spectral ratio event pair for the Japanese dataset with cross-correlation (CC)=0.90 at station KZMH.E from the Hi-net network.Dashed and solid lines denote the time windows for signal and noise estimation, respectively.eGf, empirical Green's function.(c) S-wave single spectra estimates for the windows shown in (b).(d) Fitted Swave spectrum for the target event spectrum shown in (c).(e) Fitted S-wave spectral ratio from the spectra shown in (c).

Figure S9 .
Figure S9.Single spectra and spectral ratio fitting estimates, and comparison of original and calculated magnitudes.(a) Number of S-wave spectral estimate per seismic station and distribution of single spectra and spectral ratio stress drop estimates.(b) Comparison between corner frequencies (fc) obtained from single spectra and spectral ratio fitting methods.(c) Corner frequency with respect to moment magnitude (Mw).Green circles and orange squares indicate single spectra and spectral ratio estimates, respectively.The number in brackets gives the total number of estimates.The dashed lines denote stress-drop isolines for a shear-wave velocity of 3,400 m s -1 at a depth range between 7 and 20 km.(d) Comparison between original magnitude of the Japanese Meteorological Agency magnitude (MJMA) and calculated moment magnitudes from S-wave single spectra fitting.The dashed black line indicates a 1:1 fitting line, while the orange dashed line indicates the linear regression fitting line between original and calculated magnitudes.Error bars in (b-d) represent the 95% jackknife confidence interval.

Figure S10 .
Figure S10.Comparison of stress drop values obtained from single spectra and spectral ratio fitting methods.Gray dots are individual stress-drop estimates (n = total number of estimates), squares are average stress drop values calculated at 10 km depth intervals with 50% overlap.Error bars represent one standard deviation.Note that the averages values of pre-Tohoku spectral ratio fitting estimates are statistically not representative due to the low number of estimates (n = 39).

Figure S11 .
Figure S11.Static-stress-drop values of earthquakes in the Japanese forearc before (01/01/2008-03/11/2011) and after (03/11/2011-12/31/2011) the Tohoku earthquake.(a) Hypocenters of earthquakes with estimated static stress drop along the Sendai and Iwaki profiles (see Figure 1 for traces of profiles).(b) Histograms with logarithmic scale showing the log-normal distribution of all stress drops estimates of forearc earthquakes before and after the Tohoku earthquake.The ∆ BBBB all values indicate the average and 1 standard deviation of all stress drop estimates.(c) Number of stress-drop values as function of time.(d) Mean stress drop values as function of time.The number of stress drop values in (c) and the mean stress drop in (d) are calculated for 100-day-intervals before the Tohoku earthquake and for 10-dayintervals after the Tohoku earthquake.The different time intervals used for the averaging in (d) consider differences in the total number of events, which is much lower for the time before the Tohoku earthquake.Error bars in (d) indicate the 1 standard deviation.

3 a
Friction coefficient µ' is a model input parameter and assigned to megathrust fault segments for the preseismic (pre) and postseismic (post) model steps.b The average shear stress on the megathrust segments is an output parameter and calculated for the preseismic (pre) and postseismic (post) model steps for the respective µ' values ( a ).c The stress drop is an output parameter calculated by subtracting the average shear stress values ( b ) for the postseismic model step from the shear stress values of the preseismic model steps.The values are illustrated in Figure 3.d-gStress drop values marked with d, e, f, and g have been averaged for the respective megathrust segments in Figure3for reasons of illustration.

a
Events with a magnitude (M) larger or equal to the magnitude of completeness (Mc = 3.4) only.bNumber of events in forearc areas, in which the failure tendency increases (L > 0) or decreases (L < 0) in the finite-element models.c Number of events in forearc areas, in which the deviatoric stress sdev increases (Dsdev > 0) or decreases (Dsdev < 0) in the finite-element models.

Table S1 .
Segmentation and properties of the frictional contact in the finite-element models for Japan, resultant megathrust shear stresses, and stress drop values

Table S2 .
Number of aftershocks