Influence of Subduction Zone Dynamics on Interface Shear Stress and Potential Relationship With Seismogenic Behavior

Whether tectonic convergence at subduction zones is accommodated predominantly through seismic or aseismic deformation, the former potentially generating large earthquakes, varies considerably between subduction margins. This margin‐scale variability has previously been linked to overriding plate deformation, trench migration, and their influence on the plate interface stress state. While these processes are linked to mantle‐scale dynamics, it is unclear how such dynamics influence interface stress. We systematically analyze the interface stress state in a suite of 2‐D thermo‐mechanical subduction models, where slabs display a range of morphologies that arise from diverse multiscale interactions with adjacent mantle and the overriding plate. We demonstrate that the thickness of the interface layer varies dynamically, in response to Poiseuille flow induced by slab bending or unbending, leading to associated effects on interface shear stress at typical seismogenic depth. Lower shear stress occurs when slab unbending is significant, which is commonly associated with trench retreat and draping of the slab as it impinges on the higher‐viscosity lower‐mantle. Conversely, higher shear stress is associated with limited slab unbending, which is promoted by negligible trench migration and vertically subducting slabs. We conclude that the diversity of slab dynamics may cause large variations in interface stress state between and maybe within margins. This is an additional variable that potentially controls seismogenic behavior, and we compare broad stress estimates for Circum‐Pacific margins to previous studies. Although predicted shear stress varies with observed seismogenic behavior, more detailed constraints on stress state are needed to test for correlation.

forces driving subduction (e.g., Garel et al., 2014;Holt et al., 2015). Such large-scale models have been used to infer a low interface stress for sustained subduction (10 MPa order; Duarte et al., 2015), but have not yet been used to understand stress variations between subduction zones. Here we explore this novel perspective, quantifying the interface stress state in mantle-scale geodynamic models in which subduction self-consistently evolves, with the aim of testing whether some slab dynamics are conducive to higher interface shear stress than others. We analyze a suite of numerical subduction models described by Garel et al. (2014), who modeled thermo-mechanical subduction dynamics in 2-D (Figure 1), for a broad range of subducting and overriding plate ages. We then use this analysis to estimate relative variations in the interface stress state of Circum-Pacific subduction zones thought to have contrasting seismogenic behavior.
Global variations in subduction seismogenic behavior have been correlated with contrasts in large-scale geodynamic processes, such as aseismic margins with overriding plate extension and/or slab rollback, and seismogenic margins with shortening or negligible deformation of the overriding plate and/or negligible rollback (Conrad et al., 2004;Heuret et al., 2012;Schellart & Rawlinson, 2013;Scholz & Campos, 1995Uyeda & Kanamori, 1979;Wallace et al., 2012). It is therefore of interest how these processes relate to the interface stress state. As the numerical models we analyze reproduce variations in slab rollback, but not overriding plate stress, we primarily analyze the relationship between rollback and interface stress, as well as related slab-mantle dynamics.
The Gutenberg-Richter earthquake frequency-magnitude distribution, describing the relative frequency of large and small earthquakes in a particular region, has been reproduced in models of frictional slip on faults with heterogeneous strength, geometry, and/or stress distributions (Ben-Zion, 1996;Hillers et al., 2007;Ripperger et al., 2007). Such heterogeneous faults can deform at a range of background stresses and can subsequently differ in slip styles despite having identical frictional properties. The probability of a rupture growing to a large area and therefore moment magnitude in these models is dependent on the background shear stress acting on the fault at the time of rupture initiation (Fang & Dunham, 2013;Ripperger et al., 2007).
BEALL ET AL. 10.1029/2020GC009267 2 of 20  Garel et al. (2014) to be analyzed in this study. The models reproduce a variety of trench retreat velocities and slab morphologies, depending on the initial subducting and overriding plate ages (age SP and age OP ). Squares denote initial trench position. (a and c) show the end-member models chosen in this study as examples of trenches that are stationary or rapidly retreating, while (b) is an intermediate case. Alternatively, these models also indicate that large ruptures could occur at relatively low background stress on smooth, continuous faults. Contrasts in energy radiation imply that great earthquakes may vary between these smooth (low complexity) and rough (high complexity) types, such that the role of tectonic stress may differ (Ye et al., 2018). Fault stress state may also depend on a number of factors, including pore pressure, variations in frictional strength, and previous seismicity (Brodsky et al., 2020), as well as regional tectonic stresses (Hardebeck, 2010). An understanding of how interface stress may vary between margins is therefore required to assess these competing influences on seismogenic behavior.
The influence of geodynamics on megathrust seismicity has been explored with geodynamic models incorporating the earthquake cycle (e.g., Brizzi et al., 2020;Van Dinther et al., 2013), which have the advantage of explicitly coupling both processes. However, such models are computationally expensive and cannot capture some aspects of the larger subduction system, such as dynamic driving forces (plate kinematics are typically imposed) or slab-lower-mantle interaction. Interaction between the slab and mantle is thought to be an important aspect of the subduction zone force balance (Schellart, 2004;Scholz & Campos, 1995) and is required to reproduce the diverse range of imaged slab morphologies on Earth (Garel et al., 2014;Goes et al., 2017). We aim to accurately capture the long-term interface stress state associated with this range of mantle-scale slab dynamics, with the limitation that the earthquake cycle cannot be explicitly modeled.

Modeling Background
We analyze subduction models from Garel et al. (2014), examples shown in Figure 1, which were constructed using the finite element, control-volume code Fluidity (Davies et al., 2011;Kramer et al., 2012). These models include both subducting and overriding plates, with the upper-lower mantle transition represented by a factor of 30 viscosity jump, at 660 km depth. A 5 km thick weak decoupling layer, at the interface between subducting and overriding plates, and a free-surface, facilitate trench motion. This layer is thin compared to typical geodynamic models (e.g., 15 km in Holt et al., 2017), achieved through mesh refinement, though thicker than typical fault-zones (100s m scale; Rowe et al., 2013). Typical interface stress magnitudes of ∼10 MPa (Duarte et al., 2015) are reproduced, indicating that the stress scaling is realistic. Plate and mantle rheology depend on temperature and stress. Thermal structure controls slab thickness, density, and viscosity. Accordingly, slab strength (resistance to bending and/or stretching deformation dependent on both effective viscosity and thickness) and buoyancy (the integral of density over slab thickness) evolve self-consistently with the underlying thermal state. As the models are 2-D, they are most applicable to subduction zones in which along-strike variation in slab dynamics is negligible or gradual enough that out-of-plane stress heterogeneity can be ignored.
The numerical approach, which is underpinned by anisotropic mesh optimization, maintains computational efficiency even in large domains, thereby reducing the control of boundary conditions on the resulting dynamics (the model domain is 2,900 km deep and 10,000 km wide). An initial slab shape is prescribed within the mantle and subduction then self-consistently develops through time. The weak layer is initiated across the entire slab surface and is entrained during subduction, decoupling subducting, and overriding plates (shown conceptually in Figure 2). The location of this weak decoupling layer is tracked using a volume fraction, the evolution of which is described by a linear advection equation. In order to avoid excessive numerical diffusion of the weak layer into neighboring regions, this is discretized on the control volume mesh using the minimally diffusive HyperC face-value scheme (Wilson, 2009). Anisotropic adaptive mesh optimization, with a minimum edge-length of 500 m, ensures that this weak layer is well-resolved. All materials follow a composite rheology combining diffusion and dislocation creep, which introduce temperature and strain-rate dependencies, as well as a combination of a yield stress and Peierls creep, which are stress limiting and allow slab bending to occur. The viscosity of the combined deformation mechanisms is incorporated as a harmonic mean of the effective viscosities for each. Viscosity is capped between 10 18 -10 25 Pa s. The weak layer viscosity is capped at ≤10 20 Pa s and prescribed a low friction coefficient of 0.02. Further details of the modeling setup and solution strategies are described in Davies et al. (2011) and Garel et al. (2014).
This study focuses on quantifying how mantle-scale slab dynamics influence interface shear stress. We analyze a set of 16 models that reproduce a variety of slab dynamics (with various slab morphologies and BEALL ET AL.

10.1029/2020GC009267
trench velocities), controlled by the range of subducting slab and overriding plate ages, termed age SP and age OP , which vary between 20, 30, 40, 65 and 100 Ma in different combinations. We do not analyze all of these combinations (which would be 25 models), as we exclude models that involve extreme yielding and entrainment of the overriding plate into the interface layer, because their interface thickness evolution was no longer systematic and potentially unrealistic. Model time-steps <20 Myrs are excluded, such that the prescribed initial subduction geometry has a negligible influence on our results.
Younger subducting slabs in the models are initially thinner and, therefore, more buoyant and deformable, compared to older slabs. Subduction of older slabs typically involves trench retreat, whereas younger slabs subduct with a stationary to slightly advancing trench ( Figure 1). Older overriding plates are thicker and limit rollback and slab unbending, though the influence on rollback is less pronounced in the models of Garel et al. (2014) compared to those with a fixed overriding plate (Alsaif et al., 2020). The overriding plate does not have a composition-dependent density, although it is gravitationally stable (due to its lower temperature and higher viscosity) over the evolution times of our models and can represent either oceanic or continental lithosphere, the latter represented by older plate ages.
Interface layer thickness and shear stress can be directly measured in the depth range relevant to the seismogenic zone, however it is more complex to quantify slab dynamics. The resistance of the mantle to slab rollback ( Figure 2) can play a significant role in the interface force balance (Scholz & Campos, 1995) and is quantified here by the horizontal trench retreat velocity, measured at the surface and relative to a stationary lower-mantle. This is likely to represent the force required to drive mantle flow, which is approximately proportional to rollback velocity in 3-D free-subduction models (Stegman et al., 2006), though an exact relationship has not yet been constrained for 2-D models with an overriding plate. Significant slab unbending ( Figure 2) as it enters the asthenosphere is also associated with slab rollback (Lallemand et al., 2008;Scholz & Campos, 1995) and is later shown to be a driver of interface Poiseuille flow. Slab pull is also a key part of the subduction zone force balance, though is likely to contribute less to the subduction force balance than mantle flow during rollback (Schellart, 2004).
Both slab bending moment and pull are visualized and measured directly from the slab stress state. Slab pull is quantified by calculating the average axial deviatoric normal stress within the slab at 100 km, defining the base of the slab as the 800 o C isotherm. The magnitude of slab bending torque at 100 km depth is quantified by calculating the bending moment in a slab-perpendicular cross-section (as shown in Figure S1). A clockwise torque, related to increasing curvature, is taken as positive. Bending moment is zero when there is no switch in sign, as occurs in a small number of cases when there is no axial compression.
BEALL ET AL.
10.1029/2020GC009267 4 of 20 Figure 2. Schematic of the weak layer, which is analogous to a layer of weak (or weakened) sediment and oceanic crust. The weak layer with incoming thickness L 0 is transported with a vertically constant velocity profile by the moving seafloor. The weak layer thickness at depth, L w , varies with time and depth, reflecting the volume-flux controlled by the weak layer velocity profile (shown). The interface stress at depths relevant to the seismogenic zone is hypothesized to be influenced by deeper slab dynamics, identified as mantle resistance during rollback, slab bending moment, and slab pull.

Simplified Analytic Model
The weak layer evolves freely: although its 5 km thickness is imposed as the incoming thickness (L 0 in Figure 2) and initial condition at depth, it thickens or thins in response to the surrounding force-balance. This variation in thickness influences the interface stress, which varies inverse proportionally with thickness for constant convergence velocity v c and stress below yielding. Simple shear within the interface weak layer would transport weak material at a volume-flux of 0.5v c L w , for depth-dependent weak layer thickness L w (Figure 2). At steady state the incoming sediment thickness of L 0 (transported at the surface with a flux of v c L 0 ) would then increase to L w = 2L 0 for interface segments deforming by simple shear.
Mass flow within the interface weak layer is also modified by Poiseuille flow: down-or updip flow driven by pressure gradients that cause deviations of the velocity profile within the weak layer from simple shear ( Figure 3). Superposition of simple shear with updip Poiseuille flow (due to a downdip increase in pressure, assuming lithostatic pressure has been subtracted) then generates a concave velocity profile with decreased volume-flux (Figure 3a), increasing L w and potentially decreasing interface shear stress. Superposition with downdip Poiseuille flow (resulting from a downdip decrease in pressure) generates a convex velocity profile with increased volume-flux (Figure 3b), decreasing L w and potentially increasing stress. A simplified mass balance calculation, assuming interface layer properties relevant to the models of Garel et al. (2014), indicates that a downward pressure gradient of about 24 MPa/km is required to reduce the interface layer to L w = L 0 , but only an upward gradient of 0.9 MPa/km to thicken it to L w = 3L 0 (Text S1 and Figure S2). The analytic model shows that these pressure gradients depend on the model L 0 , v c , and weak layer viscosity, however the dependence on the latter two parameters would be neglected if the interface stress during simple shear is constrained.

Results
As there are 16 models, each with contrasting subduction dynamics and evolution through time, we present our analysis in three sections that demonstrate the interplay between slab dynamics and interface stress in varying granularity. First, we compare two models chosen as end-member cases that represent negligible and significant trench retreat, and analyze their stress states in detail. Second, we compile a data set capturing multiple time-steps from every model, in order to generalize the relationship between instantaneous slab dynamics and interface stress. Finally, we demonstrate how slab dynamics and interface stress can evolve significantly through time for a given subduction model.

End-Member Comparison
The model with age SP = 100 Ma and age OP = 20 Ma generates rapid trench retreat, while the trench is quasi-stationary when age SP = 20 Ma and age OP = 20 Ma (Figures 1 and 4a). These models are taken as end-member examples. The end-member with trench retreat has an interface layer that is thicker and deforms at lower shear stress, compared to the stationary trench end-member ( Figure 4b). The largest contrast in thickness occurs at 10 km depth, where L w = 4L 0 and 0.8L 0 for the retreating and stationary trench models respectively after 30 Myrs. However, the shear stress contrast is limited at this depth as both models are at the frictional yield stress. At depths of 20-30 km, relevant to the seismogenic zone, the stationary trench model has a weak layer thickness 3× thinner and shear stress almost 2× larger than the trench retreat model ( Figure 5). This high stress in the stationary trench model is sufficient to cause yielding down to depths of 20-35 km (varying temporally), compared to about 10 km for the retreating BEALL ET AL.  trench model ( Figure S3). The stress in the trench retreat model increases with depth until the models have similar interface stress at 100 km depth, such that the stress contrast is most significant in the upper half of the interface.
Velocity profiles are extracted perpendicular to slab dip and centered at 20 km depth ( Figure 5, profile positions shown in Figure 4c). They exhibit concave and convex shapes, explained by an increasing downdip pressure gradient for the model with a thick interface layer and a decreasing pressure gradient for the thin interface layer, respectively. The corresponding effective pressure gradients required to produce such deviations from simple shear (solving for P x in Equation S1) are 15 MPa/km and −0.5 MPa/km, which are similar to the magnitudes predicted from the simplified analytical model ( Figure S2).
BEALL ET AL.

10.1029/2020GC009267
6 of 20 Downdip pressure gradients within the interface layer are influenced by stresses in the adjacent slab and overriding plate at the lithosphere scale. These stresses are analyzed as the direction of principal deviatoric compressive stress and the maximum deviatoric shear stress at each point of an interpolated grid ( Figure 4c). The total stress acting on the interface is the sum of the deviatoric stress (shown) and pressure (defined to exclude lithostatic pressure). Slab-dip-perpendicular normal deviatoric stress is effectively zero within the interface layer, due to its low viscosity. This normal stress, as well as pressure, in the adjacent slab and overriding plate material is therefore converted completely into the pressure gradients that drive Poiseuille flow within the interface layer.
The overriding plate adjacent to the interface layer in the trench retreat end-member model is in slab-perpendicular deviatoric compression at all depths (Figure 4c), varying from a maximum of 150 MPa at 30-40 km depth to a maximum of 50 MPa at <10 km depth. The larger deviatoric compressive stress at depth is responsible for the increasing downdip interface pressure. This compression at depth reflects upwards pressure from the slab (red arrow in Figure 4c), which is experiencing axial stress driving unbending as evident in deviatoric axial compression at the slab top (depths > 40 km), while unbending also results in compressional pressure at the slab top. This unbending corresponds to a reduction in slab curvature, rotating the slab anticlockwise. Though the direction of principal deviatoric stress within the slab top is orthogonal to the overriding plate, this can be reconciled with a force balance in which the slab-perpendicular deviatoric compression in the overriding plate is balanced by compressional pressure (rather than deviatoric stress) in both the interface layer and the slab. The shallow mantle wedge return flow may contribute to the overriding plate deviatoric compression, however this also occurs in models with a thick overriding plate and therefore deeper and/or limited mantle wedge flow. Slab dynamics are therefore interpreted to have the greatest influence on the interface layer pressure gradient, as further demonstrated in Section 3.2.
The overriding plate in the stationary trench end-member model is mostly in slab-dip-perpendicular deviatoric tension close to the interface (Figure 4c), contrasting the retreating trench model. The upper 10 km is in compression of similar magnitude to the deeper deviatoric tension (50-100 MPa), resulting in the large normal stress gradient responsible for downdip Poiseuille flow. The slab top is in dip-perpendicular deviatoric compression, but will also be experiencing high tensile pressure due to the slab bending away from the overriding plate, reflecting the tensional stress applied to the overriding plate (red arrow). This is compatible with lower pressure in the weak layer at increasing depth, until the onset of unbending stresses at about 70 km depth (which also corresponds to a switch to overriding plate compression). The slab is deforming BEALL ET AL.
10.1029/2020GC009267 7 of 20 by downdip stretching, corresponding to the tensional downdip deviatoric stress, at depths of 20-70 km, resulting from the slab pulling downward and away from the overriding plate. This pull contributes to the larger slab curvature relative to the retreating trench model. While the slab begins to unbend at about 40 km depth in the retreating trench end-member model, this transition to unbending occurs deeper at 70 km in the stationary trench model. This common unbending at >70 km depth may reflect the increasing similarity in interface stress between the models at these depths.

Instantaneous Interface Shear Stress and Thickness
Subduction dynamics vary through time in each of the 16 models. Causes for this include the evolving thermal and stress fields, interaction(s) between the slab and lower-mantle and, potentially, slab break-off. We quantify slab dynamics, as well as weak layer thickness and stress, for multiple time-steps from each model and treat these as instantaneous data, from which we check for general correlations.
Interface thickness (relative to L 0 ) and shear stress is analyzed for the 10-50 km depth range, covering the typical seismogenic depth range (Heuret et al., 2011;Hyndman et al., 1997), calculated as an average over a surface parallel to slab dip. These are compared to trench retreat velocity (related to mantle resistance forces), as well as slab bending moment and axial slab pull at 100 km depth ( Figure 6). The slab-mantle dynamics are related to the relatively shallow interface segment analyzed through the mechanism of slab bending and Poiseuille flow demonstrated earlier. This interaction can also be considered as a force and torque balance between the deep and shallow slab segments, though quantifying this balance is beyond the scope of our study. As some model time-steps involve anomalously low interface stress related to very low convergence velocities during periods of stalling subduction, statistics are calculated only for model steps with v c > 1.5 cm/yr (outlined points in Figure 6).
Average interface stress varies from 8 to 19 MPa and model data are generally clustered according to age SP , with both maximum stress and variability decreasing with increasing age SP . This clustering coincides with a narrow range in rollback velocity for each age SP . The highest interface shear stress occurs in models with negligible rollback (i.e., stationary trenches) and approximately decreases with increasing rollback velocity, with a correlation coefficient of r = −0.37 (Figure 6a). The interface thickness has an inverse relationship with rollback velocity, increasing with increasing rollback velocity and with a correlation r = 0.83 (Figure 6b). As stress below yielding is proportional to interface thickness for constant v c and homogeneous stress, the lower correlation in Figure 6a likely largely reflects variability in v c , however in nature plate velocities are homogenized along-strike at the plate scale. These comparisons indicate that rollback velocity, representing varying influence of mantle resistance on the interface force balance, has a first order control on interface thickness and stress range, while other dynamics are also required to explain stress variation.
The trench retreat end-member model involves significant axial deviatoric compression, likely linked to slab-lower-mantle interaction, which may reduce the axial slab pull force and therefore influence the interface stress. There is a broad correlation between interface shear and slab pull stresses (Figure 6c), with r = 0.55, where models with high interface stress generally have high slab pull stress. Low interface stress typically corresponds to low pull stress, with the slab even in net compression in some cases. Models with varying age SP overlap more than Figure 6a and high slab pull variation appears to correspond to high stress variation. There is little correlation between slab pull and interface thickness (Figure 6d), likely indicating that the pull stress influences interface stress by modulating the convergence velocity and/or stress distribution. It therefore appears that an interplay between mantle resistance during rollback and slab pull is responsible for the variation in interface stress.
Models with high interface stress generally have high bending moments, corresponding to reduced unbending as all models with v c > 1.5 cm/yr involve either unbending or negligible torque. Models with small bending moments (corresponding to high unbending) broadly have low interface stress, with a correlation of r = 0.44. An inverse correlation between slab bending and interface thickness has r = −0.73, indicating that slab bending influences interface stress by modulating the interface thickness, as expected for the Poiseuille flow process described earlier. Slab bending moment is also analyzed at 50 km depth, where there was clearly a contrast in bending between the end-member models in Figure 4c ( Figure S4). There is also a

Signi cant Unbending
Minor Unbending correlation between this shallow bending and interface stress, with r = 0.55, indicating the deep slab torque translates also to a shallow torque. The slab bending moments are similar in magnitude to the unbending component estimated (∼10 17 N m) for slab rollback by Lallemand et al. (2008). Another consequence of variation in slab unbending is its effect on the downdip length of the interface, which potentially varies inverse-proportionally with interface stress (Figure S5), such that negligible unbending may produce a shorter interface which focusses subduction forces.
While we have demonstrated relationships between interface stress and aspects of deep slab dynamics, there is scatter around the calculated linear regressions in Figure 6, which is likely to be improved in more careful force balance calculations. For example, the shear force acting at the interface changes orientation with dip. If the interface shear forces, as well as slab pull, are all resolved in the vertical direction, the correlation coefficients of Figures 6a, 6c, and 6e improve to r = −0.54, 0.68, and 0.60 respectively (Figures S6a-S6c). Variations in the strength of the overriding plate are also likely to influence the interface stress. For example, the model with age SP = 65 Ma and age OP = 20 Ma has high interface thickness and low stress, relative to other models with age SP = 65 Ma. This model involves significant overriding plate yielding (Figure 1b). The reduced mechanical coupling between the subducting and overriding plates that follows appears to promote thickening of the interface layer, without the equivalent low bending moment or rollback velocity required in other models (Figures 6b and 6f).
Average interface shear stress was also calculated for the deeper 50-100 km depth range ( Figures S6d-S6f). In this deeper depth range interface stress is higher, ranging from 16 to 45 MPa, and there do not appear to be any correlations with slab dynamics. This poor correlation may relate to the more uniform stress in the deep interface shown for the end-members in Section 3.1 and closer proximity to return mantle flow which couples the plates.

Interface Layer Time-Dependency and Slab-Mantle Interaction
Having demonstrated the instantaneous relationship between interface properties and slab dynamics in Section 3.2, we now briefly explore how these relationships evolve through time in two models, in order to illustrate the time-dependent variability that is possible for subduction zones. Variation in weak layer thickness and shear stress appears to be primarily caused by changes in slab-mantle interaction, which are particularly significant given the nonlinear rheologies utilized. Generally, slab bending moment and pull is at a maximum when young or torn slabs do not extend into the lower-mantle and/or penetrate vertically through the 660 km transition, and at a minimum when the slab is draped over the transition. These characterizations are illustrated by the following examples.
The evolution of the model with age SP = 20 Ma and age OP = 20 Ma (stationary trench end-member) is shown in Figure 7. The pressure gradient is calculated from the interface Poiseuille flow required to account for deviations from simple shear at 10 km depth (using Equation S1), where there is a minimal time delay for volume-flux variation diffusing downdip. All other data are calculated in the same manner as for Figure 6. The slab is relatively weak and prone to breaking, which occurs at 25-30 Myrs (Figure 7a). This slab break-off is associated with a simultaneous episode of large downdip-decreasing pressure gradient within the weak layer, increased slab pull, increased interface stress, and reduction in interface thickness. Both the slab pull stress and the interface stress reduce once the slab extends into the lower-mantle again (Figure 7b), while the interface thickens slowly.
The model with age SP = 40 Ma and age OP = 65 Ma involves a switch from downdip decreasing pressure gradient and relatively thin interface layer (Figure 7c), to downdip increasing pressure gradient and interface layer thickening (Figure 7d). Despite the model involving moderate trench retreat (Figure 6), the "horizontally deflected" slab morphology (identified by Garel et al., 2014) results in reduced slab-lower-mantle stress transmission and subsequently a steep slab dip in the mid-upper-mantle corresponding to reduced slab unbending. Once the slab drapes across the lower-mantle, the pressure gradient reverses and there is a slow transition to decreased slab pull (involving a switch from axial tension to compression) and interface stress, while the interface slowly thickens. This transition is effectively recorded in the hook morphology of the slab sinking through the lower-mantle.

Discussion: Estimates of Interface Stress and Comparison to Seismogenic Behavior
Seismogenic behavior (seismic vs. aseismic deformation) has been hypothesized to be influenced by slab and overriding plate dynamics (Heuret et al., 2012;Schellart & Rawlinson, 2013;Scholz & Campos, 2012) and more generally deviatoric stress (Nishikawa & Ide, 2014;Scholz, 1968). Exploration of a link between these ideas is limited by uncertainty in how interface stress varies. Interface stress has previously been inferred to be relatively homogeneous between margins (England, 2018), anomalously high at aseismic margins (Gao & Wang, 2014) or anomalously high only in regions supporting significant topography (Dielforder et al., 2020;Lamb, 2006). Our results indicate that slab dynamics may exert a considerable influence on interface stress. In this section, we draw on our model results to qualitatively estimate interface shear stress from slab dynamics for four Circum-Pacific regions, comparing to previous characterizations of stress state. We then briefly revisit qualitative correlations between stress state and seismogenic behavior for these Figure 7. Variation in interface dynamics through time for two models with high variability in interface stress through time (as in Figure 6), compared to slab morphology and strength (depicted by viscosity). Squares denote initial trench position. (Left) A period over which the slab is unsupported by the lower-mantle and dominated by downward slab pull (a) is associated with a high contribution of downdip Poiseuille flow, interface thinning, and high interface stress. Once the slab encounters the lower-mantle (b), the pressure gradient and shear stress decrease. (Right) Despite significant rollback, the model initially has a downdip Poiseuille contribution thinning the interface, related to low slab-lower-mantle interaction (c). Once the slab drapes the lower-mantle (d), the pressure gradient switches to downdip pressure increase, after which interface stress and slab pull slowly decrease and interface thickness increases. regions (Figure 9), chosen to include relatively long subduction margins, applicable to our 2-D models, that are known to have contrasting seismicity and subduction dynamics.
Subduction involving limited slab unbending (positive or weakly negative bending moment) and high slab pull at 100 km depth, generally associated with negligible trench retreat and minor slab-lower-mantle interaction, is predicted in the models (Figure 6) to be associated with a thin interface layer and, subsequently, high deviatoric shear stress within the seismogenic zone depth range (Figure 8a). Subduction with significant slab unbending and limited slab pull, promoted by rapid trench retreat and draping of the slab across the lower-mantle, is predicted to promote a thick interface layer and low shear stress (Figure 8b).
Modeled shear stress can vary between the two end-member regimes by a factor of >1.5 at 10-50 km depth ( Figure 6). This variation is larger than the range of background shear stress used by Fang and Dunham (2013) to reproduce a range of earthquake size probability distributions in models of seismic slip on rough faults (compared to a factor of 1.2 background stress variation between their models with a fault roughness amplitude to wavelength ratio of 0.01), where ruptures spanning an entire fault are relatively improbable at low background stress. Subduction zone b-values mostly range from approximately 0.9 to 1.3 (Nishikawa & Ide, 2014), which corresponds to a range in shear stress of 100 MPa order following the compilation of Scholz (2015), much greater than that modeled here. Dal Zilio et al. (2018) used models of seismic slip within collisional orogens to reproduce seismic slip catalogs with b-values ranging from 0.8 to 1.1 for models with orogen differential stress ranging from 150 to 300 MPa, related to varying convergence rate. This is a similar factor of relative stress variation as between our models, though with higher absolute magnitudes.
As our models are designed to explore mantle dynamics rather than exactly reproduce observables such as plate velocities or topography, stresses are not expected to have exact absolute magnitudes or reproduce near-surface variation. It is also difficult to directly relate the modeled stress state effective at million-year time-scales with stress variation within the earthquake cycle. We therefore take the simplified approach of characterizing regions with slab unbending and bending stresses as having relatively low and high shear stress respectively when only the influences of the modeled long-term slab dynamics are considered.
Slab bending and the resulting slab stress state can be inferred from both slab morphology, imaged with seismic tomography, and the focal mechanisms of slab earthquakes. Slab-lower-mantle interaction varies between slabs penetrating the lower-mantle with little change in morphology and slabs that stagnate and drape along the lower-mantle (Goes et al., 2017), the variation relating to the history of trench motion, as occurs in the models of Garel et al. (2014). Slab-lower-mantle interaction in geodynamic models influences the slab stress state, controlling whether the slab at intermediate depths (∼100-300 km) is bending or unbending (Figure 4; Alpert et al., 2010). Varied slab-lower-mantle interaction is indeed reflected in the orientations of principal deviatoric stresses in the slab, inferred from earthquake focal mechanisms (Goes et al., 2017;Isacks & Molnar, 1971). The focal mechanisms of slab earthquakes at "intermediate" depths tend to have slab-dip-orientated compressional axes (P-axes) when slabs are draped along the lower-mantle, as expected for slab unbending, or slab-dip-orientated tensional axes (T-axes) when slabs penetrate into the lower-mantle, as expected for bending. Alpert et al. (2010) compiled the average normalized centroid moment tensor solutions of earthquakes occurring within 50-200 km wide (along-strike) strike-perpendicular slices and 50 km depth bins starting at 100 km depth (isolating intraslab earthquakes from overriding plate and megathrust earthquakes), calculating the alignment between the downdip moment tensor component and a downdip T-axis. We assume that our measurements of slab bending and pull at 100 km depth are applicable to the intermediate depth range analyzed by Alpert et al. (2010), as intermediate depth earthquakes are generally characterized as having similar focal mechanisms, however we note any known variability. It is ambiguous whether earthquakes with slab-dip-orientated P-axes represent unbending or axial compression (bending or axial tension for T-axes), however our modeled slab unbending is characterized both by a strongly negative bending moment and decreased slab pull, so that it is unnecessary to make this distinction.
The average orientation over a depth range of 100-350 km is shown for selected slices (red and blue rectangles, Figure 9). High seismic P-wave velocity anomalies of ≥0.5% (relative to a radial reference model) at 575 km depth were extracted from the UU-P07 model (Amaru, 2007;Hall & Spakman, 2015), as an indication of stagnated slabs at the lower-mantle transition and therefore any recent trench retreat that may have resulted in slab unbending.
Regional seismogenic behavior is typically characterized by instrumentally constrained seismic coupling or records of great earthquakes (e.g., Bilek & Lay, 2018;Pacheco et al., 1993;Uyeda & Kanamori, 1979). We plot seismic coupling coefficients for subduction margin segments (trench colors, Figure 9) calculated by Heuret et al. (2011) as the ratio of the effective seismic slip rate (for earthquakes spanning  to the subduction convergence velocity. Megathrust earthquakes with M w ≥ 8 are also shown (spanning 1900-2018), compiled by van Rijsingen et al. (2018). Seismic coupling is significantly weighted toward large earthquakes, such that low seismic coupling reflects a low number of recorded large earthquakes. There are significant uncertainties in such characterization of seismogenic behavior, due to a limited instrumental record (McCaffrey, 2007) and interpretations of aseismic behavior should be considered hypotheses. Regimes of compressional, extensional, or neutral (either no deformation or strike-slip) overriding plate strain regimes characterized by Heuret et al. (2011);Heuret et al. (2012) are also plotted (labeled as C, E, and N in Figure 9).
Comparison of estimated interface stress to previous characterizations of seismogenic behavior is a test of whether geodynamic influences at the length-and time-scales studied here may influence seismogenesis. This would be superimposed on many other influences at a range of scales. Fault-zone scale heterogeneity of frictional properties, sediment thickness, or pore pressure has been associated with variability in seismogenic behavior (Pacheco et al., 1993;Scholl et al., 2015;Wallace et al., 2012). Bletery et al. (2016) argue that large earthquakes occur on flat megathrust segments (low curvature), possibly contradicting our inference that decreased slab bending is associated with low megathrust deviatoric stress. However, there is a complex dependence on what controls shallow slab geometry. For example, subduction of buoyant material promotes slab flattening, while also reducing unbending as the slab enters the asthenosphere (Van Hunen et al., 2002). Brizzi et al. (2020) showed that an increase in incoming sediment thickness in a subduction zone with imposed convergence velocity increases sediment accretion, trench retreat, and megathrust length, promoting large earthquakes. Increasing megathrust length may therefore offset decreasing shear stress during trench retreat. Margins experiencing rollback (e.g., Tonga, Ryukyu) are commonly erosive margins, rather than accretionary (Clift & Vannucchi, 2004), indicating that a complex interplay of dynamics may influence accretion.
BEALL ET AL.

Slab Bending End-Members: Peru, Chile, and Sumatra
The focal mechanisms of slab earthquakes at the margins of Peru, Chile, and Sumatra (Figures 9a and  9b) are clearly indicative of slab dip-parallel tension (consistent alignment of T-axes with slab dip; Alpert et al., 2010), consistent with slab bending or stretching, and are therefore predicted to have high interface shear stress. The high interface shear stress predicted for the Chile margin is broadly consistent with the findings of Lamb (2006) and Dielforder et al. (2020), who both calculated that high integrated interface stress supports the Andes, particularly in northern Chile. These studies predicted shear stress in southern Chile to be approximately half of that in northern Chile, though still moderate (up to approximately 50 and 100 MPa respectively at 40-50 km depth). Heuret et al. (2012) characterized the back-arc stress regime as compressive in Peru and northern Chile and neutral in Southern Chile (Figure 9a). Their neutral regime represents a combination of trench-perpendicular shortening and strike-slip, of which the latter component cannot be reproduced in 2-D models, however the shortening component is likely to be most relevant to interface stress.
The back-arc region of Sumatra is characterized as neutral (Figure 9a; Heuret et al., 2012) and is undergoing transpression, in order to accommodate oblique subduction (Bellier & Sébrier, 1995). The shear stress in Sumatra has previously been estimated to be of similar magnitude to southern Chile and therefore relatively moderate, but lower than northern Chile, corresponding to a thinner overriding plate with minor topography (Dielforder et al., 2020;Lamb, 2006). There was no obvious relationship between shear stress and overriding plate age, and therefore thickness, in our models ( Figure 6), though thick overriding plates have previously been found to encourage plate coupling (Sharples et al., 2014).
The slabs at Peru and Sumatra each penetrate or do not reach the lower-mantle transition (limited slab extent at 575 km depth in Figure 9), consistent with the high slab bending observed. The Chile slab is an exception, which flattens at 660 km while the adjacent Peru slab penetrates, consistent with higher relative trench displacement (Goes et al., 2017). This mismatch between slab morphology and bending stress may be a result of the narrow along-strike extent of draping slab or other plate coupling processes. We also note that much subduction in Peru generally occurs with a shallow dip (flat subduction; Gutscher et al., 2000), with related unbending at ∼50 km depth, before transitioning back to bending at ≥ 100 km depth (Sandiford et al., 2019), likely complicating the interface force balance.
Seismic coupling in Chile, Sumatra, and southern Peru is generally high and each of these segments have experienced many great earthquakes (Figures 9a and 9b), implying a possible correlation with the relatively high shear stresses estimated. The seismic coupling in northern Peru, however, appears to be low and there is no instrumental record of great earthquakes in the last 100 years.

Slab Unbending End-Members: Ryukyu and Kermadec-Tonga
Slab dip-parallel compression, and therefore slab unbending, is most clearly evident for the Ryukyu and Kermadec-Tonga margins (Figures 9c and 9d; Alpert et al., 2010). The Kermadec-Tonga and Ryukyu slabs are both draping across the lower-mantle (Pownall et al., 2017;van de Lagemaat et al., 2018), consistent with trench retreat indicated by plate reconstructions and back-arc basin opening (Jolivet et al., 1994;Schellart et al., 2006). The Ryukyu slab has a tear that is propagating from the south-west and does not yet appear to have completely separated the slab (Pownall et al., 2017), which may explain why the slab is still unbending. Slab unbending and rollback at the Ryukyu and Kermadec-Tonga margins is estimated to correspond to low interface shear stress.
BEALL ET AL.

10.1029/2020GC009267
15 of 20 Figure 9. Comparisons of average slab earthquake focal mechanisms (rectangles, inferring slab bending stress, red, and slab unbending, blue; Alpert et al., 2010), seismic coupling (trench coloring; Heuret et al., 2011), and draping slabs at the lower-mantle transition indicated by seismic P-wave anomalies (Amaru, 2007;Hall & Spakman, 2015). M w ≥ 8 megathrust earthquakes (1900-2018van Rijsingen et al., 2018) and trench boundaries (Bird, 2003) are also shown. Overriding plate regimes are labeled as compressional (C), extensional (E), or neutral (N) following Heuret et al. (2011), Heuret et al. (2012. Slabs that are clearly bending or stretching appear to be commonly associated with high seismic coupling (a and b; Sumatra and Chile, though only southern Peru) and slabs that are unbending are generally associated with low seismic coupling (c and d; Ryukyu and Kermadec-Tonga), while the bending style and relationship with seismicity is less clear for other margin segments.
The Ryukyu and Kermadec-Tonga margins were classified by Heuret et al. (2012) as having extensional overriding plate regimes and to be aseismic. Scholz and Campos (2012) associated rapid slab rollback at the Kermadec and southern Tonga margins with low seismic coupling. Lamb (2006) calculated that the shear stress at the Tonga margin is slightly lower or higher (within a few MPa) than southern Chile or Sumatra, depending on whether or not its low trench fill results in an unusually high friction coefficient. This differs from our estimate of anomalously low shear stress corresponding to what appears to be an ideal example of slab unbending. Our models with the most prominent slab rollback had old slabs, resulting in higher slab density, faster convergence, and lower subduction temperatures, which are all likely to promote high shear stress in a lithosphere scale force balance. However, these effects are offset by increased in-plane slab compression ( Figure 4c) and the resistance of the asthenosphere during rollback (Schellart, 2004), likely contributing to the low net shear stress. The force balance of Lamb (2006) does not include these mantle forces, which may explain the discrepancy.
Both margins are generally considered to be aseismic on account of having low seismic coupling and a paucity of instrumentally recorded great earthquakes (Bilek & Lay, 2018;Pacheco et al., 1993;Scholz & Campos, 2012), which may then be correlated with the low estimated shear stress. It has been argued that northern Tonga has moderate seismic coupling related to margin curvature near the slab edge (Scholz & Campos, 1995), the interface stress state of which could be explored in future 3-D models of slabs with finite widths.

Japan-Izu-Bonin-Mariana Trenches
The interface stress states inferred on the basis of slab bending (Figures 8 and 9c) disagree with previous estimates for the Japan-Izu-Bonin-Mariana margins. The slab appears to be unbending at 100-300 km depth at the Japan trench, conflicting with previous inference of a compressive overriding plate (Heuret et al., 2012) and high integrated interface shear stress (Dielforder et al., 2020). Slab bending appears to be occurring at the Izu-Bonin-Mariana margins, conflicting with previous estimates of overriding plate extensional regime behind the Mariana margin (Heuret et al., 2012) and low calculated interface shear stress at the Izu-Bonin margin (Lamb, 2006). The difficulty in using the models to predict the interface stress at these margins may be related to a reorganization of subduction dynamics ∼5 Ma (Hall, 2002), which potentially resulted in more complicated bending dynamics than were modeled.
The characterizations of average intermediate-depth slab stress state by Alpert et al. (2010) appear to agree with slab morphologies in both regions. The unbending slab at the Japan trench is draping across the lower mantle ( Figure 9c; Goes et al., 2017). The bending slab at the Mariana trench is penetrating into the lower-mantle, with little evidence of recent rollback (Fukao & Obayashi, 2013;Goes et al., 2017). However, in both cases, there is evidence that the overriding plate stress state is in disequilibrium with or is decoupled from mantle-scale slab dynamics. At the Japan trench, slab earthquakes at 100-150 km depth have focal mechanisms with downdip T-axes (Alpert et al., 2010) potentially indicative of bending or stretching at the 100 km depth we analyzed, from which our model would instead infer high interface stress. Yang et al. (2018) reproduced this pattern in geodynamic models, coinciding with the initial stalling of Japan trench retreat. At the Mariana trench, back-arc spreading is still occurring, controlled by motion of the overriding plate (Heuret & Lallemand, 2005) and interaction with the Ryukyu subduction zone, lowering the geodynamic coupling between the overriding and subducting plate in geodynamic models of the region (associated with low pressure within the plate interface; Holt et al., 2017;Holt et al., 2018).
It is therefore plausible that the intermediate-depth slab bending dynamics are not a dominant control on the interface stress state at the Japan and Izu-Bonin-Mariana margins and previous inferences of low interface stress may result from overriding plate dynamics not modeled here. The Japan trench has hosted many great earthquakes and appears to have high seismic coupling (Figure 9c), which has been argued to correlate with overriding plate compression (Heuret et al., 2012). Overriding plate extension has likewise been associated with the lack of great earthquakes and low seismic coupling at the Mariana trench. Future modeling incorporating overriding plate forcing is required to test any correlations between geodynamic influences on interface stress and seismogenic behavior.

Along-Strike Variation in Slab Dynamics
While we have focused on long margin segments most applicable to our 2-D analysis, there is potential for future studies to explore the extent to which along-strike variation in subduction dynamics results in contrasts in interface stress. For example, unlike the Ryukyu slab, the slab subducting at the Nankai trench does not extend below 100 km in parts (Pownall et al., 2017;Wu et al., 2016) and experiences bending or stretching at 50-100 km depth (Bailey et al., 2012). There is also a contrast between slab bending at the Kermadec and Hikurangi margins (Figure 9d), potentially corresponding to a reduction in recent slab rollback . In both cases, there is evidence of an along-strike increase in interseismic coupling in regions of increased slab bending, which have been associated with a switch in overriding plate deformation style (Wallace & Beavan, 2010;Wallace et al., 2009).
At the Sumatra-Java margins, it is also unclear to what degree the inferred slab bending and high interface stress at Sumatra extends to Java. Heuret et al. (2012) classified both Sumatra and Java as having similar overriding plate regimes, however they found that Java has relatively thinner trench sediment fill that may discourage the occurrence of great earthquakes. Raghuram et al. (2018) used forebulge flexure models to argue that high plate coupling occurs in Sumatra, where the slab is short, due to horizontal forcing resulting from the adjacent subduction of the longer Java slab (P-wave anomalies in Figure 9a). This may be analogous to imposing a constant v c in our models with young age SP and episodic stalling. While Java appears to have hosted fewer great earthquakes and to have lower seismic coupling than Sumatra, the seismogenic behavior remains poorly constrained (Scholz & Campos, 2012). Our 2-D models are unable to capture alongstrike variations in slab dynamics, however there is potential for slab dynamics to influence transitions in interface stress state in the regions discussed, requiring future exploration using 3-D models. There is also potential to use 3-D models to explore the stress state of narrow margins with high trench curvature, such as the Sandwich and Antilles margins, which appear to have low seismic coupling (Heuret et al., 2012).

Conclusion
The relationship between subduction dynamics and interface shear stress has been analyzed for data compiled from numerical thermo-mechanical subduction models by Garel et al. (2014). We find that the interface layer shear stress is broadly proportional to slab bending and slab pull, and inversely proportional to trench retreat velocity. It follows that higher shear stress is predicted at seismogenic depths for slabs that are dominated by vertical slab pull and limited unbending, while lower shear stress is predicted for slabs with significant unbending and reduced slab pull, typically associated with trench retreat. These relationships were used to qualitatively estimate relative interface stress using previous characterizations of slab bending for Circum-Pacific regions and compared to previous estimates of stress and seismogenic behavior.
At the Tonga-Kermadec, Sumatra, southern Peru, Chile, and Ryukyu margins, where axial slab stress state is most clearly characterized, estimates of high or low interface stress associated with bending or unbending broadly agree with previous estimates and correlate broadly with aseismic or seismic characterizations respectively (with the exception of northern Peru). Other margins analyzed are anticorrelated or ambiguous. The predictions of interface stress from slab bending style is difficult to reconcile with previous stress characterizations at the Japan-Izu-Bonin-Marianas margins, which we interpret to be related to recent reorganization of subduction dynamics and overriding plate forcing not modeled. Our findings indicate that there is potential for large-scale subduction dynamics to produce considerable variation in interface stress between margins thought to have contrasting seismogenic behavior. Further modeling is required to quantify how this long length-and time-scale stress influence interacts with the earthquake cycle, as well as other critical influences on fault stress, such as fluid pressure, trench sediment thickness and heterogeneity of frictional fault-zone properties.