The link between great earthquakes and the subduction of oceanic fracture zones

Abstract. Giant subduction earthquakes are known to occur in areas not previously identified as prone to high seismic risk. This highlights the need to better identify subduction zone segments potentially dominated by relatively long (up to 1000 yr and more) recurrence times of giant earthquakes. We construct a model for the geometry of subduction coupling zones and combine it with global geophysical data sets to demonstrate that the occurrence of great (magnitude ≥ 8) subduction earthquakes is strongly biased towards regions associated with intersections of oceanic fracture zones and subduction zones. We use a computational recommendation technology, a type of information filtering system technique widely used in searching, sorting, classifying, and filtering very large, statistically skewed data sets on the Internet, to demonstrate a robust association and rule out a random effect. Fracture zone–subduction zone intersection regions, representing only 25% of the global subduction coupling zone, are linked with 13 of the 15 largest (magnitude Mw ≥ 8.6) and half of the 50 largest (magnitude Mw ≥ 8.4) earthquakes. In contrast, subducting volcanic ridges and chains are only biased towards smaller earthquakes (magnitude


Introduction
Earthquake supercycles (Sieh et al., 2008) occur on timescales of up to or beyond 1000 yr (Gutscher and Westbrook, 2009), defying prediction using traditional methods (Stein et al., 2012).For instance, no instrumentally recorded great (moment magnitude M w ≥ 8) subduction zone earthquake has occurred along the Cascadia margin (Fig. 1), but there is evidence for 13 events in the last 7500 yr with average repeat times of ∼ 600 yr, including a magnitude 9 event on 26 January 1700 (Goldfinger et al., 2003).There are many other regions with a history of great subduction earthquakes and "supercycle" recurrence (Gutscher and Westbrook, 2009).These regions are not adequately represented in traditional earthquake hazard maps, leading to a failure to predict locations of giant earthquakes based on these maps (Stein et al., 2011(Stein et al., , 2012)).Digital earthquake catalogues combined with the characteristics of regional fault systems do not allow reliable differentiation of regional risk levels if earthquake cycles are up to an order of magnitude longer than the ∼ 100 yr time span covered by these catalogues.An alternative method to forecast long-term seismicity is based on the global strain rate map (Bird et al., 2010), but regional differentiation between high-risk and low-risk areas for great earthquakes is poor.This problem has given rise to the use of probabilistic methods such as Monte Carlo methodologies (Parsons, 2008) to fit wide ranges of distribution parameters to short paleoseismic series.Lay and Kanamori (1981) developed a conceptual model in which major subduction zone earthquakes are driven by strong coupling between the downgoing and overriding plates, driven by the subduction of asperities, i.e. aseismic ridges on the downgoing plate which cause strong coupling at the plate interface.Here we combine  (Amante et al., 2009), on a Robinson projected map: subduction coupling zones (blue bands) (see text for coupling zone model description), oceanic fracture zones (dark gray) and oceanic volcanic chains and aseismic ridges (pink) (Matthews et al., 2011), intersection points of fracture zones with subduction zones (yellow squares), intersection points of the volcanic chains and ridges with subduction zones (green squares), largest 15 instrumentally recorded earthquakes (M w ≥ 8.6) (red stars), largest 50 earthquakes (M w ≥ 8.4) (light blue circles), and all other significant earthquakes (small beige circles) (NGDC/WDC, 2011).Inv FZ, Investigator Fracture Zone, Krus FZ, Krusenstern Fracture Zone, Phil SZ, Philippine Subduction Zone, Mar SZ, Marianas Subduction Zone, SW-P SZ, Southwest Pacific Subduction Zone, Jap SZ, Japan Subduction Zone, Kam SZ, Kamchatka Japan Subduction Zone, Aleu SZ, Aleutian Subduction Zone, Cas SZ, Cascadia Subduction Zone, C-Am SZ, Central America Subduction Zone, S-Am SZ, South America Subduction Zone, L-An SZ, Lesser Antilles Subduction Zone.
this conceptual approach with a quantitative analysis involving an integrated set of global digital geophysical data sets to investigate spatial associations between significant earthquakes and different types of subducting asperities.The effect of aseismic ridge and seamount subduction on seismic coupling and earthquake rupture behavior and overriding plate deformation has been investigated at many localities (Das and Watts, 2009).A detailed study of the tec-tonic setting along the Japan Trench (Mochizuki et al., 2008) led to the conclusion that subducting seamounts are associated with weak interplate coupling.This observation has not been tested globally, but casts doubt on the idea that volcanic edifices on ocean crust are the most obvious candidates for barriers that locally inhibit faulting for long periods of time, leading to great earthquake supercycles.Oceanic fracture zones represent another form of subducting asperities that are quite different from volcanic edifices in that they are often accompanied by elevated, continuous ridges that represent uplifted edges of normal ocean floor (Sandwell and Schubert, 1982).Their effect on earthquake rupture has been investigated regionally, for instance along South America (Contreras-Reyes and Carrizo, 2011;Robinson et al., 2006;Carena, 2011), Alaska (Das and Kostrov, 1990), Sumatra (Ammon et al., 2005) and the Solomon Islands (Taylor et al., 2008), but not globally.A recent global digital fracture zone data set based on vertical gravity gradients derived from satellite altimetry data (Matthews et al., 2011) reveals a total of 59 fracture zone-subduction zone intersections.It is striking to observe that many of these are in close proximity to great earthquake epicentres (Fig. 1), raising the question as to whether this observation is supported by a statistically robust association, and ultimately a physical link.This data set includes the Kashima Fracture Zone, whose landward extension straddles the location of the 11 March 2011 Tohoku-Oki earthquake (Fig. 1).This fracture zone is well expressed in offsets of marine magnetic anomalies (Nakanishi et al., 1992) and appears as a linear feature in the vertical gravity gradient derived from satellite altimetry (Matthews et al., 2011).It is characterised by a trough bounded by two ridges that are up to 2 km above the surrounding seafloor (Nakanishi, 1993), similar to many other major fracture zones (Sandwell and Schubert, 1982).It may be expected that the subduction of prominent fracture zone ridges affects long-term seismic coupling and thus an elevated probability of seismicity.Fracture zones are associated with ridges as much as 3 km above the surrounding abyssal seafloor (Bonatti, 1978;Sandwell and Schubert, 1982;Bonatti et al., 2005;Collete, 1986;Wessel and Haxby, 1990), potentially leading to enhanced coupling between the downgoing and overriding plate, and sustained for long periods of time.

Overview
Our analysis assesses the spatial association between shallow subduction-based earthquakes and the location of intersections between fracture zones or volcanic ridges/chains and subduction zones, computed as a function of earthquake magnitude.The methodology is broken down as follows: (1) a significant earthquakes catalog is filtered to include only those events constrained within the coupling zone defined at the intersection between subducting plates and the overriding lithosphere; (2) a spatial data set is derived in which existing global digital fracture zone and volcanic ridge/chain data sets are intersected with the subduction coupling zones to form a basis for the association analysis; (3) we assess the association strengths as a function of earthquake magnitude using a methodology called "Top-N" analysis, well suited for analysing skewed earthquake magnitude distributions and (4) the sensitivity of the associations to the arbitrary case (or assessing the Null hypothesis, i.e. the case in which there is no association) is computed; and (5) the computed associations are expressed spatially, indicating regions fitting the hypothesis.Finally, we discuss our approach to assess the effect of convergence rates.

Analysing shallow subduction-based earthquakes in subduction coupling zones
The majority of known mega-thrust earthquakes are known to occur along subduction zones at relatively shallow depths at the coupling interface between the overriding and downgoing plates.We construct a subduction coupling zone model by combining the Slab1.03-dimensional global subduction zone model (Hayes et al., 2012) with lithospheric thickness model TC1 (Artemieva, 2006) , 1998).The coupling zone of remaining slabs not covered in either model was constructed with an extent of 150 km from their surface expression following Bird (2003) and using slab dip angles from Lallemand et al. (2005).The reasoning behind this spatial partitioning is that it forms a constrained physical boundary in which megathrust earthquakes are expected to originate, and is not constrained by known spatial extents of pre-recorded ruptures.This is due to the long periodicities of larger events leading to poor representivity in earthquake catalogs.
The NGDC (NGDC/WDC, 2011) significant earthquakes catalog is used in this study, consisting of a monolithic catalog of events skewed towards high magnitude earthquakes and including the most recent events.This catalogue is well suited for our analysis as it is up-to-date and complete, and we are less concerned with relatively small epicentre location errors (Engdahl et al., 1998) because our targeted associations occur on larger spatial scales (of the order of ∼ 100 km).We only used earthquakes in NGDC's for post-1900 events, reducing the entire data set from 5539 to 3157 earthquakes.Earthquake magnitudes are determined using the moment magnitude scale (McCalpin, 2009) for 761 events, whereas for the remaining events we use the maximum magnitude available, considering that the magnitudes of older events, determined from obsolete scales, are generally underestimated for large magnitudes (see Appendix A).These measures include the Richter scale, which underestimates large earthquake magnitudes; surface wave magnitude, with moderate improvements over the Richter scale; and body-wave magnitude, which is accurate only for smaller www.solid-earth.net/3/447/2012/Solid Earth, 3, 447-465, 2012 R. D. M üller and T. C. W. Landgrebe: Great earthquakes and oceanic fracture zones magnitude events.Some very old records use an intensity scale, which we regard as too inaccurate for our purposes.Subsequent filtering and isolation of events originating from the subduction coupling zone reduces the data set to 1486 observations.

Intersections between fracture zones and volcanic ridges/chains with subduction zones
The analysis relies on the identification of both fracture zones and volcanic chains/aseismic ridges that occur in the vicinity of subduction zones.Fracture zone-subduction zone intersections using fracture zone identification from Matthews et al. (2011) were flagged automatically, while a combination of bathymetry and gravity anomaly data were used to assess fracture zone locations within close proximity to subduction zones, taking into account that sediments on the downgoing plate seaward of the trench may partly obscure bathymetric expressions of fracture zones.This resulted in a total of 59 identified intersection points (Fig. 1).Volcanic chains and aseismic volcanic ridges have been compiled based on Coffin and Eldholm (1994), and subduction zone intersections were computed as in the fracture-zone case.Features on the seafloor in the proximity of subduction zones were classified to be either in the process of being subducted or not.A total of 14 locations were identified (Fig. 1).
The data set selection in this study thus comprises large, well-defined bathymetric features in the vicinity of the various subduction trenches.The reasoning is that the very large, consistent bathymetric anomalies may be playing a substantial role in increasing subduction coupling and/or creating the necessary conditions for temporary "locking" within the coupling zone.Smaller, less well-defined features identified in geophysical data are excluded, also including minor neartrench fractures due to uplift.
We combine the fracture zone and volcanic chain/aseismic ridge subduction boundary intersections with the filtered earthquake database, and consider the spatial domains formed by our CouplingZone1.0model.We project fracture zone-subduction zone (FZ-SZ) intersections onto the subduction coupling zone along the axis of the fracture zone, resulting in linear traces spanning the width of the coupling zone.These coupling zone intersections form the basis for the association calculations, with the regions adjacent to these fracture zone traces within the subduction coupling zone analysed as a function of the perpendicular distance away from them.In this way, a data-driven association analysis is undertaken in which earthquake associations can be computed adaptively in terms of the estimated 3-dimensional nature of both the subducting and over-riding plates, and taking the broad orientation of the subducting features into account.Selected widths are used to form zones centred on the linear traces, in which the regions on either side of the lines are used to investigate proximal earthquakes.The approach taken here is to assess the association sensitivity to a variation in these widths.

Computing coupling-zone spatial associations via Top-N analysis
The analysis presented in this paper investigates magnitude relationships between earthquake locations in the vicinity of fracture zone and volcanic chain/ridge intersections with the coupling zone.The structure of the bathymetric features was used to project their extensions into the nearby coupling zone, maintaining the same azimuth as in their oldest geophysical expressions seaward of a given trench.The resultant intersection is bounded by the width of the coupling zone.This results in linear spatial features that serve as a reference for undertaking the analysis of associated earthquakes for a range of proximities, allowing the sensitivity of the association to be quantified.In Fig. 7 fracture-zone intersections with subduction zones are shown for a 100 km buffer region, demonstrating how the spatial associations undertaken for this study are computed.The buffer regions are progressively increased in size to trade-off the strength of the association with the specificity of the targeted area, calibrated against the entire coupling zone area.
We apply a type of information filtering system technique called "Top-N analysis", which is widely used in searching, sorting, classifying, and filtering very large data sets (Cremonesi et al., 2010) to investigate the association of significant earthquakes as a function of magnitude with subduction coupling zone segments with and without fracture zone or volcanic ridge/chain intersections.This approach suits the magnitude distribution of earthquakes, where larger events are more infrequent.The Top-N analysis progressively assesses how strongly particular sets of subduction zone segments are associated with sets of sorted earthquakes in a given magnitude range.As the total number of Top-N earthquakes, from the largest event down to a given cutoff magnitude, is increased by progressively including smallermagnitude events, the so-called recall is computed, defined as the number of Top-N earthquakes associated with the target coupling zone regions divided by N. The resultant statistical measure represents an intuitive description of the effectiveness of a given target data set (fracture zone-subduction zone (FZ-SZ) intersections, volcanic ridge/seamount intersections, or randomly chosen subduction coupling zone segments) in accounting for the location of significant earthquakes on record.
The analysis is described formally as follows: the significant earthquakes data set is denoted E, consisting of a list of latitude (θ), longitude (λ) and magnitude (M) tuples such that E = [e 1 , e 2 , . . ., e N e ] for a data set size of N e , and the i-th element of E is denoted e i = (θ e i , λ e i , M e i ).In this analysis E is then re-ordered by sorting it in descending order in terms of magnitude, resulting in E s .Target locations are projected into the coupling zones as described previously, resulting in lines traversing the coupling zone.The list of N t projected target lines pairs is defined as: (1) The Top-N methodology involves computing the ratio of the highest-N earthquakes within a specified region of interest (ROI), specified by a buffer distance d ROI in kilometres around L t , denoted the Top-N recall.When studying the recall for a particular N, the performance evaluation score is defined as recall-N, and is calculated by summing the number of N sorted earthquakes (i.e.[ e 1 , e 2 , . . ., e N ]) associated within the ROI regions.Formally this procedure involves creating a binary vector A of length N, defined as A = [a 1 , a 2 , ..., a N ].The i-th item in A is determined as follows: where CZ is the coupling zone polygon geometry, and F considers the association within a thresholded buffer around L j bounded by CZ.
The function G(L j , CZ, d ROI ) creates a buffer polygon of width d ROI km around L j , clipped by CZ.The target earthquake e s(i) can then be tested for association via a standard point-in-polygon test, denoted inpolygon.Now that the binary vector A has been computed, the recall-N score can be determined via

Sensitivity analysis: ruling out random effects
The study objective is to evaluate the nature of significant earthquakes in the vicinity of target locations along subduction zones.It is important to compare the nature of these associations with the "random-association" case (which we call the arbitrary case) in order to ascertain whether the resultant associations could occur at random, and to calibrate the strength of the association (i.e. the "sensitivity" and "specificity" characteristics).The approach taken is to carry out repeated analyses along subduction zones in which arbitrary target locations are specified.These arbitrary target location selections are repeated 100 times in a stochastic fashion to ascertain respective statistical variability, allowing for a thorough ruling out of an association occurring at random (the Null hypothesis pertaining to the topic of this study).Subduction coupling zone partitions are generated perpendicular to the coupling zone axis at a spacing of 20 km, providing a set of 2634 partitions (Fig. 6).Each partition element is extended along strike of a given subduction zone using widths ranging from 50 to 400 km to enable a sensitivity analysis with regard to spatial proximity of earthquakes and FZ-SZ intersections.
The arbitrary case methodology comprises the extraction of a number of points randomly sampled along subduction zones, drawing the same number of random samples as there are target locations (i.e.59 fracture-zone intersection locations, and 14 intersections pertaining to major volcanic chains/ridges, respectively).Thus, a repeated selection of 59 and 14 virtual target regions is undertaken, followed by the same Top-N association analysis.The process is repeated 100 times (exceeding this amount resulted in consistent statistics), with the statistics from various runs summarised via median, and 20th/80th percentile error bars.The arbitrary case methodology thus repeatedly simulates a similar scenario to the targeted case, providing a relative measure by which the statistical significance of the associations can be ascertained for different proximities from the target zones.

Spatial hazard-zone construction
In order to relate the FZ-SZ intersection regions to a spatially meaningful interpretation, it is important to establish a baseline/reference spatial zone that can be compared with.The spatial zone comprises the regions in which all subduction earthquakes are known to occur, i.e. the coupling zone defined earlier.This allows a direct comparison to be made between the regions in which the associations discussed in this paper are strong, and the entire reference region.This approach also has application to long-term earthquake hazard risk assessment.The baseline coupling zone surface area is computed to be ∼ 1.088×10 7 km 2 , which is about 2.1 % of the Earth's surface area (∼ 51.095 × 10 7 km 2 ).The approach evaluates the proportion of earthquakes associated with particular subduction coupling zone sub-regions (e.g. 150 km wide centred on FZ-SZ intersections).These regions of interest form the FZ-SZ target zones, which can then be compared to the baseline area to compute its "specificity", i.e. to what extent the baseline area has been reduced.This is then assessed together with the association strength ("sensitivity") to obtain a spatially meaningful assessment of the association strength for different buffer-widths across the coupling zones.The computed fraction of the coupling zone area for chosen buffer widths 50, 100, 150, 200, 250, 300, and 400 km, respectively, is presented in Table 1, having been computed taking into account the fact that neighbouring regions can overlap somewhat.

Geodetic plate convergence rates
Present-day inter-plate convergence velocities from the Global Strain Rate Map Project (Kreemer et al., 2003) are used to estimate the trench perpendicular convergence rates and azimuths at subduction zones to assess the role of these parameters in modulating the effect of subduction asperities on generating significant earthquakes.Relative convergence speed and azimuths within 7.5 degrees (in both latitude and longitude dimensions) of target locations are extracted in order to provide sufficient spatial information to accurately define the boundary between the two plates (the grids have a 1 degree resolution).The extracted boundary is subsequently used to compute smoothed convergence velocities.

Intersections between fracture zones and volcanic ridges/chains with subduction zones
The Top-N analysis is initially carried out using a subduction coupling zone segment width of 150 km following the observation that topographic anomalies associated with fracture zones are rarely wider than this (Fig. 8), and accounting for spatial uncertainties in the projected location of asperities in the coupling zone.The analysis reveals a significant relationship between FZ-SZ intersections and large earthquakes (Fig. 2a).For earthquakes with magnitudes larger than 8.0, the Top-N recall for FZ-SZ intersections diverges sharply from the arbitrary case (Fig. 2a), suggesting that there is a relationship between great earthquakes and the subduction of topographic anomalies associated with FZ-SZ intersections, and that this correlation becomes more pronounced with increasing earthquake magnitudes.For the 50 largest earthquakes (M w ≥ 8.4), 50 % are associated with FZ-SZ intersections as opposed to 25 % of randomly selected subduction zone corridors.For the 15 largest events (M w ≥ 8.6), 13 of which are associated with FZ-SZ intersections, the difference rises to about 70 % (Fig. 2b).If fracture zones did not play a special role in driving great earthquakes, we would expect only 3 of the 15 largest subduction earthquakes on record to be associated with coupling zone regions centred on subducting fracture zones.In contrast, the volcanic ridge/seamount chain subduction zone intersections analysed here display a far weaker association, and only for magnitudes less than 8, as compared with randomly selected subduction zone locations (Fig. 2c, d).We note that this data set is relatively small (14 intersections), so the associations do have substantial uncertainties.The robustness of our results is tested via a sensitivity analysis with respect to the width of the coupling zone segments (Fig. 2e, Appendix D).It demonstrates (1) that the 15 largest megathrust earthquakes on record (M w ≥ 8.6) are significantly biased towards FZ-SZ intersection regions and (2) that an intersection corridor width of 150 km represents a threshold.As the coupling segment width adjacent to FZ-SZ intersections is widened from 50 to 150 km, there is a steep increase in the association of the top 15 earthquakes with FZ-SZ intersection coupling zones, whereas a further increase in coupling zone width does not capture any additional events (Fig. 2e).A similar association is still visible for the top 50 (M w ≥ 8.4) earthquakes, with an inflection point at a corridor width of 150 km (Fig. 2e), leading to the conclusion that great earthquakes with magnitudes larger than 8.4 are preferentially associated with subduction coupling zone regions 150 km wide centred on FZ-SZ intersections.The initial Top-N analysis (Fig. 2a, b) suggests that this relationship holds for all great earthquakes (M w ≥ 8.0), but for events with magnitudes between 8.0 and 8.4 the association is not as clearly dependent on the fracture zone corridor width.Previous studies have suggested that subduction earthquake size distributions may depend on convergence rates (McCaffrey, 1994;Gutscher and Westbrook, 2009;Ruff and Kanamori, 1983).Our analysis reveals that all subduction earthquakes with magnitudes of 7 and more at FZ-SZ intersections are associated with relatively high mean convergence rates of about 65 mm a −1 , whereas great earthquakes at FZ-SZ intersections with magnitudes over 8.5 are distinguished by their link with relatively shallow mean slab dip angles (measured directly below a given earthquake epicentre) of around 20 • , but smaller significant earthquakes occur over wide ranges of slab dips, with relatively shallow mean dips of 23-27 • (Fig. 3a).These observations lead to the conclusion that very shallow slab dips (∼ 20 • ) combined with subducting fracture zones give rise to particularly large great earthquakes.

Long-term great subduction earthquake hazard map
Applying the Top-N association analysis in the context of the spatial hazard map shows that 87 % of the 15 largest, half of the largest 50 and 44 % of the largest 100 earthquakes are associated with FZ-SZ intersections, an area restricted roughly to 32 % of the subduction coupling zone (Fig. 5 and Table 1).There is a long-standing controversy over whether the present state of seismological science inhibits reliable differentiation of the risk level in particular geographic areas along subduction zones, or whether the feasibility of assessing at least the long-term regional earthquake potential is promising in principle, based on long-term earthquake activity and the plate tectonic setting (Sykes et al., 1999).McCaffrey (2007,2008) goes as far as to suggest that given a long enough subduction trench and sufficient time, giant earthquakes may well happen at any subduction zone.Our analysis provides overwhelming evidence that FZ-SZ intersections are associated with a significantly elevated probability of long-term great earthquakes as compared to the remainder of subduction segments, confirming that knowledge of the tectonic setting can help differentiate long-term regional earthquake risk (Sykes et al., 1999).Our results therefore provide a way to objectively test long-term earthquake hazard maps, a need discussed extensively after the 2011 Tohoku-Oki event that was not predicted by previous hazard maps (Stein et al., 2011(Stein et al., , 2012;;Avouac, 2011;Geller, 2011).A local subducting fracture zone alone may not be sufficiently noteworthy to suspect a link with seismic hazards without having recognised a strong global link between subducting fracture zones and great earthquakes, as demonstrated here.This connection provides critical additional information for seismologists to pinpoint particular tectonic environments that are more prone to strong seismic coupling and great earthquake supercycles than average subduction settings.Candidate locations (Fig. 5) for such earthquake supercycles will need to be scrutinized in greater detail, including their historical earthquake records and geological data such as co-seismically displaced coral reefs (Taylor, 2011;Sieh et al., 2008) and prehistoric tsunami deposits (Satake et al., 2003).We identify a total of 25 candidate regions, located along the Java, Japan, Aleutian, Central and South American, Scotia, Lesser Antilles and Cascadia trenches (Fig. 5).
Our results suggest the possibility that the Tohoku-Oki 2011 giant earthquake cycle, and the related distributed set of asperities (Tajima and Kennett, 2012), may be related to strong coupling due to the obliquely subducting Kashima Fracture Zone (Fig. 1), whose offshore location is well mapped, mainly based on ship marine magnetic and seismic reflection data (Nakanishi, 1993;Nakanishi et al., 1992).Its association with subduction earthquake cycles has not been investigated.It is, however, conceivable that a very obliquely subducting, fracture zone ridge with varying elevation along the strike, as is common for large fracture zone ridges (Tucholke and Schouten, 1988;Cande et al., 1995;Croon et al., 2008), causes subduction zone segmentation as observed by Tajima and Kennett (2012).We regard this as more likely than a subducting seamount to be responsible for the Tohoku-Oki 2011 event, as suggested by Duan (2012), given that there is no observational evidence for a subducting seamount close to its epicentre, including no sign of any seamount chain intersecing the subduction zone in the vicinity of the epicentre.The closest seamount chain is subducting about 200 km to the south of the Tohoku-Oki 2011 epicentre (Duan, 2012).The 300 km long asperity recently mapped by Hashimoto et al. (2012) may possibly reflect the subducting Kashima fracture zone.
The Cascadia margin is of particular interest in this context, because there is circumstantial evidence that the main area of slip of the great 1700 Cascadia earthquake was centred on the intersection of the southern Cascadia margin with the Blanco Fracture Zone (Fig. 1).Satake et al. (2003) compare 6 Cascadia rupture models, which they rank based on paleoseismological evidence and comparisons between modelled tsunami heights in Japan.Amongst alternative displacement models for the northern, central and southern Cascadia subduction zone, the southern model, which extends 440 km southward from central Oregon to northern California with an average slip of 21 m, is the most effective in terms of accounting for tsunami run-up in Japan (Satake et al., 2003).Even though this model does not account for the entire pattern of coastal subsidence following this event (Leonard et al., 2010), these results raise the possibility that subduction of a fracture zone ridge may have played a role in triggering the 1700 event.The Blanco Ridge, associated with the large-offset (350 km) Blanco Fracture Zone, displays relief of up to 1 km (Embley and Wilson, 1992).Its extension to the Cascadia margin has been interpreted as continuing in a relatively straight fashion to the southern Oregon margin, where the Blanco slide has been mapped (Goldfinger et al., 2000), and this is confirmed by satellite gravity data (Sandwell and Schubert, 1982).The massive failure of the southern Oregon slope, where the Blanco fracture zone intersects the subduction zone, is expressed by consecutive slide events, and has been interpreted as the result of subduction of a linear ridge, as opposed to subduction of a nearby seamount province, as seamount subduction elsewhere suggests that the typical upper plate expression of a subducted seamount is a relatively narrow deformation trail, unlike the large slope failures of southern Oregon (Goldfinger et al., 2000).This linear ridge would in all likelihood correspond to the Blanco Ridge, in which case the 1700 Cascadia megathrust event would fit our global observations and conceptual model.
There are no FZ-SZ intersections along the Philippines, Marianas, Tonga-Kermadec or any of the Southwest Pacific trenches (Fig. 5).This does not rule out great subduction earthquakes along any of the latter subduction zones, but events in these regions would not be due to an association with a fracture zone intersection.Our analysis shows that data mining techniques primarily developed for extracting useful knowledge hidden in enormous volumes of electronic data on the Internet have great potential for the analysis of large, multidimensional and statistically skewed sets of geological and geophysical data.Recognising the connection between the subduction of fracture zones and great earthquakes globally has the potential to transform long-term earthquake hazard map generation.

Great earthquake rupture propagation and fracture zone ridges
We review rupture propagation associated with FZ-SZ intersections during three well-mapped great earthquake events to evaluate the role of subducting fracture zone ridges in the time-dependence of stress release during great earthquakes.
During the 2004 Sumatra-Andaman earthquake (M w 9.1) the rupture was initiated on the seismically imaged subducting Simeulue Ridge, which is 60 km wide, elevated by 1 km on its western flank, and up to 3 km along its eastern flank (Franke et al., 2008).The ridge is a buried extension of the so-called 96 • Fracture Zone (Kopp et al., 2008) (Fig. 4a) and defines a major segment boundary for the great 2004 Sumatra-Andaman earthquake associated with a step in the slab east of the ridge (Franke et al., 2008).The 2004 event was initiated on the western flank of the subducting Simeulue Ridge (Fig. 4a), and simultaneously climbed eastwards onto the crest of the ridge, while also propagating to the northwest and southwest onto a double fracture zone system called the 94 • and 93 • fracture zones (Kopp et al., 2008), then crossing these and initiating a separate rupture along an unnamed fracture zone west of the 93 • fracture zones (Fig. 4a).The majority of slip distribution during the event is associated with the four fracture zones involved in the rupture, suggesting that Sumatra-Andaman Earthquake (Robinson, 2007) and (b) the 2001 M w = 8.4 Peru earthquake (Robinson, 2007) depicting the modelled rupture process through time for 5 m slip contours, and (c) the 1986 M w = 8.0 Andreanoff islands earthquake (Das and Kostrov, 1990) showing the modelled rupture process through time for a momentmagnitude contour of 3 × 10 20 Nm.Earthquake epicentres are shown as red stars, overlain over vertical gravity gradient maps (Sandwell and Smith, 2009), with fracture zone locations (Matthews et al., 2011) shown as yellow dashed lines, with their interpreted coupling zone extensions outlined in faint yellow dashes.Coloured polygons illustrate the progression of the rupture process relative to the inception at the epicentre, providing a visualisation of the role that the fracture zones played during the event.Inset plots show ship bathymetry (top) and gravity anomalies (bottom) along magenta profiles crossing key fracture zones, with red arrows indicating fracture zone location.The seismically imaged Simeulue ridge (Franke et al., 2008), associated with the subducting 96 • Fracture Zone (Kopp et al., 2008), is outlined as black dashed line.the associated fracture zone ridges represent anomalously coupled regions due to their elevation associated with a shallow slab dip (Supplement Spreadsheet 1).The 2001 Peru earthquake (M w 8.4) nucleated near the top of the Nazca fracture zone lateral ramp, briefly stalling on the Nazca ramp before stepping down the ramp, breaking the adjacent flat segment and finally stopping at the Iquique ridge ramp (Carena, 2011;Robinson et al., 2006) (Fig. 4b).The fracture zone increased plate coupling between the two sides of the fault, resulting in a heterogeneous rupture, with the main stress release focussed on the Nazca FZ-SZ intersection (2006) (Fig. 4b).There are alternative slip models for this event, but we did not find published models that show either slip or energy release as a function of time after the event.
The 1986 Andreanof Islands (Aleutians) earthquake (M w 8.0) rupture history is available as computed moment release through time (Das and Kostrov, 1990).It illustrates that this event was initiated on a subduction segment bounded by the Adak and Amlia fracture zones (Lu and Wyss, 1996), before propagating west and climbing onto a ridge adjacent to the Adak fracture zone (Fig. 4c) where the main moment release was concentrated.This event represents an example where most of the elastic energy was dissipated by the rupture propagating onto the strongly coupled Adak fracture zone ridge, with a peak-trough topographic elevation of over 1000 m (Fig. 8), preventing the rupture from propagating far into the next segment.

Physical model for fracture zone-subduction zone interaction
Global analyses for how fracture zones and other aseismic ridges may relate to the occurrence of large earthquakes were first carried out by Mogi (1969), who suggested that great shallow earthquakes preferentially occur on local depressions.In another global analysis, Kelleher and McCann (1976) concluded that where bathymetric rises or irregularities interact with active trenches the largest shallow earthquakes are generally smaller and less frequent than events along adjoining segments of the plate boundary.These conclusions are contradictory to our results; however, it needs to be kept in mind that global analyses of all subduction zone asperities will likely not reflect a single physical mechanism.
Studies focussed on the role of fracture zones in seismic hazard generation have largely focussed on their regional role in segmenting subduction zones by providing barriers, separating adjacent subduction segments from each other, thus containing ruptures to individual segments (Bilek, 2010).Lu and Wyss (1996) analysed the segmentation of the Aleutian plate boundary and found that the regional segmentation boundaries are controlled by subducting fracture zone, which they interpreted as zones of weakness across which stress may not be transmitted fully.Based on an equivalent study of four segments and their barriers along the Central Chile subduction zone, Metois et al. (2012) concluded that the segments are characterised by higher coupling and separated by narrow areas of lower coupling.If this model held true globally, then we would observe precisely the opposite of what we find, namely great earthquakes should be biased towards subduction segments between fracture zones and the segments centred on fracture zones should be biased against large seismic events.
These three great earthquake rupture histories, together with our global analysis, lead to the following model for why subducting fracture zones are frequently associated with great seismic events.Fracture zones in cross-section form either large steps or represent valleys bounded by topographic ridges; these two morphologies are known as "Pacific-type" and "Atlantic-type" fracture zones (Matthews et al., 2011).Pacific-type large-offset fracture zones display elevation steps in cross-section due to age-offsets at transform faults, where a given fracture zone originates, and these steps are enhanced by flexure due to differential cooling and thermal lithospheric contraction and associated flexural bending, producing flexural ridges (Sandwell and Schubert, 1982) whose internal structure is identical to that of normal ocean crust.In contrast, Atlantic-type small-offset fracture zones are characterised by asymmetric valleys with a prominent ridge on their old side (Collette, 1986), or with ridges on both sides (Matthews et al., 2011).Both types are found in all ocean basins, and mixed types are observed (Matthews et al., 2011).In both categories fracture zone ridges reflect uplifted edges of normal ocean crust, characterised by pronounced relief, lateral continuity and structural integrity, giving rise to strong, persistent subduction interface coupling and stress concentration maintained over long periods of time.
To investigate the nature and elevation of subducting fracture zone ridges we extract ∼ 300 km long bathymetry and gravity profiles orthogonal to subducting fracture zones about 100 km away from a given trench (see Fig. 8 and Appendix E for selected profiles).This distance is required as fracture zone topography is often completely covered with sediments close to the trench (Franke et al., 2008;Robinson, 2007).We use an established method to measure fracture zone topography (Matthews et al., 2011) that captures peaktrough fracture zone elevations, which are found to range from 200 to 1200 m (Fig. 3b), with a mean elevation of about 500 m, and mean gravity anomaly of ∼ 60 mGal (Supplement Spreadsheet 1).The mean age offset of the fracture zones is about 5 million yr.Many fracture zone ridges are simply a function of flexure resulting from differential lithospheric cooling of juxtaposed lithosphere of different ages (Sandwell and Schubert, 1982), with younger oceanic crust being more elevated relative to older crust on the opposite side, and their internal structure is therefore that of normal ocean crust (Bonatti et al., 2005).However, uplifted fracture zone ridges can also be generated via extensional and compressional periods experienced by transform faults (Bonatti et al., 2005;Pockalny et al., 1996).As a consequence we do not find a simple correlation between oceanic age offsets and fracture zone topography (Supplement Spreadsheet 1).After an earthquake nucleates on a given subduction segment the rupture will propagate; upon encountering a nearby fracture zone step or ridge, the rupture is required to step up or down, or to stop.The stress driving the rupture to propagate up, breaking an anomalously coupled zone, will depend on the fracture zone ridge height as well as on its angle relative to the coupling interface.If most of the seismic event's elastic energy has been dissipated, the rupture will be stopped, with a fracture zone step/ridge acting as a barrier, but if not, it will step over the barrier, triggering a large event if the seismogenic zone at the FZ-SZ intersection has been locked for a long period of time.Whether a subducting fracture zone acts as a barrier or leads to the generation of a great earthquake thus depends on the tradeoff between the magnitude of the yield shear stress when a rupture reaches the barrier and the elastic energy carried by the tip of the propagating rupture.This mechanism provides a physical model to account for our observed bias of great earthquakes towards FZ-SZ intersections.Our results and interpretations support the idea of cascading earthquake nucleation from small patches into larger earthquakes (Parsons and Velasco, 2011) under particular circumstances.The physics of such dynamic triggering is clearly complex and likely involves dependencies between subducting fracture zone ridges, faults and ramps (height, length and orientation), surface-wave stresses, pore fluids and aseismic transient slip, such that a larger nucleation area might require a greater amplitude and/or duration of stress change (Parsons and Velasco, 2011).
In contrast, subduction of volcanic chains/ridges and seamounts influences deformation of the overriding plate (Bilek et al., 2003;Scholz and Small, 1997), but seamounts are known to result in weak seismic coupling (Mochizuki et al., 2008;Singh et al., 2011), and this is confirmed by our observations (though we acknowledge we are limited in confidence due to a limited number of identified volcanic chain/ridge intersections).Seismic studies of seamounts have revealed their complex internal layering structure, and their separation from the underlying ocean crust by a decollement surface, both contributing to their disintegration and shearing off during subduction (Das and Watts, 2009), and thus preventing them from initiating strong local coupling at the subduction interface.As a consequence the seismogenic behaviour of subducting seamounts is controlled by the development of an adjacent fracture network during subduction (Wang and Bilek, 2011).The complex structure and heterogeneous stresses of this network provide a favourable condition for aseismic creep and small earthquakes, but an unfavourable condition for the generation and propagation of large ruptures (Wang and Bilek, 2011).Our results support this model overwhelmingly, as subducting volcanic ridges and chains are found to be associated with an excess of relatively small seismic events (magnitude < 8), but not great earthquakes, suggesting that subducting volcanic edifices are not able store significant stress, as stress is released either in small earthquakes, or aseismically.

Conclusions
We have investigated well-defined fracture-zones and volcanic chains/ridges and seamounts that are in the process of being subducted, and their associations with the locations of significant earthquakes.Linking subducting bathymetric features with large earthquakes is a particularly interesting question, since if recurrence-cycles are indeed several hundreds of years in duration, then hazard predictions methods involving the use of existing data may be ineffective at predicting many of the more devastating events.The analysis involved investigating the associations as a function of eathquake magnitude, making use of a data-mining methodology called "Top-N" analysis for coping with the statistically skewed nature of earthquake magnitude distributions.The analysis investigated the associations in the coupling zones between subducting and overriding plates using recent 3-dimensional data sets, and utilised a comprehensive sensitivity analysis of computation of the arbitrary case to thoroughly test the null hypothesis, i.e. that the associations could occur at random.The results in this paper revealed a striking association between fracture-zone/subduction-zone intersections, and the very largest subduction-events on record, with the effect diminishing as magnitudes decrease below moment magnitude 8.0.A similar effect was not demonstrated by volcanic chains/ridges and seamounts, but the conclusions are limited by a small data set size.The analysis went on to demonstrate that the most important fracture zones generally have very large physical offsets, which could play an important role in influencing the subduction coupling.Three modelled rupture events were then discussed, with the spatiotemporal nature of the rupture propogations correlated with the approximate physical locations of the fracture zones.This analysis could have important implications for the field.Figures B1, B2, B3, B4, B5 and B6 present a number of regional maps detailing the tectonic settings where fracturezone/subduction-zone (FZ-SZ) interactions are prominent.These include: the East Indian ocean region along the Java-Sunda trench; the Japanese region; the Aleutian and Alaskan regions considering the geological settings proximal to the Kamchatka and Aleutian subduction-zones, together with a number of prototypical fracture-zone intersection locations associated with large earthquakes; the Central American and South American regions depict numerous target locations along the Andes and Central America.The various plots are superimposed on the Sandwell and Smith vertical-gravity gradient grid (Sandwell and Smith, 2009), as was used in locating and digitising fracture-zones.

Appendix C Constructing a global lithosphere-subduction coupling zone
The interface between the overriding lithosphere and the down-going slab is the coupling zone in which shallow mega-thrust earthquakes occur.In this study the coupling zone is isolated in 3-dimensions, allowing earthquake catalogs to be filtered out.This enables us to compute slab dip in the coupling zone, allowing associations to be assessed comparing flat-slab to steep-slab subduction settings.
In this global study we have made use of the NOAA significant earthquakes catalog to investigate relationships between subducting fracture zone and volcanic ridges/chains with subduction thrust-type earthquakes.It is important to isolate only those earthquakes occurring in the coupling zone at the interface between the overriding and downgoing plates.The approach taken is to use recent global models and data sets to localise these 3-dimensional structures, capitalising on the recent publication of the Slab 1.0 3dimensional subduction model (Hayes et al., 2012), as well as a number of other data sets, described below.Surfaceexpressions of subduction zones form the upper boundary of the coupling zone, defined from the global plate boundary dataset (Bird, 2003).Corresponding three-dimensional models are extracted via the Slab 1.0 3-dimensional global subduction zone model (Hayes et al., 2012).Since the Slab 1.0 model did not cover all regions defined by the global plate boundary data set, missing regions were modelled using the Regionalised Upper Mantle (RUM) seismic model (Gudmundsson and Sambridge, 1998).A few remaining regions were assumed to extend 150 km from their surface expressions, according to global average estimates of slab geometry pertaining to the upper regions of subducting slabs (up to 125 km in depth with 30-50 • subduction dips), as discussed in Lallemand et al. (2005).Intersecting models of the Lithosphere-Asthenosphere Boundary (LAB) with recent 3dimensional representations of subduction zones allows coupling zone regions to be defined.The LAB is defined using a combination of continental and oceanic models.The continental model makes use of the global thermal model TC1 (Artemieva, 2006), in which the 1300 degree isotherm is a standard proxy for the LAB.The oceanic model makes use of the approach in Rychert and Shearer (2009) to estimate thicknesses in regions not represented by the TC1 data set.A nearest-neighbour interpolation is used to make use of the nearest available measurements.Thus constructing the overall coupling zone involved the integration of 5 global data sets.The coupling zone model is referred to as the "CouplingZone1.0"data set.The intersection of the Slab 1.0 model with the LAB model is derived by calculating the 3dimensional intersection between the two models.In Fig. C1 the approach taken to intersect the continental model with the 3-dimensional slab model is illustrated in an example region.As listed in Supplement Spreadsheet 1, 11 fracture zones were found to be associated with 13 of the largest 15 earthquakes in the filtered catalogue.In Fig. 8 bathymetric profiles are shown attributed to these 11 fracture zones.In Fig. E1 the respective gravity anomaly values for the 11 highlighted fracture zones are shown, complementing the results in Fig. 3b of the main text.A good correlation between the two results can be observed, which is expected.
Fig. 2. (a)Evaluation of coupling-zone filtered earthquake associations with intersections between 59 fracture zones and subduction zones, sorted by magnitude N and using a buffer of 150 km on either side of the projected fracture zone.The sorted event magnitudes are plotted against the recall, the number of Top-N earthquakes associated with the fracture zone intersection regions divided by N ; this is compared with the same number (59) of randomly generated subduction zone intersection segments drawn from the global coupling one, repeated 100 times for the entire filtered earthquake data set, shown as median values with the 20th/80th percentile error-bars; (b) same as (a) but for the top 100 (M w ≥ 8.1) significant earthquakes;(c, d) same as (a, b) for volcanic ridge/chain intersections with subduction zones; (e) top-N great earthquake analysis in proximity of subduction zones illustrating the baseline-normalized recall (sensitivity) traded off against the reduction in hazard surface area (% of original hazard area, all subduction earthquakes), called the hazard specificity.Note the sharp inflection point for the top 15 earthquakes at 150 km, illustrating that 87 %, i.e. 13 out of 15 largest earthquakes, all occurred within 150 km of a FZ-SZ intersection.

Fig. 4 .
Fig. 4. Three great earthquake rupture case studies for (a) the 2004 M w = 9.1

Fig. 5 .
Fig. 5. ETOPO1 global relief map(Robinson projection)  with oceanic fracture zones shown as black lines.Global subduction coupling zone (solid blue regions) is overlain with fracture zonesubduction zone intersection regions (solid red bands), significantly more prone to great earthquakes than the blue bands, based on our analysis.Intersections between volcanic chains/ridges (solid yellow bands) are biased towards an above-average occurrence of smaller earthquakes, but this result is more uncertain due to the relatively small number of these intersections.Overlain are the largest 15 instrumentally recorded earthquakes (Mw = 8.6) (white filled stars) and the largest 50 earthquakes (Mw = 8.4) (while filled circles).Subduction zone labels as in Fig.1.

Fig. 6 .
Fig. 6.Illustrating the approach taken for computing arbitrary associations for the Northwest Pacific region.Yellow lines along the coupling zone define a collection of arbitrary spatial regions (at approximately 20 km spacing) that are used to construct experiments involving repeated unbiased sampling at random.The Slab 1.0 3dimensional subduction zone model is shown in blue, coastlines in white, and fracture zones in black.

Fig. 7 .
Fig. 7. Selected region in South America to demonstrate the definition of fracture-zone subduction-zone regions of interest (filled blue regions), constrained by the underlying coupling zone (filled yellow regions).Green contours correspond to depth contours of the Slab 1.0 3-dimensional subduction model, with black lines defining fracture zones.Coastlines are shown for context (grey).The fracturezone intersection regions are shown for a 100 km buffer situation.

Fig. B1 .
Fig. B1.Regional tectonic settings and key data for the East Indian Ocean region.The following colour schemes are used, differing slightly to those in Fig. 1 of the main paper for enhanced background contrast: subduction zones (blue bands), intersection points of fracture zones with subduction zones (yellow squares), intersection points of the volcanic chains and ridges with subduction zones (green squares), largest 25 earthquakes (red stars), earthquakes magnitude 8.3 and above (light blue circles), all other significant earthquakes (small brown circles).The Sandwell and Smith vertical-gravity gradient grid (Sandwell and Smith, 2009) is shown in the background.Conic equidistant projections are used throughout.

Fig. B4 .
Fig. B4.Tectonic setting and key data for the Aleutian and Kamchatka regions.See Fig. B1 for legend.

Fig. B5 .
Fig. B5.Tectonic setting and key data for the Alaskan region.See Fig. B1 for legend.

Fig. B6 .
Fig. B6.Tectonic setting and key data for the Central American region.See Fig. B1 for legend.

Fig. C1 .
Fig.C1.Illustration of the intersection between a lithosphere thickness model of the overriding plate with a 3-dimensional slab model (coloured contours).Green colours reflect the subduction coupling zone, here conservatively taken as the entire intersection region between the overriding and downgoing plates, and red colours reflect the slab surface below the lithosphere-asthenosphere boundary of the overriding plate.The black points define the edge of the coupling zone.
Fig. E1.Gravity anomalies associated with the fracture zones corresponding to those shown in Fig. 3b.

Table 1 .
Calculated areas and fraction of coupling zone for chosen buffer widths around fracture-zone intersections.

Landgrebe: Great earthquakes and oceanic fracture zones
(Kreemer et al., 2003)rouped in five magnitude categories plotted as a function of median slope of the subduction coupling zone(Hayes et al., 2012)versus the median (trench perpendicular) geodetic convergence rates(Kreemer et al., 2003), with error-bars depicting 20th/80th percentiles.(b)Elevation of fracture zone ridges relative to the adjacent ocean floor for the 13 largest earthquakes on record associated with FZ-SZ intersections amongst the top 15 events (see Figs.1 and 2).Bathymetry profiles ∼ 300 km in length were extracted ∼ 100 km seaward of the trench, perpendicular to fracture zones.Fracture zone ridge heights were computed by measuring the total height between ridge crests and fracture zone valleys.