Imbalanced Moment Release Within Subducting Plates During Initial Bending and Unbending

Internal deformation within the downgoing plate in subduction zones to accommodate the bending of the plate as it starts to subduct is reflected in widespread intraplate seismicity. This seismicity, extending from the outer rise and outer trench slope down to intermediate depths within the slab, exhibits a wide variety of focal mechanism types, often indicative of the accumulation and recovery of down‐dip curvature through the combination of normal‐faulting and thrust‐faulting earthquakes. In the idealized case, where the bulk internal deformation is recovered and slabs descend as a straight plate into the deeper mantle, we might expect the seismic moment released in both extension and compression to balance. However, a number of factors may complicate this: the thermal, compositional, and rheological evolution of the slab as it subducts, changes in the proportion of deformation accommodated seismically, and whether the slab undergoes any permanent deformation. These factors may result in intraslab deformation and seismicity unrelated to the slab curvature. Here, we assess earthquake moment release in intraslab settings around the world, focusing on those subduction systems with relatively simple slab geometries. Consideration of several regions around the western Pacific and eastern Indian Ocean indicates that substantially more deformation is accommodated seismically during bending than during unbending, and that in both settings, significantly more moment release reflects down‐dip extension than down‐dip compression. This suggests that, although the location of seismicity is clearly related to changes in slab curvature, there is a component of permanent, unrecovered down‐dip extension in many subducting slabs.

The implications of intraslab seismicity remain unclear. There are longstanding questions regarding the mechanisms that permit seismicity to occur at the depths and pressures seen within slabs (Hacker et al., 2003;Isacks & Molnar, 1971;Kirby, 1987). But there are also questions regarding the deformation modes that these earthquakes represent. While smaller intraslab earthquakes may be related directly to the release and passage of free fluids (e.g., Halpaap et al., 2019), larger-magnitude seismicity is expected be the result of the plate-scale stress field in the source region, and therefore provides a vital insight in to the orientation the stress field and the geodynamics of the slab, even if the rheological conditions allowing brittle deformation are related to local mineralogy and the occurrence of metamorphic transitions within the slab (e.g., Hacker et al., 2003;Isacks & Molnar, 1971;Peacock, 2001).
Initial work analyzing the intraslab stress field through intermediate-depth earthquake focal mechanisms (e.g., Alpert et al., 2012;Isacks & Molnar, 1969, 1971 suggested that intraslab deformation was dominated by the influence of axial plate stresses (i.e., slab pull; ridge push; tractions on the edges of slabs, and lower mantle resistance). In contrast, more recent work, benefiting from better resolution in earthquake locations, suggests that changes in curvature (bending/flexure/slab buckling) provide a first-order explanation for the location and orientation of intraslab seismicity in many global subduction settings Myhill, 2013;Sandiford et al., 2019Sandiford et al., , 2020Willemann & Davies, 1982), at least at intermediate depths (deep-focus earthquakes may be more complex). This particularly applies to slabs at shallow depth within the mantle, before interactions with the mid-mantle transition zone at ∼660 km become important. While increasing evidence from the polarization of double seismic zones (e.g., Bloch et al., 2018;Igarashi et al., 2001;Kita et al., 2010) and the correspondence between slab seismicity and slab curvature (e.g., Myhill, 2013;Sandiford et al., 2020) suggest that stresses relating to the bending and buckling of the slab dominate in driving intraslab seismicity, the degree to which this interacts and overprints the regional stress (i.e., the net component of the in-plane stress within the plate) field remains uncertain. There is also an important distinction between elastic bending of the plate (where mechanical work is stored as elastic energy with little seismic expression), and the plate's eventual failure by faulting and/ or buckling. In the latter case, mechanical work is dissipated through the generation of heat, seismic energy, or the latent heat required to drive changes in plate structure (e.g., grain size reduction). Only in this latter case can intraslab seismicity occur, where stresses are too great to be supported elastically.
Here, we test a simple hypothesis, illustrated in Figure 1: in regions where slab morphology is relatively simple during the initial stages of subduction, with the development and recovery of a single predominant down-dip bend, does the seismic moment released in the development of the bend (increase in curvature) match that released during the unbending process (decrease in curvature)? In assessing this, we consider intraplate seismicity beneath the outer rise and in the near-trench environment and intraslab seismicity to be one family of events-all potentially relating to the development and recovery of slab curvature.
Our hypothesis, and the discussion section that follows, relies on the assumption that intra-plate seismicity (reflecting the brittle failure of the plate in discrete, seismogenic events) is representative of the stress state at the location at which the earthquake occurs, which in turn reflects the accumulation of strain. We further assume that the majority of the stress within the plate is supported by the sections which deform through brittle failure (e.g., Forsyth, 1980;Jackson et al., 2008;Wiens & Stein, 1983), and that sections of the plate that are instead deforming through other processes, particularly viscous creep, support only a small proportion of the stresses applied to the plate. This assumption allows us to disregard the contribution that the viscous part of the plate will make to the distribution of bending stresses with depth into the slab, and to focus on the seismogenic deformation. We assume that the majority of strain within the deforming volumes is accommodated during earthquakes and that aseismic slip on faults is low. Given the typically lower activation stresses for aseismic creep compared to the yield stresses involved in seismic slip, were significant aseismic slip occurring, we would expect deformation within the volume to concentrate into these regions, and largely deactivate other potentially seismogenic fault planes. That we still see widespread intraslab seismicity suggests this is not occurring, and that our assumption is reasonable. However, we are aware that other readers may find these assumptions simplistic. Our analysis assumes that the bending strain is dominated by the contribution of the advective component; intuitively this amounts to the curvature changes produced as the slab moves over a fixed template corresponding to its shape (e.g., Kawakatsu, 1986;Sandiford et al., 2020). In practice, this means we assess the seismic moment release in terms of the estimated down-dip curvature gradient of the slab at the source location, as described in further detail in Section 4. Lastly, we assume that the seismogenic cross section of the slab does not narrow significantly over Figure 1. Schematic representation of initial intraplate deformation and seismicity. (a) Shows the internal deformation field in the slab, highlighting regions of down-dip strain rate related to slab bending, with blue indicating down-dip extension, and red down-dip compression. Strain with the viscously deforming region (shaded gray) is not shown. Dashed line indicate the neutral fiber between bending regions. Dotted line indicates the point of maximum down-dip curvature. Mechanisms are shown in cross-section, and indicate the style of seismicity expected to accompany each region of strain rate (precise orientations will vary). Inset panels show two conceptual rheological models (under our assumptions) for the stress state in regions in increasing curvature, with constant (left) and depth-dependent (right) yield stresses. (b) It is as in (a), but flattened into a reference frame relative to the top of the plate. (c) Shows schematic curves for the accumulation of slab dip (black), curvature (gray, solid), and the down-dip rate of change in curvature (gray, dashed).

Seismicity Catalogs
Most subduction zones are host to prolific seismic activity associated with various elements of the subduction process (seismicity on the main subduction interface, within the downgoing plate, and within the overriding plate). To isolate a subset of the earthquake catalog that is associated only with changes in the downdip curvature of the descending plate, we require accurate information on both the location of the earthquake and the style of deformation that the earthquake represents (exemplified in the moment tensor). In this study, we draw on two seismicity catalogs, each with different strengths. We rely on the gCMT catalog (Dziewonski et al., 1981;Ekström et al., 2012) for information regarding earthquake moment tensors and earthquake magnitudes. However, the relatively long period seismic data (>40 s, for both body-and surface-waves) upon which the gCMT catalog is based precludes the precise determination of earthquake locations, leading to relatively large uncertainties, particularly in source depth, particularly for smaller earthquakes. Instead, we draw on the ISC-EHB catalog (Engdahl et al., 2020) for locations. This catalog, derived from the inversion of phase arrival times for a multitude of phases, offers the most comprehensive and accurate routinely calculated global catalog of earthquake locations available. In combining the two catalogs, we are limited in both the magnitude range we can use, typically determined by the completeness of the gCMT catalog, and the time duration to the catalogs. Here, we use catalogs extending from 1970 to the end of 2016, the current end date for the ISC-EHB.
To combine these two catalogs, we first take regional subsets of both catalogs, 1° greater in every direction than a given study area, allowing a buffer zone around the edge of our study area such that we avoid the problem of losing events due to small mismatches between catalog locations around the periphery of a regional study area. The origin time of each catalog entry is converted to decimal years. We then test every event in the gCMT catalog against entries in the ISC-EHB catalog, looking for entries that match within given tolerance in terms of origin time (<3 × 10 −7 yrs, or about 10 s), inter-catalog event separation (<50 km), and magnitude (within two magnitude units). Tolerance thresholds are adapted slightly for regional variability, particularly in the location accuracy of larger events in the gCMT catalog.
To pass through to our final catalog, we require an entry in the gCMT catalog to be associated with a single entry in the ISC-EHB catalog. In cases where two matches exist in the ISC-EHB catalog to one gCMT entry, we discard the result-this limitation particularly affects events within aftershock sequences, especially after major interface seismicity (e.g., near Honshu shortly after the 2011 M w 9.0 Tohoku-Oki earthquake). Discarded results, both those with no matches, and those with multiple potential matches, are checked manually to ensure that no major (large-magnitude) earthquakes are excluded. In general, this approach usually yields matches for >95% of earthquakes in the gCMT catalog.
For each correlated event, we then take the source time and hypocentral location from the ISC-EHB catalog, and the moment tensor and moment magnitude from the gCMT catalog, and use this combined catalog for the rest of the study.

Region Selection and Structure
We focus our study on 13 regions, selected on the basis of a relatively simple slab morphology characterized by a single major bend, followed by a relatively straight section of slab at intermediate depth, based on the slab models of Hayes et al. (2012Hayes et al. ( , 2018; hereafter referred to as "Slab1" and "Slab2," respectively). These 13 regions include 6 regions from the NW Pacific (Aleutians, Kuril-Kamchatka, Honshu, Ryukyu, Bonin, and Marianas), 3 from the SW Pacific (The Solomon Islands, New Britain, and the New Hebrides), 2 sections of the Tonga-Kermadec-Hikurangi subduction system (Tonga-Kermadec and Kermadec-Hikurangi), and 2 from Indonesia and the eastern Indian Ocean (Sumatra and Java-Sumba). In Figures 2-5 we show processing results from four example regions (Honshu, Kuril-Kamchatka, Tonga-Fiji-Kermadec, and New Britain), and include processing results from the other nine regions in Supporting Information (Figures S2-S10 in Supporting Information S1). Seismicity catalogs considered in each of these 13 regions are provides in Supporting Information.
In assessing the relationship between seismicity and slab geometry, we rely on the Slab1 and Slab2 models for the surface of the subducted plate. We note the slight circularity in using models for the slab geometry that are 5 of 18 partially derived using earthquake locations (particularly for Slab1) to interpret the geometrical context of the same seismicity. As Slab2 also uses numerous other constraints (e.g., seismic tomography, reflection, and refraction data), and as we are interested in the changes in slab morphology, rather than its actual location in space, we consider this circularity to be only a minor concern with this model, and prefer this more recent compilation of slab models where possible. In a limited number of cases, where the regional slab geometry from Slab2 clearly deviates from the location of nearby seismicity (Tonga-Fiji-Kermadec and Bonin), we revert to using the older Slab1 model. White-green contours show slab depth from Slab2, and our digitized trench line (dashed). Yellow lines show the along-strike extent of the the study area. Circles show earthquakes in our combined catalog, scaled according to earthquake magnitude. Circle outlines are shaded by slab depth. Blue points are down-dip extensional earthquakes updip of the first zero-crossing in curvature. Red points are down-dip compressional earthquakes updip of the first zero-crossing in curvature. Gray earthquakes are those excluded from our catalog, as discussed in Section 2.4. (b) Slab cross section. Earthquakes are shown by the orientation of the T-axis of the focal mechanism, shaded as in (a). The slab model is shown as a histogram, based on the discrete profiles used in our reprojection scheme, spaced at 10 km. (c) Slab-relative cross section. Earthquakes are plotted as function of down-dip distance and depthinto-slab. T-axes are rotated to be slab relative. (d) Histogram of slab dip as a function of distance. Red line shows the mean slab dip. The yellow line shows mean slab dip, converted into distance along the slab surface. (e) As in (d), but for slab curvature. Purple lines show the maximum in mean curvature in slab-relative distance (dashed line), and the first zero crossing in mean curvature in slab relative distance (solid line). (f) It is as in (d), but for rate-of-change of slab curvature. (g) Combined mechanism from the summation of all down-dip extensional faulting updip of the first zero-crossing in slab curvature, and a stereonet showing the distribution of P/N/T axes for those earthquakes. Both are expressed in a slab-relative coordinate system. (h) It is as in (g), but for all down-dip compressional earthquakes. 6 of 18

Earthquake Reprojection
To isolate a set of earthquakes limited to intraslab deformation, and to interpret these in the context of the local slab geometry, we need to reproject our earthquake catalog into a slab-relative reference frame ( Figure 1b).
For each region, we draw on the plate boundaries of Bird (2003), or, where these visibility deviate from the bathymetric trench, our own determination of the local trench line from available bathymetric data. In relating earthquakes to changes in slab geometry, we reproject our earthquake dataset into a slab-relative reference frame, and for the determination of metrics relating to slab geometry involved in this (dip, curvature, and rate-of-change in curvature), we again rely on the chosen regional slab model.
We start by merging the relevant slab model for each region with a flat bathymetry at its up-dip extent, to extend it out onto the oceanic plate. Following the approach of Sandiford et al. (2020), we then determine a set of trench-perpendicular azimuths at 20 km intervals along the plate boundary, and extract the slab geometry along each trench-perpendicular profile. The study area and seismicity catalog, is then limited by the two profile lines at either end of the selected region (yellow lines on Figures 2-5 and S2-S11 in Supporting Information S1). For Figures 2-5 show four example regional study areas, for Honshu, the Kuril-Kamchatka Arc, Tonga-Fiji-Kermadec, and New Britain, and cross sections through the slab model and earthquake dataset, with the horizontal axis being distance from the trench. For every earthquake in our combined earthquake catalog, we associate the earthquake with a profile line based on the closest profile to the earthquake location, when both are projected to the Earth's surface. We then calculate the earthquake location in terms of down-dip distance along the slab surface profile and perpendicular distance from the closest approach to the slab surface (i.e., distance along the slab surface and depth into the slab; shown in Figures 2-5c).

Isolating Intraslab Seismicity
Finally, we apply a set of filters aimed at isolating the subset of seismicity which is associated with the down-dip deformation of the slab. In order to be retained, earthquakes are required to pass all of the following filtering stages: 1. Earthquakes are required to have a minimal double couple component to their moment tensor. Here, we limit this by requiring the percentage double couple (γ) to be >70%. We follow Jackson et al. (2002) in defining this as: where λ i are the eigenvalues of the moment tensor. 10.1029/2021JB023658 9 of 18 2. We exclude all earthquakes where the depth from the ISC-EHB catalog is >15 km above the local depth of the slab surface from Slab2 (or Slab1 for the Tonga-Fiji-Kermadec and Bonin regions). For completeness in removing shallow earthquakes within the over-riding plate, we also remove earthquakes landward from the trench at distances where no slab model is present. 3. To remove earthquakes associated with motion on the subduction interface, we follow Sandiford et al. (2020) in determining a similarity coefficient (χ) between the moment tensor of each earthquake, and a predicted subduction interface moment tensor based on the orientation of the local slab surface. We define this similarity as: where ‖M‖ is the norm of the moment tensor and : is the tensor double dot product. Here, we predict the interface moment tensor based on the local slab geometry at the location of each earthquake, and assume that interface events are pure dip-slip compressional deformation. Earthquakes with a χ value above a threshold value are deemed too similar to the expected interface deformation, and excluded. The threshold value of χ is determined manually for each region, based on the χ/frequency distribution, where the interface is marked by the onset of a rapid increase in the number of earthquakes per χ increment as χ approaches 1. We tested an approach where the interface moment tensor was calculated using the slab model geometry and a rake value based on the regional plate motions (allowing a non-dip-slip component to the moment tensor), however this was found to produce less clear χ distributions in cases where subduction is oblique, due to the common partitioning of slip between trench-normal deformation on the interface, and trench-parallel deformation within the forearc (McCaffrey et al., 2000). 4. To remove earthquakes associated with strike-slip deformation within the slab, we exclude all earthquakes with a null axis within 45° of the direction normal to the slab surface at its closest location. 5. To remove earthquakes predominantly representing along-strike deformation within the slab, we exclude earthquakes with a null axis within 35° of the down-dip direction on the local slab surface. 6. Earthquakes are required to be located up-dip of the first point where the slab curvature on their associated slab profile returns to zero, down-dip of the trench. This step is aimed as isolating our dataset from complexities in the slab geometry beyond the first initial bend and the recovery of this curvature. 7. In certain cases, most notably for Honshu during the aftershock sequence of the 2011 Tohoku-Oki earthquake, we add an additional filter, designed to remove earthquakes occurring during the horizontal extension of the forearc region of the over-riding plate. We exclude earthquakes arc-wards of the trench, between 12.5 and 60 km depth, with depths shallower than the local slab depth, and with a T-axis, that is, closer to vertical than the P-axis Figures 2-5a, 5c shows the impact of applying all of these filters, with all earthquakes excluded by these steps shown in gray, and those retained shown in blue (downdip extension) and red (downdip compression).

Regional Examples
Figures 2-5 show four regional examples of the processing approach described, for northern Honshu, the Kuril-Kamchatka Arc, Tonga-Fiji-Kermadec, and New Britain, respectively. Figure 2 shows perhaps the most straightforward example of an initial bend developed through the outer rise and outer trench slope, with pervasive downdip extensional faulting at shallow depth to 30-40 km arc-wards of the trench, which matches the peak in curvature of the slab. This is then followed by a rapid transition to down-dip compressional faulting at shallow depths within the downgoing slab as the curvature is recovered, with some limited down-dip extensional seismicity at greater depth on the outside of the unbend. Moderate magnitude seismicity is sparse in the depth range ∼200 to 400 km, where the curvature change appears to be minimal. Even in this simple case, however, it is important to note that although the slab geometry and dip are relatively consistent across the study region (Figures 2b and 2c), as we consider curvature and rate-of-change in curvature, increasingly complex variability emerges, emphasizing the need to consider the localized slab geometry along each profile at the location of each earthquake. Figure 3 also shows a more complex case for the Kuril-Kamchatka Arc. As demonstrated in Figures 3a and 3b, the overall slab geometry here is also comparatively simple, dominated in all cases by a single major bend, before straightening out at ∼200 km depth. However, increased variability in the along-strike geometry of the slab (Figures 3d-3f) again emphasizes the need to consider each earthquake in its local geometrical context. Once again, we see the initial development of plate curvature expressed through shallow down-dip extensional seismicity, starting seawards of the trench beneath the outer rise. This is accompanied by limited deeper compressional seismicity on the inside of the bend beneath the outer rise (Craig, Copley, & Jackson, 2014). Around the section of peak curvature, expected to mark the transition from predominant bending to predominant unbending, we see the sudden cessation in down-dip extensional seismicity and the onset of widespread down-dip compressional seismicity, accompanied by minor down-dip extensional faulting at greater depth into the slab (Figure 3c). During unbending, as the plate returns to zero-curvature, seismicity is separated into shallow down-dip compression and deeper down-dip extensional zones, as is typical of western Pacific margin slabs (Sandiford et al., 2020).
Matching detailed studies of the depth extent of seismicity in outer rise regions around the world (Chapple & Forsyth, 1979;Craig, Copley, & Jackson, 2014), we expect the complete separation between extensional and compressional seismicity, and, in the absence of any exterior variation in the stress field, a consistent depth of the separation between the two. The slight overlap seen on Figure 3c suggests that either the resolution in earthquake locations available from the ISC-EHB catalog, the resolution of the slab model, or a combination of the two, slightly obscures the finer details of the transition between deformational regimes with depth into the slab.
In Figure 4 we show another example, which highlights a number of remaining limitations. This region-the northern end of the Tonga-Fiji-Kermadec subduction system-is one of the most active intraslab settings in the world, with high productivity rates of both intermediate depth and deep-focus seismicity. The first limitation is the reliance on existing slab models. In this case, when using Slab2 ( Figure S1 in Supporting Information S1) there is a clear spatial deviation between the intraslab seismicity and the slab surface, with the slab model consistently underestimating slab dip (or the seismicity being consistently mislocated), leading to a progressive increase in depth-within-slab with distance along slab. In the example shown in Figure 4, we instead use Slab1 (Hayes et al., 2012) as the slab model. In this case, we make this choice based on a clear divergence between the location of intraslab seismicity (from the ISC-EHB catalog), and the slab location in Slab2 (shown in Figure S1 in Supporting Information S1). Comparing Figures 4 and S1 in Supporting Information S1, the overall slab shape, the variation in slab curvature, and the location of inflection points in slab curvature are broadly similar, but while the slab depth (integrated dip) in the case of Slab2 is consistently less than would be predicted by the location of seismicity, Slab1 (more directly constrained by seismicity) produces a slab top that is much more consistent with the location of intraslab seismicity, as highlighted by the two slab relative cross sections. A similar deviation is only seen for two other regions (Bonin and Kermadec), when Slab2 similarly under-predicts the depth of the slab. In these cases, we again revert to Slab1. The second limitation, connected to the first, is the reliance on slab models in limiting the down-dip extent of the seismicity considered. In both Figures 3a and 4a, the summation includes small groups of events that are substantially further down-dip than the cut-off limit imposed for most of the subduction zone (at 150°E, 50.5°N and −178°E, −22 to −20°N, respectively). This is again a result of method described above, where the depth range over which the slab model indicates unbending (i.e., where the downdip curvature remains positive, prior to returning to zero) can vary substantially. This is in keeping with the approach described above, and technically correct, but again highlights a potential problem in cases where the slab geometry is poorly defined or complex. This potential uncertainty is notably absent for the simple slab geometry under Honshu (Figure 2). In contrast, the increasingly complex geometries of the Kuril-Kamchatka (Figures 3d-3f) and Tonga-Kermadec (Figures 4d-4f) cause such issues to arise, although the moment contribution from such earthquakes is comparatively small.
Despite these issues, Figure 4 shows a broadly similar pattern of seismicity to Figure 3, with the clear separation of shallow down-dip extension and deeper down-dip compression in the outer rise, switching to shallow downdip compression and deeper down-dip extension at intermediate depths. While absolute separation is again not imaged using our combined earthquake catalog, more detailed studies of the outer rise region along Tonga-Kermadec have shown that this is the case seaward of the trench (Craig, Copley, & Jackson, 2014;Lay et al., 2013;Todd & Lay, 2013), and we see no reason why it would not persist at greater depths. We suspect that the lack of an absolute separation between down-dip extension and down-dip compression in this region is a conflation of increased uncertainties in the earthquake locations, and increased uncertainty in the slab models, both stemming from the paucity of near-field data to constrain earthquake locations and the regional velocity structure.
The final example we show (New Britain; Figure 5) illustrates the case of comparatively simple slab geometry, with a single dominant bend, and a pattern of seismicity illuminating the full deformation distribution. Shallow down-dip extension and deeper down-dip compression beneath the outer rise give way to a complex region in which no clear spatial separation between mechanism types is seen, before separating into shallow down-dip compression and deeper down-dip extension as the slab starts to unbend.
For all four regional examples, we calculate composite moment tensors by summing all down-dip extensional (g) and down-dip compressional (h) earthquakes, after rotation into a slab-relative reference frame, along with stereonet plots of the P, N, and T axes for the relevant earthquake selections. In each of the three cases shown, the composite mechanisms for down-dip extensional and down-dip compressional seismicity have similar orientations, supportive of the concept that faults initiated during the initial development of the bend are reactivated at intermediate depth with the opposite sense of motion (Chen et al., 2004;Ranero et al., 2005).
While we show only 4 of our 13 regional studies in this manuscript, we include similar plots for the other regional studies in Supporting Information (Figures S2-S10 in Supporting Information S1).

Moment Summation
Our intention in this study is to consider the deformation accommodated seismically during the initial flexural cycle (bending and unbending) as the downgoing plate enters the subduction zone. We have therefore limited the seismicity catalog considered to only earthquakes located up-dip of the first zero-crossing in the curvature of each slab profile, in effect, the point where the slab geometry first returns to being "straight" (shown by the solid purple lines of Figures 2-5c, 5e). The total flexural deformation between the oceanward extent of deformation and this point should sum to zero (when considering both strain accommodated both seismically and aseismically). We then further subdivide the seismicity catalog, using two methods to assign earthquakes to regions of "bending" and regions of "unbending." In the first, we determine the point of maximum curvature on each profile, and assign earthquakes up-dip of this point as "bending," and down-dip as "unbending." In the second, we instead assign earthquakes based on the gradient in curvature (Figures 2-5f); earthquakes located where the calculated gradient in curvature is positive are designated to be related "bending," and those corresponding to negative values designated to be related to "unbending." We also consider further separation into earthquakes related to down-dip extension, and those related to down-dip compression, based on the orientation of the P and T axes with respect to the slab dip vector. Note that, while we illustrate these divisions in Figures 2-5 using the mean slab profiles, the point of peak curvature, and the gradient of curvature are calculated for each profile, with each earthquake using the values from their closest associated profile.
In a number of cases (e.g., Aleutians, Figure S2 in Supporting Information S1; Solomons, Figure S6 in Supporting Information S1) where the determination of the deeper slab geometry is limited by a paucity of data, we note that the slab models used default to a flat slab with zero dip at depth. This could, in theory, lead to an under-estimation of the first down-dip return to zero curvature, and lead to the exclusion of some earthquakes that would otherwise be classified as "unbending." However, the lack of a well-constrained slab model in these regions is partly due to a lack of earthquakes, and hence we do not consider this to be a problem for our analysis, both in these specific regions, and in general.
In Figure 6, we show the results for the moment summation of all 13 regions considered here. In each region, we sum earthquakes by their relation to regions of bending and unbending (in the case of Figure 6, on the basis of the rate of change in curvature), and also separate by mechanism type. We also highlight in each case the contribution of the largest earthquake within each grouping, shown by the white bars. In many cases, the moment release is dominated by one major event (e.g., the M w 8.3 1977 Sumba earthquake for the Java region; the M w 7.9 2014 Rat Islands earthquake for the Aleutian Arc)-a problem which limits the extent to which the overall deformation state of any given slab, when taken in isolation, can be assessed from available earthquake catalogs, which may not be entirely representative of the long term deformation pattern.
In Figure S12 in Supporting Information S1, we show an equivalent set of summations to those in Figure 6, but instead defining the separation into bending and unbending based on whether earthquakes are updip of the point of maximum curvature on the closest slab profile to the earthquake location, or between the point of maximum curvature and the first subsequent zero-curvature point. As these figures show, this different method of separation makes little difference for the majority of regions. The only major changes are for the Ryukyu and New Britain subduction zones, where the different definition changes the moment balance between unbending and bending because the largest earthquake in the population changes from bending to unbending (or the inverse).
Most individual regions we consider are subject to sampling bias, given both the dominance of a small number of individual events, and that, with only a 46-yr catalog duration, we are looking at a relatively small portion of the seismic cycle in such regions. Hence, instead of further interrogating individual regions, we combine our regional catalogs into a single composite catalog. Individual regions differ in a number of crucial ways, including different rates of subduction and degree of curvature, leading to differing strain rates; differing seismogenic thicknesses; and differing geodynamic settings. Of the 13 regions we consider, and shown on Figure 6, we exclude one, Sumatra ( Figure S9 in Supporting Information S1), from further consideration. Sumatra shows a notably different pattern in intraplate seismicity, with the shift from seismicity evidencing bending to seismicity suggesting unbending occurring slightly seawards of the trench (Craig & Copley, 2018)-significantly further up-dip than in any other subduction zone, and in a manner that does not match the long-term strain implied by the geometry of the plate interface (Singh et al., 2012). This notably different behavior of the intraslab seismicity may be a result either of a dynamically evolving overriding plate, affecting the near-trench intraslab stress field (Craig & Copley, 2018), or a consequence of the diffuse intraplate deformation occurring seawards of the trench within the Indian Ocean to accommodate the differential motion of India with respect to Australia (e.g., Geersen et al., 2015;Wiens et al., 1985). We therefore exclude it from our further discussion, and from Figure 7.  Figure 7 shows the summation of moments from all 12 remaining regions from Figure 6. Again, the white bars show the contribution of the largest event to each bin (with the largest single event being the M w 8.3 1977 Sumba earthquake). The upper panel shows the results when earthquakes are divided based on the rate of change in curvature at the location of each earthquake, the lower panel when they are separated based on whether they are up-dip or down-dip of the maximum curvature point on their associated profile. Again, the different mechanism of defining the separation from bending to unbending makes little difference to the overall moment balance. Two observations stand out, regardless of which separation approach is used.
First, substantially more seismic moment occurs up-dip of the point of maximum curvature than occurs between that point and the first zero-crossing point, suggesting that more seismic moment is released in the initial bending process than in the recovery of the same total curvature. Given that, as slabs descend into the mantle, they will heat up, their seismogenic thickness decreases, and increasing amounts of strain will be accommodated through ductile deformation, this seems reasonable. However, the majority of the cases considered here are old, cold subducting plates, where the slab temperature and rheology are unlikely to evolve rapidly during the initial phases of subduction. This variation may instead be explained by the requirement for deeper earthquakes to occur in regions where the evolving mineralogical composition of the subducting plate leads to the release of fluid, permitting seismogenic failure to occur.
Second, moment release through earthquakes accommodating down-dip extension is significantly higher than that released in earthquakes accommodating down-dip compression. While it remains subject to sampling problems related to the dominance of the largest individual earthquake, this latter trend also holds for all of the individual regions studied ( Figure 6), with the exception of Bonin. This asymmetry is most pronounced in the initial bending region, in line with expectations that the supported stresses likely increase with depth, as confining pressure and hence yield stress both increase, and with the increasing contribution of the ductile lithosphere in supporting a small proportion of the total stress.

Geodynamic Implications
Inferring geodynamic processes (fundamentally dependent on stress and strain) from moment release can be a complex process, often with multiple non-unique implications (Chapple & Forsyth, 1979;Mueller et al., 1996). Following Kostrov (1974), moment release from a population of earthquakes can be related to volumetric strain (ϵ ij ) using: Where V is the volume under consideration and μ is the shear modulus, under the assumption that all strain is accommodated seismically. In the context of subducting slabs, both the shear modulus and the volume (in effect, the seismogenic cross section of the slab) under consideration will vary with depth, as the rheology of the slab evolves as it descends into the Earth's interior and is subject to increased pressures and temperatures. Exactly how these parameters vary will be different in different slabs, depending on their geometry and the thermo-chemical structure of the incoming plate. However, broadly speaking, we expect the shear modulus to increase slightly with increasing down-dip distance from the trench, and the volume to decrease slightly with distance from the trench, as the plate heats up. In the majority of cases considered here, however, where slabs are usually old and cold at the point of subduction, these effects are likely to be minimal over the depth range we consider.
The narrow vertical separation of regions in horizontal extension and compression in outer rises (Craig, Copley, & Jackson, 2014) suggests that the core of a bending tectonic plate which deforms elastically over geological timescales, and is unbroken in earthquakes, is relatively narrow (≲5 km). However, the extent of this elastic core will vary with the elastic limits of the material. While estimates of the effective coefficient of friction on intraplate oceanic faults suggest that it is low (Craig, Copley, & Middleton, 2014), the elastic core width will still be expected to increase with increasing depth and confining pressure. This suggests that the elastic core would widen to some degree with increasing depth, again potentially decreasing the amount of strain, and related moment release, that is accommodated though permanent seismogenic deformation. This may play a role in reducing the seismic moment release down-dip of the inflection point in curvature, and may address the decrease in moment release between the accumulation and recovery of curvature. The depth of the transition from down-dip compression to down-dip extension will also change between regions of bending and unbending, as the bending-related stresses are superimposed on any in-plane stresses arising from additional plate-driving and resisting forces (e.g., slab-pull and ridge push). In the regions considered here, with the exception of south of Ryukyu, and potentially northernmost Tonga, the depth of the transition from normal to thrust faulting beneath the outer rise region is more than half the depth of the deepest earthquakes there (Craig, Copley, & Jackson, 2014). In the perhaps simplistic view that bending strain is linearly proportional to plate-perpendicular distance from the elastic core, and following Equation 3, this matches with a greater predicted moment release in the bending region through down-dip extension than through down-dip compression-consistent with the results shown in Figure 7. In the absence of in-plane stresses, unbending should result in the complete reversal of accumulated strain, and would be matched by a mirrored moment release. To explain the greater extensional moment release through unbending as well as bending implies the addition of in-plane stresses which, averaged over all regions considered, are down-dip extensional.
Thus, the dominance of moment release through down-dip extensional earthquakes over down-dip compression throughout the initial curvature cycle indicates that, although the pattern and along-dip distribution of seismicity is strongly related to changes in the plate curvature (Sandiford et al., 2020), a portion of the deformation that takes place in these regions is unrelated to the bending of the slab, but is instead unrecovered down-dip extension of the whole slab. That this permanent deformation seems to take place in regions of increased change in curvature suggests that the additional stresses associated with slab bending are necessary to push the slab beyond its elastic limit and into the regime of brittle failure, and that in-plane forces (e.g., slab pull), while modulating the depth of the transition between bending stress regimes, are insufficient alone to produce widespread intraslab seismicity. This pattern is consistent with the observation that mature oceanic lithosphere is generally able to support the stresses transmitted from subducted slabs without undergoing significant deformation (i.e., a necking instability).
The conclusion-that apparent bending-related seismicity also accumulates a component of long-term bulk strain-is supported by seismological and bathymetric evidence for the orientation of active fault planes in both the outer rise and the regions of intraslab seismicity (Masson, 1991;Myhill & Warren, 2012;Ranero et al., 2005;Warren et al., 2007), wherein an overall preferential dip-direction for the active fault plane is consistent the accumulation of down-dip strain (i.e., slab thinning). While our study has deliberately focused on earthquakes at intermediate depth (and shallower) within the slab, that the depth-dependence in focal preferential mechanism orientation observed by Warren et al. (2007) in the Tonga subduction zone, where earthquakes shallower than 300 km show a preferential plane orientation but those below do not, raises the possibility that this accumulation of bulk strain via the modulation of bending-induced seismicity only occurs during the initial stages of subduction. Figure 8 shows a simple conceptual model for the accumulation of permanent down-dip deformation through the recovered bending cycle through the superposition of a down-dip in-plane stress, moving the neutral fiber down (in the bending case) and up (in the unbending case). The change in the depth of neutral fiber that results from the addition of an in-plane stress leads to the bending strain field being different to the unbending strain field, and allows the accumulation of permanent bulk strain. In the idealized example shown, the strain resulting purely Figure 8. Simplified sketch illustrating how strain imbalance may be achieved through addition of an in-plane stress, superimposed on bending stresses. Red indicates down-dip compression, blue indicates down-dip extension. Crosshatched areas show strain accommodated through permanent (i.e., non-elastic) deformation likely to be seismogenic if conditions allow. Gray shaded area shows the elastic core separating areas of potentially brittle failure. Yellow dashed line shows the elastic limit on strain, assumed to be depth-independent in this conceptual model. from the in-plane stress could be accommodated elastically, in the absence of any bending stresses, and therefore would not be expected to produce seismicity in a straight section of slab. The sensitivity of the neutral plane depth is itself dependent on the strength profile of the subducted lithosphere-the lower the yield stress of the plate, the greater the sensitivity of the neutral plane depth will be to variation in the in-plane stress, and the more profound the effect on the accumulation of permanent deformation will be.
Despite the improved resolution in earthquake locations available from the ISC-EHB catalog, we are typically unable to image accurately the separation between upper and lower seismic zones through bending regions. To do so requires either detailed teleseismic analysis not yet routinely undertaken (e.g., Craig, 2019;Florez & Prieto, 2019), or high-quality local seismic data not globally available (e.g., Sippl et al., 2018;Wei et al., 2017). However, our results suggest that we would, in general, see a upward-migration in the depth within the plate of the elastic core upon moving from initial bending to initial unbending (note that, due to the summation approach used here, this may not be the case in all subduction regions considered, but we expect this to be the case in the majority of regions). Again, we emphasize that the results shown in Figure 7, on which this interpretation is based, are drawn from the combination of multiple regions, in which the stress state of the whole slab will vary. The scenario we describe appears to be the average case, but we do not expect this to necessarily apply in all subduction zones when considered individually (e.g., not the case for Tonga-Kermadec; see Figure 6). This discussion has focused on the seismic expression of intraplate strain, but this is only part of the total intraplate strain, and significant strain is accommodated within the plate through ductile deformation of material at higher temperatures. However, the material undergoing ductile deformation at high temperatures is unlikely to support significant stresses in comparison to material deforming brittlely at lower temperatures, and thus we consider the contribution that ductilely deforming sections of the slab make to the overall support of the intraplate stress field to be negligible. Ductile strain is excluded from Figure 8, which focuses on deformation in areas potentially capable of producing earthquakes. In the future, as our understanding of the rheology of subducting plates, and our observational seismic catalog, develop, it may be possible to construct geodynamic models for the subduction process that can determine the ways in which stress is supported in the slab through brittle, semi-brittle, and viscous processes, and to more directly relate modeled estimates of volumetric strain to the distribution of seismicity.
The simple rheological models considered in Figures 1a and 8, and that underpin our discussion, assume a deformation regime in which stresses within the plate are predominantly supported on faults that are both brittle and seismogenic. One complexity we have not considered is the potential for a transitional regime at the base of brittle section of the subducting plate, where strain is accommodated by stable sliding on velocity-strengthening material, rather than through rapid, velocity-weakening seismogenic failure. Such a region would remain capable of supporting significant stress, but would be undetectable using our seismicity catalogs. Were such a region to be present at the base of the stress-supporting plate, it would lead to our seismicity catalogs systematically under-representing down-dip compressional deformation in regions of increasing curvature, and under-representing down-dip extensional deformation in regions of decreasing curvature. This would impact on our interpretation of Figures 6 and 7, as the moment balances observed would not be fully representative of the supported stresses, or accumulated strains. While it would alter the details of the conceptual model presented in Figure 8, it would not alter the overall implications.

Conclusions
Seismicity within subducting slabs is predominantly concentrated in areas of the slab where the down-dip slab curvature is rapidly changing, suggesting that slab bending stresses exert a strong first-order control on the location of intraslab earthquakes. In this study, we have studied the intraslab seismicity associated with a number of subduction zones around the western Pacific and eastern Indian Ocean, where down-dip slab geometries are relatively simple, characterized by the development and recovery of a single major down-dip bend. The shallow intraslab seismicity associated with this first bend and subsequent unbending demonstrates that there is significantly more moment released in down-dip tensional intraslab seismicity than in down-dip compressional earthquakes, in the initial regions of bending and unbending. This imbalance in moment release across the cycle of the initial slab bending indicates that, overall, slabs are undergoing a degree of permanent down-dip extension, consistent with models of (mild) slab necking. We propose a model wherein this accumulation of permanent down-dip strain arises from the variation in the depth within the slab of the neutral fiber associated with each bend, resulting from the superposition of down-dip intraslab stresses on the bending stresses. The down-dip stress alone is usually insufficient to produce brittle failure in the slab, but can act to modulate the depth of the transition between bending and unbending, leading to the accumulation of permanent strain from low intraslab stresses.