Plate bending earthquakes and the strength distribution of the lithosphere

SUMMARY

This study investigates the dynamics and constitutive behaviour of the oceanic lithosphere as it bends and yields during subduction. Two main observational constraints are considered: the maximum bending moment that can be supported by the lithosphere, and the inferred neutral plane depth in bending. We particularly focus on regions of old lithosphere where the ‘apparent’ neutral plane depth is about 30 km. We use subduction modelling approaches to investigate these flexural characteristics. We reassess bending moment estimates from a range of previous studies, and show a significant convergence towards what we call the ‘intermediate’ range of lithosphere strength: weaker than some classical models predict, but stronger than recent inferences at seamounts. We consider the non-uniqueness that arises due to the trade-offs in strength as well background (tectonic) stress state. We outline this problem with several end-member models, which differ in regard to relative strength in the brittle and ductile regimes. We evaluate the consistency of these models in terms of a range of constraints, primarily the seismic expression of the outer rise. We show that a 30 km neutral plane depth implies that net slab pull is not greater than about 2 TN m⁻¹. In contrast, models with low brittle strength imply that regions with a 30 km neutral plane depth are under moderate net axial compression. Under these conditions, reverse faulting is predicted beneath the neutral plane at depths >30 km. We show that moderate variations in background stress have a large impact on the predicted anelastic dissipation. We suggest brittle reverse faulting is a marginal phenomenon which may be inhibited by moderate changes in background stress.

1 INTRODUCTION

The depth distribution of strength in the oceanic lithosphere is thought to arise mainly from the interaction between the normal-stress sensitivity of faults in the brittle domain and the temperature dependence of intracrystalline creep. An important aspect of this strength distribution is the relatively low flexural strength relative to axial strength, because the strain developed in bending has an inverse character to the strength profile. The influential view that subduction bending is a major resisting force in mantle convection dynamics (e.g. Conrad & Hager 1999) has therefore receded somewhat (Capitanio et al. 2009; Leng & Zhong 2010; Buffett & Becker 2012). Nevertheless, there is still significant debate about the role and interplay of various deformation mechanisms, the maximum stress levels in bending, and the interaction of bending related stress with other sources of stress—usually anticipated to be the background component associated with the driving/resisting forces of plate motion.

We start by clarifying some assumptions and, in particular, what we mean by strength. We assume that the response of the lithosphere in subduction-related loading can be adequately described by simple bending of a plate in plane strain, and that the stress state that develops is dominated by this bending, although it may be modified by background stresses. This set of assumptions seems justified because, while the oceanic lithosphere away from subduction trenches exhibits diverse focal mechanisms (normal, reverse and strike slip), the ∼100 km region near trenches is overwhelmingly dominated by normal faulting in the upper 15–35 km, with subhorizontal tension axes oriented perpendicular to the trench, consistent with a state of effective tension induced by bending (Stein & Pelayo 1991; Craig et al. 2014a). Note that we use the term ‘effective tension’ to describe the deviatoric stress state associated with (Andersonian) normal-faulting; all principal stresses are of course compressive in the depth range that is relevant to the strength of the lithosphere.

Under typical values of subduction-hinge curvature (e.g. 5× 10⁻⁶ m⁻²), elastic moduli, and laboratory-constrained deformation mechanisms, the lithosphere is expected to undergo comprehensive yielding (Goetze & Evans 1979). In this situation the bending moment saturates, such that it does not increase with additional strain (bending). By lithospheric strength, we hence mean the saturation moment of a given rheological model. Saturation moment depends not only on strain (i.e. curvature) but on the strain rate. Hence is is only valid to talk about the saturation moment of a given strength model, when the strain rate model is specified. The assumptions we make in estimating saturation moment, along with other aspects of the modelling, are discussed in Section 4.1.

In this study we focus particularly on lithosphere of 80–120 Myr age, as this region is well represented in the current global distribution of subducting plate by age, and is also where the inferred neutral plane depth can be estimated with most confidence. To aid discussion, we characterize the saturation moment of different rheological models as being ‘weak’, ‘intermediate’ and ‘strong’. The proposed strength of the oceanic lithosphere in previous studies spans this range (Chapple & Forsyth 1979; McNutt & Menard 1982; Levitt & Sandwell 1995; Garcia et al. 2019; Bellas & Zhong 2021). These ranges of saturation moment are shown in Fig. 1, and are based on a simple ‘family’ of strength models discussed in Section 6. The regions are defined as follows: for 100 Myr-old lithosphere, weak models have saturation moments (<1); intermediate models (1–2) and strong models (>3). These ranges are measured in units of N-m m⁻¹ (N) × 10¹⁷. In the following we use (N) to describe the bending moment per unit-length of trench.

For lithosphere in the 80–120 Myr range, outer rise seismicity exhibits an apparent transition from normal- to reverse-faulting mechanisms at a depth of ∼ 3 km (e.g. Fig. 3). It is usually (although not always) assumed that this observation reflects a stress transition across the neutral plane associated with bending (Chapple & Forsyth 1979). The thermal thickness of plates in the 80–120 Myr range is 100–130 km (Parsons & Sclater 1977; Richards et al. 2020), so that the inferred neutral plane lies at close to a quarter the plate depth. This in turn suggests that plate strength lies mainly in the upper half of the plate.

Classical yield strength envelopes (YSEs) based on experimental constitutive laws generally predict satisfactory neutral plane depths (30–40 km for ∼100 Myr lithosphere), under the assumptions that the background stress is insignificant, the friction coefficient is relatively high (i.e. comparable to Byerlee’s Law), and when the exponential stress dependence of olivine (referred to here as LTP for low-temperature plasticity) is taken into account (e.g. Byerlee 1978; Goetze & Evans 1979; McNutt & Menard 1982; Hunter & Watts 2016). We note that while this classical model predicts satisfactory neutral plane depths, it typically predicts that LTP will function as the dominant deformation mechanism beneath the neutral plane, with the brittle (frictional) strength in compression several times higher than the steady-state flow stress (Goetze & Evans 1979). This prediction conflicts with observations of deeper seismicity below the neutral plane indicating trench-perpendicular compression.

However, many of these standard assumptions might be questioned. For instance, if net axial force (see 2.4) due to slab pull is significant (i.e. much greater than typical ridge push magnitudes) and/or if outer rise normal faults have low friction coefficient (<0.3; e.g. Craig et al. 2014b) the ∼30 km apparent neutral plane depth is difficult to reconcile, unless differential stress levels in the ∼30–60 km region are also significantly reduced relative to typical LTP models (e.g. Goetze & Evans 1979; Mei et al. 2010). The net effect is a much weaker strength model, and would be broadly consistent with conclusions that have emerged from recent studies of seamount flexure (Zhong & Watts 2013; Bellas et al. 2022).

The torque that acts on the plate due to the flexural topography, referred to as the ‘observed moment’ by previous studies, provides a rheology-independent constraint on the bending moment (Goetze & Evans 1979; McNutt & Menard 1982). The mathematical background is provided in Section S1. The bending moment and observed moment are defined in eqs (S5) and (S9). Together, we refer to bending moment and neutral plane estimates as the ‘flexural characteristics’. There is a complimentary aspect to these constraints as the neutral plane is sensitive to the relative (depth) distribution of stress in the lithosphere, while the observed moment has sensitivity to magnitude of those stresses (albeit via the first moment integral). However, the constraints are not sufficient to uniquely determine the strength distribution with depth, and the non-uniqueness is a key focus of this study.

After preliminary sections that cover the background (Section 2), data (Section 3), an overview of modelling approaches (Section 4) and modelling insights (Section 5), the study consists of two main pieces of analysis. In Section 6, we consider a range of existing data sets that constrain the maximum bending moment in subducting lithosphere. We argue for a reinterpretation of the observed moment data of Levitt & Sandwell (1995), in terms of a weaker strength model than originally proposed. A key conclusion is that when interpreted through a ‘common lens’ of the maximum bending moment, a range of studies have converged on a similar ‘intermediate’ strength for old oceanic lithosphere (Chapple & Forsyth 1979; Levitt & Sandwell 1995; Hunter & Watts 2016; Garcia et al. 2019). In Section 7, we consider the less well-constrained part of the problem, which is the relative distribution of the strength with depth. We characterize a parameter space consisting of the relative strength of the brittle and ductile parts of the lithosphere, along with the background stress and explore the implications of different parts of the parameter space.

2 BACKGROUND AND PREVIOUS WORK

2.1 Flexural strength of the oceanic lithosphere

There is an extensive body of literature that has used the deflection of the plate surface due to loading to infer the mechanical properties of the oceanic lithosphere. The primary environs are subduction zone outer rises and seamounts. The intention here is not to attempt a summary, but to recall for the reader that: (1) studies of subduction zones have inferred a large range of strength models (from what we call weak, through to strong; see Fig. 1) and (2) that studies of seamount flexure have increasingly made the case for weak or intermediate models.

Flexural studies of subduction zones offer different perspectives on acceptable models of lithosphere strength, as underscored by the work of Chapple & Forsyth (1979), McNutt & Menard (1982), Levitt & Sandwell (1995), Hunter & Watts (2016) and Garcia et al. (2019). In each case these studies address the inverse problem of generating optimal models of plate strength in terms of either: (1) direct fitting of forward models with observed topography or gravity profiles (Chapple & Forsyth 1979; Hunter & Watts 2016; Garcia et al. 2019) or (2) comparison between the ‘observed’ moment associated with flexural topography, and the predicted saturation moment of different strength models (Goetze & Evans 1979; McNutt & Menard 1982; Levitt & Sandwell 1995). We note that the observed moment is typically calculated around the first zero crossing of the flexural topography, as this eliminates the effect of the net axial force term in the moment balance (see Section S1).

Chapple & Forsyth (1979) solve the thin-plate equations using an iterative numerical approach that incorporates the moment variation with curvature for different ad hoc depth-strength profiles. The profiles are calibrated so that they approximately match the inferred neutral plane constraint, based on seismological data available at the time. Their optimal model, which is developed with reference to old lithosphere (100–140 Myr) is shown in Fig. 1. It produces a saturation moment of  1.8 × 10¹⁷ N, and lies within what we have called the intermediate strength range.

Garcia et al. (2019) solve the thin-plate equations for 3-D subduction segments with variable loading. Like Chapple & Forsyth (1979), the constitutive model includes anelastic deformation (yielding). In this case, they incorporate a moment-curvature relationship underpinned by a classical strength model (olivine creep and brittle failure, see Section 4.2). They vary the friction coefficient to produce a strong (μ = 0.6) and a weak (μ = 0.3) model (both assuming a hydrostatic pore pressure). They find that the ‘weak’ model can better fit data across Pacific basin subduction zones. This model produces a saturation moment of 1.3 × 10¹⁷ N for 100 Myr lithosphere, and lies at the lower end of our intermediate strength range. We note the ‘weak’ and ‘strong’ models described in Garcia et al. (2019) both lie within what we call the intermediate range (suggesting that these authors had already narrowed the focus of their investigation to that strength region).

Hunter & Watts (2016) optimize elastic deflection models for circum-Pacific subduction zones, based on both constant-thickness and variable thickness plate models. To compare these results to laboratory constrains, YSEs are computed for various strength models, then converted to equivalent elastic plate thicknesses, and compared with the optimized values from the inversion. They argue that their inversions are consistent with several combinations of friction coefficient and low-temperature plasticity, including the parameters of Mei et al. (2010), with friction co-efficients of either 0.6 or 0.3 (in both cases assuming hydrostatic pore pressure). These combinations can be shown to describe strength models that lie in the intermediate range (e.g. Fig. 1).

An alternative way to interpret the results of Hunter & Watts (2016), is to consider the bending moment predicted by the elastic flexure models at the first zero crossing. This subjects their results to an equivalent analysis as is undertaken by Goetze & Evans (1979), McNutt & Menard (1982) and Levitt & Sandwell (1995), in which an elastic plane solution is used to fit the flexural topography, and then the moment at the first zero crossing is estimated. Data from fig. 10 of Hunter & Watts (2016) was re-digitized, and the moment values at the first zero crossing are shown with red crosses in Fig. 1. These estimates reflect inversion for averaged flexural deflections (’trench-bin’) typically aggregated across several thousand kilometres. In contrast, the blue points in Fig. 1 represent estimates from individual trench profiles [from fig. 1 Levitt & Sandwell (1995)].

The studies of McNutt & Menard (1982) and Levitt & Sandwell (1995) use similar approaches: estimating the observed moment based on uniform elastic plate fitting. However, they arrived at quite different conclusions: weak and strong models respectively (following our definition of strength). There are two reasons for this. First, the estimates derived in McNutt & Menard (1982) mainly lie within the weak strength region (green points in Fig 1), while the larger data set of Levitt & Sandwell (1995) contains moment values across weak, intermediate and strong regions (red points in Fig 1). We discuss potential reasons for this discrepancy in Section 6.3. The second issue relates to the interpretation of scatter in the data. Levitt & Sandwell (1995) argue that an acceptable lithospheric strength model should provide a ceiling to the observed moment data, which implies models in the strong region; their proposed model is shown with a dashed blue line in Fig. 1. In contrast, McNutt & Menard (1982) sought a strength model that would provide a ‘best fit’ to the data. This issue of interpreting scatter in the observed moment data is the focus of Section 6.

In summary, previous studies addressing lithospheric strength in subduction zones have advocated models throughout what we define as the weak, intermediate and strong regions (as shown in Fig. 1). However, we must be attentive to the assumptions, ambiguity, and interpretative aspects of some of these conclusions. For instance, the basis on which Levitt & Sandwell (1995) advocate a strong lithosphere model is tied to the interpretation of the scatter in the data.

Meanwhile, a number of studies have concluded that observations of lithosphere flexure around Hawaii and seamounts more generally, require that the lithosphere is substantially weaker than inferred at subduction zones Zhong & Watts (2013), Pleus et al. (2020), Bellas & Zhong (2021) and Bellas et al. (2022). We do not consider seamount flexure directly in this study, however our analysis has implications for the the purported ‘conundrum’ of strong plates at subduction zones and weak plates around seamounts.

2.2 The neutral plane as a constraint on strength

For lithosphere in the 80–120 Myr range, outer rise seismicity exhibits an apparent normal-reverse mechanism transition depth of ∼30 km (e.g. Fig. 3, see Section 3). Note that all references to depths refer to depth beneath the surface of the plate (seafloor). It is usually assumed that this observation reflects a transition across a neutral plane in the bending of the plate as subduction proceeds. A contrasting view is discussed the following section.

Chapple & Forsyth (1979) noted that ‘the distribution, mechanisms, and depths of the earthquakes associated with the bending of the lithosphere provide the information needed to reduce greatly the non-uniqueness inherent in modelling topography alone’. In Chapple & Forsyth (1979), this process involved testing ad hoc strength-depth models to match the neutral plane constraint as well as flexural-topographic profiles from old subducting lithosphere. They noted that the apparent neutral plane depth for lithosphere in the 100 Myr range is ∼30 km, an observation that we discuss in the next section. The strength model derived by Chapple & Forsyth (1979) is shown with a dashed cyan-coloured line in Fig. 1, and lies in the range we describe as intermediate strength.

In exploring different ways to weaken classical strength models, McNutt & Menard (1982) argued against frictional-strength reduction because the resulting model produced a neutral plane depth of ∼40 km for 100 Myr lithosphere. Their preferred strength model instead achieves strength reduction by appealing to significantly reduced creep strength (LTP) through modification of the activation energy in an existing flow law. Interestingly, this produces the opposite problem: a neutral plane that is shallower (20–25 km for 100 Myr lithosphere) than the apparent depth inferred from seismicity.

In reality, the background stress contribution also effects the depth of the neutral plane. If the typical stress state of lithosphere prior to subduction (and its variability) was well constrained, the neutral plane would potentially greatly reduce non-uniqueness, as Chapple & Forsyth (1979) propose. Unfortunately, because it is not, the problem remains under-constrained. We address this issue in Section 7.

2.3 The outer rise normal-reverse transition

Outer rise normal-faulting earthquakes, where subhorizontal tension axes are oriented perpendicular to the trench, have long been associated with plate bending, (e.g. Stauder 1968). It should be noted that the seismicity is not concentrated at the bathymetric rise/fore-bulge, but rather occurs throughout the domain of plate curvature increase, which may extend beneath the forearc (Sandiford et al. 2020). In this study, the term ‘outer rise’ is used to refer to the domain of curvature increase in subducting plates, and the associated dynamics.

Extensional earthquakes in the outer rise region are generally located at systematically shallower depth than less frequent reverse events (e.g. see Fig. 3). These reverse-mechanism earthquakes have usually been attributed to compressional bending stresses beneath a neutral plane (Chapple & Forsyth 1979; Craig et al. 2014a). The distribution of outer rise reverse events is non-uniform, with few or no occurrences along some subduction zones (such as eastern Sunda and the Aleutians). Most activity has been along the Kurile, Japan, Philippine, Solomon Islands, Tonga, Kermadec and Chile zones (Craig et al. 2014a). Even in those regions there is notable clustering, with a concentration of events in central-to-southern Honshu, central-to-northern Kuriles, northern Ryuku, southern Philippines and central Chile. Tonga-Kermadec represent the only region where outer rise compressional events are common across several thousand kilometres (the reader is referred to maps in Craig et al. (2014a); Ye et al. (2021).

The study that has most directly challenged the idea of a reverse faulting as a consequence of plate bending induced shortening/compression is Mueller et al. (1996b). They argued that the purported association is compromised by both laboratory and seismologic observations. Of course, this perspective would invalidate attempts to infer a neutral plane on the basis of the normal-reverse seismicity transition. The argument can be summarized in the following points:

1. Brittle failure is an essential condition for seismicity. Reverse-faulting earthquakes, occurring at a depth of 30–50 km, are inconsistent with Coulomb-type frictional failure unless the effective friction coefficients drastically reduced (μ < 0.1). But the flexural topography of the outer rise requires that friction coefficients cannot be this low.

2. There is insufficient evidence for material instability, that is fault development with unstable stress drop, in the semi-brittle or plastic regimes, as was suggested by Chapple & Forsyth (1979).

3. Attributing both thrust- and normal-faulting outer rise events to the same mechanism cannot directly account for the important observation that the former occur much less frequently than the latter.

Instead, Mueller et al. (1996b) propose that reverse-faulting earthquake mechanisms in subduction outer rises are associated with regions of anomalous compression, which develop in response to strongly localized subduction-resisting forces, such as interplate asperities or buoyant subducted crust. They argue that many such earthquakes may be shallower than previously reported, thus obviating the incompatibility with Coulomb-type frictional strength considerations.

We do not have the scope to address each of these points, but will instead focus on the seismological data (with the benefit of more than 20 yr’ additional accumulation). In Section 3, we review the data from a range of sources. In our assessment, these now provide a much more compelling case that in many regions the outer rise exhibits a transition from normal faulting to reverse faulting, in close proximity.

2.4 Background sources of stress in the bending lithosphere

Although flexure can clearly be modified by the background stress state, in the modelling of subduction zone flexure, there seems to be little consensus as to the overall pattern of such contributions. Some early studies include a net background component essentially as a free parameter, and advocated large compressive forces (Hanks 1971). Parsons & Molnar (1976) showed that applied moments, which had been neglected in some previous studies, could produce equivalent deflection without requiring large axial compression. In many subsequent studies, the background stress is ignored, or plays a subordinate role (Turcotte et al. 1978; Chapple & Forsyth 1979; Hunter & Watts 2016; Garcia et al. 2019).

Different models for the sources of background stress, imply differences in terms of reference state, stress magnitude, and spatial and temporal variability. A useful dichotomy is articulated in Mueller et al. (1996b), who state that while ‘driving forces are broadly distributed and... resisting forces [may be] relatively localized’. Mueller et al. (1996b) proposed that where outer rise seismicity exhibits reverse faulting mechanisms, the stress state is dominated by localized resisting forces, such as asperities in the interplate zone.

On the other hand, if the most important driving force for oceanic plates is net slab pull (of up to perhaps 10 TN m⁻¹; Conrad & Lithgow-Bertelloni 2002; Faccenna et al. 2012), a state of effective tension seems the logical a priori choice for the vast majority of the global subduction system; it would seem strange not to see the ‘signal’ of such strong background stress in some aspect of the flexural behaviour, for instance in the depth to the neutral plane.

However, the relative partitioning of plate driving forces remains debated (Becker & O’Connell 2001; Conrad & Lithgow-Bertelloni 2002; Sandiford et al. 2005; Copley et al. 2010; Ghosh & Holt 2012; Husson 2012; Stotz et al. 2018). If the mantle has an active driving role beneath the plates (Stotz et al. 2018; Liu et al. 2021), or if gravitational potential energy changes across the trench-outer rise are as large as some studies suggest (Bessat et al. 2020), a compressional rather than (effective) tensional background stress state may prevail. The numerical model we discuss in Section 5 provides some important insights into this question.

Flexure of the lithosphere is commonly modelled as a cylindrical bending problem (i.e. 2-D plain strain). Yet the background deviatoric stress state is clearly 3-D; under Andersonian assumptions the order of two remaining principal stresses is variable, as is the angular relationship to the trench or the plane of bending. In the 2-D modelling of subduction flexure, the background stress is usually incorporated in terms of a prescribed net force, with either a (effective) tensional or compressional sense. In this study we refer to this simplification (of the 3-D stress) as the net axial force denoted as F_net, and defined in eq. (S4). Throughout the paper we refer to the net axial force magnitude as being ‘low’ (≤3 TN m⁻¹), ‘moderate’ (∼5 TN m⁻¹) and ‘high’ (≥10 TN m⁻¹).

3 DATA

3.1 Outer rise seismicity

Fig. 3(b) shows well-constrained centroid depths of outer rise earthquakes, after Craig et al. (2014a) and sources therein. Although there has historically been some controversy about the veracity of the systematic depth offset (e.g. Mueller et al. 1996b), this pattern is clearly evident for lithosphere >80 Myr, as shown in Fig. 3(b). These data also reaffirm the conclusion of Chapple & Forsyth (1979), that the apparent normal-reverse transition depth for 80–120 Myr lithosphere is at about 30 km. The distribution becomes far more ambiguous, however, for lithosphere younger than 80 Myr, although we note that the low proportion of subducting lithosphere in the 60–80 Myr range (see Fig. 3a) is responsible for producing the substantial gap in data in this age range.

The conclusions drawn from the global distribution of moderate-magnitude outer rise earthquakes (e.g. Craig et al. 2014a, Fig. 3) are generally consistent with more localized insights gained from finite-fault modelling of individual events, doublet occurrences and microseismcity (e.g. Obana et al. 2012; Todd & Lay 2013; Ye et al. 2021). Selected examples are summarized here.

Finite-fault modelling of the M_w 7.5 2020 Kurile earthquake (Ye et al. 2021), confirms that rupture was limited to depths beneath the anticipated neutral plane depth estimated from hypocentres for lithosphere of that age (e.g. Fig. 3). The reported aftershocks of the 2020 Kurile event are concentrated at depths of 27–45 km, consistent with the expectation that these events lie between the neutral plane and an isotherm close to 600 °C (Emmerson & McKenzie 2007).

A number of instances of normal-reverse doublets provide further evidence that the stress state in the outer rise transitions from effective tension to compression with depth. These include the 2012, Honshu Japan doublet, the 2011 Kermadec doublet, the 2006–2009 Kurile doublet, and the 2012 Phillippine trench events [see Ye et al. (2021) Fig. 7 for summary]. Aftershocks recorded with OBS after the 2012 Honshu doublet, showed normal-faulting seismicity down to about 30 km, relative to the plate surface (see Obana et al. 2014, Fig. 2).

Following the 2011 Tohoku-Oki earthquake, a permanent seafloor seismic network was installed offshore Honshu (S-net). While S-net is not used to routinely determine earthquake hypocentres, a recent study us S-net data to relocate small events in the outer rise (Zhao et al. 2022, fig. 11). The relocated hypocentres suggest clustering into shallower and deeper groups, with a separation between the two occurring about 30 km relative to the plate surface. The deeper cluster corresponds to events with reverse focal mechanisms (personal communication with author). This is broadly consistent with the 30–35 km neutral plane depth inferred based on OBS data obtained following the 2011 Tohoku-Oki earthquake (Kubota et al. 2019). S-net data has the capacity to significantly improve our future understanding of outer-rise seismic expression, including its temporal and spatial variation.

In summary, we think that a normal-reverse transition with systematic depth offset, is a consistent phenomenon in global outer rise seismicity in older plates, albeit unobserved in some regions. We think that the simplest explanation for this seismic expression relates to stress/strain conditions across the neutral plane of bending.

3.1.1 The apparent depth of the normal-reverse transition

The data and studies summarized in the previous section suggest that the apparent normal-reverse transition depth for 80–120 Myr lithosphere is around 30 km, although there is some regional variability. A key question is to what extent this apparent transition is an intrinsic property of the strength distribution of the lithosphere, or alternatively exhibits a bias due to regions where the background stress is anomalous and effects the neutral plane depth. End-member perspectives can be found in Chapple & Forsyth (1979) and Mueller et al. (1996b). In the former, the apparent neutral plane depth is assumed to be the intrinsic depth for old lithosphere (i.e. the depth produced in regions of isotropic background stress). In the latter study, reverse seismicity only occurs where the background stress is anomalously compressional, and represents a localized deviation from the background trend.

It is logical to think that if there are significant global variations in the background stress in the subducting lithosphere, more reverse seismicity will occur in the regions with the highest compression (or at least the smallest extensional component). This is because, relative to a ‘average’ background stress state, a more compressional state will result in a shallower neutral plane, increasing the volume that exists between the neutral plane and the brittle–ductile transition (i.e. the volume of rock where brittle compressional rupture can occur). As such, it is seems plausible that the apparent normal-reverse seismicity transition may biased by regions where the background stress is more compressional (which may simply mean less tensional). In our view, however, the argument that outer rise reverse seismicity is restricted to regions of anomalous resistance, where flexure is overwhelmed by net compression (Mueller et al. 1996b), is inconsistent with a range of subsequent detailed observations (Todd & Lay 2013; Obana et al. 2014; Zhao et al. 2022).

3.2 Observed moment data

The observed moment is the torque around a given point, due to the load of the flexural (uncompensated) topography seaward of that point. The derivation of this relationship is outlined in Section S1. and is independent of the details of the particular rheological model under consideration.

Fig. 1 shows moment data from several previous studies. In the case of McNutt & Menard (1982); Levitt & Sandwell (1995) these are the observed moment ‘sensu stricto’, and reflect the value at the first zero crossing. In the case of Garcia et al. (2019), they represent the tabulated bending moment at maximum curvature from the flexural inversion (which utilizes a yielding rheology). Typically these values were within 5–10 per cent of the saturation moment for the given rheology and assumptions. In the case of Hunter & Watts (2016), the moment data were redigitized from the bending moment profiles in the original study, estimated at the first zero crossing. These data reflect the variable elastic plate inversion (A2) from Hunter & Watts (2016).

Our understanding of rock deformation mechanisms suggests that bending moments should be close to saturated at curvature values encountered in the trench slope (Goetze & Evans 1979). Ideally, moment data would reflect the saturation moment of the lithosphere, in which case different data sets could be compared without additional correction for curvature (which often has significant associated uncertainty). As shown in the Fig. S7, for the majority (∼70 per cent) of the profiles in the Levitt & Sandwell (1995) data, the curvature estimates imply that the plates have reached >75 per cent moment saturation by the first zero crossing (see Section S5). This is also consistent with the results of the numerical model discussed in Section 5. Hence, most of the data in Levitt & Sandwell (1995) should exhibit less than 25 per cent moment deficit; the fact that the scatter is significantly larger than these considerations suggest, is a key issue in the current study.

4 METHODOLOGY

4.1 Subduction models—overview

Two modelling approaches are used to make quantitative predictions related to the bending of the lithosphere during subduction. These are referred to as (1) YSEs and (2) numerical subduction modelling. The results from the numerical subduction modelling inform certain aspects of YSE construction (e.g. strain rate patterns), as well as the interpretation of moment data (e.g. moment deficit). YSEs are used to produce ‘forward models’ of flexural behaviour, in order to compare with observations. In both modelling approaches, the constitutive behaviour follows a ‘classical strength model’, which combines brittle faulting with both low and high temperature plasticity (as described in Section 4.2). Fig. 2 shows an example of the YSE approach. The implementation details are described in detail in Section S4.

The numerical subduction model was developed with the ASPECT code (Bangerth et al. 2020). The model simulates the evolution of an idealized 2-D subduction system, or the one-sided sinking of a thermal boundary layer due to thermal buoyancy. The implementation details are described in detail in the Section S2. All parameters for the numerical model are provided in Table S1. Fig. 5 shows stress and deformation patterns in the numerical subduction model, described in Section 5.

In recent studies, the applicability of laboratory-derived flow lows, particularly those describing low-temperature plasticity, has been an important issue (Hunter & Watts 2016; Bellas et al. 2022). This question is strongly dependent on the assumptions made in the model, in particular the temperature and strain rate. YSEs are often constructed under the assumption of a constant strain rate (e.g. Hunter & Watts 2016; Garcia et al. 2019). However, we assume that strain rates vary with distance from the neutral plane, a choice that is consistent with deformation patterns in the numerical model (and described in Section S5). This assumption means that for given flow law parameters, we compute higher creep stress than those of Hunter & Watts (2016) and Garcia et al. (2019). Fig. S6 highlights these differences.

The YSEs are based on a plate cooling model, with parameters derived from Parsons & Sclater (1977) (and identical to Hunter & Watts 2016; Garcia et al. 2019). For 100 Myr lithosphere, the inferred depth limit for well-constrained earthquakes is ∼50 km depth, and corresponds to a temperature of ∼700 °C. The Parsons & Sclater (1977) model remains relatively consistent with more recent approaches, incorporating the joint inversion of plate subsidence and heat flow, and based on temperature and pressure dependent thermal properties (Richards et al. 2018). The temperature profile is warmer than assumed in Bellas et al. (2022), due to the higher diffusivity in that study [1.0, rather than 0.8 (⁠|$\times 10^{-6}\, \rm {m^{2}\,s^{-1}}$|⁠)]. This higher diffusivity is equivalent to a 25 per cent increase in apparent age (i.e. |$\sqrt{\kappa t} =$| constant). For typical strength models discussed in this paper, this would produces a ∼25 per cent increase in saturation moment. A comparison of temperature models is shown in Fig. S8.

4.2 Constitutive behaviour

4.2.1 The classical strength model

The classical strength model of the oceanic lithosphere combines elastic deformation, brittle deformation (frictional sliding on optimally oriented faults) and separate expressions for low temperature (aka LTP or Peierls) and high temperature (aka power law or dislocation) creep (e.g. Goetze & Evans 1979). Despite the success of this model, remaining sources of uncertainty include: (1) the significant variability in published experimentally derived equations for LTP; (2) significant variations in inferred friction coefficient based on different geophysical studies and (3) strength properties in the semi-brittle regime.

Semi-brittle behaviour is thought to operate throughout a significant depth, within the strongest part of the lithosphere (Ohnaka 2013). A defining characteristic of semi-brittle behaviour is that strength has both strong pressure and temperature sensitivity. It is often assumed that semi-brittle deformation will lead to a zone in which typical differential stresses are smaller that the peak stress given by the intersection of the brittle and ductile stresses (Ohnaka 1992, 2013; Molnar 2020). The paucity of constraints on stress–strain rate–temperature relationships in this domain means that the semi-brittle regime is routinely neglected in YSE modelling, which is also the case here. However, we note that the typical depth range of reverse faulting (typically 30–50 km for ∼100 Myr lithosphere) overlaps the anticipated semi-brittle domain (Molnar 2020).

Moreover, the brittle strength of rocks is not restricted to the frictional strength of pre-existing faults. Other important processes include the saturation of frictional strength at high confining pressure, the temperature dependence of the fracture strength, changes in microcrack behaviour and thermal shear instability (non-Coulombic faulting, Stesky et al. 1974; Shimada et al. 1983; Renshaw & Schulson 2004; Ohnaka 2013; Molnar 2020). Such processes may even lead to negative brittle strength-depth gradients (Ohnaka 2013). It is plausible that outer rise reverse faulting at ∼25–50 km, may involve some of these processes.

4.2.2 Frictional sliding

The Coulomb strength criterion states that frictional strength on sliding surfaces is linearly related to the effective normal stress :

where μ is the friction coefficient. For faults of optimal orientation, the strength can be expressed in terms of the limiting value of the differential stress σ_xx = (σ₁ − σ₃) as a function of depth (y):

where ρ is density, g is gravity and λ is the pore pressure factor. The sign convention used here is that positive differential stress is extensional; the upper signs (±) in eq. 2 apply to normal faults and the lower signs for reverse faults.

Predictions of Coulomb theory, incorporating standard laboratory-derived frictional coefficients (e.g. 0.6–0.8), and hydrostatic pore-pressure, are generally consistent with measured differential stress in the shallow upper crust (<1 km) (Townend & Zoback 2000). Certain rock types and mineralogies (e.g. shales, clays and serpentines) can exhibit much lower friction coefficients (e.g. 0.1–0.3, Ewy et al. 2003).

4.2.3 Intracrystalline creep

At high temperatures (⁠|$T \gt 0.7 \, T_{\textrm {melt}}$|⁠) and low stresses (<100–200 MPa), olivine creep exhibits a power-law relationship between strain rate and differential stress, known as dislocation creep. These conditions are relevant to the lower lithosphere and asthenosphere. The form of the power-law creep is:

where |$\dot{\varepsilon }$|⁠, Δσ, T, P, R, are the strain rate, differential stress, temperature (K), pressure (Pa) and gas constant; A is a pre-exponential factor, E is the activation energy, V is the activation volume.

At lower temperatures (⁠|$T \lt 0.6 \, T_{\textrm {melt}}$|⁠) and higher stresses, relevant for the shallow mantle lithosphere, a stronger dependence of the strain rate on stress is observed, which is usually described by an exponential law, known as low-temperature plasticity (LTP), or Peierls creep.

The form of the LTP law generally used is:

where several additional parameters appear: σ_P is the Peierls stress constant, and p, q are additional parameters that depend on the geometry of the impediments to dislocations (Mei et al. 2010).

5 INSIGHTS FROM NUMERICAL MODELS

The numerical subduction modelling approach was briefly introduced in Section 4.1, while the implementation is described in detail in Section S2. This approach, in which stress and strain (rates) emerge dynamically, provides important and complimentary insights alongside the kinematic-based YSEs. This section provides a summary of those insights, while a more detailed description of the analysis is contained in Section S3.

5.1 Contributions to subducting plate topography

Fig. 4 shows the subducting plate topography from the numerical model, once the stress state, geometry and velocity of the subduction hinge has reached quasi-steady state (5 Myr after start of simulation, same time step as shown in Fig. 5). We note that the numerical model has a free-surface, and that the model topography (solid red line, Fig. 4) represents the deformed uppermost elements.

The analysis of model topography is described in detail in Section S3, and demonstrates that the subducting plate topography has three primary contributions. The first is the isostatic adjustment related to the density of the cooling plate, accounting for about ∼2.5 km of topographic variation from the ridge to the trench. The second is the dynamic topography, mainly associated with anomalous ‘dynamic pressure’ variation in the asthenosphere, producing about ∼0.5 km of topography. The (positive) pressure anomaly in the asthenosphere is consistent with the shallow return flow to the ridge (Schubert et al. 1978), as well as the accommodation slab rollback; similar levels of topographic response are observed in 3-D models (Holt 2022).

The third contribution, is the flexural topography, that is the topography related to gradients in shear stress resultant (shown in detail in Fig. S4). The black line in Fig. 4 shows the residual (flexural) topography, once the isostatic and dynamic pressure contributions are removed.

5.2 Stress and strain rate patterns

Fig. 5 shows stress and deformation patterns in the numerical model. The scalar fields show downdip stretching/shortening, and downdip (effective) tension/compression, as described in the figure caption. In the subduction hinge there are two main regions of deformation: in the subducting plate immediately seaward of the trench (the ‘outer rise’), where curvature is increasing, and in the unbending zone, where curvature reduces.

Stress profiles at four locations are shown in the right hand panel of Fig. 5. The blue line labelled x₀ is the first zero crossing, based on analysis of the flexural component of the topography. The black line is located at the point of maximum bending moment [labelled M_b(max)]. The green line shows location of partial slab unbending; here the stress profile has a zig-zag character due to elastic effects (cf. Engdahl & Scholz 1977; Mueller et al. 1996a). The red line is located where the stress pattern has undergone full reversal compared with the outer-rise. Note that when the deformation is dominated by analastic processes (yielding), unbending is a stress-inducing process.

Although the model is clearly dominated by bending-related stress, it also provides important insights in terms of ‘background’ stress patterns. In particular, the model exhibits systematic variations in the net axial force (F_net), due to the gravitational potential energy associated with the flexural topography. To demonstrate, we consider the 2-D depth-integrated horizontal force balance (e.g. Fleitout & Froidevaux 1983; Ghosh et al. 2009):

The overbars represent depth integration. The term on the left hand side of eq. (5) is approximately equal to F_net, the net axial force, as discussed in Section S1. The first term on the right hand side is the GPE and the second term is the basal shear stress. Hence, eq. (5) states that gradients in net axial force will be related to gradients in the gravitational potential energy (GPE), as well as the shear stress acting on the base of the plate. In the flexural domain, horizontal gradients of the vertical shear stress resultant cannot be neglected, and it is not valid to assume that the vertical stress (σ_yy) is the lithostatic pressure (P_L).

Fig. 6 shows the variations of the terms in eq. (5), based on the numerical model results (at the same time step as shown in Fig. 4). The values were calculated by interpolating and numerically evaluating the integrals from the model ridge height to a depth of 150 km. Arbitrary constants are added to the basal shear stress and GPE terms to facilitate comparison with the net axial force (F_net shown in black).

Toward the right-hand side of the region shown in Fig. 6, at ∼1000 km from the trench, gradients in the GPE are negligible, and the variations in F_net is balanced by the integrated basal shear stress. Near to the trench, F_net varies much more dramatically, and is closely correlated with the change in GPE. Even though the flexural topography is supported by the strength of the lithosphere, gradients in the vertically integrated stress (the GPE) are still important. In our model, the implied contribution of this ‘trench GPE’ is about ∼2.0–2.5 TN m⁻¹. It is interesting to note that this estimate of trench GPE has a very similar magnitude to the long-wavelength isostatic GPE attributable to the mid-ocean ridges relative to old oceanic lithosphere (Coblentz et al. 2015).

The trench GPE values we infer are about a third of the magnitude (∼7 TN m⁻¹) reported by Bessat et al. (2020), by applying the density moment integral to similar numerical subduction models. However, our results suggest that the density moment would significantly overestimate the trench GPE, due the non-isostatic nature of topography. The key insight from our model is that the net axial force varies in a way that is consistent with trench GPE, this leads to a ‘less tensional’ background stress state near the trench, where the net axial force almost disappears.

Strain rate patterns in the numerical model are shown in Fig. S5. Within the strong part of the plate, strain rates show a linear increase away from the neutral plane, in agreement with the simple bending approximation. The bending rates approximately double between the first zero crossing (x₀) and the point of maximum bending moment [M_b(max)], indicating that both curvature and curvature gradient are increasing between these two locations. When numerical model convergence velocities (14.5 cm yr⁻¹) are re-scaled to a reference value of 10 cm yr⁻¹, we find that strain rates of ∼10⁻¹⁵ s⁻¹ are typical of the strong part of the lithosphere undergoing ductile creep (40–50 km).

5.3 Moment deficit

An important question related to the interpretation of observed moment data is the potential amount of deficit between the bending moment at the first zero crossing, and the saturation moment (e.g. Levitt & Sandwell 1995). As shown in Fig. 5, the bending moment at x₀ is about 85 per cent of the maximum bending moment [at M_b(max)]. The moment deficit has been assumed to relate primarily from the potential increase in curvature between the two locations (Levitt & Sandwell 1995; Garcia et al. 2019). Using YSEs calibrated to the same conditions as the numerical model, and knowing the strain rate change between the two points, we find that about one-third of the change is due to strain rate increase, whereas the remaining two-third is due to the increase in curvature and the reduction of the elastic core.

6 INTERPRETATION AND COMPARISON OF BENDING MOMENT DATA

6.1 Defining lithosphere strength

In this section the YSE approach is used to develop a simple ‘family’ of strength models with varying saturation moment. These strength models are used to: (1) define regions in age-moment space (the weak, intermediate and strong regions in Fig. 1) and (2) help discuss and interpret different existing age-moment data sets (which is the topic of the next section).

The family of strength models is created by generating a strong ‘reference’ model, and then simply scaling this model to achieve weaker variations. The reference model is constructed using the ‘Set 1’ flow law parameters given in Table 2. We then choose the friction coefficient so that a 30 km neutral plane is achieved for zero net axial force (or an isotropic background stress) for 100 Myr lithosphere. In this exercise we assume zero pore fluid pressure for simplicity. The required friction coefficient to generate the 30 km neutral plane is unrealistically high (μ = 1.5). However, the aim at this stage is simply to provide an upper limit for subsequent analysis. This reference model is shown by the outermost model in the left hand panel of Fig. 2. The upper bound to the ‘strong’ region (shown in red) in Fig. 1, represents the variation of strength with age for the reference model.

Having constructed this reference model, we simply scale the YSE (or the age-moment relationship) to achieve an arbitrary strength (saturation moment). Based on this scaling approach, we describe three strength regions, which are defined by the value of the saturation moment for 100 Myr lithosphere: weak models have saturation moments (<1); intermediate models (1–2) and strong models (>3), measured in units of N-m m⁻¹ (N) × 10¹⁷. These regions are shown in green, tan and red, in Fig. 1 and in Fig. 2. We also use this scaling approach as a basis for fitting observed moment data (as we describe in the next section). The solid blue line in Fig. 1 is the result of such an approach, that is it is simply a scaled version of the original model discussed above.

The solid black and red lines in Fig. 1 and the right hand panel of Fig. 2, show two different ‘intermediate’ strength models from recent subduction-flexure studies. The black line shows the preferred model from Garcia et al. (2019). The red line shows one of the models that was determined to be consistent with flexural inversion from the study of Hunter & Watts (2016) (they find that several combinations of flow law parameters and friction coefficients are viable). We use an identical thermal model to these studies. We note that the age variation of strength models is rather insensitive to the choice of parameters. The solid blue, black, and red lines in Fig. 1 have virtually identical moment-age variation, although having significant differences in strength-depth patterns, particularly in the LTP domain (as shown in Fig. 2).

6.2 Bending moment estimates: interpreting scatter

The object of this section is to bring together estimates from a range of previous studies for the bending moment of old lithosphere; ideally these estimates represent locations close to moment saturation, an assumption that is assessed in the analysis. A key focus is how the data set of Levitt & Sandwell (1995) is interpreted in terms of competing strength models.

The observed moment data of Levitt & Sandwell (1995) show an age dependency, as well a significant component of scatter (e.g. Fig. 1). We want to interpret this data in terms of competing lithospheric strength models, which in turn requires a model for the scatter. We expect there to be three major contributions: (1) real variations in the saturation moment for lithosphere of given age (for instance, due to variations in the net axial force, thermal perturbations or strain rate variation(s); (2) a deficit between the bending moment at the point being measured (e.g. the first zero crossing) and the saturation moment (‘moment deficit’) and (3) estimation errors. The moment deficit effect can only lead to underestimates. The other types of effects could in principle lead to random or systematic bias in the scatter.

Levitt & Sandwell (1995) assumed that the scatter primarily derives from a moment deficit at the first zero crossing. Hence, they appealed to a strength model in which the age-moment curve provides a ceiling to the observed moment data. Models in the ‘strong’ range tend to satisfy this interpretation. Their preferred model has a saturation moment of about 2.5 × 10¹⁷ N, for 100 Myr lithosphere, and is shown with the dashed blue line in Fig. 1. Strength models in the ‘strong’ range are consistent with relatively high friction coefficients, comparable to those determined by experiments (Byerlee 1978); depending on assumptions regarding the thermal model and strain rate, they are also consistent with a variety of LTP laws.

However, the interpretation of Levitt & Sandwell (1995) implies that in many cases, at the first zero crossing, the plate is considerably ‘undersaturated’. For instance, many of the observed moment values of <1 × 10¹⁷ N for lithosphere >100 Myr, would imply bending moments of about ∼1/3 of the saturation moment. As shown in Fig S7, for the curvature estimates corresponding to the moment data, all profiles are expected to have reached >65 per cent moment saturation, with the majority (>70 per cent) being >75 per cent saturated. This distribution is consistent with the behaviour observed in the numerical model, where the bending moment at the first zero-crossing (x₀)is about 85 per cent of the maximum bending moment. Additionally, Fig. S7 shows that many of points that show the highest deficit have high curvature. Indeed, the trend of the moment deficit with curvature is opposite to what would be expected. Hence, we argue: that (1) where the lithosphere has a moment deficit at the first zero-crossing point, this deficit will typically be less that 25 per cent of the saturation moment and (2) moment deficit is inconsistent with the trends of the scatter in the data of Levitt & Sandwell (1995). This implies that estimation errors also play a significant role.

If we assume that the combined effect (both real effects and error) is to produce random (or at least not heavily biased) scatter, then we require a rheological model that minimizes the misfit. For the family of strength models described in the previous section, the misfit is minimized (RMS) with a model that has a friction coefficient of μ = 0.25. This model is shown with the solid blue line in Fig. 1. It lies in the lower half of what we call the intermediate strength region, and exhibits a typical saturation moment close to 1.5 × 10¹⁷ N, across the range of old lithosphere (100–150 Myr).

6.3 Bending moment estimates: comparison of different data sets

In the previous section we suggested that a combination of effects (both real effects and error) could lead to to more-or-less random scatter in the Levitt & Sandwell (1995) data set. Additional studies provide an important test of this interpretation.

As shown in Fig. 1, there is a considerable consistency between the studies of Chapple & Forsyth (1979), Hunter & Watts (2016) Garcia et al. (2019) and our reinterpretation of the Levitt & Sandwell (1995). These studies each infer strength models in the mid-intermediate strength range, and lie within about ∼20 per cent of the best-fitting strength model we developed for the Levitt & Sandwell (1995) data. In terms of comparing these data sets, there are two additional features that are important to highlight.

First, we note that moment estimates from Hunter & Watts (2016) (red circles Fig. 1) represent ‘trench-bin’ averages: each of these points contains information averaged over a typical distance of several thousand kilometres. On the other hand, the data points from Levitt & Sandwell (1995) represent flexural inversion over single bathymetry/gravity profiles. The fact that an average (or a best fit) interpretation of the Levitt & Sandwell (1995) data, leads to very similar moment estimates to the ‘trench-bin’ results of (Hunter & Watts 2016), helps clarify the connection between these two studies.

Secondly, we note that the results (moment estimates) of the aforementioned studies are generated under a range of different modelling approaches. For instance, these inversions used a range of forward flexural modelling techniques: constant-thickness elastic plate (Levitt & Sandwell 1995), variable thickness elastic plates (Hunter & Watts 2016), non linear flexure with ad-hoc yield limits (Chapple & Forsyth 1979) and non-linear flexure based on classical yield strength model (Garcia et al. 2019). These modelling methodologies are summarized in Section 2. This suggest that a range of different approaches are able to consistently determine that the bending moment of the lithosphere saturates in the ‘intermediate’ strength range.

What is also clear in the compilation of data in Fig. 1, is that only the data set of Levitt & Sandwell (1995) has been used to advocate a model in the ‘strong’ region. However that conclusion was tied to a specific interpretation of the scatter, which has proven to be inconsistent with results from additional studies (Hunter & Watts 2016; Garcia et al. 2019). Of the data compiled in Fig. 1, the only results are not consistent with intermediate strength are those of McNutt & Menard (1982) (green points). These estimates mainly lie in upper part of the weak region. As in Levitt & Sandwell (1995), these estimates represent results from individual profiles, but unlike Levitt & Sandwell (1995) they are based on modelling bathymetry data, rather than joint inversion that includes gravity. This is expected to make them less reliable. However, as we noted the Levitt & Sandwell (1995) data also contains points (outliers) in the weak strength region. These do not appear in trench-averaged approaches (e.g. Hunter & Watts 2016). For these reasons we think it is appropriate to disregard the estimates of McNutt & Menard (1982). However, to the extent that there is valid information in the McNutt & Menard (1982) data, it would imply an even weaker strength model.

The model we generated to fit the Levitt & Sandwell (1995) data has a low effective friction coefficient (μ = 0.25, assuming zero pore fluid pressure). We note that a range of recent studies have converged on similar values for the frictional strength (Garcia et al. 2019; Pleus et al. 2020; Bellas et al. 2022), as shown in Fig. 2. The models of (Hunter & Watts 2016; Garcia et al. 2019) both assumed hydrostatic pore pressure, which offsets the slight higher friction coefficient (μ = 0.3), relative to the best-fitting model shown in Fig. 1 (μ = 0.25). While this convergence towards low frictional strength is worth highlighting, there are important trade-offs between strength in the brittle and ductile regimes. As we discuss in the following section, this means that the frictional strength is not as well constrained as it might appear on the basis of the results presented in this section.

7 STRENGTH AND NON-UNIQUENESS

7.1 Relative strength and background stress

Different combinations of brittle and ductile strength can be combined to produce a strength model with the same saturation moment (e.g. an intermediate strength model). The depth of the neutral plane seems to offer the most direct means of trying to constrain the relative strength in each regime, as was pointed out in Chapple & Forsyth (1979). However, this is hampered by the additional influence of background stress. This leaves us with two main sets of observational constraints: those which reflect the integrated plate strength (e.g. the saturation moment) and the neutral plane depth. However, there are at least three key ‘variables’: brittle strength, ductile strength and the net axial force (factors like the convergence rate, orientation of seafloor fabric, out-of plane stress, are expected to provide additional controls). We can think of these key variables as a parameter space. Together they lead to predictions of the saturation moment as well as the neutral plane depth (as well as additional predictions that can be compared to the seismic expression). Despite the non-uniqueness, it may be still be useful to try to investigate different regions of the parameter space, particularly in reference to additional constraints.

Strength models proposed in previous studies span a wide range within the parameter space. We highlight four such models to illustrate the problem. These are shown in Fig. 7, evaluated for 100 Myr lithosphere. We refer to these as models A, B, C, D, following the labels on the Figure axes. Three of these these models (A, B, C) utilize variations of the LTP formulation originally described in Goetze & Evans (1979), the (‘Set 2’ parameters shown in Table 2). The ductile strength in the HTP/LTP regions is varied through changes the activation energies, as described in the figure caption. In all cases the LTP activation energy lies in the uncertainty range specified by the original study (100–130 kcal mole⁻¹, or 418–544 kJ mole⁻¹). In this sense, all of the models within Fig. 7 have acceptable ductile strength, within the range of uncertainty given in Goetze & Evans (1979).

Each of the models shown in Fig. 7 have intermediate strength, but the relative strength in the brittle and ductile regions is different and hence the intrinsic neutral plane depth varies. Therefore, the models imply differences in terms of the background stress required to facilitate regions with a 30 km neutral plane (which produce the apparent neutral plane depth shown in Fig. 3). Each model is shown with different levels of net axial force, as labelled in the legend shown on Fig. 7(b).

The maximum seismogenic limit for each of the models is shown with T_s. This estimate is based on different assumed mechanisms for brittle failure in compression. The estimate is highly speculative for Models C and B and is based on a ‘stress amplification’ model described in the following section.

We will try to evaluate the consistency of these different strength models in light of additional constraints. The constraints primarily relate to the seismic expression of outer rise bending, as inferred from the distribution of well-constrained centroid depths (e.g. Fig. 3), and following a similar logic to Chapple & Forsyth (1979). The key features we refer to are the depth limit of normal faulting (∼30 km); the maximum depth of limit of (reverse) earthquake nucleation (∼50 km); the capacity to predict a neutral plane depth of ∼30 km within assumed variations of background stress. These constraints are relevant to older plates, particularly the range 80–120 Myr and are shown schematically with the blue and red bars in Fig. 7. The general pattern of strength increase with age leads us to expect that the neutral plane depth would increase by about 5 km for very old lithosphere (∼150 Myr), as shown in the top panel of Fig. 1.

We differ from Chapple & Forsyth (1979) in that we do not assume that these features are characteristic of an isotropic background stress state; rather we argue that a consistent strength model should be capable of generating these features within reasonable variations in background stress. Several other constraints are also referred to: independent constraints on the frictional strength of outer rise normal faults; the capacity to explain the anomalously shallow neutral plane (∼25 km) along Tonga–Kermadec (e.g. Fig. 3); as well as a priori assumptions about the background stress state. Although Fig. 7 shows specific strength models, the discussion here is intended to represent, somewhat more broadly, the regions of the parameter space represented by these models.

7.2 Brittle failure in compression

Before we evaluate the strength models, we need to address the issue of the mechanism of brittle failure in compression. Model A is the only model that can directly explain this brittle failure in terms of frictional sliding. In model D, the brittle behaviour in the lower part of the seismogenic zone (>20 km) is assumed to reflect processes other than frictional sliding, and the magnitude of the brittle strength is selected try to match the bending moment inferences (following Chapple & Forsyth 1979).

Models B and C do not provide an explicit mechanism for brittle failure beneath the neutral plane, where the strength is instead limited by LTP. This is a common feature of most classical strength models, although it has seldom been viewed as a critical limitation. Several studies have taken a different view, and insisted that a consistent strength model needs to account for brittle failure in the depth range where outer-rise compressional earthquakes are observed (Chapple & Forsyth 1979; Mueller et al. 1996b). We tend to agree with this perspective. Even if there is a significant depth range in which the constitutive behaviour is semi-brittle, it does not seem satisfactory to simply assume that the stress in a semi-brittle domain should be equivalent to the LTP strength—a rheological process that is itself unrepresentative of brittle failure, even if it may coexist with brittle failure mechanisms transiently.

One possibility was recently suggested by Toffol et al. (2022), termed a ‘stress amplification’ model. Although their study is focused on intermediate depth seismicity, it also has relevance for reverse faulting in the outer rise. The model proposes that weak domains of hydrated peridotite can create large stress amplifications which are capable of exceeding the frictional strength, even at large confining pressure. An interesting feature of this model is that although it can produce large stress amplifications, the ‘flexural characteristics’ (neutral plane depth, bending moment) appear to be insensitive to the perturbations (see e.g. Toffol et al. 2022, fig. 5). This is unsurprising given flexural strength depends on the integrated stress state (resultants, moments) both across the plate and laterally. Stress amplification can therefore provide a solution as to why brittle behaviour might occur within a rock mass where the macroscopic strength (e.g. the saturation moment) is determined by the LTP strength.

Stress amplification requires that strain rates can increase above the bulk plate-scale value around weak inclusions. We do not know what strain rates are appropriate, but we note that the amplification magnitudes reported in Toffol et al. (2022) vary from about four to five times. To characterize a seismogenic thickness (T_s) for model B and C, we assume that stress amplification can increase the LTP stress locally by five times, and define the T_s is the intersection of the perturbed stress and the frictional strength. As in Toffol et al. (2022), we have neglected the fact that fracture strength will place an upper limit on frictional strength, with 1 GPa being a typical limit strength based on experiments in igneous rocks in the 500–600 °C range (Stesky et al. 1974; Shimada et al. 1983; Ohnaka 2013).

We adopt this stress amplification model as a tentative solution to brittle failure in compression for Models C and B. For model A, T_s(1) refers to the seismogenic thickness based on the background stress state, although T_s(2) refers to the seismogenic thickness based on a stress amplification model.

7.3 Evaluating strength models

Model A is very similar to the preferred model proposed by Garcia et al. (2019). The friction coefficient is low, consistent with the inferred dip angles of newly formed outer rise faults (Craig et al. 2014b), and a relatively high pore fluid pressure is assumed (λ = 0.36). The intrinsic neutral plane depth is about 37 km. A compressional net axial force of −5 TN m⁻¹ results in a neutral plane depth of 32 km (close enough to be considered consistent with the apparent neutral plane depth). Model A provides an explanation for brittle failure beneath the neutral plane in terms of frictional sliding, although the maximum depth of frictional behaviour is about 40 km [T_m(1)] and hence additional mechanisms of brittle failure may be required to facilitate earthquakes to ∼50 km. A stress amplification model could explain these depths [T_m(2) ∼50]. In regions of of moderate slab pull (5 TN m⁻¹), model A predicts normal faulting to 42 km, significantly deeper than the 30 km limit that is typically reported; this would be even deeper if the strength model was evaluated for lithosphere older than 100 Myr. Meanwhile, model A implies that regions that exhibit a ∼30 km neutral plane are in moderate net axial compression (−5 TN m⁻¹). This could be considered problematic in light of a priori assumptions about tectonic stress in subducting plates. Finally, even with moderate compression, the neutral plane (32 km) is significant deeper than inferred across Tonga–Kermadec (∼25 km).

Model B is based on a strength model proposed by McNutt & Menard (1982). In this model the HTP/LTP strength is reduced by changing the activation energy, while the frictional strength remains relatively high (note also that there is no pore pressure in this model). The intrinsic neutral plane depth is about 26 km. A (effective) tensional net axial force of 5 TN m⁻¹ results in a neutral plane of about 30 km. Hence, model B predicts a 30 km neutral plane under conditions that are consistent with inferences of moderate net slab pull (Bird et al. 2008). Clearly much higher estimates are also present in the literature (e.g. Conrad & Lithgow-Bertelloni 2002). The seismogenic depth (based on the stress amplification model) in model B is very low T_s = 35 km. This is difficult to reconcile with the observation of earthquake hypocentres to at least 50 km in lithosphere of 80–120 Ma. Secondly, the brittle strength (μ = 0.6, λ = 0.0) implies dry faults, which seems unlikely given geophysical inferences of outer rise hydration (Ranero et al. 2003), as well as the preponderance for in situ intraplate stress measurements to conform to hydrostatic strength expectations (Zoback & Townend 2001).

Model C is similar to the one presented in the study of Goetze & Evans (1979) (see caption for details). It balances strength in the brittle and ductile domains so that the intrinsic neutral plane depth is ∼30 km. The assumption is hydrostatic Byerlee strength. Within the assumed range of background stress variations, model C remains consistent with the observed maximum depth of outer rise normal faults. The friction coefficient is typical of laboratory experiments, in situ stress measurements, and some observations of outer rise fault dip angles (Obana et al. 2023). However, it is significantly higher than the seismological constraints of Craig et al. (2014b). For Model C, the proxy for the seismogenic thickness (T_s) is 41 km, which is smaller that the maximum depth extent of outer rise seismicity, but not as dramatic as the disagreement in model B. We note a trade-off is emerging here: models which more accurately reproduce the maximum depth of normal faulting, tend to also have a T_s that is too shallow. Allowing for moderate compressional stress, model C predicts a neutral plane of 27 km, fairly close to the inferred depth along Tonga-Kermadec.

Model D is a variation on the ad hoc yield-strength model of Chapple & Forsyth (1979). It is therefore automatically consistent with some of key constraints, such as the maximum depth of seismicity (T_s is 50 km by construction). In the original study, the differential stress is truncated at 100/600 MPa, whereas in Model D, the truncation is 100/400 MPa. This provides a bending moment similar to the other models shown in Fig. 7. The model was originally presented and discussed in the context of an isotropic background stress state. It is useful, however, to consider the consistency of the model in relation to assumed variations in background stress. Overall, model D has the greatest sensitivity in terms of the variation of neutral plane depth as function of background stress. Under moderate net axial compression, model D can explain the 25 km neutral plane depth inferred at Tonga-Kermadec. However, like model A, under moderate net axial effective tension, the predicted depth of normal faulting (37 km) becomes deeper that is typically observed.

In our view Model B, based on McNutt & Menard (1982) is the least consistent in view of additional constraints. Fundamentally, the ductile strength is so low that seismic rupture beyond about 35 km seems unlikely. An implication of rejecting model B, is that regions that give rise to the apparent 30 km neutral plane depth cannot be supporting a net axial force of ∼ 5 TN m⁻¹ (or higher). We expand on this argument in the following section.

None of the other models (A, C, D) is completely consistent with additional constraints. However, rather than advocate a preferred model, we think it sufficient to summarize the predictions of each models, and the areas of (in)consistency, under the assumptions we have outlined. Our main intention here is to highlight the constraints that can be brought to bear on different strength models, and the challenges that remain in successfully explaining the deformation mechanisms, brittle behaviour, and relative strength properties (Chapple & Forsyth 1979).

Of the models shown in Fig. 7, model A is unique in that it offers a direct explanation for outer rise reverse faulting in terms of frictional sliding. Previous studies have argued that a low friction coefficient would generate an inconsistent strength model and therefore could not be invoked to explain brittle failure in compression (Chapple & Forsyth 1979; Mueller et al. 1996b; McNutt & Menard 1982). However, this mechanism is not able to explain earthquake nucleation to the full extent of the observed depth range ∼50 km.

We invoked the stress amplification model of Toffol et al. (2022) to try to quantify a relevant seismogenic depth in models where the macroscopic strength is governed by LTP (e.g. B and C). Even under rather lenient assumptions (stress amplification of a factor of 5, the maximum reported in Toffol et al. 2022) the seismogenic thickness remains too shallow. Interestingly, the stress amplification model works more successfully for model A, which is again due to the low friction coefficient assumed in that model.

Overall, the presence of deeper compressional seismicity, and particularly the capacity for earthquake ruptures to nucleate (most clearly illustrated by the smaller-scale seismicity shown on Fig. 3) remains difficult to reconcile with existing rheological formulations for any of the models shown. Resolving this will require both additional observational data, to better understand the spatial distribution of this deeper compressional seismicity, and additional work on the nature of seismicity close to the brittle to ductile transition (e.g. Aharonov & Scholz 2019; Toffol et al. 2022).

8 DISCUSSION

8.1 Constraints on net slab pull

In the previous section we suggested that regions of old lithosphere which exhibit a 30 km neutral plane cannot be supporting ‘typical’ levels of net slab pull (e.g. 5 TN m⁻¹). This reasoning was based on the inconsistent properties of a strength model that could give rise to a 30 km neutral plane depth under this background stress state. Now we consider constraints on net slab pull more generally.

The net axial force is given by the integral of the differential stress across the thickness of the lithosphere (e.g. eq. S4). In simple concave-down bending, net slab pull must be transmitted entirely through the part of the lithosphere experiencing effective tension (i.e. above the neutral plane). In this region, the differential stress is expected to be controlled by fault frictional strength. Integration of eq. (2), to the depth of the neutral plane provides and upper limit on the net slab pull, which is therefore limited by the square of that depth. Observations of neutral plane depths, in conjunction with estimates of the brittle strength, imply limits on net slab pull. Assuming a friction coefficient of 0.6, along with hydrostatic pore pressure (λ = 0.3), with zero cohesion, leads to a maximum F_net = 6.9 and 12.3 TN m⁻¹ for 30 and 40 km neutral planes, respectively.

Of course the above estimates are upper limits, as they factor in no compressional stress beneath the neutral plane. Observations of seismicity, as well as laboratory estimates of olivine creep, provide compelling indication for significant strength to at least 600–700 °C, or ∼50 km (for 100 Myr lithosphere). A simple assumption we can make is that the ratio of net force in the extensional and compressional parts of the envelope scales with ratio of depth above and below the neutral plane, to 50 km. In this case, estimates of the maximum net slab pull reduce to 2.3 and 9.2 TN m⁻¹ for 30 and 40 km neutral planes in ∼100 Ma lithosphere.

If the 30 km neutral plane is assumed to be globally representative for old lithosphere (e.g. Chapple & Forsyth 1979), we could rule out net slab pull being greater than a few TN m⁻¹. However, if the apparent neutral plane depth is biased by regions of anomalous compression, then slab pull may be focused along other regions where the neutral plane is deeper (although likely to be unconstrained due to the absence of reverse seismicity). Nevertheless, this simple calculation shows the geodynamic relevance of constraining the neutral plane depth, and therefore the importance of earthquake hypocentre accuracy (Craig et al. 2014a).

8.2 Intraplate stress variations: potential implications for outer rise seismicity

A characteristic of outer rise reverse seismicity is the uneven distribution of events relative to the shallower normal faulting (Mueller et al. 1996b; Craig et al. 2014b; Ye et al. 2021). Global compilations reveal clustering in certain regions as well as the absence of moderate-size compressional events along significant (>2000 km) lengths of trench (as discussed in Section 2.3). We now consider the potential role of intraplate stress variations in controlling such variability. To explore this idea, we follow Zoback et al. (2002) and adopt a value of ±5 TN m⁻¹ to represent the typical magnitude of the integrated intraplate stress variations.

Although the mechanism of brittle reverse seismicity remains opaque (as discussed in the previous section), it is instructive to consider general dynamic interactions between the net axial force and bending. To demonstrate, we consider model C (e.g. in Fig. 7c). Increasing the net force by +5 TN m⁻¹ results in an increase the neutral plane depth by about 4.5 km. This has the effect of reduces the maximum differential stress that is developed beneath the neutral plane. Furthermore, the integrated volumetric strain strain rate beneath the neutral plane will be reduced. This is because in simple-bending, the axial strain rate depends on distance from the neutral plane.

A useful way to quantify these dual effects is to consider the analastic dissipation rate, as was proposed by Lorinczi & Houseman (2009) for understanding the distribution of seismic moment. Fig. 8 shows the predicted dissipation for model C, with the same prescribed variations in net axial force. The strain rate is based on a simple bending (rate) model, being proportional to distance from the neutral plane. The results demonstrate that the dissipation beneath the neutral plane is particularly sensitive to the neutral plane depth, and hence the background stress.

In this respect we think that the typical magnitude of intraplate stress variations, producing ∼5–10 km increase in the neutral plane depth (assuming 100 Myr lithosphere), may be sufficient to explain the absence of seismicity along some margins. In this view, outer rise reverse faulting is a marginal process: it probably requires a neutral-to-moderately compressional background stress, and may be largely suppressed with only moderate changes in the background stress. In all of these situations, however, the stress state remains dominated by bending, while the background stress plays a modifying role (e.g. Fig. 7). This presents a counterpoint to the argument of Mueller et al. (1996b), that links outer rise compressional seismicity to anomalous, spatially localized subduction resistance, where axial forces dominate.

Tonga-Kermadec provides some important observations in light of the framework described above. This region stands out globally as having the most numerous and broadly distributed pattern of outer reverse faulting (see Ye et al. 2021, Fig. S5). It also stands out as having an unusually shallow neutral plane depth; Northern Tonga is particularly atypical, but across the entire system an apparent neutral plane depth of 25–30 km seems a reasonable inference (e.g. Fig. 3). These two factors support the argument above: a more compressional background stress state results in a shallower neutral plane, and much higher dissipation in the compressional part of the plate. Of course, this inference about the background stress state is consistent with the longstanding idea that the Tonga-Kermadec slab is anomalously compressional (Isacks & Molnar 1971; Gurnis et al. 2000).

8.3 Energy dissipation in bending

Conrad & Hager (1999) suggested that bending of the lithosphere during subduction may constitute a major component of dissipation in mantle convection, being around 60 per cent of a slab’s potential energy release. Several studies have presented a counter arguments, based on kinematic models incorporating more sophisticated rheological models (Buffett & Becker 2012), numerical modelling based on multilayer strength models (Capitanio et al. 2009) and convective scaling arguments (Leng & Zhong 2010).

Buffett & Becker (2012) estimated much lower bending dissipation values: ∼10 per cent of the slab’s potential energy. This conclusion was developed for a strength model with a saturation moment of  4.5 × 10¹⁷ N (estimated for old lithosphere). In this study we have argued that a strength range of 1–2 × 10¹⁷ N is compatible with data for old lithosphere (>100 Myr). As the dissipation is linearly dependent on the saturation moment (for a given subduction hinge geometry), our study suggests further reduction in terms of the importance of bending dissipation for mantle convection.

9 CONCLUSIONS

The current study has a close affinity to the work of Chapple & Forsyth (1979); they combined flexural modelling with seismological constraints in order to develop a consistent (albeit ad hoc) model for the yield strength distribution of old oceanic lithosphere. The resulting strength model remains broadly consistent with a range of additional studies and data (as shown in Fig. 1).

We survey a range of historical seismicity data, including well-constrained centroid depths of outer rise earthquakes. This data increasingly supports the view that normal faulting transitions to reverse faulting within numerous segments of the global subduction zone outer rise system, consistent with a stress state dominated by flexure [as was originally interpreted by Chapple & Forsyth 1979), but questioned by other studies (e.g. Mueller et al.1996b)].

We highlight the considerable consistency between the studies of Chapple & Forsyth (1979), Levitt & Sandwell (1995), Hunter & Watts (2016) and Garcia et al. (2019), regarding the strength of old oceanic lithosphere. All of these studies can be interpreted as favouring ‘intermediate’ strength models of the lithosphere. While there may still be a discrepancy in the strength inferred in seamount versus subduction-related flexure, the problem could be less acute than proposed by Bellas et al. (2022).

The way in which different deformation mechanisms interact to produce this ‘intermediate’ strength is less clear. In trying to constrain the relative strength distribution with depth, the neutral plane constitutes an important constraint. In contrast to Chapple & Forsyth (1979), we do not assume that the ‘apparent’ neutral plane depth (30 km for 80–120 Myr lithosphere) is equivalent to the ‘intrinsic’ neutral plane depth. That is, we allow for the fact that the apparent neutral plane depth could be biased by regions due to systematic (non-isotropic) background stress deviations.

We characterize this problem as a parameter space, broadly defined as the relative strength in the brittle (extensional) regime, the strength in the ductile region, as well as the net axial force. In doing so do we need to make certain a priori assumptions about what might constitute reasonable variability in the background stress. We consider a range of end member models, based on previous studies, each with intermediate strength characteristics, but different relative strength characteristics. The problem is inherently non-unique, but additional constraints are of some value in assessing such models.

Evaluating these strength models in light of additional constraints leads to several conclusions: (1) none of the investigated models are entirely consistent with additional constraints, if the assumed variations of net axial force (±5 TN m⁻¹) are applicable to the global subduction system; (2) in our view only one of these models can be ruled out. This end-member is characterized by high frictional strength, and very low ductile strength, as proposed by McNutt & Menard (1982); (3) in regions of ∼100 Myr lithosphere, which exhibit a 30 km neutral plane, net slab pull should be no greater than about 1–2 TN m⁻¹ and (4) models characterized by very low frictional strength (e.g. Garcia et al. 2019) are considered viable. However, they imply that regions exhibiting a ∼30 km neutral plane are under moderate net axial compression (∼−5 TN m⁻¹). Under these conditions, frictional reverse faulting is predicted beneath the neutral plane at depths of ∼30–40 km.

In Section 8, we address the global variability in outer rise seismicity (particularly the presence/absence and clustering of reverse seismicity). We again consider this variability in terms of the impact of background stress state. We show that the magnitude of the anelastic dissipation beneath the neutral plane is particularly sensitive to the neutral plane depth, and hence the background stress. We suggest that outer rise reverse faulting may be a marginal phenomenon, with moderate variations in the background stress being sufficient to inhibit brittle behaviour. In the framework proposed here, background stress plays a modifying role, with stress patterns still dominated by flexure (cf. Mueller et al. 1996b).

The study integrates numerical modelling to support the analysis and data interpretation. Some key insights from the numerical model are: (1) the subducting plate is about 85 per cent moment saturated at the first zero crossing in the outer rise; (2) in the strong part of the lithosphere, strain rates are well approximated by the simple bending model (proportional to distance from the neutral plane); (3) representative strain rates in the strong ductile part of the lithosphere are |$\mathcal {O}(10^{-15})\,\, \mathrm{s}^{-1}$| for subduction velocities of ∼10 cm yr⁻¹ and (4) the net axial force varies throughout the trench-outer rise domain, and is closely correlated with the estimated contribution of the GPE. This leads to negligible net axial force beneath the trench.

SUPPORTING INFORMATION

Figure S1. Numerical model domain. Model height is 2900 km, and width is 11 600 km (aspect ratio of 4). Coloured contours show shape of slab at 0, 5, 10 Myr. Scalar field is temperature at 5 Myr. Convergence velocity at 5 Myr is shown with arrows, and is dominated by rollback. The texture in the main panel is a line integral convolution showing the direction of flow in the mantle. Inset panels show the adaptive mesh refinement at 5 Myr.

Figure S2. Snapshot of model dynamics in quasi-steady state (5 Myr after start of simulation). Left-hand panels show subduction hinge, right-hand panels show larger part of the model domain: (a) Strain rate magnitude and (b) differential stress magnitude. Both of these scalar values are proportional to the square root of the second invariant (J2) of the deviatoric tensors; (c) Deformation map. Colours show the dominant deformation mechanism. Elastic deformation is defined as locations where the anelastic mechanisms contribute <25 per cent of the total strain rate.

Figure S3. Benchmark for isoviscous simple shear, adapted from Farrington et al. (2014). Blue line shows analytical solution, the black lines show results from ASPECT.

Figure S4. The left-hand side shows model topography (top panel) and cumulative moment (bottom panel) based on either the flexural topography or the vertical shear stress resultant. The right-hand side shows the same data at different scale. In the bottom panels the bending moment evaluated at the first zero crossing is shown with the black dashed line: this is the value that the estimates should recover. The flexural component of the model topography is shown in the black line, the best-fitting (L2) uniform-thickness elastic plate solution is shown in the blue line, as well as the pseudo-topography (the gradient of the vertical shear stress resultant) in the dashed red line. In the lower panels the horizontal dashed black line shows the bending moment calculated from the model at x₀.

Figure S5. Strain rates in the bending lithosphere. Solid lines shows the axial strain rate in the numerical model at the first zero crossing (x₀) in blue and the maximum bending moment in black [M_b(max)]. These locations are shown in Fig. 5 of the main manuscript. The strain rates have been scaled to a representative convergence rate of 10 cm yr⁻¹, using the proportionality between bending rate and convergence velocity. The stair-step nature of the strain rates is due to the post-processing, in which strain rates are averaged across elements. The red dashed lines show the strain rate assumptions we use in YSEs in this study. The black line shows the much lower (constant) strain rates used in Garcia et al. (2019). Although Hunter & Watts (2016) nominally consider bending rates in their YSE models, they implement these as a constant strain rate, with a typical magnitude of 10⁻¹⁶ s⁻¹ (as far as we could determine based on values shown in the study).

Figure S6. Comparison between the strain rate assumptions used in this study, and those of (Hunter & Watts 2016). All YSEs have the same creep parameters (Set 1, Table 1). Two different friction coefficients are shown. Both resulting strength models were considered consistent with the optimized elastic plate thickness values in Hunter & Watts (2016). The difference between the red and black lines is due to the strain rate model we adopt in this study, which is proportional to distance from the neutral plane, and has significant higher values that the constant rate (10⁻¹⁶ s⁻¹) assumed in Hunter & Watts (2016). The result of our assumptions is a stronger model.

Figure S7. Top panel shows the distribution of curvature values (shown as the reciprocal—radius of curvature, or ROC) from the profiles of Levitt & Sandwell (1995). Black line in top panel shows the evolution of moment as a function of curvature, for a strong lithosphere model. Bottom panel shows the predicted moment deficit, which is the predicted moment, minus the observed moment, using the rheological model as shown in top panel. The trend of the moment deficit should follow the dashed line, if moment deficit is a primary cause of the scatter, which is not the case.

Figure S8. Temperature model comparison for 100 Myr lithosphere. Black dashed lined shows plate model parameters from Parsons & Sclater (1977). The red line shows the model from Richards et al. (2018), which is a numerical-based model based on joint inversion for subsidence and heat flow. The blue line shows the temperature model used in Bellas et al. (2022). The cooler temperature profile is primarily due to the higher diffusivity assumed (1.0, rather than 0.8 × 10⁻⁶ m² s⁻¹). The grey regions shows 600–800 °C, where the temperature dependence of strength is largest.

Table S1. Parameters used in the numerical model.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. DS was supported by the Australian Research Council (Discovery grant DP150102887). TJC was supported in this work by the Royal Society under URF\R1\180088. TJC was also supported through COMET, the UK Natural Environment Research Council’s Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics. The authors thank the developers of open source software packages that enabled this research (in particular ASPECT, The Geodynamic World Builder, SciPy, Numpy, Matplotlib, PyVista) Bangerth et al. (2020), Hunter (2007), Fraters et al. (2019), Virtanen et al. (2020), Harris et al. (2020) and Sullivan & Kaszynski (2019). The authors thank the reviewers Fan Zhang and Shijie Zhong for their constructive reviews which have helped to improve the paper. We also thank Fern Storey of GJI for their assistance.

DATA AVAILABILITY

Input files and a description of code modifications to reproduce the numerical model can be found at https://github.com/dansand/subduction_GJI2022. Data plotted in Fig. 1 were sourced from tables in previous studies (Levitt & Sandwell 1995; Garcia et al. 2019; Craig et al. 2014a), except for data of Hunter & Watts (2016), which were redigitized by the authors. Data in Fig. 3 are from Craig et al. (2014a). Microseismicity data are from Hino et al. (2009), Gamage et al. (2009), Obana et al. (2012), (Obana et al. 2014) and Obana et al. (2018) for Honshu, and Zhu et al. (2019) and (Chen et al. 2022) for the Marianas (shallow and deeper regions, respectively).

References