Revisiting interseismic deformation in Nankai: focusing on slip-deficit accumulation in the ETS zone and comparison with Cascadia

Various stress-releasing phenomena, such as episodic tremor and slip (ETS) and low-frequency earthquakes, occur at the downdip seismogenic zone in southwest Japan. However, it is unclear how much net stress and slip deficit accumulate at these depths during the interseismic phase. Here, we perform both elastic and viscoelastic earthquake-cycle forward models and reassess the locking state in Nankai from a synthesized perspective with the aid of geodetic modeling results. Our results suggest that the overestimation of the locking depth due to ignoring Earth’s viscoelasticity is much smaller (less than 5 km) in this early interseismic subduction zone compared to that (~ 10 km) of late-inter-seismic margins. Considering viscoelastic modeling results and other physical arguments, the preferred steady-state viscosities for the continental and oceanic mantle are 5 × 10 19 Pa s and 10 20 Pa s, respectively. We find a clear trade-off between the full locking depth and the width of the transition zone when explaining both horizontal and vertical geodetic data, demanding other data to further resolve this inherent ambiguity. Unlike in Cascadia, partial megathrust locking in Nankai likely penetrates into the ETS zone, leaving no intervening gap between the shallow megathrust, where hosts large earthquakes, and the ETS zone. Assuming locking extends into the downdip of the ETS zone (i.e., 40 km), we propose a preferred viscoelastic locking model with a full locking depth of 18 km and a broad transition zone spanning a 22-km depth range. In this model, the downdip half portion of the transition zone corresponds to the ETS zone, which can accumulate certain slip deficit in a largely creeping and partially locked state. However, most of the accumulated slip deficit in the ETS zone may be accommodated aseismically simultaneously by stress-releasing phenomena, leaving limited to no budget to release during future megathrust earthquakes. We suggest that precise documentation of total slip during slow slip events, along with refinement of viscoelastic locking models, will provide new insights into the net slip budget available in the ETS zone. This will help assess the potential of future coseismic and/or postseismic slip penetrating into the ETS zone in Nankai, Cascadia and other subduction zones.


Graphical Abstract 1 Introduction
Downdip seismogenic limit of a megathrust interface determines the maximum landward extent of coseismic rupture and is a key parameter for understanding subduction-zone dynamics.Traditionally, it has been proposed that downdip seismogenic limit is primarily controlled by either isotherm of ~ 450 ℃ or fore-arc Moho intersection, whichever shallower (e.g., Hyndman et al. 1997;Oleskevich et al. 1999;Hyndman 2013).In some circumstances, additional factors (e.g., stress state, elastic material properties and structures of upper and lower plates) may also influence the seismogenic extent (e.g., Wallace et al. 2009;Lay et al. 2012;Saffer and Wallace 2015).In the early twenty-first century, various stress-releasing events, such as episodic tremor and short-term (weeks) slow slip (ETS), long-term (months or longer) slow slip events (SSE), low-frequency earthquakes (LFE) and very low-frequency earthquakes (VLFE), were discovered in a ~ 10-km depth band (hereafter called the ETS zone for simplicity) along the strike direction at the downdip extension of the seismogenic zone in Nankai and Cascadia (Dragert et al. 2001;Obara 2002;Rogers and Dragert 2003;Shelly et al. 2007).Despite tremendous advances in documenting and modeling the kinematic and mechanical behaviors of the ETS zone since their discovery, particularly in Cascadia and Nankai (Hirose and Obara 2005;Schwartz and Rokosky 2007;Rubinstein et al. 2007;Maeda and Obara 2009;Wech et al. 2009;Holtkamp and Brudzinski 2010;Bartlow et al. 2011;Nishimura et al. 2013), it is still under heated debate about their momentduration scaling law (Ide et al. 2007;Michel et al. 2019;Dal Zilio et al. 2020;Ide and Beroza 2023), underlying faulting mechanisms and physics (Liu and Rice 2007;Audet and Bürgmann 2014;Gao and Wang 2017;Bürgmann 2018;Im et al. 2020;Behr and Bürgmann 2021;Ando et al. 2023;Perfettini and Molinari 2023) and their interplay with megathrust events at shallower depth (Ito et al. 2007;Segall and Bradley 2012;Obara and Kato 2016;Uchida et al. 2016).
Spatially, the ETS zone has been proposed to be isolated at depth to some degree and separated by an intervening gap from the shallower seismogenic zone in many margins (e.g., Obara et al. 2011;Husker et al. 2012;Takagi et al. 2016) and particularly evidenced in Cascadia (Wech et al. 2009;Hyndman 2013;Hyndman et al. 2015).This separation along depth is proposed to be mainly due to rheologic variations resulting from the interplay between temperature and high-fluid pressure (Gao and Wang 2017) and may occasionally be filled by long-term SSE (e.g., Cavalié et al. 2013;Takagi et al. 2016;Rousset et al. 2019;Nuyen and Schmidt 2021) or even penetrated through by some energetic deep ETSs (Kano et al. 2019;Bartlow 2020) during the interseismic phase of the great megathrust earthquake cycle at shallow depths.Despite being rarely observed or tightly constrained, many studies of global subduction zones have also suggested that the ETS zone (and the intervening gap) is prone to creep aseismically as afterslip under the coseismic stress perturbations following shallow great megathrust earthquakes, as in Mexico, Ecuador, Alaska, Nankai and Hyuganada (off Kyushu, Southwest Japan) cases (Yarai and Ozawa 2013;Gualandi et al. 2017;Rolandone et al. 2018;Takagi et al. 2019;Sherrill & Johnson 2021;Okada & Nishimura 2023).Such slip/stress-release behaviors of the ETS zone and the intervening gap during these earthquake-cycle phases indicate that slip-deficit or tectonic stress has already been accumulated in the fault system as a background process during the long-term interseismic phase of megathrust earthquakes and/or the quietness time between aseismic events (i.e., inter-SSE period), considering the balance between the sum of seismic and aseismic slip and the long-term fault slip (i.e., plate convergence rate) (Avouac 2015).Therefore, if the slip-deficit or locking degree is high enough at the ETS zone, it should result in observable crust deformation, which could be, in turn, inverted for the slip-deficit rate with the back-slip method (Savage 1983).Such slipdeficit rate estimates and comparisons between different subduction zones, such as Nankai and Cascadia, should provide kinematic constraints for inferring the dynamics of the global ETS zone as a whole.
Two difficulties, however, lie in determining the locking at the ETS zone: first, geodetic data generally have poor resolution at the ETS depth (30-40 km and deeper) (e.g., Wang & Tréhu 2016); second, long-term interseismic deformation involves viscoelastic mantle relaxation (e.g., Wang et al. 2012;Li et al. 2020;Li and Chen 2023).Despite the resolution issue, previous studies using elastic inversions suggest that the locking degree at ETS depth may differ during the inter-SSE and the long-term interseismic phase at multiple subduction zones (Radiguet et al. 2012(Radiguet et al. , 2016;;McCaffrey 2014;Li et al. 2016;Xie et al. 2020;Saux et al. 2022) and may also vary on a decadal timescale in Japan subduction zones (Nishimura et al. 2004;Loveless and Meade 2016).These results together suggest that there may indeed exist certain short-wavelength signals from the relatively narrow ETS zone buried in the long-wavelength strong signals from the shallow large-scale megathrust locking.
On the other hand, Li et al. (2015), through both synthetic simulations and a case study of surface deformation in North Chile, suggested that ignoring viscoelastic relaxation during the late interseismic phase in a purely elastic Earth model would lead to an overestimation of downdip locking depth, i.e., the similar depth of the ETS, during this phase.This conclusion is later confirmed by many case studies of other margins (e.g., Li et al. 2018;Itoh et al. 2021;Diao et al. 2022).Given the elapsed time of last great megathrust earthquake and the average recurrence interval (Ishibashi 2004), Nankai in southwest Japan is expected currently during the earlyto-mid-interseismic phase and features a short-wavelength viscoelastic deformation pattern (Li et al. 2020;Li and Chen 2023).It is, however, unclear how viscoelastic relaxation affects estimation of the locking depth in this phase.Previous locking studies of Nankai commonly assume either a purely elastic Earth (e.g., Ozawa et al. 1999;Miyazaki and Heki 2001;Tabei et al. 2007;Liu et al. 2010;Yoshioka and Matsuoka 2013;Nishimura et al. 2018), or a vertically layered viscoelastic Earth (e.g., Ito and Hashimoto 2004;Noda et al. 2018), in which ignoring an elastic slab may also lead to a determination of more back-slip at the ETS depth (Li and Chen 2022).Overall, it is, therefore, not well explored the downdip seismogenic limit and slip-deficit accumulation at the ETS zone in Nankai using contemporary deformation and the framework of viscoelastic Earth with an elastic slab.Quantifying these aspects during the interseismic phase in Nankai and comparing them with Cascadia would help understand the overall slip budget at the ETS zone, the potential of coseismic slip penetrating into the ETS zone, and afterslip during the postseismic phase at the ETS zone (e.g., Dessa et al. 2009;Jiang and Lapusta 2017;Sherrill and Johnson 2021).Li et al. (2020) focused on the long-term (i.e., the time scale of the earthquake cycle) mantle rheology and the time-dependent interseismic deformation.They used a series of two-dimensional viscoelastic finite-element models to explore the earthquake-cycle deformation in Nankai with the century-long leveling data (Thatcher 1984) and contemporary Global Navigation Satellite System (GNSS) velocity field (Nishimura et al. 2018) across the Shikoku Island.They found that locking duration (i.e., the elapsed time since the fault is relocked) and locking depth both affect the time-dependent interseismic deformation pattern, and the long-term continental mantle viscosity should be in the order of 10 19 Pa s.In their models, the transition zone (i.e., from fully locked to fully creeping) downdip of the seismogenic zone was fixed to be 10 km width.Here, we employ similar models and additionally purely elastic models to explore the kinematic state at depth of the ETS zone using the same GNSS data.Unlike Li et al. (2020), we treat the transition zone as an unknown and investigate its existence, depth and width.In this way, we gain insights into the slip budget at the ETS zone.As an additional step forward from Li et al. (2020), we explore the impacts of subduction-zone rheology structure, particularly the viscoelastic oceanic mantle and the elastic structure of the entire model domain, on inland surface deformation.With the aid of other observations (e.g., depth of the ETS zone), we extensively discuss the performances of elastic and viscoelastic models in terms of fitting the horizontal and vertical GNSS data, the associated implications for subduction dynamics in Nankai, and the similarities and differences between Nankai and Cascadia, which aids in understanding of dynamics of the ETS zone.

Geodetic data and finite-element models
We consider a 100-km-wide, margin-normal swath profile representative of the Nankai subduction zone to explore the downdip seismogenic limit in southwest Japan (Fig. 1a).This profile passes through the main source area of the 1946 Nankai Earthquake (Sagiya and Thatcher 1999;Ishibashi 2004;Sherrill and Johnson 2021) and the ETS belt beneath the Shikoku Island (e.g., Obara et al. 2010;Nishimura et al. 2013) (Fig. 1a).The locking estimation in SSE-occurring subduction zones is found to be sensitive to the time period during which surface velocities are measured (e.g., Saux et al. 2022;Maubant et al. 2023).Here, we consider the inland velocity field estimated during a relatively stable interseismic period from 2005 to 2009 with seasonal signals removed (Nishimura et al. 2018).This period is free of large earthquakes (especially the 2011 Mw 9.0 Tohoku-Oki earthquake) and large transient events (Nishimura et al. 2018).Small transient events were later detected in this time period from the same working group using improved time-series analysis techniques (e.g., Fujita et al. 2019;Okada et al. 2022).The seaward motions of these events may be partly removed as seasonal signals due to their low magnitudes but may still result in lower overall interseismic rates to some extent (e.g., Frank 2016;El Yousfi et al. 2023).Thus, the stable locking in this study may be slightly underestimated in terms of locking depth and/or locking degree.These velocities are defined with respect to the Amurian Plate reference frame (Nishimura et al. 2018).The landward horizontal velocity field decays with trench distance from ~ 50 mm/year near the Pacific coast to ~ 4 mm/year near the Sea of Japan coast (Fig. 1a), suggesting a typical interseismic deformation pattern due to certain megathrust locking.Moving away from the trench, the vertical velocity field features subsidence near the Pacific coast (6 mm/year maximum), a distributed uplift mainly on Shikoku Island (6 mm/year maximum), and a secondary subsidence zone on Honshu Island (1-2 mm/year) (Fig. 1a).
We build finite-element models using the open-source code PyLith (Aagaard et al. 2013) and perform both elastic and viscoelastic modeling to explain the GNSS velocity field.To compare the GNSS data within the relatively far reference frame and limit boundary effects, we set our model domain to 2000 km wide and 500 km deep.Our model is composed of four domains: oceanic and continental plates and oceanic and continental mantle (Fig. 1b).Our model thus includes the first-order control of the slab on the viscous mantle flow (i.e., the viscous flow cannot penetrate the slab) during both postseismic and interseismic phases (e.g., Pollitz et al. 2008;Johnson and Tebo 2018;Li and Chen 2022).Both the oceanic and continental plates are set to be 30 km thick (Li et al. 2020).Curved slab geometry and megathrust geometry are incorporated in the models based on the local 3D seismic velocity model of Hirose et al. (2008).The rigidity values of the elastic plates and mantle are assumed to be 48 and 64 GPa respectively, and the Poisson ratio is assumed to be 0.25 throughout the model domain.Displacements perpendicular to the lateral and bottom model boundaries are fixed to zero, while displacements parallel to these boundaries are not constrained.The top boundary is set free.
To investigate the locking state of the downdip portion of the megathrust, we define two free parameters to characterize the megathrust: the downdip depth D of full locking (i.e., locking degree = 1) from the trench and the depth range W of the transition zone (which may overlap the ETS zone) from the downdip of full locking depth (Fig. 1b).At further depths, the megathrust interface is imaged to be greatly thickened (e.g., Nedimović et al. 2003) and is thought to be fully creeping (e.g., Wada and Wang 2009).The locking of the transition zone in the dip direction is assumed to linearly decease to zero.In the purely elastic models, the mantle is also assumed to be elastic.In the viscoelastic models, we define two additional free parameters: the Maxwell viscosity of oceanic mantle η o and the Maxwell viscosity of continental mantle η c (Fig. 1b).Since smoothing and regularization of con- ventional inversion schemes potentially dampen out the true locking of the ETS zone, we in this study mainly rely on forward simulations of surface deformation in fitting the GNSS data.Megathrust locking is simulated by kinematically prescribing a back-slip rate distribution (Savage 1983) that is equal to the locking ratio distribution multiplied by the profile-parallel convergence rate (63.3 mm/ year, Nishimura et al. 2018).In the viscoelastic models, given mantle viscosity, the earthquake recurrence time is important to the model results (Li and Chen 2023) and is assumed to be 150 years according to the study of historical earthquakes (e.g., Ishibashi 2004).For simplicity, we assume that the slip deficit accumulated over each earthquake cycle is entirely released by the coseismic rupture (Li et al. 2020).We spin up our viscoelastic models into a steady state by running them over ten earthquake cycles (e.g., Hetland and Hager 2006) and take the velocities at 60 years after the last earthquake (the same as the time that the GNSS velocities were measured) to compare with the GNSS data.We evaluate the goodness of fit by the root mean square (RMS) misfit and consider not only the horizontal deformation as conventionally done, but also the vertical deformation, which is sensitive to and may be diagnostic of viscoelastic mantle structure and kinematics on the megathrust (e.g., Li and Chen 2022).
Beside the GNSS constraints, we use the prior knowledge of the downdip ETS zone depth of 30-40 km (e.g., Obara

Sensitivity of Nankai rheology structure
We first use the GNSS data to explore the sensitivity of subduction-zone rheology with both elastic and viscoelastic models, with no transition zone (i.e., W = 0 km) (Figs. 2 and 3).A new finding compared to Li et al ( 2020) is that elastic models exhibit an overall good fit to both horizontal and vertical GNSS data (Fig. 2a, andb), with an optimal full locking depth of ~ 30 km (i.e., the locking depth with minimum misfit in Fig. 2c).This depth is consistent with previous elastic studies (e. they require less fault back-slip to explain the data, and 3D elastic models may lead to full locking depth deeper than 30 km.However, since the ETS zone (i.e., 30-40 km depth) is expected to be partly locked or even almost creeping, as the average slip rate of slow slip events is lower than one third of the plate convergence rate (e.g., Okada et al. 2022;Nishimura et al. 2013), the full locking should, to some extent, be shallower than 30 km depth.We interpret this overestimation of purely elastic models as partly due to ignoring interseismic viscoelastic relaxation (e.g., Li et al. 2015Li et al. , 2018)), despite Nankai being thought to be in the early interseismic phase (Li et al. 2020).Elastic models may overestimate the locking depth by ~ 10 km in the later interseismic period time due to stronger viscoelastic relaxation (e.g., Li et al. 2015Li et al. , 2018)).The difference between elastic and viscoelastic models was not explored in Li et al. (2020), and the 3D viscoelastic effect and its time-dependency are issues we are working on.Viscoelastic models without a transition zone prefer 1-3 km shallower full locking depth than the elastic models (Fig. 3a) to fit the horizontal GNSS data, confirming the previous findings (e.g., Li et al. 2015Li et al. , 2018)).This result is not surprising due to the short-wavelength viscoelastic deformation in the early interseismic phase (Li et al. 2020;Li and Chen 2023), consistent with the shortwavelength elastic model predictions (Li et al. 2015).Despite the optimal elastic model shows a slightly lower RMS misfit (comparing the gray curve to color-coded curves in Fig. 3a), the misfit difference is below the horizontal observation errors (~ 1 mm/year).Because incorporating the Earth's viscoelasticity is thought to be more physical (Wang et al. 2012;Li et al. 2015), we believe the obtained viscoelastic locking depth to be more realistic.
A relatively lower oceanic mantle viscosity (10 19 Pa s) apparently results in overall larger horizontal and vertical RMS misfit (blue curves in Fig. 3a, b), consistent with previous notes of high oceanic mantle viscosity (10 20 Pa s) (Wang et al. 2012).However, a relatively lower continental mantle viscosity (10 19 Pa s) also results in overall larger horizontal and vertical RMS misfit (longdashed curves in Fig. 3a, b, and Figure S1), suggesting that the continental mantle viscosity may be larger than the previously preferred interseismic value of 10 19 Pa s in Nankai and other margins (e.g., Wang et al. 2012;Johnson and Tebo 2018;Li et al. 2020).Judging from both the full locking depth of minimum horizontal and vertical RMS misfit, its spatial relation with the ETS zone (i.e., full locking depth < 30 km), and the fact that dehydration of the subducting slab causes a lower mantle wedge viscosity than that of the oceanic mantle (Wang et al. 2012), we prefer 5 × 10 19 Pa s for the continental mantle and 10 20 Pa s for the oceanic mantle (green solid curve in Fig. 3a, b).Despite contrasting tectonic settings, these values in Nankai are similar to the values obtained from North Chile (Li et al. 2015).In general, the viscoelastic models prefer a deeper locking depth from the vertical data (Fig. 3b) than the horizontal data (Fig. 3a).This result is the same for the elastic modeling results (Fig. 2c).Since inclusion of a cold nose in the model would make the modeling results more similar to those of elastic models (Fig. 2), we think that a cold nose is unlikely the solution to this inconstancy.

Inference for locking state at the Nankai ETS zone
We further explore the existence and width of the transition zone (i.e., W) using various viscoelastic structures (Figure S2) and illustrate the main findings with the preferred viscoelastic structure (i.e., η c = 5 × 10 19 Pa s and η o = 10 20 Pa s) (Fig. 4).We find that horizontal and vertical GNSS data generally favor low and high W, respectively (Figs. 4a, b and S2).We also find a clear trade-off between D and W in explaining both horizontal and vertical GNSS data (Fig. 4a, andb).That is, a relatively lower D combined with a relatively higher W can equally explain the data (e.g., dark blue and dark green color in Fig. 4a, andb, respectively).This trade-off is also clearly seen in the results of other viscoelastic models and even the elastic model (Figures S2 and S3), indicating inherent ambiguity in determining both D and W with GNSS data alone.This result thus suggests revisiting previous kinematic (locking) models that assume a constant transition zone.Overall, models with locking penetrating into the ETS zone (i.e., 30 < D + W < 40 km) exhibit a good fit to both horizontal and vertical GNSS (e.g., Fig. 4a, andb), suggesting certain locking and hence slip-deficit accumulation at the ETS zone.This finding is consistent with the optimal value (i.e., D + W = 37.5 km) determined in Li et al. (2020) with the assumption of W = 10 km, η c = 10 19 Pa s and η o = 10 20 Pa s, and also a recent advanced inversion work of Sherrill et al., (2024).In comparison, models with D + W < 30 km (e.g., below the green line in Fig. 4a, andb), which indicate the existence of an intervening gap (e.g., Hyndman 2013; Gao and Wang 2017), perform less well in fitting the horizontal (Figs.4a and S2a) and particularly vertical GNSS data (Figs.4b and S2b).From the horizontal RMS misfit, models tend to yield limited or no transition zone (i.e., dark blue part below the green line in Figs.4a and S2a).We therefore argue for limited or no intervening gap in Nankai, suggesting the long-term SSEs fully bridge the ETS zone and the shallow megathrust zone.Given the depth range of the ETS zone (i.e., 30-40 km) and the thermomechanical structure in Nankai (e.g., Hyndman et al. 1995;Gao and Wang 2017), we think that the models with D + W > 40 km are less physical (white shadings in Figs. 4 and S2).
Assuming megathrust locking extends to the downdip of the ETS zone (i.e., D + W = 40 ± 2 km), we further discuss the preferred D and W with the specific results of the preferred viscoelastic structure (magenta dashed lines in Fig. 4a, b).Results show that horizontal and vertical GNSS data prefer consistent D (18 ± 2 km) and W (22 ± 2 km) (i.e., yellow shadings in Fig. 4c, d, respec- tively).From the misfit curves, the preference for these values is clearer from the horizontal data (red curves in Fig. 4c, d), suggesting that horizontal data is more diagnostic of the locking state of downdip seismogenic zone than the vertical data (blue curves in Fig. 4c, d).Compared with the results of Li et al. (2020) (i.e., D = 27.5 km and W = 10 km), this preferred model (i.e., D = 18 km and W = 22 km) has a relatively shallow full locking depth and a relatively wide transition zone, half overlapping with the ETS zone.If assuming a shallower megathrust locking, the preferred D increases and the preferred W decreases (Figure S4), due to the trade-off between these two parameters.The viscoelastic structure also influences the determinations of the locking depth (Li et al. 2020) and W but reasonable structures show similar values of D and W (lower right panels in Figure S2).This suggests that the main earthquake source might be shallow (< 20 km) with possibly large tsunami potential.Because we assume linear decay of locking in the transition zone, the maximum locking degree of the preferred model in the ETS zone (Fig. 4) is ~ 0.5, i.e., ~ 30 mm/year slip-deficit accumulation.This value is almost identical to the long-term average slip rate of the ETS zone (e.g., Okada et al. 2022;Nishimura et al. 2013), suggesting limited final slip-deficit accumulated through the earthquake cycle and hence possibly limited to no coseismic slip and/or afterslip penetrating into the ETS zone.Since this deduction relies on several above-mentioned assumptions, it suffers from uncertainties in model parameters and requires validation from earthquake source documentation of future events.Traditional asperity models at subduction zones (e.g., Scholz 1998;Lay et al. 2012) suggest that the downdip seismogenic zone is composed of small-scale asperities and hence has a limited capacity for accumulating slip deficit.From this perspective, our preferred model is compatible with traditional asperity models.

Comparison with Cascadia
Nankai in southwest Japan and Cascadia in northwest North America are two of the most well-studied and well-instrumented seismic gaps at subduction zones.These two margins have long been found sharing many similar features, such as oblique plate convergence, relatively young and warm subducting slabs, low stress states, and ample ETS events (e.g., Wang and He 1999;Wang 2000;Dragert et al. 2001;Hirose & Obara 2005).However, when synthesized with thermal models (Hyndman et al. 1995;Hyndman 2013), the kinematic locking state obtained in Nankai in this study contrasts with that of Cascadia (e.g., Schmalzle et al. 2014;Li et al. 2018), particularly in terms of the ratio of D + W to H (i.e., the elastic upper-plate thickness) (Fig. 5).Despite the suppression of SSE signals in the geodetic data (McCaffrey et al. 2013;Nishimura et al. 2018) and the inclusion of interseismic viscoelastic effects, the associated viscoelastic locking models (this study and Li et al. 2018) may still underestimate the locking depth to some extent due to some residual or unresolved SSE signals (e.g., Frank 2016;Saux et al. 2022;El Yousfi et al. 2023).Considering the thermal structures in the two margins (Gao & Wang 2017), their contrasting ratios of D + W to H are still thought to be validated (Fig. 5).Nankai and Cascadia fit perfectly into the two end-member models, where the downdip seismogenic limit is primarily controlled by the fore-arc Moho intersection and the isotherm of ~ 450 ℃, respectively (e.g., Hyndman et al. 1997;Oleskevich et al. 1999;Hyndman 2013).Given the ETS-favoring rheological and petrological conditions at the fore-arc Moho intersection (Gao & Wang 2017) and similar Moho depths in the two margins, their contrasting locking depths result in a clear intervening gap in Cascadia and no such a gap in Nankai (Fig. 5).Multiple other lines of evidence also support the existence of an intervening gap in Cascadia, including thermal structure and paleo-seismic coastal subsidence (Hyndman 2013), while limited or no similar evidence of the gap has been reported in Nankai.How the presence of intervening gap influences the earthquake and long-term geodynamics requires further investigations.
The kinematic representations in Nankai and Cascadia (Fig. 5) cannot be directly used to infer earthquake dynamics or frictional behavior, particularly at the ETS depths (e.g., Wang and Dixon 2004).One important reason is the stress shadowing effects on both updip and downdip portions of the seismogenic asperity (Herman and Govers 2020;Lindsey et al. 2021).Further theoretical investigations of kinematic locking with laboratoryderived laws (e.g., the rate and state friction law) have proven insightful for understanding earthquake and tsunami hazards (e.g., Baranes et al. 2018;Yang et al. 2019;Ramos et al. 2021;Melgar et al. 2022).However, it is still debated whether the brittle-ductile transition zone participates in the coseismic rupture or creeps aseismically after the earthquake (e.g., Bilham et al. 2017;Jiang and Lapusta 2017;Rolandone et al. 2018).In other words, it is largely unclear if there is a persistent spatial physical separation, or segmentation, between seismic and aseismic regions along both fault dip and strike directions.Likely, both short-term dynamic factors and long-term tectonic factors play roles.A comprehensive comparative study of different subduction zones (even collision zones), together with long-term across-time-scale earthquake-cycle simulations (e.g., Julve et al. 2024), is hence recommended for the future.
Despite recent tremendous progress on documenting and modeling SSE in the two margins, such contrasting long-term seismogenic behaviors (Fig. 5) tend to be underappreciated.Mechanically, this ratio (i.e., (D + W)/H) determines the tectonic stresses introduced into the viscous domain and hence varies modeled viscoelastic earthquake-cycle deformation (e.g., Hashima and Sato 2017;Li and Chen 2023).Further validation of these modeling results and refinement of the first-order physical processes throughout the entire earthquake cycle require a long-term continuous monitoring strategy and comparison studies of margins that are in different earthquake-cycle phases (Wang et al. 2012;Li and Chen 2023).

Conclusion
We perform a series of 2D elastic and viscoelastic forward locking models to fit the horizontal and vertical GNSS data in the southwest Japan subduction zone and pay particular interest in the existence, width, and locking state of the transition zone at the downdip seismogenic zone.Based on the modeling results, we draw the following conclusions.We confirm the previous finding that elastic model tends to overestimate the locking depth (Li et al. 2015) but the overestimation is much small due to the early interseismic phase in Nankai (Johnson and Tebo 2018;Li et al. 2020).The inland GNSS data is sensitive to the lower bounds of oceanic and continental mantle viscosity (i.e., both > 10 19 Pa s).Based on the modeling results of models, we propose preferred values for η c and η o as 5 × 10 19 Pa s and 10 20 Pa s, respectively.Assuming the locking extends Fig. 5 Seismogenic limits and their spatial relations with the ETS zones and thermal structures in (a) Nankai and (b) Cascadia, compared with the elastic upper plate and viscoelastic continental mantle along depth (sidebars on the right).The seismogenic states in the two margins are synthesized from both viscoelastic locking models (this study and Li et al. 2018) and thermal models (Hyndman et al. 1995;Hyndman 2013), representing long-term averaged states over multiple SSE cycles.These thermal models suggest that the megathrust is fully locked in the temperature range of 150-350 ℃ and partially locked or creeping in the temperature ranges of 0-150 ℃ and 350-450 ℃ in the updip and downdip portions of the fully locked zones, respectively.The partial locking near the trench in Nankai is also supported by seafloor geodetic observations (Yokota et al. 2016).The intervening concept in Cascadia is originated from Hyndman (2013).L-SSE, S-SSE and tremor depths are mainly from Dragert et al. (2001), Obara (2002), Rogers and Dragert (2003), Obara (2011) and Nuyen and Schmidt (2021).Grey contours are thermal structures from Gao and Wang (2017).H represents the elastic thickness of continental plate.Red line in the inset plots show the location of two-dimensional cross section.Figure is after Gao and Wang (2017) into the downdip of the ETS zone (i.e., 40 km), we propose a preferred model with a full locking depth of 20 km and a broad transition zone spanning a 20-km depth range.In this model, downdip half of the transition zone overlaps with the ETS zone.Therefore, unlike the Cascadia case (Hyndman 2013), there is likely limited or no intervening gap, and partial locking penetrating into the ETS zone, introducing more tectonic stresses into the viscous domain in Nankai (Fig. 5).Due to SSEs and other stress-releasing events, there may be neglectable slip-deficit accumulation at the ETS zone, indicating limited or no coseismic slip and afterslip may propagate into the ETS zone in Nankai.This interpretation is contrary to some geodetic modeling and numerical studies.However, a precise quantification of the slip budget for coseismic slip and afterslip demands a refinement of the interseismic locking model in terms of slab geometry and rheology structure, along with a continuous documentation of both short-term and long-term slow slip events.To deepen the understanding of kinematic and dynamic behaviors of the ETS zone over the earthquake cycle, comparison studies of the Nankai, Cascadia and other subduction zones are required.

Fig. 1 a
Fig. 1 a Seismotectonic settings and horizontal and vertical GNSS velocities at Nankai.The GNSS velocities are from Nishimura et al. (2018).Red shading is the preferred coseismic slip distribution of 1946 and 1944 megathrust earthquakes from Sherrill and Johnson (2021).Yellow shading represents regions of short-term slow slip events during 1996-2012 (Nishimura et al. 2013), and dark green dots indicate tremors observed during 2006-2009(Obara et al. 2010).Blue shading is the source area of 1968 Mw 7.5 Hyuga-nada earthquake(Yagi and Kikuchi 2003).The gray dashed lines are depth contours of the plate interface at 10-km intervals, fromHirose et al. (2008).Dashed green lines show the swath profile of the used geodetic data and solid green line shows the location of megathrust geometry used for finite-element models.b Illustration of the four explored free parameters of the modeled subduction zone system.The subduction structure is simplified for illustration purpose and the realistic slab and megathrust geometry has been considered in the finite-element models(Li et al. 2020)

Fig. 2 Fig. 3
Fig. 2 Performance of elastic models without a transition zone (i.e., W = 0 km) on fitting (a) horizontal and (b) vertical GNSS data (black dots with error bars).c RMS misfit as a function of full locking depth

Fig. 4
Figure 4c & 4d Okada et al. 2022;Nishimura et al. 2013)himura et al. 2013;Obara & Kato 2016;Takagi et al. 2019;Yabe et al. 2023)2023)to assess the GNSS-constrained parameter spaces of D and W. In particular, the long-term (over a decade) average slip of SSEs can be in a rate as high as one third of the plate convergence rate (e.g.,Okada et al. 2022;Nishimura et al. 2013), indicating the partial locking in the ETS zone and D + W > 30 km.