The role of elasticity in slab bending

Previous studies showed that plate rheology exerts a dominant control on the shape and velocity of subducting plates. Here, we perform a systematic investigation of the role of elasticity in slab bending, using fully dynamic 2‐D models where an elastic, viscoelastic, or viscoelastoplastic plate subducts freely into a purely viscous mantle. We derive a scaling relationship between the bending radius of viscoelastic slabs and the Deborah number, De, which is the ratio of Maxwell time over deformation time. We show that De controls the ratio of elastically stored energy over viscously dissipated energy and find that at De>10−2 , substantially less energy is required to bend a viscoelastic slab to the same shape as a purely viscous slab with the same intrinsic viscosity. Elastically stored energy at higher De favors retreating modes of subduction via unbending, while trench advance only occurs for some cases with De<10−2 . We estimate the apparent Deborah numbers of natural subduction zones and find values ranging from 10−3 to > 1, where most zones have low De<10−2 , but a few young plates have De > 0.1. Slabs with De<10−2 either have very low viscosities or they may be yielding, in which case our De estimates may be underestimated by up to an order of magnitude, potentially pointing towards a significant role of elasticity in ∼60% of the subduction zones. In support of such a role of elasticity in subduction, we find that increasing De correlates with increasing proportion of larger seismic events in both instrumental and historic catalogues.


Introduction
The bending of subducting plates at the trench controls how potential energy supplied by the sinking plate is partitioned between plate deformation and motion into the mantle [Conrad and Hager, 1999;Buffett, 2006;Capitanio et al., 2007]. Depending on slab rheology and external constraints (e.g., upper plate or mantle forcing), slab bending may consume a substantial part of the available potential energy and slow down subduction. The partitioning between bending dissipation and kinetic energy determines not only the velocity, but also the style of subduction and its expressions at the surface, including retreating versus advancing trench motion, compressional versus extensional stress in the upper plate, and whether the plate sinks into or stalls above the lower mantle Goes et al., 2008;Schellart, 2009;Stegman et al., 2010;Ribe, 2010]. This leads to the question whether elasticity, which allows for nondissipative deformation, could affect the energetics of bending enough to influence subduction dynamics.
Many papers have studied how a freely subducting plate (driven solely by its own buoyancy and resisted by mantle viscosity, without interaction from an upper plate or the large-scale mantle-plate system) bends under its own weight as it sinks into the mantle. The diversity of rheologies and methods used in these studies leads to a wide range of possible slab dynamics. Numerical simulations [Ribe, 2010;Li and Ribe, 2012] and laboratory experiments [Schellart, 2008] where the lithosphere is treated as an isoviscous plate show that the contrast in viscosity between the lithosphere and mantle, together with the ratio of slab thickness over depth of the upper mantle control the mode of subduction, while the plate's negative buoyancy does not affect slab morphology and subduction mode. In contrast, in models with viscoplastic rheologies Stegman et al., 2010], plate buoyancy does affect the slab geometry. Ribe [2010] argues that reason for this is that the extent of yielding is a function of time available for bending and this time is set by the slab-density-controlled velocity of sinking. In viscoelastic models [Capitanio and Morra, 2012], there also appears to be an effect of density on subduction style, but the actual role of elasticity in slab bending and its effect on overall subduction dynamics has yet to be studied systematically. This is the aim of the present study.
Studies of lithospheric surface deformation often do treat the plates as viscoelastic, but, in overall subduction dynamics and mantle convection, elasticity is generally assumed to play a negligible role. The main argument for this is that mantle processes occur on time scales longer than the estimated Maxwell relaxation time for mantle rocks, s Maxwell 5g=E, the ratio of the viscosity of the material, g, over its Young's modulus, E. The dimensionless number that describes the relative importance of elastic and viscous deformation is called the Deborah Number, De, (first defined by Reiner [1964]) and equals the ratio of the Maxwell relaxation time and the time scale of deformation, s obs . De ( 1 implies viscous deformation dominates. For De values approaching 1 both viscous and elastic processes play a role, while for De ) 1, deformation is dominantly elastic. For mantle properties, Maxwell time equals a few hundred to a few thousand years, i.e., it is much less than mantle deformation time scales which are on the order of ten thousands to millions of years, and De ( 1. But for lithospheric viscosities in the range of 10 23 to 10 25 Pa s, the higher end of the Maxwelltime range overlaps with tectonic deformation time scales.
It is debated whether this implies that elasticity is important in subduction. Although some subduction models have included elasticity in their rheology [e.g., Shemenda, 1992;Hassani et al., 1997;Toussaint et al., 2004;Sobolev and Babeyko, 2005;Bonnardot et al., 2008], they did not analyze how elasticity affected their models' behavior. The initiation of subduction seems facilitated by the effects of elasticity because these allow for the rapid release of stored energy to contribute to localized weakening [Regenauer-Lieb et al., 2001;Regenauer-Lieb et al., 2012]. In contrast, mantle convection models that include elasticity show that it does not affect the pattern of convection or Rayleigh-Taylor instability growth for Earth's range of parameters [Moresi et al., 2002;Muhlhaus and Regenauer-Lieb, 2005;Kaus and Becker, 2007]. However, these studies do indicate that elasticity may enhance slab roll back and decrease energy dissipation in subducting slabs. Using a regional free-subduction model, Funiciello et al. [2003a] showed that elasticity is a necessary ingredient to obtain an acceptable range of solutions comparable with the natural range of trench profiles, trench retreat velocities, and dip angles. In their models, slab bending at the trench as well as the bending forced by the interaction with the 660 km discontinuity introduce local strain rates that are so high that De > 0.5 locally. It might be expected that the high contribution of elasticity indicated by these local De would affect the overall style of subduction. Yet, Capitanio et al. [2007], using a similar viscoelastic model but spanning a wider range of downgoing-plate densities, did not find a clear effect of elasticity on overall subduction and trench velocities.
Furthermore, many studies have estimated the effective plate viscosity during bending in natural subduction zones to be only about two orders of magnitude above mantle viscosity, i.e., around 10 23 Pa s [Zhong and Gurnis, 1994;Conrad and Hager, 1999;Buffett and Rowley, 2006;Wu et al., 2008;Funiciello et al., 2008;Di Giuseppe et al., 2009;Capitanio et al., 2009;Schellart, 2009;Ribe, 2010;Loiselet et al., 2010;Capitanio and Morra, 2012;Buffett and Becker, 2012;Alisic et al., 2012], implying a low effective De. The relative weakness of the slab in bending has been interpreted by some as the result of yielding that weakens the slab in the high-stress bending region [Buffett and Heuret, 2011;Alisic et al., 2012] and interrupts the elastic transfer of energy from the subducted slab to the plate at the surface, and hence makes the elastic part of the deformation unimportant for subduction dynamics [Billen and Gurnis, 2005;Arredondo and Billen, 2012].
Many of these studies used plate curvature or velocities observed in natural subduction zones to estimate slab viscosities [e.g., Conrad and Hager, 1999;Buffett and Rowley, 2006;Wu et al., 2008;Capitanio et al., 2009;Loiselet et al., 2010;Buffett and Becker, 2012]. When the bending energy required to achieve observed curvature and velocities is estimated under the assumption of a viscous plate rheology, high values are obtained ranging from 60% to over 100% of the gravitational energy [Conrad and Hager, 1999;Buffett and Heuret, 2011], for plate viscosities that are a factor of 10 to 100 times that of the upper mantle. Such dissipation estimates may be lowered if yielding or nonlinear temperature and stress-dependent rheology including pseudoplasticity are considered [Buffett and Heuret, 2011;Rose and Korenaga, 2011]. Elastic storage of energy could potentially further decrease the energy required for bending. This stored energy would be available for unbending below the trench, as well as for release in earthquakes.
In this paper, we will analyze whether elasticity plays a role in controlling overall subduction behavior in natural subduction zones. To do this, we first perform a systematic assessment of the importance of elasticity in slab bending geometry using a set of basic viscoelastic free subduction models with a wide range of rheological parameters. These analyses lead us to a geometrical scaling law and allow us to evaluate how elasticity affects the style of subduction. Next, we assess the importance of elasticity in bending energetics Geochemistry, Geophysics, Geosystems by comparing elastic, viscoelastic and viscoelasto-plastic cases. Finally, we use our geometric scaling law to estimate apparent Deborah numbers for natural subduction zones from their present-day subduction morphologies and velocities. We then compare these De estimates with the ratio of the number of large over number of small earthquakes observed in different subduction zones, the so-called b-value. B-values have been proposed to decrease with stress [e.g., Schorlemmer et al., 2005;Narteau et al., 2009], and hence may be expected to decrease with the ability of the slab to store elastic energy. We find that there is a correlation in that lower b-values correlate with higher apparent De, showing that elasticity may indeed be important in the subduction process.

Model Description
We performed a large set of subduction models to obtain our geometrical scaling analysis and evaluate bending energetics. The models are basic in that no upper plate is included, and plate rheology is either uniformly elastic, viscoelastic or viscoelasto-plastic. However, they are advanced in that they do include a fully dynamic slab and trench, a free top surface, and assume slabs are of finite width with viscous mantle resistance including the effect of flow around the slab's sides. Below we describe the general model setup (subsection 2.1), the implementation of plate rheology (subsection 2.2), other model parameters (subsection 2.3) and how we measure the Deborah number in our simulations (subsection 2.4).

Model Setup
We carried out numerical simulations of a plate freely subducting into an infinite-depth isoviscous mantle, as illustrated in Figure 1. This type of model was originally set up by Funiciello et al. [2003a], and updated by Capitanio et al. [2007], and has been further validated and applied in a wide range of studies Goes et al., 2008;Capitanio et al., 2009Capitanio et al., , 2010Goes et al., 2011;Capitanio and Morra, 2012]. As in those models, the plate is driven by its own negative buoyancy, has a free top surface and is subject to viscous mantle drag and isostatic support. Its trailing edge is kept free in all models studied here. Subduction is initiated via a downward force at the slab tip until it reaches a depth of 100 km. Subsequent subduction is self-sustained. Different than in the models from Funiciello et al. [2003a] and Capitanio et al. [2007], we do not include an upper-lower mantle interface, allowing us to isolate the rheological effects on bending behavior at the trench.
The plate is modeled in 2-D (i.e., neglecting along-strike variations in strain), while the mantle drag response is obtained using the analytical solution of a rigid 3-D ellipsoid sinking in a viscous medium [Capitanio et al., 2007]. This formulation takes into account that mantle drag in 3-D is affected by mantle flow around the plate's sides, where the width of our plate is set to 1000 km. This is a crucial aspect of 3-D dimensional subduction that facilitates trench motion, which is intimately related to bending energetics [e.g., Funiciello et al., 2003a;Capitanio et al., 2007]. The mantle drag is applied to the plate via a set of dashpots, and the mechanical energy conservation equations are solved inside the oceanic plate using a Lagrangian-Eulerian finite element software, ABAQUS [Hibbitt and Sorensen, 2000]. The oceanic plate of length L and thickness h is modelled by a Lagrangian grid of 250 or 500 (for longer plates) by 20 nodes, with constant node spacing along the plate and a refinement of vertical spacing towards the plates top. The four-node elements are bilinear in velocity and constant in pressure. The subduction behavior obtained with these models is

Winkler foundation
Mantle drag s=0 x y Figure 1. Schematic of the model setup and definition of parameters. The slab of thickness h, density q L , Young's modulus E, viscosity g L , and in some cases yield stress r Y is freely sinking into an infinite mantle of density q M and viscosity g M , achieving a subduction velocity U p , deflection w and slab dip h. Slab deformation is modeled in 2-D, but the analytically calculated mantle drag that is applied along the slab's sides and tip is for 3-D flow around a finite-width plate. The plate is supported in the top few kilometers by a Winkler foundation of thickness h w . The top surface and the trailing edge are free.
In all our models, after a transient phase due to the application of the initiation force, the deflection evolves linearly with the buoyancy force meaning that a quasi steady-state is reached and slab geometry remains constant in time until the entire plate has subducted. We measured geometric parameters, plate velocities, and energy as a function of time during this quasi steady state.

Constitutive Relations
We tested different rheologies for the subducting plate ranging from purely elastic to viscoelasto-plastic cases. To keep the models simple, we take all rheological properties as constant across the entire plate thickness. In the purely elastic cases, the behavior is determined by the Young's modulus E and the Poisson ratio m, according to Hooke's law. In the models that also include viscous behavior, Maxwell linear viscoelasticity is used, such that the total strain rate, _ e, is the sum of the creep strain rate, _ e V , and the elastic strain rate, _ e E : where g L is the lithospheric viscosity, r is the deviatoric stress tensor, and _ r is the stress rate tensor. After characterizing the effect of elasticity by comparing elastic and viscoelastic cases, we will analyze a set of viscoelastoplastic cases. In these models, we use perfect plasticity based on a Von-Mises yield stress criterion, i.e., a yield stress, r Y, that is constant and independent of strain rate.

Model Parameters
The variables used in this paper and their values are all defined in Table 1. The plate density (q L ) is greater than the mantle density (q M ) and Dq5q M 2q L is the driving force in the system. We systematically varied the plate's properties: density contrast from that appropriate for very young and very old oceanic plates (10-120 kg/m 3 ), viscosity between 10 and 10 5 times mantle viscosity, thickness between 40 and 120 km, and Young's modulus from a mantle/lithosphere value of around 200 GPa to values two orders of magnitude lower and a factor of 10 higher. The Poisson ratio m is set to a constant value of 0.3 in all models. In the viscoelasto-plastic simulations, yield stress was varied from very high values (1 GPa) to values that would lead to slab break-off (50-75 MPa depending on the other rheological parameters of the plate) [Evans and Goetze, 1979;Billen and Hirth, 2007].
In all models, an isostatic restoring force, R w , is implemented via a nondissipative Winkler foundation. It acts to support the plate where the surface is not submerged into the mantle deeper than a depth h w . R w depends on density of the infill material at the trench (usually water), q i , and depth of the top of the plate, z, as: R w 5ðq M 2q i Þgz. Near the trench, the switch from support to no support once the slab dives below h w leads to an additional bending moment. Although this force is implemented to account for the isostatic restoring force associated with the trench depth and infill, we changed the strength of the Winkler foundation, by varying h w , over a large range, from 15 to 70% of the slab pull. This gives us a simple way of investigating how forces applied at the trench, which in nature would include forces imposed by the upper plate, may affect bending behavior.

Deborah Number Definition
The Deborah number is the dimensionless number that indicates whether a viscoelastic system behaves as an elastic solid or as a viscous fluid. It is the ratio of the Maxwell time (g L =E) and an observation time characteristic of the deformation process of interest. Funiciello et al. [2003a] used a local strain rate to show that elasticity could play a role in bending. Here, we are interested in the overall bending process, and choose the time that a particle of the slab spends in the bending region as our observation time, s obs . This time is equal to the length of the bending portion of the plate, l b , divided by the subduction velocity U p . Our Deborah number is thus expressed as: We follow Ribe [2010] to determine l b . Let's define K as curvature, i.e., the change in slab dip in s-direction, where s denotes the curvilinear coordinate along the slab measured along the plate's mid surface. _ K ðsÞ is the time derivative of K, i.e., the rate of change of the plate's curvature, also called the curling rate. The bending length l b is defined as the length of the portion of the plate's midsurface where _ K differs from zero (typically a few percents) [Ribe, 2010].
The curling rate can be derived from the model velocities and slab shape using equations (14) and (15) in Ribe [2010]. Figure 2a gives a more visual definition of l b , and shows how _ K goes through a negative peak during bending at the trench and a positive peak during unbending at larger depth before returning to zero. The thus defined length scale l b is fully dynamic. Ribe [2010] used l b to describe the bending response of a purely viscous sheet. In viscous plates, the bending length keeps evolving with time. A difference in our viscoelastic cases is, as shown in Figure 2b, that the bending length reaches a constant quasi steady state value. This is one visible effect of elastic unbending.
Since l b is a difficult parameter to measure in natural subduction zones, as one would need the dynamic evolution of the entire curved profile, we correlated it with the usual geometrical parameters monitored in subduction zones. Previous studies have used the approximation l b 52 R min [Buffett, 2006;Buffett and Heuret, 2011], where R min is the minimum radius of curvature. However, Ribe [2010] notes this is an incomplete description. When we measure the bending length during the quasi steady state phase of our 54 viscoelastic cases, we find it is better described by the following relationship ( Figure 2c): where h is the maximum slab dip expressed in radians. This expression, essentially the geometrical relation for arc length, captures that bending length increases with the radius of curvature, but also depends how far the plate bends along that curvature, i.e., what dip is achieved.

Bending Geometry
The first step in our analyses is to determine how bending geometry varies as a function of rheology. We characterize bending in terms of the minimum radius of curvature, following a range of previous studies [e.g., Wu et al., 2008;Buffett and Heuret, 2011;Rose and Korenaga, 2011]. We derive a scaling relation between the minimum radius of curvature for viscoelastic cases, R VE , relative to the minimum radius for equivalent purely elastic cases, R E , (i.e., cases with no viscosity, but with the same plate thickness, elastic parameters, and density contrast) and the Deborah number (subsection 3.3). This relation will allow us to estimate apparent Deborah numbers from observed minimum radii of curvature (subsection 5.1), as the equivalent elastic R E can be calculated from plate thickness and elastic parameters (subsection 3.2), based on well established bending theory (subsection 3.1). Furthermore, bending at the trench sets the stage for the subsequent evolution of slab morphology during its transition through the upper mantle, and we will assess whether elasticity affects the overall modes of subduction (subsection 3.4).

Bending Theory
For the analysis of the purely elastic cases, we can refer to existing analytical solutions that describe the bending behavior of laterally infinite elastic beams or beams clamped at their trailing edge [Caldwell et al., 1976;Turcotte and Schubert, 2002]. Caldwell et al. [1976] showed that the deflection, w E , of a semi-infinite elastic beam subject to a vertical end force and a given bending moment is: where x is the horizontal coordinate, A is a constant depending on the boundary and initial conditions and a E is the flexural length defined as: where k is the hydrostatic restoring force per unit deflection, k5 q M 2q i ð Þ g and D is the flexural rigidity, which is determined by the elastic parameters E and m, and slab thickness h. The beam's curvature equals the second derivative of the deflection, which can easily be calculated and minimized in order to find the minimum radius of curvature along the profile: where C 0 is a constant with dimensions m 21 depending on boundary and initial conditions.
It is worth noting that in the case of a purely viscous plate of uniform Newtonian viscosity [de Bremaecker, 1977], the expressions are similar, but the instantaneous flexural length is now expressed: In a similar manner to the elastic case, we can show that the minimum radius of curvature in that case is proportional to a 2 V . Thus, in both elastic and viscous cases, bending is fully described by the flexural length and the boundary conditions, and hence this parameter is also expected to play a significant role in Maxwell viscoelastic cases, where bending occurs through a superposition of elastic and viscous deformation.
Our models are more realistic than the analytical problem in a few respects: (1) the mantle traction forces are taken into account, and (2) the isostatic restoring force contributes only until the depth where the slab is fully submerged in the mantle. Since we varied the isostatic restoring force by changing h w rather than changing ðq M 2q i Þ, we recast the flexural parameter as follows: where h ref w is a reference thickness for the Winkler foundation, which we take equal to 6 km. Increasing h w leads to a higher restoring force and a lower value of a m . In the following, we will always use a m as the flexural length, and we will show that bending in purely elastic models is only a function of this flexural length.

Elastic Bending
A set of 25 purely elastic simulations was carried out where we systematically varied the main parameters of our system (h, Dq, E, h w ). Figure 3a shows the resulting shape of the slab for different values of a m , when the same buoyancy force, Dqghl, is applied. We determined the minimum radius of curvature, R E , from the profile of curvature for each of our purely elastic cases, and derive a scaling law for how R E increases with a m (Figure 3b): where C 1 51:1744310 26 m 21:35 and n 5 2.35 as determined by fitting. This scaling law is in good agreement with the analytical solution from Caldwell et al. [1976], confirming that bending geometry is mainly controlled by plate rigidity, also in our more dynamic system. Note that Capitanio and Morra [2012] also found that the flexural profile of models like ours can, at shallow depths, be fit with the analytical solution to the flexural equations.
Interestingly, the bending radius has little sensitivity to the slab-mantle density difference. This is a consequence of the dynamic system, where the trench can migrate as slab pull evolves, allowing the slab to maintain a bending geometry controlled by its elastic properties. However, as predicted by the analytical solutions, the isostatic restoring force is important as it contributes to the bending moment: increasing this force increases bending, i.e., decreases bending radius.

Viscoelastic Bending
A second series of 54 simulations was performed using a uniform viscoelastic rheology, where we systematically vary plate viscosity g L , in addition to h, Dq, E, and h w . Figure 4 shows how the stress state and plate geometry differ for otherwise equivalent elastic, viscoelastic (De $ 10 22 ) and essentially viscous (De $ 10 24 ) plates. Here the Maxwell time for the viscoelastic case is 0.53 Myr, while for the viscous end member it is only 5.2 kyr. In all three cases, the same buoyancy force is applied meaning that the same volume of slab (i.e., length measured from the trench to the slab's tip 3 thickness, which evolves as well) has subducted (although in the elastic case only at a very small angle and to a shallow depth). In the elastic plate, bending introduces symmetric patterns of high stress at the top and bottom of the plate, while in the viscous case, due to relaxation effects, there is an offset between the high stress region at the plate's bottom downdip and the high stress region at the top of the plate. The viscoelastic case displays a mixture of the two end members' stress distributions. Viscous relaxation also increases deflection compared to the elastic case for the same slab pull and isostatic restoring force.
We calculate the Deborah number from the bending length and plate velocity measured in every simulation. Then, we compared the radius of curvature obtained in viscoelastic simulations, R VE , with the one measured in an equivalent purely elastic case, R E , and plotted this ratio as a function of the Deborah number (Figure 5). At high Deborah numbers (De > 10 21 ), the radius of curvature, R VE , is identical to the one observed in purely elastic simulations meaning that the deformation is essentially elastic. At low Deborah numbers (De < 10 23 ), R VE is much less than R E and deformation is Figure 5. The relative importance of elasticity in slab bending geometry, illustrated by the variation of the ratio of viscoelastic minimum radius of curvature, R VE , over the minimum radius of curvature of the equivalent elastic case (same slab properties), R E , as a function of Deborah number. Bending is essentially viscous for De 10 23 , where R VE ( R E , and largely elastic for De > 10 21 , where R VE % R E , while both deformation mechanisms play a role in between these regimes (domains delimited by arrows at the top). The dashed line corresponds to the best fit scaling law (equation (15)). largely achieved viscously. In between these two regimes, a transitional regime takes place where both mechanisms are important.

Geochemistry, Geophysics, Geosystems
To help understand the shape of the distribution in Figure 5, we can go back to the linear viscoelastic constitutive law, and cast it in terms of the De parameters: where s max is the Maxwell time, and we allow the observation time s obs to be scaled by a constant C 2 . According to Buffett [2006] and under the assumption of pure bending, _ e ss 52zU p @K @s : where z is the vertical coordinate, and @K @s is again the variation of curvature along the slab coordinate s. Following Ribe [2010], the variation in curvature is inversely proportional to bending length times radius of curvature: where C 3 is a constant. Using the expression for observation time and equations (10), (11), and (12), we obtain: Using e E ss 5r=E for a purely elastic rheology and the fact that, under the assumption of pure bending, e E ss 52z=R E (R E could be seen as an equivalent elastic radius of curvature) we obtain: We searched for the best fit of that expression to our simulations. We allowed the exponent to vary as well, to be able to fit the lower Deborah number cases where viscous stretching might occur and hence the assumption of pure bending may not fully hold. We obtained the following scaling law, illustrated by the dashed line in Figure 5: This simple derivation gives us a good explanation of the behavior observed in our models. Geometrical law 15 gives us the possibility to evaluate an apparent Deborah number for natural subduction zones (see section 5) since R VE is measured from Benioff zone profiles, and R E can be calculated using equation (9), by assuming the plate's thickness and elastic properties.

Subduction Mode
Next, we illustrate that the effect of elasticity on bending of viscoelastic slabs indeed affects large scale dynamics. We do this by considering the different modes of subduction that are produced for different Deborah numbers and the extent of forcing at the trench (in our models represented by the large range of  R w and thereby a m ). Several studies [Funiciello et al., 2003b;Bellahsen et al., 2005;Schellart, 2009;Ribe, 2010;Stegman et al., 2010;Capitanio and Morra, 2012;Garel et al., 2014] characterized the mode of subduction selected by a slab as a function of plate rheology and the physical parameters in the system. They found a range of subduction modes (illustrated in the inset of Figure 6): (1) a viscous-beam type retreating mode, termed ''strong retreat'' (SR); (2) a ''weak retreat'' (WR) mode where the slab sinks at a large dip (but < 90 ) and flattens when resistance is encountered at the base of the upper mantle; (3) a fold-andretreat (FR) mode where the slab sinks almost vertically and piles up on top of a higher-viscosity lower mantle while the trench first advances and then retreats and (4) an ''advance'' (A) mode where the slab is overturned (angle > 90 ) when it reaches the base of the upper mantle and the trench subsequently advances. These studies all considered purely viscous or viscoplastic slabs and found that high viscosity slabs either strongly retreat or advance depending on forcing.
As we did not implement any change in the mantle at the 660 km depth, we identified the mode of subduction by the angle of the tip of the slab when it reaches that particular depth, as was done by Ribe [2010], who discusses that this condition is what controls the subsequent evolution of the subduction mode and accompanying trench motion. Figure 6 displays the thus predicted mode of subduction, as a function of the Deborah number and the dimensionless elastic flexural length, a m =h, where the latter is a proxy for the external forcing at the trench, with lower values corresponding to stronger forcing.
We find that at higher De numbers, the stored elastic energy facilitates unbending once the slab is beyond the trench. This means that retreating modes become the favored way of subduction. Folding and advancing modes only occur at low De < 3 Á 10 23 , where the slab behaves essentially viscously. It is worth noting that the advancing mode is only triggered in our simulations with a very high Winkler foundation force, unrealistic when scaled back to Earth. Elasticity requires a very strong additional bending moment in order to achieve a slab roll over. Using the same kind of models, Capitanio et al. [2010] obtained such a mode by adding a strong ridge push force at the free end of the down-going plate. This further indicates that this mode, which in seismic tomography is only seen below India [ Van der Voo et al., 1999], is not easily realized. Weakly retreating slabs are common across all De with some amount of external forcing. Do note that throughout the Cenozoic, most trench segments do indeed retreat [Goes et al., 2011;Sdrolias and M€ uller, 2006].

Bending Energetics
In the previous section, we have shown that elasticity can significantly affect bending geometry as well as the large-scale style of subduction if De exceeds 10 22 . Another way to understand the role of elasticity in bending is by considering how it affects the energetics.

Elastically Stored Energy
An important implication of elasticity resides in the amount of energy that is stored rather than dissipated during bending. In Figure 7, we display measurements of the maximum elastic energy stored in our numerical slabs during subduction, EE max , compared to the total internal energy change in the slab (i.e., elastically Figure 7. The relative importance of elasticity in deformational energy for viscoelastic slabs, illustrated by the variation of the maximum elastic strain energy, EE max , scaled to the total strain energy, E total , in the system, as a function of the Deborah number. For De within the range where bending geometry appears controlled by both viscous and elastic deformation (Figure 5), the energy stored elastically ranges from 0.1 up to 10% of the total. The dashed line is the best fit scaling law from equation (17). stored 1 viscously dissipated), E tot , as a function of the Deborah number. The proportion of strain energy stored elastically increases with the Deborah number, reaching 10% for De > 10 21 and 100% for De > 1, when the slab behaves as a purely elastic body.

Geochemistry, Geophysics, Geosystems
We can derive a scaling law from the definition of elastic and total strain energy changes: Using Maxwell's linear viscoelasticity, and Hooke's law of elasticity, and again allowing for constant scaling, by C 4 , of s obs , we find: A fit of this relation to the simulations yields C 4 50:7, and describes the results very well (dashed line in Figure 7). Buffett and Heuret [2011] found that if bending only occurs via viscous deformation, the resulting dissipation estimated from observed radii of curvature can exceed the gravitational energy available from subduction. They invoked a plastic rheology to resolve this paradox. However, as discussed above, at higher De, elastically stored energy can also significantly affect the overall energetics of the subduction.

Viscous Energy Dissipation in Bending
In Figure 8, we compared the energy required for bending one of our viscoelastic plates, denoted as U B , with that required to bend an equivalent (i.e. same viscosity, velocity, thickness, and radius of curvature) viscous plate. To estimate the latter, we use the analytical expression for the energy dissipated viscously in bending based on geometry and under the assumption of pure bending derived by Buffett [2006] and Buffett and : The result is striking. For cases with low Deborah number (De < 10 22 ), the viscous dissipation given by equation (18) is indeed a good estimate of the dissipated energy in bending as measured in our models. However, at higher De, this relation strongly overestimates the energy dissipated in the bending. For De around 0.1, the viscoelastic slab only requires 10% of the energy required to bend an equivalent purely viscous slab to the same radius of curvature. This implies that studies that estimate effective viscosity from bending radius may be underestimating slab viscosity while overestimating how much energy is required to subduct stronger slabs. The additional role of plasticity will be assessed below.

Energetics With Plasticity
We ran 27 simulations with the additional effect of plastic yielding. Do note that yielding only becomes important at higher Deborah numbers where slabs can sustain sufficiently high stress. For example, for a yield stress of 600 GPa and a Young's modulus of 200 GPa, plastic yielding only affects slabs with De > 5 Á 10 22 (supporting information Figure S1). Plastic yielding leads to local weakening of the slab at its top and base where bending stresses are highest (supporting information Figure S2). This results in higher bending than if no yielding is possible, i.e. smaller l b than for slabs with otherwise the same properties but no plasticity. However, yielding does not strongly impede the propagation of stress, as the core of the plate remains unaffected, thus leading to less reduction of U p than of bending length. The net result of this is a smaller observation time, and hence higher De.
To test how the Deborah numbers relate to elastic energy storage, we measured the maximum elastic energy stored scaled by internal energy changes for all viscoelasto-plastic cases, and plot EE max =E tot as a function of De (Figure 9). Although it does not fully describe the points, the relation that fit the viscoelastic cases (equation (17)) gives a pretty good prediction for the viscoelasto-plastic cases. This illustrates that even with yielding, the Deborah number as estimated according to equation (2) provides a measure of the proportion of stored over dissipated energy in bending.
However, when we infer the Deborah number from plate geometry using scaling law 15, this provides an estimate, which we called the apparent Deborah number, De app , that is systematically lower, by as much as an order of magnitude, than the De number estimated from deformation time scale. This implies that based on the analysis of slab geometry, we would underestimate the effect of elasticity if applied to a slab that undergoes yielding. Other nonlinear, temperature or stress-dependent, rheologies expected to govern lithospheric strength would also lead to local weakening [e.g., Billen and Hirth, 2007;Rose and Korenaga, 2011]. Hence, we would expect, similarly to the effects of plasticity documented here, that such rheologies give rise to lower dissipation relative to a uniform viscosity slab of the same shape [Rose and Korenaga, 2011], and hence De app probably gives a lower estimate of actual De. Keeping these effects of plasticity in mind, we will now apply our scaling relations to natural subduction zones to evaluate if elasticity may be important there.

Elasticity in Natural Subduction Zones?
Having established what the possible consequences of elasticity on slab dynamics could be, we next estimate apparent Deborah number for natural subduction zones (subsection 5.1). Although real subduction zones are expected to be substantially more complex than our models, including additional forcing from the upper plate and mantle and lateral variability in plate buoyancy and strength, our analysis should be able to provide a first-order assessment of whether the effect of elasticity could be significant in subduction zone bending on Earth. We further test the De app estimates by comparing them with a measure of earthquake activity, as seismicity is an expected expression of elastic energy storage and release (subsection 5.2).

Apparent Deborah Numbers From Subduction Parameters
Using the scaling laws derived from our viscoelastic simulations (equations (9) and (15)), we can now estimate the apparent Deborah number of natural subduction zones from the observed range of radii of curvature, which we equate to R VE . We use Rose and Korenaga [2011]'s compilation of R min which is a slightly extended version of the compilation by Wu et al. [2008]. The maximum De app is capped at 1, as the scaling law breaks down for R VE =R E > 1:047. From equations (2) and (3), we can also estimate the apparent Figure 9. Maximum elastic energy, EE max , scaled to the total strain energy, E total , in the system as a function of the Deborah number, De, for all viscoelastoplastic simulations. Filled circles are measured De, while open circles display the apparent Deborah numbers estimated from the geometric scaling law that we derived for our viscoelastic models (equation (15)). The dashed black line is the energy scaling law previously derived for viscoelastic simulations (see Figure 7). Both De estimates converge for low De where the slabs are too weak to yield, but the apparent Deborah number underestimates the actual De and hence the importance of elasticity if plastic yielding occurs in the slab.
Geochemistry, Geophysics, Geosystems The equivalent elastic radii of curvature R E are calculated using a constant Young's modulus of 70 GPa [Afonso et al., 2005], and a constant restoring force assuming a trench of 6 km depth filled with water. If we underestimate E, or the restoring force, by a factor of 2, we would overestimate, respectively, underestimate, De app by a factor of $2 (equations (8), (9), and (15)). As plate thickness, we use thermal thickness, estimated from age at the trench, assuming conductive thickening up to a maximum thickness of h 5 100 km at t 5 80 Ma, i:e:; hðtÞ5minð11:2 ffiffi t p ; 100Þ where t is in Ma, following Ribe [2010]. There are two common ways of defining elastic plate thickness. (1) Elastic seismic waves image lithospheric thicknesses. Below the oceans, seismic thickness closely follows thermal thickness [Nishimura and Forsyth, 1989;Maggi et al., 2006].
(2) Elastic plate flexural models that are fit to observed bathymetric and gravity profiles yield effective elastic thickness, T e . At trenches, T e tends to be only about half as large as plate thermal thickness [e.g., Levitt and Sandwell, 1995]. In principle, and as assumed in our models, the whole thickness h of a viscoelastic plate has the potential to behave elastically. Depending on the forcing the plate is subjected to, it may not, and T e less than h would be an expression of partially nonelastic plate behavior on tectonic time scales. Our De app scaling laws use h, and provide an alternative way of estimating the contribution of elasticity, one that can be linked to bending energetics. Note that overestimating h by a factor of 2 would lead to an overestimate of R E and hence an underestimate of De app , by a factor $8. Thus, the most significant uncertainties in R E (in h and h w ) would both lead to an underestimate of De app .
About 85% of the transects yield De app below 10 22 (Figure10a). At face value, this would imply that elasticity is of minor importance in most subduction zones (Figures 7 and 8). The majority of our subduction transects yield apparent plate viscosities between one and three orders of magnitude higher than mantle viscosity (Figure 10b). This range is similar to that inferred for subducting plate viscosity during bending from a large number of studies using a wide range of methods, all assuming plates are purely viscous or viscoplastic [Zhong and Gurnis, 1994;Conrad and Hager, 1999;Buffett and Rowley, 2006;Billen and Hirth, 2007;Wu et al., 2008;Di Giuseppe et al., 2009;Ribe, 2010;Loiselet et al., 2010;Buffett and Heuret, 2011;Alisic et al., 2012], again indicating that subduction bending may be quite well described as a viscous process.
However, plastic yielding has been proposed to occur in a number of or even most zones [e.g., Levitt and Sandwell, 1995;Billen and Gurnis, 2005;Buffett and Heuret, 2011;Arredondo and Billen, 2012;Alisic et al., 2012]. If plasticity indeed plays a role, then any De app > 10 23 may underestimate actual De by up to an order of magnitude, and elasticity would play a role in slab bending at about 60% of the subduction transects (Figure 10a).
Our De app estimates span a large range from 10 24 to elastic values of 1, and c app ranging from 1 to 10 5 (Figure 10). In a few zones, Nankai, Cascadia, the southern tip of Chile, central America and eastern Alaska, high to very high values indicate significant elastic behavior. Apart from Nankai, which is near a triple junction, Figure 10. (a) Apparent Deborah numbers for natural subduction zones, estimated from equation (15) that relates bending radius and De for viscoelastic slabs, for each zone as well as a histogram. We obtain a wide range of De estimates ranging from purely viscous to purely elastic behavior. If yielding occurs, these De could be underestimated by up to an order of magnitude pointing toward a higher importance of elasticity. (b) Effective viscous contrast for natural subduction zones estimated from equations (2) and (3). The histogram on the righthand side displays a mean value around 100, consistent with other studies. Subduction zone labels are taken from Lallemand et al. [2005].
Geochemistry, Geophysics, Geosystems  (Figure 11b), a clear decreasing trend emerges. This, surprisingly, indicates that younger plates may behave more elastically than older plates. This trend is largely the result of the age-dependent elastic thickness, which we based on the observation of an age trend in both seismic and T e estimates of elastic thickness [e.g., Nishimura and Forsyth, 1989;Levitt and Sandwell, 1995;Maggi et al., 2006]. Age-dependent h leads to R E values that increase with age. By contrast, R VE does not have any clear trend with age (Figure 11a), although the highest values do occur for some of the youngest plates. Buffett and Heuret [2011] interpreted the lack of an age trend in R VE as evidence of yielding, which they required at all of their old subducting plates and some intermediate aged ones. Do note that there is a lot of scatter in R VE and De app , and certainly not all young plates are characterized by high De. However, some young and intermediate aged plates have bending radii that approach elastic R E , while all the oldest plates bend much more strongly than expected from R E .

Deborah Numbers and Seismicity
At higher Deborah numbers, the plates store more elastic energy and hence should have higher seismic potential. Furthermore, higher De could increase plate coupling by increasing normal stress across the contact between the lower plate and upper plate. Indeed stress transfer between the downgoing plate and the plate interface has been previously invoked to explain observed temporal correlations between subduction interface earthquakes and seismicity in the outer rise and in the down-dip Benioff zone [Astiz et al., 1988;Lay et al., 1989]. So we might expect that higher De is associated with higher seismicity rates and/or larger earthquake magnitudes.
To test this, we compare our De app estimates with the seismic b-value, which is the slope of the earthquake magnitude-frequency distribution. Low b-values correspond to distributions with high numbers of large events compared to the number of small earthquakes. We might expect that higher elastic energy storage provides a more favorable condition for larger events. Furthermore, some previous studies have related the b-value to the level of stresses in the lithosphere [Schorlemmer et al., 2005;Narteau et al., 2009]: the higher the stresses, the lower the b-values. Our models showed that higher De slabs reach higher maximum stress levels (supporting information Figure S1). Thus, from both arguments, we might expect higher De to correlate with lower b-value.
For this analysis, we defined 24 boxes around the Pacific, Sandwich, and Antilles subduction zones (see Figure 12 and supporting information Table S1). To enclose a reasonable number of earthquakes, these boxes each comprise several of the profiles from Lallemand et al. [2005] and extend 1 seaward and 4 trenchward. Their along-strike lengths are defined in such a way that the average De app of a particular box remains of the same magnitude as the De app of each of the individual profiles contained in that box. Measuring b-values carries substantial uncertainties, as seismicity rates tend to vary significantly during the seismic cycle between large events and seismicity catalogues are on the short side relative to the repeat time of large earthquakes. We therefore use two different measures of the relative rate of small and large events, one based on the recent instrumental Global CMT Catalogue [Dziewonski et al., 1981;Ekstr€ om et al., 2012], and one based on the historical catalog from the past century [Engdahl and Villaseor, 2002;Heuret et al., 2011].
The Global CMT catalogue spans about 30 years and appears to provide a relatively stable rate of events with magnitudes exceeding 5.5 since 1982 [e.g., Main et al., 2008;Kagan, 2010]. In this catalog from 1982 to 2010, we counted the number of earthquakes located inside each box and we fitted a power law trend to the cumulative number of events as a function of moment magnitude between a range of 5.5 and 7. For larger magnitudes, the length of this catalogue is too short to estimate recurrence times and hence earthquake rates for individual subduction zones. Only boxes with at least 30 events are considered. This fit gives us a b-value that we plot as a function of the average De app of the box. The b-values clearly decrease with Deborah number (Figure 13a).
To confirm the trend, we also calculated b-values using all the earthquakes contained in the boxes with a Deborah number of the same magnitude. This averages out the effect of local tectonics thus highlighting the influence of rheology on the b-value. The b-values thus determined are better constrained, since the number of earthquakes per Deborah number is larger. The resulting points (red squares in Figure 13a) show an even stronger correlation (correlation coefficient R50.95) and hence confirm the trend from the b-values in the individual boxes.
For the historic catalogue, we looked at the average magnitude per event, M MME , as a proxy for b-value, as the historic catalogue suffers from strong variations in catalogue completeness with magnitude. The long catalog does give a reasonably good estimate for the total long-term moment rate, which is dominated by the largest events. We use the long-term moments from the 1900-2007 catalog compilation of Heuret et al. [2011] (the moments corresponding to their M MRR are summed to obtain moment for each of our zones). We tested two estimates of the number of events per zone: (1) a number of events with magnitude greater than 5.5, N long !5:5 , which we scaled up from the CMT catalogue, as this catalog probably provides the most reliable event rate, (2) the number of events with M ! 7 from the long-term catalog, N long !7 , as Pacheco  Table S1). Light red lines represent plate boundaries.
Geochemistry, Geophysics, Geosystems and Sykes [1992] claim their catalogue is quite complete for events above M7, although they only include events deemed to be main shocks. M MME is the average moment per event, converted to moment magnitude. This average magnitude should be high for regions with low b-value and low for regions with a high b-value.
We plot M MME as a function of the Deborah number in Figure 13b. We indeed get a trend of increasing M MME with increasing De, which is better defined when using N long !5:5 (R50.62), than with N long !7 (R50.38). These trends are the result of well-defined negative correlations between De app and N long !5:5 and De app and N long !7 , and a scattered positive correlation between M long 0 and De app . Given the large uncertainties in the long-term catalogue, much better trends are probably not to be expected. Hence the M MME trend seems to further support a correlation between De and b-values.
There are various caveats to these b-value analyses, as the number of events per zone, ranging from 30 to 1000, is on the low side, we did not attempt to select events based on style of deformation and there has been discussion about whether variable b-values are statistically resolvable from the global earthquake catalogue and whether apparent variability is not just due to the complex nature of earthquake patterns [Kagan, 2010]. More work is required to test this further. However, we consider it significant that the trends between b-values and De app are shown both by the shorter and longer catalogues, and the type of correlation we find is as physically expected due to the role of elasticity.

Discussion
Many previous studies investigated correlations between seismic activity (e.g., seismicity rates, M max or moment release rates) with subduction parameters. Few simple correlations are found, showing that the interplay between driving forces and frictional properties that governs earthquake potential is complex [e.g., Pacheco et al., 1993;Heuret et al., 2011]. A correlation with downgoing plate age or thermal parameter, which governs slab pull as well as plate rheology, has been suggested [Ruff and Kanamori, 1980;Kanamori, 1986], but is not statistically strong [Pacheco et al., 1993;Heuret et al., 2011]. In addition, there are correlations with overriding plate stress state and motion [Scholz and Campos, 1995;Conrad et al., 2004;Heuret et al., 2011]. Both trends may be understood considering that the downgoing-upper plate force balance controls normal stress and thereby locking and seismic coupling across the fault interface [Scholz and Campos, 1995]. Conrad et al. [2004] furthermore found that seismic coupling (the ratio of seismic over total interplate strain) correlated with the fraction of upper mantle slab pull they required to reproduce plate motions in a global dynamic model, such that subduction segments with stronger coupling had lower slab pull fractions. and either (green circles) the number of events with M w ! 5:5 from the CMT catalogue scaled to 107 years, or (black circles) the number of events M w ! 7 from the historic catalogue [Engdahl and Villaseñor, 2002]. Dashed lines are the best linear fits. Higher estimates of De appear to correlate with a stronger contribution of large earthquakes to the total seismicity rate.
Geochemistry, Geophysics, Geosystems The subduction segments with the strongest seismic coupling estimates, Nankai, southern Chile, Eastern Aleutians/Alaska, and Mexico-Colima [Scholz and Campos, 1995], are also the segments where we find the highest De app , while those with the weakest seismic coupling, Southern Marianas, Kermadec, part of the Kuriles and Java-Timor, are also those where we find the lowest De app . We would interpret this as follows: high De app are the result of a balance of upper and lower plate forcing which leads to (a) high normal stress and strong seismic coupling [Scholz and Campos, 1995;Heuret et al., 2011], and (b) relatively mild and hence significantly elastic bending, (where R VE is not much less than half R E , the radius of curvature expected for an elastic plate with age-dependent thickness subsection 5.2). Garel et al. [2014] showed that variations in upper plate forcing lead to a larger range of subduction styles for younger, lower slab pull, plates than older ones. This stronger sensitivity of young plates to outside forcing could explain the large range of R VE =R E for young relative to old plates.
However, while in terms of bending the lower De app zones could be considered weaker, Conrad et al. [2004] found that the low-coupling zones required high slab pull transmission. This further argues for the apparent bending weakness of older plates to be a consequence of plastic yielding rather than low effective viscosity [Buffett and Heuret, 2011], as this would still allow force transmission through the remaining strong core. By contrast, the strongly coupled, high De app plates may be dissipating a more significant amount of their potential energy in their interaction with the upper plate, as found in some of the dynamic models by Arcay et al. [2008], thus resulting in a lower fraction of transmitted slab pull to the unsubducted plate. This again emphasizes that the interplay between subduction forces and rheology that governs both large scale dynamics and earthquake potential is complex.

Conclusions
Although often neglected, the role of elasticity during the bending of subducting slabs into the trench is uncertain and potentially important in controlling plate velocities and subduction styles. Using a set of freesubduction models with elastic, viscoelastic, and viscoelasto-plastic rheologies, we analyzed the role of elasticity in plate bending at the trench. Results are characterized as a function of the Deborah number, De, the ratio of Maxwell over deformation time, in our case defined as the time the plate takes to move through the bend.
Our model slabs behave essentially viscous for De less than 10 23 , and largely elastic for De above 0.1. For De $ 0:1, slabs store about 10% of the energy in the system elastically, and they require 90% less energy to achieve a given bending radius than purely viscous slabs with same viscosity, thickness and buoyancy force. At De exceeding 10 22 , elasticity facilitates unbending to the extent that subduction in retreating modes (where the slab flattens in the transition zone), becomes dominant over modes where the trench advances or stagnates (which would cause folding slabs).
From the viscoelastic models and elastic models, we derived a geometrical scaling relation that allows us to estimate apparent Deborah numbers from comparing measured subducting-plate minimum radii of curvature with radii expected for age-dependent elastic thickness. Models including plasticity show that these De app can underestimate the actual De by up to an order of magnitude if local weakening occurs. Using the bending radius compilation by Rose and Korenaga [2011], we obtain De app less than 10 22 for most subduction zones. However, if a significant number of slabs plastically yield during bending, as previously suggested [e.g., Levitt and Sandwell, 1995;Billen and Gurnis, 2005;Buffett and Heuret, 2011;Buffett and Becker, 2012;Arredondo and Billen, 2012], the actual De may exceed 10 22 in as much as 60% of the transects we evaluated, and elasticity is more important than often considered.
Our analyses return a large range of apparent De. The largest range, including the highest De, are found where relatively young plates are subducting, and the lowest De for the oldest plates. This can be explained if older plates have a higher chance of reaching their yield stress due to higher slab pull forcing. For young plates, with lower slab pull, external forces can lead to a wider range of bending radii and hence De. Do note that if low apparent De are due to localized weakening, older plates may appear weak in bending, but could still transmit slab pull efficiently through an intact strong core.
A further confirmation that elasticity in slab bending may be more important than generally considered, is given by the systematic decrease in b-values (a measure of the relative rate of small versus large earthquakes) with increasing De app that we found in instrumental and, with more scatter, in historic catalogues.