An explanation for the age independence of oceanic elastic thickness estimates from flexural profiles at subduction zones, and implications for continental rheology

Most properties of oceanic lithosphere are widely observed to be dependent on the age of the plate, such as water depth, heat flow, and seismogenic thickness. However, estimates of the ‘effective elastic thickness' of oceanic lithosphere based on the deflection of the plate as it enters a subduction zone show little correlation with the age of the incoming lithosphere. This paradox requires reconciliation if we are to gain a full understanding of the structure, rheology, and behaviour of oceanic lithosphere. Here, we show that the permanent deformation of the plate due to outer-rise faulting, combined with uncertainties in the yield stress of the lithosphere, the in-plane forces transmitted through subduction zones, and the levels of noise in bathymetric and gravity data, prevents simple elastic plate modelling from accurately capturing the underlying rheological structure of the incoming plate. The age-independent estimates of effective elastic thickness obtained by purely elastic plate modelling are therefore not likely to represent the true rheology of the plate, and hence are not expected to correspond to the plate age. Similar effects may apply to estimates of elastic thickness from continental forelands, with implications for our understanding of continental rheology.


Introduction
Constraining the mechanical properties of the lithosphere is important for understanding how it supports and transmits stress, and also the controls on the location and characteristics of deformation. The influence of composition and layering on lithospheric strength in the continents remains a controversial topic (e.g. Jackson et al., 2008;Burov, 2010). However, oceanic lithosphere should represent a case in which there is relatively little variation in composition, either vertically or laterally. The rheological structure of the plate is believed to be composed of a single strong layer at the top of the plate, underlain by much weaker material beneath. Numerous studies have tried to estimate the 'effective elastic thickness' (T e ) of oceanic lithosphere (the thickness of the layer within the plate capable of supporting elastic stresses over geological timescales -hereafter referred to as the elastic thickness) from the deflection of the incoming oceanic plate into subduction zones, as observed in either bathymetric or gravitational data (e.g., McNutt, 1984;McQueen and Lambeck, 1989;Levitt and Sandwell, 1995;Bry and White, 2007; Contreras-Reyes * Corresponding author. E-mail address: craig@geologie.ens.fr (T.J. Craig). and Osses, 2010, see also Fig. 1). The loading of the incoming oceanic plate by the overriding forearc, along with the negative buoyancy of the downgoing oceanic lithosphere, causes the plate to bend down into the subduction zone. The rheological properties of the plate control the width and amplitude of the trench itself, along with the wavelength and amplitude of bathymetric features further seawards including the outer rise (see Fig. 2).
Many of the observable properties of oceanic lithosphere show a general dependence on the age of the lithosphere, and are therefore likely to be related to the thermal structure of the plate. The subsidence of the seafloor, the surface heatflow, the seismological structure of the plate, and its seismogenic thickness, all show a general correlation to the plate age (Parsons and Sclater, 1977;Wiens and Stein, 1983). The seismogenic thickness is expected to be a similar indicator of the mechanical strength of the plate to the elastic thickness (at least in terms of general trends), as both are strongly dependent on the thermal structure of the lithosphere, regardless of where within the lithosphere the plate-driving forces are supported. Given the direct correspondence seen between plate age and seismogenic thickness (Wiens and Stein, 1983;McKenzie et al., 2005;Craig et al., 2014), a similar age dependence might be expected from estimates of the elastic thickness.   1. Elastic thickness measurements from published studies using trench flexure, and thermal structure as a function of plate age. Diamonds are from McNutt (1984), squares are from McAdoo et al. (1985), triangles are from McQueen and Lambeck (1989), stars are from Levitt and Sandwell (1995), circles are from Bry and White (2007), hexagons are from Contreras-Reyes and Osses (2010), and octagons from Chang et al. (2012). Black points are those estimates derived from modelling profiles stacked along-strike for a wide region. Grey points are those derived from modelling individual bathymetry or gravity profiles. Error bounds, where estimated are shown, with small black triangles indicating that the maximum values are unconstrained. All error bars are large enough to be seen -for points with no visible error bar, no error estimate was given. Isotherms are calculated using the model of McKenzie et al. (2005) for a plate 106 km thick and a mantle potential temperature of 1315 • C. The vertical axis corresponds to depth within the plate for the thermal structure, and elastic thickness for the modelled elastic thickness points.

Previous observations of elastic thickness
A compilation of studies that have estimated the elastic thickness of oceanic lithosphere, based on the deflection of the plate at oceanic trenches observed using either bathymetric or gravitational profiles, is shown in Fig. 1 (see figure caption for the sources of these estimates). As this compilation demonstrates, individual estimates of the elastic thickness show no direct correspondence with the age of the plate where the measurements were made. For example, some estimates of the elastic thickness from trenches where the downgoing plate age is ∼150 Ma are as low as 15 km, whilst some based on trenches where the downgoing plate age is ∼25 Ma are as high as 30 km. This lack of age dependence holds true even when comparing estimates derived using the same data, and the same technique, from the same study. This paper seeks to understand why individual estimates of the elastic thickness from oceanic trenches show such little dependence on plate age, in contrast to other indicators of plate structure and rheology. Such an explanation is important for understanding the way in which the lithosphere supports applied stresses, and how these are reflected in the deformation of the plate.
Previous estimates of the elastic thickness in the oceans have relied upon modelling the deflection of the plate under an applied load -either a seamount or at an oceanic trench. However, in order to make an estimate of the mechanical properties of the plate, a rheology must be assumed. Bathymetric and gravitational profiles have been shown to be insufficient to distinguish between different rheological models for oceanic lithosphere, with different combinations of elastic, elastic-plastic and viscous responses all able to fit the deflection data equally well (Chapple and Forsyth, 1979;Forsyth, 1980). The occurrence of earthquakes associated with the bending of the plate at the outer rise (e.g. Stauder, 1968;Chapple and Forsyth, 1979;Craig et al., 2014) indicates that the process of plate flexure in this setting it not purely elastic, but involves the brittle failure of the plate resulting in permanent deformation. This brittle failure of the plate in earthquakes is commonly treated in modelling studies as the bulk plastic, rather than elastic, deformation of the plate at points exceeding a certain limiting yield stress (e.g., Chapple and Forsyth, 1979;Cattin et al., 2001). Whilst not strictly correct, as this approach distributes the deformation throughout the plastic region, rather than localising it onto discrete planes, it is thought to be an adequate approximation when considering the large-scale deformation.
Simple elastic modelling remains the most straightforward and widespread approach that is used in assessing plate rheology from deflection profiles. Whilst such models do not account for the observed permanent deformation, it is important to investigate how to interpret the results of purely elastic inversions of a plate that in reality breaks in earthquakes. In particular, the range in elastic thickness estimates, their typically low values relative to the seismogenic thickness (which approximately follows the 600 • C isotherm; Fig. 1; McKenzie et al., 2005), and the lack of a direct correlation between elastic thickness and plate age, are features that it is important to understand if we are to fully comprehend the rheology of the oceanic lithosphere.
Here, we investigate whether the simplifying assumption that the plate behaves elastically, when parts of it actually undergo permanent deformation, can be responsible for the independence of elastic thickness estimates on age. We also study the influence of variations in the in-plane force transmitted through the downgoing plate (essentially stresses within the plate that result from the application of far-field plate driving forces from outside the bending region) on the recovered elastic thickness. Here, we invert synthetic elastic-plastic deflection profiles to retrieve their apparent effective elastic thickness assuming a purely elastic rheology. The results of this comparison are used to aid the interpretation of the results of previously published elastic thickness inversions of real data.

Modelling
The mathematics of both elastic and elastic-plastic bending are well known, and described in detail elsewhere (e.g., McAdoo et al., 1978;Chapple and Forsyth, 1979). The approach taken here is to calculate synthetic profiles for the elastic-plastic case, using a range of rheological properties and applied in-plane stresses, and then to invert the resultant synthetic bathymetric profiles using the purely elastic approach of Bry and White (2007). In this study, we have only worked on synthetic bathymetric profiles. However, the effects investigated will apply to both bathymetric and gravitational data. For simplicity, a simple rheological model in which the yield stress is constant with depth is used in all elastic-plastic models shown here. These models are constructed for the region seaward of the trench, where the deflections are small compared with the horizontal extent of the bending region. The assumption of a constant yield stress simplifies the calculations, allowing the elastic solution to be used at low values of the bending moment (i.e. far from the trench) as stresses are below the yield stress throughout the plate along a vertical section. Use of a rheology with a depth-dependent yield stress would alter the actual values of the parameters discussed below, but would not change the overall trends, which are the main focus of this study.
The approach taken here follows previous studies (e.g., Chapple and Forsyth, 1979) by making three assumptions about the plate. First, in both elastic and elastic-plastic models, the plate is treated as a thin layer overlying an inviscid halfspace, with the underlying halfspace contributing nothing to the strength of the plate. Second, the coordinate system used is based on horizontal distance, rather then arc-length along the plate. This simplification only becomes problematic when the plate dip or curvature is no longer small, and so is of little relevance to the models shown here. Third, the thin-layer approximation is used. This approximation treats the plate as a plane sheet, rather than a spherical cap, and assumes that the normal stress perpendicular to the layer remains zero throughout. Whilst this approximation changes the details of the stress distribution within the plate, it has no effect on the support of the bending moment within the plastic layer.
The model used is depicted in Fig. 2. The deflection of the bending plate is defined in terms of the bending moment (M(x)) which varies along the length of a bending profile, and the inplane force (T ), which is the horizontal force applied to the end of the plate, and is constant along a profile. The in-plane force essentially arises from far-field forces applied from outside the bending region (e.g., slab pull, ridge push). The deflection (w(x)) of a plate undergoing elastic-plastic behaviour is given by the coupled second-order differential equations where M(x) is the bending moment as a function of distance from the trench, T is the in-plane stress, ρ is the density contrast between the plate and overlying water, and E and ν are Young's modulus and Poission's ratio (McAdoo et al., 1978). The rheological parameters of the model enter the formulation in defining the depths and horizontal normal stresses (σ xx (z)) at the top and base of the elastic core (z 1 and z 2 respectively). These stresses are determined by satisfying the expressions for the in-plane force (T ) and bending moment (M), integrated over the thickness of the plate such that In the case of a constant yield stress, the deflection of the plate far seaward from the trench will follow the elastic solution whilst M is small and the whole thickness of plate is able to support the applied stresses elastically, rather than undergoing plastic deformation. As such, the elastic solution, w 0 is a deflection constant used to scale the deflection of the elastic model, and z M is the plate thickness, can be used to calculate the deflection of the plate up until the point when the bathymetric gradient dips down towards the trench (McAdoo et al., 1978).
The co-ordinate system used is such that w = 0 at x = 0, with x = 0 corresponding to the first zero-deflection point seaward of the trench. Bathymetric profiles are then calculated by integrating Eqs.
(1) and (2) trenchwards from the top of the forebulge using a 4th order Runge-Kutta scheme (Press et al., 2007). Stress profiles are monitored at each spatial step to determine the correct expressions for the deflection and its derivatives based on whether plastic yielding is occurring at the top, base, or both, of the plate. Once the maximum moment point has been reached ( Fig. 2), where dM dx = 0, the plate ceases to bend and fail plastically, and the deflection is instead described by elastic unbending (McAdoo et al., 1978), given by Boundary conditions in each case are determined by requiring the elastic-plastic solution to match the elastic solution at the top of the forebulge, and by fixing the magnitude of the deflection at the point of transition from bending to unbending (the maximum moment point), to values chosen to reproduce typical trench deflections (0.5-2.5 km).
Having generated synthetic bathymetric profiles for plates undergoing elastic-plastic deformation, these profiles are then inverted using the approach of Bry and White (2007) to determine the effective elastic thickness under the assumption that oceanic lithosphere can be represented as a uniform elastic beam (e.g. Watts, 2001). In this approach, the synthetic bathymetric profile is fitted with an elastic profile calculated using where M T and V T are the bending moment and point load applied at the trench, and c and d are constants to compensate for long-wavelength regional tilting and vertical shifts in the data (Bry and White, 2007). D and α are defined previously in Eq. (6). The best-fitting elastic profile is calculated by minimising the misfit function, H , defined as where N is the total number of points along the inversion profile, h n is the 'observed' height at the nth point of the profile, w(x n ) is the height of the elastic beam at the nth point of the profile, and σ n the standard deviation of the 'observed' value (Bry and White, 2007). In the case of our synthetic profiles, the standard deviation is undefined, as we are not averaging several data profiles along strike of a trench, and so is instead replaced by an arbitrary constant. We use the same misfit function as the most recent global review of elastic thickness estimates from trench deflection profiles (Bry and White, 2007) to facilitate the comparison between our synthetic study, and its conclusions, and the trends seen when using real data. This technique has been used to analyse extensive datasets, both bathymetric and gravitational, and so allows comparison between the synthetic approach taken here and the results obtained from real data. Bry and White (2007) investigated the influence of profile length on the results of elastic inversions, determining that the results varied with the length of profile used, with shorter profiles typically resulting in lower values of elastic thickness. This likely results from shorter profiles lacking sufficient length to accurately constrain the long-wavelength tilt due to the variation of seafloor depth with plate age along a profile (c in Eq. (8); Bry and White, 2007). The synthetic elastic-plastic profiles determined here are therefore inverted using 900 km lengths, and the outcome compared only to the results from inversion of long-profile datasets of Bry and White (2007), except when investigating the influence of this effect (see Section 5).

The influence of yield stress
An example set of elastic-plastic profiles, and their corresponding best-fit elastic profiles, are shown in Fig. 3. In cases where substantial parts of the plate are yielding plastically, the plate curvature is much higher than in the purely elastic case, for a plate of the same overall thickness. This results in the effective elastic thickness that is recovered from an elastic inversion of an elastic-plastic profile being smaller than the true thickness of the initial elastic-plastic plate (e.g., Fig. 3e), in order to replicate the same near-trench curvature. Fig. 4 presents the relationship between yield stress and the recovered elastic thickness for a range of model results for plates with thicknesses of 50 and 20 km, yielding plastically in regions where a constant yield stress in exceeded. A range of yield stresses were used for each of nine different values of the in-plane force, chosen to span the probable range of this force (5 ×10 12 to −5×10 12 N m −1 ). Increases in the in-plane force, which widen the lateral zone over which plastic yielding occurs, increase the difference between elastic and elastic-plastic profiles if all other parameters are held constant. This difference results in a lower recovered elastic thickness for a given yield stress with increasing in-plane force in either extension or compression. When the yield stress is high, the entire plate behaves elastically, and the recovered elastic thickness is that of the elastic-plastic plate used.
These synthetic tests imply that the elastic thickness estimated by inverting bathymetric or gravitational profiles across the trench using a purely elastic model will underestimate the actual thickness of the plate capable of supporting significant stress on geological timescales in cases where plastic (or brittle) yielding occurs. Therefore, the variability seen in measurements of the elastic thickness, and their lack of direct correspondence to the plate age or seismogenic structure, are likely to partly be a result of the assumption that the plate is purely elastic. Additionally, the recovered elastic thickness depends on the in-plane force, which will vary depending on forces transmitted from the slab downdip of the subduction zone, forces arising from the plate seaward of the trench, and any forces transmitted across the subduction interface from the overriding plate. Lateral variations in these quantities (e.g., due to plate age and slab geometry) will therefore lead to apparent variations in the elastic thickness estimated using purely elastic models. Fig. 3f also shows that the misfit curve for elastic thickness has a much better defined, sharper minimum for synthetic models with small recovered elastic thicknesses, due to the increased curvature of the bending region resulting from increased levels of plastic yielding (due to either low yield stresses or high in-plane forces). Notably, for those values of estimated elastic thickness where the elastic thickness is small (see Figs. 1 and 3), the error bars -effectively indicating the width of the misfit minimum -are small, suggesting the estimate is well constrained, although the synthetic tests presented here indicate that it may actually be far different from the true thickness of the plate capable of supporting significant flexural stresses.
The estimates of apparent elastic thickness presented in Fig. 4 do not show error bars. However, it is worth noting that for cases where little plastic yielding has occurred, the misfit curve is typically quite flat (see Fig. 3f), indicating that the elastic thickness recovered, although the minimum is often close to the thickness of the input elastic-plastic plate, can encompass a wide range of values within only a small tolerance, or amount of noise.

The influence of deflection and profile length
The suites of deflection profiles used in calculating the estimates of elastic thickness used in Fig. 4 are all calculated for a fixed magnitude of the deflection at the maximum moment point. As Fig. 5a shows, the value chosen for this deflection significantly affects the recovered elastic thickness for a given yield stress. Greater deflections are reflected in much lower recovered elastic thicknesses for a given yield stress (due to greater plate deflection resulting in higher bending stresses and more plastic yielding). However, the position of the maximum moment point relative to the trench in nature is unknown, although it is likely in-plane tensional force of 2.5 × 10 12 N m −1 (in black). Also shown is the best-fit elastic plate model (dashed red). The numbers on each panel give the yield stress used in calculating the elastic-plastic deflection profiles, and the best fit purely elastic plate thickness. (f) shows the misfit vs. recovered elastic thicknesses for the five example profiles shown in (a)-(e), with each coloured line representing a different yield stress. Note that the misfit shown in this and subsequent figures is dimensionless (see Bry and White, 2007). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) to be between the trench and the top of the outer rise (Chapple and Forsyth, 1979). Because of this uncertainty, the relationship between the deflection at this maximum moment point and the deflection of the trench is not known. Therefore, due to uncertainties in the real world in the location of the maximum moment point in terms of its updip distance from the trench, or from the top of the forebulge (if observable), it is not possible to constrain the actual average yield stress of real plates from estimates of the elastic thickness made using either purely elastic or elastic-plastic models. However, on the basis that the trench deflection is typically 500-2500 m below the level of the undeflected plate, and that the trench is likely to be 0-50 km from the maximum moment point, then to reproduce the lowest (∼15 km) estimates of elastic thickness for lithosphere older than 125 Ma (Fig. 1), the average yield stress must be less than ∼200 MPa (for a constant yield stress model). At higher yield stresses, the recovered elastic thickness would be expected to be closer to the mechanical plate thickness, and little of the deformation would be accommodated by mechanisms other than elastic flexure, in contrast to the seis-micity seen globally in such locations (Chapple and Forsyth, 1979;Christensen and Ruff, 1988;Craig et al., 2014).
In the synthetic cases presented here, the deflection profile is not limited by the leading edge of the forearc as it is when modelling observed data. Deflection profiles here are limited by using profiles up to an arbitrary distance 'landwards' of the maximum moment point. Fig. 5b presents a set of tests to determine the influence that the choice of this distance has on the recovered elastic thickness, showing the same suite of profiles inverted using distance of 50, 30, 10, 0 and −10 km from the maximum moment point. The choice of this value can have a significant effect on the recovered elastic thickness for a given yield stress if there is significant plastic failure of the plate.
Similarly, the overall length of the profile seaward of the maximum-moment point can influence the recovered elastic thickness when modelling data (see Fig. 24 of Bry and White, 2007). Fig. 5c presents a set of tests, taking the same suite of profiles, and inverting using overall profile lengths of 300-900 km. In this noise-free case, the use of profiles that are too short prevents the accurate determination of the regional tilt (which should be zero in the cases shown), and instead part of the flexural signal is accommodated by allowing the plate to tilt, resulting in the recovered elastic thickness failing to accurately represent the thickness of the input elastic-plastic plate, and in some cases over-estimating this value.
In summary, without knowing the precise location of the maximum moment point in nature, we cannot use bathymetric profiles in conjunction with the distribution of seismicity to constrain the mechanical properties of the plate.

The influence of noise
The synthetic profiles used so far in this study are noise-free, which allows the difference between elastic and elastic-plastic profiles to be easily observed. When handling real data, this is not the case. Bathymetric 'noise' arises from a number of sources. At short lengthscales, faulting in the incoming plate, and the presence of seamounts and other volcanic features, results in localised variation. At longer wavelengths, regional tilting may result from variations in the age of the incoming plate along the length of the profile, or from dynamic topography of the oceanic plate. Whilst this tilt can be accounted for by inverting for a gradient across the profile, the use of profiles long enough to determine this tilt can result in regional dynamic topography being variable along the profile, and so add to the noise. As shown by Fig. 5c, profiles that are too short also risk modelling part of the flexural signal as re- To test the effect that real noise has on the inversion of bathymetric profiles, we add 'noise' to our synthetic profiles, determined by extracting bathymetric profiles from standard ocean floor away from any subduction influences from the SRTM30PLUS model of Becker et al. (2009) (see Fig. 6). Profiles are averaged over a width of 50 km, and then the mean of the profile is removed from the data before it is added to the synthetic elastic-plastic profiles.
Adding real bathymetry to our calculated profiles demonstrates that typical noise levels are such that elastic and elastic-plastic profiles become indistinguishable in nature. Fig. 6 shows the same set of profiles as Fig. 3, with noise added from the bathymetric tracks shown in Fig. 6a from the eastern Pacific -an area of seafloor unaffected by widespread intraplate volcanism, and with limited gravity anomalies at long wavelengths (>400 km), indicative of low dynamic topography. Each profile is then inverted for an effective elastic thickness, shown by the adjacent misfit curves. One of the critical features for distinguishing between elastic and elastic-plastic profiles is the amplitude and wavelength of the outer rise. However, the amplitude of seafloor noise is such that the addition of noise to our models makes sufficiently accurate observation of the outer rise difficult, even when taking the bathymetry from a relatively unaltered region of seafloor with little regional tilt.
The examples shown in Fig. 6 demonstrate that the addition of real noise flattens the misfit curve at elastic thicknesses greater than 20 km for all cases except those with the highest near-trench curvature. In the cases with lower curvature, purely elastic inversions result in acceptable fits to the data for a range of elastic thicknesses that encompasses both the mechanical thickness of the initial elastic-plastic plate used to generate the profile, and its best-fit elastic solution in the noiseless case.
Tests from sea floor with a greater regional tilt either into or away from trench, taken from the Indian and Southern Oceans (blue and red profiles, Fig. 7) fail to find a well-defined misfit minimum for an elastic plate until the curvature at the trench becomes large (i.e., the plate becomes weak). For plates where the near-trench deflection is not significantly greater (in the example shown, a factor of ∼4) than the magnitude of the regional tilt across the profile, the recovered elastic thickness can in fact overestimate the elastic-plastic thickness of the input plate (e.g., red lines of Fig. 7d, f, h).
Significant alteration to the bathymetry by volcanism or widespread tectonic deformation (e.g., green profiles, Fig. 7) prevents the accurate determination of any regional tilting, even with long profile lengths. The inclusion of substantial non-flexural signals in the profiles used can lead to false minima in the misfit curve, which do not represent the plate rheology, but the numerical best fit of the elastic plate model to the combination of flexural and non-flexural signals.
It is important to note that, whilst the difference between the synthetic elastic-plastic profiles and the best-fit elastic profile for the cases shown in Fig. 3 is clearly apparent, the addition of real noise to these datasets makes distinguishing between the two rheological cases extremely difficult using only real bathymetric data. However, the presence of earthquakes associated with the bending of the plate supports our choice of an elastic-plastic rheology for this work. Additionally, noise sufficiently perturbs the bathymetric profile that the misfit function for the elastic thicknesses becomes significantly flatter, and makes the upper estimates effectively unconstrained in many cases, particularly in cases where the maximum gradient and overall deflection are relatively low.

Oceanic plate rheology
The lack of correlation between effective elastic thickness estimates and plate age in the oceans is probably a result of the assumptions made in modelling the plate as a simple elastic beam, the influence of lateral variations in the in-plane force, and the effect of making assumptions about the moment distribution of the elastic plate (e.g. the 'plate break', Jackson et al., 2008). The constant-yield-stress model used in calculating elastic-plastic profiles here is only one of a number of possible formulations for the Fig. 7. Noise tests using bathymetry from the Indian and Southern Oceans. Caption as in Fig. 6. The three tracks are chosen to be indicative of regional tilting into and away from the trench (red and blue) and volcanic addition to the profile (green). (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) failure envelope, and may not correctly represent the true rheology of the lithosphere. However, the issues highlighted and the overall trends discussed here will apply to any realistic envelope for the failure stress of the plate. That the plate does indeed exceed this limit, whatever it may be, is demonstrated by the occurrence of seismicity associated with the bending of the plate at the outer rise, highlighting that the deformation of the plate is not purely elastic. Fig. 3 shows that bathymetric profiles with the most deflection give the best-constrained minima, and result in low apparent elastic thickness values. This likely explains why in Fig. 1 the lowest error bars, often based on a change in the misfit relative to the minimum, are on small estimates of the elastic thickness, and for larger values the upper limit is typically unconstrained. As demonstrated in the previous section, the consideration of real bathymetric noise makes accurate determination of the misfit minimum difficult, particularly when features unrelated to the plate bending are included in the study profile.
Similarly, for plates with relatively low curvature, the lack of constraint on the point at which the elastic plate ends and where the load is applied, often termed the 'plate break' in the elastic model, and on the magnitude of the applied load, results in the minima in the misfit curve being unconstrained at high values of the elastic thickness (see Jackson et al., 2008 for a discussion of this effect). This tradeoff between the point on the plate where the boundary condition is assumed, the magnitude of the load, and the recovered elastic thickness, is probably the explanation for the points in Fig. 1 that indicate an elastic thickness significantly greater than would be expected given the depth limit on earthquake occurrence of ∼600 • C (McKenzie et al., 2005). As Jackson et al. (2008) note, whether the elastic inversion scheme used allows the point of plate-break to be free, or fixes it a priori, can lead to significant variations in the recovered elastic thickness, and in the estimated error. Points coloured black in Fig. 1 all allow this parameter to vary during inversion. It should also be noted that, whilst most of the effects discussed here result in the underestimation of the mechanical plate thickness when modelled using an elastic rheology, the use of profiles that are too short to capture the longest wavelengths of the flexural signal (Fig. 5c) and the influence of significant non-flexural bathymetry (see tests using noise from the eastern Indian Ocean; Fig. 7), can lead to the overestimation of the mechanical plate thickness, in addition to the problems discussed by Jackson et al. (2008).
The occurrence of outer rise earthquakes tells us that the elastic limit of the plate has been exceeded, and the mechanism of these earthquakes indicates the orientation of the principle stresses at the point of failure (e.g., Chapple and Forsyth, 1979). However, seismological studies currently lack the ability to further constrain the shape of the yield stress envelope, as the magnitude of the stress drop in most earthquakes, and whether this represents the complete or partial release of stress, remain uncertain. Hence, even when combining the shape of the plate deflection with seismicity observations, limitations in these data, and remaining uncertainties in the in-plane force and the applied bending moment, make the determination of the plate rheology difficult.
Given the substantial evidence for the accommodation of plate bending by mechanisms other than elastic flexure, along with the limitations of the modelling approach, it is not surprising that estimates of an effective elastic thickness of the incoming plate derived from purely elastic modelling of trench deflection profiles do not correspond to a simple age-dependent trend. Various factors unaccounted for in the simple uniform elastic modelling approach can lead to the best-fit elastic model either over-or under-estimating the actual thickness of the layer within the incoming plate capable of supporting stresses on geological timescales.
Estimates for the elastic thickness of oceanic plates based on seamount loading show a general increase in elastic thickness with age at the time of loading, although the exact form of this increase is difficult to determine given uncertainties in the age of loading, as well as in the elastic thickness determined (Watts et al., 2006). The effect of yielding, and of any applied in-plane stresses, in such settings is uncertain, and the assessment and modelling of deformation in these settings is beyond the scope of our study.

Continental rheology
The specific case investigated here is that of oceanic lithosphere bending into subduction zones. However, similar issues may affect estimates of the elastic thickness from modelling the shape of forelands of mountain ranges in continental settings. By their very nature, the in-plane forces in such regions are large, with estimates of the buoyancy force exerted by mountain ranges on their forelands reaching to 5 × 10 12 N m −1 (e.g., Molnar and Lyon-Caen, 1988;Copley et al., 2010). In addition, whilst the large-scale rheological structure of oceanic lithosphere is generally free from significant lateral variations in compositional structure, the same is not true in continental settings.
Convergence rates in continental settings are typically far slower than rates of subduction. This probably explains the relatively low levels of seismicity associated with continental forelands compared with oceanic trenches, as the bending strain rates in foreland basins are much lower. However, in a number of places, seismicity related to the bending of the continental lithosphere has been observed, displaying a similar pattern to that seen at the outer rise of subduction zones (Chapple and Forsyth, 1979;Christensen and Ruff, 1988;Craig et al., 2014), with shallow normal-faulting earthquakes and deeper thrust-faulting earthquakes. In the foreland of northern India, normal-faulting earthquakes are observed down to ∼20 km depth, with deeper thrustfaulting earthquakes from ∼30 km to ∼50 km (Molnar and Tapponnier, 1978;Chen and Molnar, 1996;Chen and Kao, 1996;Maggi et al., 2000;Copley et al., 2011). Beneath the western Tarim Basin, shallow normal faults are observed in the basin, along with clustered deeper seismicity of unknown mechanism, which may be related to plate flexure (Xu et al., 2006;Sloan et al., 2011). Around the South Caspian Basin, microseismicity suggests the flexure of the lowlands east of the Caspian Sea under the northern Alborz mountains, with extensional seismicity down to ∼20 km, and compressional seismicity in the lower crust down to ∼40 km (Nemati et al., 2013). South of the Alborz mountains of northern Iran, a single poorly-constrained normal-faulting event near Tehran may indicate brittle failure in response to flexure under the mountain range (Jackson et al., 2002).
Outer-rise seismicity is typically dominated by normal-faulting earthquakes (Christensen and Ruff, 1988;Craig et al., 2014). This may be the result either of the in-plane force being mainly extensional, and probably derived from slab pull, or of a rheology in which the yield stress increases with increasing depth and confing pressure (resulting in a deeper neutral fibre, and the thickness of the plate breaking in plate-parallel extension being greater than that breaking in plate-parallel compression). In contrast, continental foreland seismicity in the underlying plate is mainly thrustfaulting, and probably reflects the fact that the principal in-plane force is instead compressional, dominated by buoyancy forces related to crustal thickness contrasts between the mountains and their forelands. Whilst such buoyancy forces can also be transmitted across the subduction zone megathrust, they are typically smaller, and the stress state in the downgoing lithosphere at subduction zones will be dominated by the combination of the bend-ing stresses with the in-plane forces derived from the slab (see Craig et al., 2014 for a discussion). Whilst the depth of the transition between extension and compression in most oceanic cases is deeper than, or within error of, half the seismogenic thickness of the plate (Craig et al., 2014), in continental settings, this transition is observed to be much shallower, with over half the seismogenic part of the plate in horizontal compression (e.g., northern India; Copley et al., 2011). This observation is again consistent with the large compressive force exerted between the mountain ranges and their forelands.
Our results presented here imply that estimates of effective elastic thickness in continental settings made using flexural profiles may not accurately represent the mechanical strength of the plate outside of the flexing region. The assumption that the plate is elastic, in cases where permanent deformation, as evidenced by bending-related seismicity, is occurring, may lead to underestimation of the true mechanical thickness of the incoming plate. Additionally, in cases where the mechanical thickness is large, the upper bound will be unconstrained, as the magnitude of the curvature is insufficient to constrain the width of the bending region resulting in a flat misfit function for large values of the effective elastic thickness. Similar effects also arise from the uncertainty in the location of the point at which the elastic plate ends (Jackson et al., 2008). Comparisons between estimates for the elastic thickness from profiles in flexural regions (rather than estimates made in the spectral domain) and other proxies for plate strength (e.g., the seismogenic thickness) are hence not necessarily expected, nor required, to show any form of direct correspondence. However, a striking feature of some recent elastic thickness estimates is that they appear to track the seismogenic thickness (Maggi et al., 2000), with the elastic thickness being the slightly lower of the two. This correspondence may imply that the permanent deformation of the continental lithosphere in these regions is relatively minor, and that the effects we discuss in relation to the oceanic lithosphere may similarly be relatively minor on the continents.

Conclusions
That estimates of the elastic thickness of oceanic lithosphere derived from the bending of the downgoing plate at subduction zones show no clear dependence on plate age is not surprising. As the synthetic tests shown in Figs. 3, 4, and 5 demonstrate, inverting bathymetric profiles using purely elastic models will result in estimates of the elastic thickness that are significantly different from the mechanical thickness of the input elastic-plastic plate in cases where permanent deformation is occurring. For such permanent deformation to occur, the yield stress which limits the stress that can be supported elastically must be exceeded. For this to occur, and to reproduce typical trench deflections (∼500-2500 m), then either the average yield stress must be low (<200 MPa) or the in-plane stresses must be high enough that the yield stress is still exceeded (see Fig. 4). The latter possibility would result in the stress field being dominated by far-field in-plane effects, rather than local bending stresses, and would be unlikely to result in the failure of the plate in both shallow extension and deeper compression, as is often seen at the outer-rise seismicity. The in-plane force will be variable around the world, depending on a number of factors relating to the downgoing slab, the oceanic plate seaward of the trench, and the transmission of stresses across the subduction interface, and this variability will result in a range of elastic thickness estimates recovered from oceanic lithosphere of similar ages. Due to the unknown yield stress and the influence of the in-plane force, models of the plate deflection at trenches, either elastic or elastic-plastic, are not able to independently constrain the strength of the incoming oceanic lithosphere.