## Mega-earthquakes go the flat way

Megathrust faults in subduction zones cause large and damaging earthquakes. Bletery *et al.* argue that certain geometric features of the subduction zones relate to earthquake size. The key parameter is the curvature of the megathrust. Larger earthquakes occur where the subducting slab is flatter, providing a rough metric for estimating where mega-earthquakes may occur in the future.

*Science*, this issue p. 1027

## Abstract

The 2004 Sumatra-Andaman and 2011 Tohoku-Oki earthquakes highlighted gaps in our understanding of mega-earthquake rupture processes and the factors controlling their global distribution: A fast convergence rate and young buoyant lithosphere are not required to produce mega-earthquakes. We calculated the curvature along the major subduction zones of the world, showing that mega-earthquakes preferentially rupture flat (low-curvature) interfaces. A simplified analytic model demonstrates that heterogeneity in shear strength increases with curvature. Shear strength on flat megathrusts is more homogeneous, and hence more likely to be exceeded simultaneously over large areas, than on highly curved faults.

Past mega-earthquakes, such as the magnitude (*M*) 9.6 Chile earthquake in 1960 and the *M* 9.3 Alaska earthquake in 1964, occurred in areas where the subducting lithosphere was relatively young (and buoyant) and the plate convergence rate was relatively high (*1*). These observations led some authors to hypothesize that maximum earthquake size is controlled by these two geological parameters (*2*, *3*). The development of space-based geodesy enabled refined measurements of plate motion that challenged the role of convergence rate (*4*–*6*). Additionally, the moment magnitude (*M*_{w}) 9.0 Tohoku-Oki earthquake (*7*) ruptured lithosphere that is over 120 million years old (*8*), ruling out lithospheric age as the dominant control. Weak correlations appear in recent data sets among a variety of parameters, including forearc structure (*9*, *10*); age, density, and buoyancy of the slab (*6*); upper plate motion (*11*); upper plate strain (*12*); long-term trench migration (*11*); trench sediment thickness (*12*); and width of the seismogenic zone (*11*, *13*, *14*). However, aside from dip angle and seismogenic width [coefficient of correlation *r* = 0.61 (*14*)], none of these parameters correlate strongly (|*r*| < 0.5) with the maximum earthquake magnitude recorded in the different zones (*11*). Stronger correlations are found with fault properties such as thrust interface dip angle (*11*) and apparent friction derived from heat flow measurements (*15*).

Today, forecasts of potential large earthquakes focus mainly on imaging the slip deficit with respect to the relative plate motion, because such deficits are thought to result in stress loading that is ultimately released in earthquakes [e.g., (*16*)]. We followed a complementary approach by analyzing large-scale geometrical features of subduction faults and assessing their possible influence on the physical conditions favoring large earthquake ruptures. Seismic moment *M*_{0} ∝ Δσ*S*^{3/2} (*17*) depends mainly on the surface area *S* over which an earthquake ruptures, because the stress drop Δσ (the difference in stress before and after an earthquake) is roughly constant among earthquakes over a large range of magnitudes (*18*). Here we demonstrate that slab interface curvature exerts a leading-order control on the spatial extent of potential ruptures in subduction zones and hence on the magnitude of the largest earthquakes.

Data sets from active- and passive-source seismology have been combined to constrain the slab1.0 model for the geometry of the world’s major subduction zones (*14*) [details and limitations of the slab1.0 model are discussed in the supplementary materials (*19*)]. We computed maps of the along-dip interface curvature *K _{s}* =

*d*θ/

*ds*, hereafter referred to as curvature [where θ is the dip angle and

*s*is the tangent to the interface pointing in the down-dip direction (

*19*)], using the slab1.0 model (

*14*,

*20*) (Fig. 1). Because

*K*is almost always positive (Fig. 1), we hereafter interchangeably refer to low-

_{s}*K*fault regions as flat, low-curvature, or planar. Similarly, we computed the along-strike gradient of the dip angle

_{s}*K*=

_{t}*d*θ/

*dt*(fig. S1). Comparison with a catalog of historical events (

*19*) reveals a tendency for mega-earthquakes to occur on relatively flat (low-

*K*) megathrusts (Fig. 1). The curvature is particularly small in the Japan-Kuriles-Kamchatka, Alaska-Aleutians, Sumatra-Java, South America, and Cascadia subduction zones, which are known to produce

_{s}*M*≥ 9.0 earthquakes. Subduction zones with large curvatures, such as the Philippines, Solomon Islands, Izu-Bonin, Santa Cruz–Vanuatu-Loyalty, and Tonga-Kermadec zones, lack historic mega-earthquakes. At a smaller scale, we also observe that the down-dip limit of mega-earthquakes offshore of Sumatra—and to a lesser extent, in the Aleutians, Alaska, Nankai, and Kamchatka—coincides with an abrupt change in the slope, suggesting that the conditions required to generate large earthquakes are directly related to the local curvature along the megathrust interface.

To further explore the role of subduction geometry, we calculated the average dip angle (Fig. 2A), the average curvature (Fig. 2B), and the average along-strike gradient of dip angle (fig. S2) from the trench to 60 km depth, and we compared these values with the magnitude of the largest megathrust earthquake *M*_{max} recorded in each subduction zone. We found that *M*_{max} anticorrelates (*r* = –0.72) with (Fig. 2A) and anticorrelates even more strongly (*r* = –0.80) with (Fig. 2B). The anticorrelation with is weaker (*r* = –0.64) but still significant (fig. S2). We suggest, on the basis of our observations, that planar slab interfaces are prone to hosting larger earthquakes. Multiple factors, including the subducting (*21*) or overriding (*22*) plate thickness and viscosity contrast (*23*), have been invoked as controls on subduction curvature; these studies (*21*–*23*) provide a useful introduction to this field of research.

The average curvature within mega-earthquake slip contours is consistent with the linear regression in Fig. 2B, with observed values always smaller than 3.76 × 10^{–6} m^{–1} (fig. S5B). Given the distribution of *K _{s}* across the different megathrusts of the world, the likelihood that all known mega-earthquakes occurred by chance on

*K*≤ 3.76 × 10

_{s}^{–6}m

^{–1}areas is 0.8% (

*19*), allowing us to affirm with >99% confidence that earthquake magnitude is related to megathrust curvature. The distribution of the average dip angle within mega-rupture contours is much broader (fig. S5A), confirming the stronger relationship of earthquake size with curvature than with dip angle (

*19*). Nevertheless, known mega-earthquakes do not exhibit a strong tendency to rupture the lower-curvature portions of their host subduction zones (

*19*). This suggests that over the long term, mega-earthquakes might rupture any portion of a low-curvature subduction thrust, such as those of the Cascadia, South America, Alaska-Aleutians, Sumatra-Java, Japan-Kuriles-Kamchatka, Ryukyu-Nankai, and Central America subduction zones. Large earthquakes are preferentially hosted by megathrusts with low dip angles in part as a consequence of larger down-dip seismogenic extents making ruptures possible over wider fault areas (

*11*,

*13*,

*14*). However, the stronger relationship of earthquake size with curvature than with dip angle noted above (Fig. 2B) warrants a more thorough description of fault loading.

The present understanding of seismic ruptures can be framed in terms of the asperity model: Earthquakes are caused by the sudden failure of locked macroscopic asperities that accumulate stress during interseismic periods (*24*). In this context, the extent of an earthquake has a leading-order dependence on the size of the rupturing asperity and the neighboring region that is pushed to failure by coseismic stress change. Recent efforts aimed at characterizing the nature of asperities and their immediate surroundings have mostly focused on potential variations in friction along faults [e.g, (*25*–*28*)]. The capacity of a fault segment to accumulate stress is bounded by the Coulomb failure criterion (1)where τ^{c} is the critical shear stress required to initiate a rupture (hereafter referred to as shear strength), μ is the coefficient of friction, *p *is the pore pressure, is the normal stress when the rupture initiates, and cohesion is neglected. The Coulomb failure criterion is met when the accumulated shear stress reaches the frictional shear strength and an earthquake initiates. From Eq. 1, one can define an asperity as a fault area characterized by a particularly high friction coefficient μ. Alternatively, one can also explain the time-dependent behavior of asperities by appealing to variations in pore pressure *p* that accompany fluid migration [e.g., (*29*)]. We explore here a different hypothesis, which is that the locations of asperities can be determined by analyzing variations in normal stress σ_{n} that are governed by large-scale geometrical characteristics of megathrust interfaces. Assuming constant friction μ and hydrostatic pore pressure *p* = ρ_{w}*gh* (with ρ_{w} being the water density, *g *the acceleration of gravity, and *h* the local interface depth), we illustrate the effect of variations in normal stress σ_{n} on variations in shear strength [this method is reminiscent of the consideration of normal force in (*16*) but invokes different causative factors, i.e., curvature instead of plate tectonic forces].

We considered a simple two-dimensional (2D) fault model: a fault interface subjected to the pressure of the upper plate mass and an unknown horizontal tectonic stress. We also treated the convergence direction as aligned with the horizontal projection of the local along-dip direction (i.e., pure dip-slip convergence). This latter simplification is justified by the dominantly dip-slip mechanisms of most megathrust earthquakes [e.g., figures 14, 15, and 16 of (*30*)]. Even though subduction zones present, to varying degrees, oblique convergences, this obliquity is usually accommodated for the most part by large strike-slip faults in the back arc. One notable exception is the Solomon Islands subduction zone, which presents a very strong oblique convergence and hosts megathrust earthquakes that typically exhibit large strike-slip components. Our idealized description of the convergence direction allows us to treat the 3D problem with a 2D model (fig. S6). Within this framework, the shear strength—defined as the critical shear stress required to initiate slip at a given location—can be expressed as a function of the crustal density ρ, the friction μ, the depth *h*, the dip angle θ, and the angle ψ between the maximum principal component of stress and the horizontal (2)τ^{c} increases with θ () (*19*)meaning that large dip angles imply greater shear strength. If the dominant control of earthquake magnitude were the amplitude of shear strength, we might expect a positive correlation between *M*_{max} and θ. The negative correlation that we observe (Fig. 2A) points to other factors.

Because the magnitude of an earthquake is primarily controlled by the area of rupture, homogeneous distributions of shear strength over large fault areas may favor the occurrence of mega-earthquakes. In this framework, the shear-strength gradient *d*τ^{c}/*ds* is a critical parameter. Treating μ, ρ, ρ_{w}, and *g* as constants, spatial variations of τ^{c} obtained by differentiation of Eq. 2 satisfy (3)where ρ > ρ_{w} and the functions *A*_{μ}(θ,ψ), *B*_{μ}(θ,ψ), and *C*_{μ}(θ,ψ) are positive for μ = 0.6 and values of θ and ψ typically encountered in subduction zones (*19*). Therefore, |*d*τ^{c}/*ds|* increases with *K _{s}* [a similar argument can be made to show that |

*d*τ

^{c}/

*dt|*increases, on average, with |

*K*| (

_{t}*19*)]. Shear-strength heterogeneity increases with curvature: The flatter the megathrust interface (lower

*K*), the more homogeneous the shear-strength distribution. Stress accumulation along faults is complex, and heterogeneities in driving stress can result from the effect of stress concentrations inherited from past events or changes in coupling associated with variations in the friction coefficient or pore fluid pressure. However, dynamic stress perturbations induced by earthquake propagation are more likely to overstep the resistance to failure over broad areas if shear strength is homogeneously distributed along a large fault surface. In contrast, barriers to earthquake propagation are more likely for heterogeneous shear-strength distributions, because the dynamic stress perturbations required to reach local values of shear strength will tend to be larger in some regions. The idea that geometrical heterogeneity may limit earthquake propagation has been previously invoked (

_{s}*15*,

*31*) and is consistent with the observation that many large earthquake ruptures terminate near subducting seamounts or other structural heterogeneities (

*31*,

*32*).

During an earthquake, fault slip will propagate as long as the dynamic stress perturbation induced by the earlier phase of the rupture, δτ, is larger than the difference between shear strength and shear stress in the surrounding area, Δτ. Hence, one way to generate mega-earthquakes is to initiate a large dynamic stress perturbation by the failure of a strongly locked asperity so that δτ > Δτ over a large surrounding fault area. In this scenario, mega-earthquakes are characterized by the rupture of one or several asperities where shear strength is particularly high. Our finding that large earthquakes rupture flat megathrusts suggests that another way to generate mega-earthquakes is if Δτ is small in a large area surrounding rupture nucleation, and thus even relatively small dynamic stress perturbations δτ can continue to propagate rupture. The 2004 Sumatra-Andaman earthquake, a unilaterally propagating 1600-km-long rupture without evidence of a localized high stress drop (*33*), exemplifies this second scenario, because it is unlikely that the energy released during the early phases of rupture could have initiated slip on fault portions as far as 1600 km away if the entire ruptured area was not already close to failure. In this situation, mega-earthquakes are not characterized by the rupture of asperities in regions where strength is much greater than in the surrounding area but, on the contrary, by the absence of strong variations in strength. The limiting case of homogeneous null shear strength results in constant creep and the absence of earthquakes of any magnitude. Hence, although the amplitude of τ^{c} is an important parameter, the observation that mega-earthquakes preferentially occur on flat subduction zones (Fig. 2) seems to indicate that the homogeneity of the shear-strength distribution is among the most critical factors in enabling the generation of mega-earthquakes.

It has recently been proposed that any subduction zone may produce a *M* ≥ 9 earthquake (*34*). Our results suggest that earthquake size is limited by curvature and, if our interpretation is correct, mega-earthquakes might be physically incapable of rupturing highly curved subduction zones with large shear-strength gradients, such as the Philippines, the Solomon Islands, Scotia arc, and Santa Cruz–Vanuatu-Loyalty; heterogeneous shear strength creates natural barriers that stop earthquake propagation. It is possible that heterogeneous shear stress could build to match a heterogeneous shear-strength distribution and allow rare large events, but such scenarios are expected to be less likely, and therefore more infrequent, than on megathrusts with more homogeneous shear strength. As indicated by Fig. 2B, some subduction zones, such as Izu-Bonin or Central America, may yet host earthquakes larger than previously recorded. Cascadia is noteworthy because the last mega-earthquake that it hosted is thought to have ruptured the entire subduction fault and therefore reached the maximum possible magnitude for this subduction zone (*35*). At a smaller scale, some areas, such as Peru, Java, or the large low-to-negative *K _{s}* fault region in Central America (Fig. 1), which includes the Guerrero gap but extends much farther, show favorable features for a possible very large rupture—implying in the Guerrero case the potential for an even larger earthquake than that inferred from the extent of the gap. Moreover, large flat portions of subduction faults may sometimes rupture as one mega-earthquake and sometimes as several smaller earthquakes. Such behavior is documented, for instance, in Nankai (

*36*), where two nearby asperities sometimes break as two separate consecutive earthquakes (as in 1854 and 1944–1946) and sometimes as one larger earthquake (1707; blue contour in the Ryukyu-Nankai box). This suggests that high-curvature regions of generally flat subduction zones may act as barriers to rupture most of the time and still sometimes be overcome by coseismic stress changes. Planar fault areas may thus host moderate-sized earthquakes for a long time and still have the potential to generate mega-earthquakes in the future.

It is possible that frictional properties and geometry are related. Such mechanical relationships are well described in accretionary wedges (*37*) and are certainly expected, though more difficult to constrain, along the deeper portions of megathrusts. Slab curvature has also been proposed to promote higher permeability (*38*–*41*) and thus to have an impact on the pore fluid pressure as well. Nevertheless, we have shown that spatial variations in friction and pore fluid pressure are not required to explain the gross distribution of mega-earthquakes. The observation that mega-earthquakes preferentially rupture flat megathrusts (Figs. 1 and 2B) is consistent with the inference that shear strength tends to be more homogeneously distributed along such subduction interfaces (Eq. 3), facilitating synchronized failure over large areas. This implies that the critical feature at play in the generation of mega-earthquakes is not the amplitude of shear strength but its spatial variations. Thus, the absence of asperities on large faults may counterintuitively be a source of higher hazard. Though our study focused on subduction earthquakes, flatness may favor large earthquakes on long strike-slip faults as well.

## Supplementary Materials

www.sciencemag.org/content/354/6315/1027/suppl/DC1

Materials and Methods

Figs. S1 to S8

Tables S1 to S2

## References and Notes

**Acknowledgments:**We thank R. Bürgmann for his valuable comments on the manuscript. The slab1.0 model is available online at earthquake.usgs.gov/data/slab/. The U.S. Geological Survey catalog of historical earthquakes (

*M*≥ 8.0 since 1900) that we used in this study can be found at http://earthquake.usgs.gov/earthquakes/world/historical.php. For the subduction zones that have not experienced any

*M*≥ 8.0 earthquakes since 1900, this catalog was complemented by the Global Centroid Moment Tensor catalog (www.globalcmt.org/). We used Generic Mapping Tools to compute the distributions of dip-angle gradients (gmt.soest.hawaii.edu/). This work was supported by NSF grant EAR-1520238, ANR project TO-EOS, and the French Ministry of Research and Education.