Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief

An ensemble averaging technique is used to remove the long‐wavelength topography and gravity field from subduction zones. >200 residual bathymetric and gravimetric anomalies are interpreted within fore arcs, many of which are attributed to the tectonic structure of the subducting plate. The residual‐gravimetric expression of subducting fracture zones extends >200 km landward of the trench axis. The bathymetric expression of subducting seamounts with height ≥1 km and area ≥500 km2 (N=36), and aseismic ridges (N>10), is largest near the trench (within 70 km) and above shallow subducting slab depths (SLAB1.0 <17 km). Subducting seamounts are similar in wavelength, amplitude, and morphology to unsubducted seamounts. Morphology, spatial distributions, and reduced levels of seismicity are considered inconsistent with mechanical models proposing wholesale decapitation, and the association of subducting seamounts with large‐earthquakes. Subducting aseismic ridges are associated with uplift and steepening of the outer fore arc, a gradual reduction in residual bathymetric expression across the inner fore arc, and a local increase in the width and elevation of the volcanic‐arc/orogen. These contrasting expressions reflect the influence of margin‐normal variations in rigidity on where and how the upper plate deforms, both to accommodate subducting relief and in response to stresses transmitted across the plate interface. The outer fore arc and arc have lower rigidity due to fracturing and thermal weakening, respectively. Similar associations with complex earthquakes and fault creep suggest aseismic ridge subduction may also be accommodated by the development and evolution of a broad fracture network, the geometrical strength of which may exceed the locking strength of a smooth fault.

The main objective of this study is to place new regional constraints on where subducting topographic features are likely present within seismogenic zones, which is important under either seismotectonic scenario for improving understanding of how bathymetric features subduct and evaluating seismic hazard.
Along strike, there are over 50,000 km of subduction fore arcs on Earth [Bird, 2003]. Most seismogenic zones lie underwater making them challenging and expensive regions to study, while Hayes, 1966;Hayes and Ewing, 1970;Watts and Talwani, 1974;Watts and Talwani, 1975] and are similarly present and well resolved by satellite-derived solutions and altimetry Smith and Sandwell, 1997;Sandwell et al., 2014].
The fore-arc region of subduction zones separates the peaks of these paired anomalies and is thus characterized by high topographic (mean 3.5 ) and gravimetric (mean 2 mGal km 21 ) gradients. These gradients can overwhelm lower-amplitude shorter-wavelength signals, thereby masking considerable detail in the structure of the fore arc, arc, and back arc. The objective of this paper is to isolate and remove the average trench-normal profile. This procedure reduces the trench-normal dynamic range, removes the steep gradients associated with the seaward and landward trench-slopes, and significantly increases the resolvability of short-wavelength lower-amplitude structure.
The process by which this is achieved is illustrated in Figure 2. For each subduction margin, the average trench-normal profile is calculated by sampling regional grids of bathymetry and free-air gravity anomaly along trench-perpendicular profiles extending >600 km up and downdip of trench axis and spaced 25 km along strike. The trench axis is defined as the deepest point using the GEBCO08 bathymetric grid [IOC, 2003]. A Fast-Fourier-Transform is used to calculate the frequency spectrum of each profile. The average frequency spectrum is then calculated from the Fourier coefficients at each frequency before an inverse transform is applied to calculate the ensemble average trench-normal profile for each data set and margin. The linearity of the Fourier transform and its inverse means spectral averages should be identical to arithmetic or stacked averages, however, small differences resulting from finite precision of the Fourier transform and associated downsampling are present for small ensembles (< 30 profiles) (supporting information Figure 1). Maintaining the geometry of the trench, grids of each average profile are created and subtracted from the original data sets to produce grids of residual bathymetry and residual free-air gravity anomaly. As opposed to methodologies involving regridding of residuals calculated along profiles, the subtraction of an Ensemble average profiles of the Airy-Heiskanen isostatic gravity anomaly are also calculated for each margin and are dashed in gray in Figure 3b. There is a close agreement between the ensemble free-air and isostatic gravity anomalies and so the long-wavelength component of both anomaly data sets inform us of the nonisostatic, dynamic, contribution to Earth's gravity field in subduction zones.  Figure 2. Figures illustrating the ensemble averaging and grid processing methodology as applied at the Tonga-Kermadec subduction zone. Regional grids of (a-c) bathymetry [IOC, 2003] and (d-f) free-air gravity anomaly [Sandwell et al., 2014] are sampled along trench-normal profiles (a and d). The spectral average is calculated from each ensemble of profiles (b and e inset). Maintaining the geometry of the trench, grids of the average profile are constructed (b and e), and subtracted from the original data sets to reveal (c) residual bathymetry and (f) residual free-air gravity anomaly.

Seismicity
Structural observations from residual gravity and topography anomalies are compared with seismological and geodetic inferences of megathrust slip-behavior. Many of the regions considered in this study ( Figure  1) are characterized by a lack of large magnitude (Mw8) earthquakes and the Global Centroid Moment Tensor (GCMT) [Dziewonski et al., 1981;Ekstr€ om et al., 2012], Engdahl, van der Hilst, and Buland (EHB) [Engdahl et al., 1998], and International Seismological Center (ISC) (http://www.isc.ac.uk/ accessed January 2015) earthquake bulletins are integrated with geodetic inferences of plate boundary locking. A large magnitude earthquake catalog compiled for this study containing the epicenters, seismic asperities, and slip distributions for over 150 of the largest (Mw 7-9.5) subduction zone earthquakes is presented and analyzed in part 2 .
We are interested here in earthquakes occurring on the megathrust and extract events of all magnitudes from the GCMT catalog (period 1976 to July 2013) with thrust faulting mechanisms (T-axis plunge 45 , Paxis plunge 45 ), centroid depth 65 km [Heuret et al., 2011;Hayes et al., 2012], and with strike and dip within 20 and 15 , respectively, of the subducting SLAB1.0 geometry [Hayes et al., 2012]. The full duration ) and magnitude range (mb 3) of the EHB catalog is downloaded and reported event times are used to relocate GCMT centroids to their EHB hypocenters. Unless stated otherwise, earthquake density is calculated as the number of events with depth 65 km within a 10 km search radius using the EHB catalog for the period 1960-2009 and the reviewed ISC catalog (Mw 3.0) for the period 2010-2012.

Global Ensemble Average Profiles
Ensemble average profiles calculated for topography, free-air gravity anomaly, and Airy-Heiskanen isostatic gravity anomaly (dashed) for each margin are displayed in Figure 3. We also display the global ensemble average profile that is calculated from, and shown overlying, the stack of average profiles for each margin segment normalized across distances 300-600 km seaward of the trench axis ( Figure 3) [Hatherton, 1969;Hayes and Ewing, 1970].
Seaward of the trench, ensemble profiles reveal a several hundred kilometer wide belt of free-air gravity anomalies of 150-60 mGal that mark the Outer Gravity High [Watts and Talwani, 1974]. This gravity anomaly is accompanied by a regional rise in topography of several hundred meters. The Outer Gravity High and steep gravity gradient associated with the seaward wall of the trench are almost entirely fit by the computed gravity effect of topography assuming a water/mantle density contrast of 2270 kg m 3 [Watts and Talwani, 1975]; and have been modeled as the flexure of an elastic plate [Hanks, 1971;Watts and Talwani, 1975;Parsons and Molnar, 1976;McAdoo et al., 1978;Bodine and Watts, 1979;Levitt and Sandwell, 1995;Bry and White, 2007]. Flexure is modulated by the effective elastic thickness (T e ) of the subducting lithosphere, which new analysis suggests is strongly linked with its age [Hunter et al., 2013].
The negative free-air gravity anomaly in the vicinity of the trench, a significant part of which cannot be explained by isostatic models, is not always coincident with the trench axis and is sometimes displaced over the landward trench slope due to displacement of crustal materials by low-density sediments. Analogous to the seaward trench slope, the gravity gradient landward of the trench primarily reflects the topography of the fore arc; however, the controls on this gradient are more complex. Elastic plate theory [Davies, 1981], shearing of a viscous accretionary wedge [Emerman and Turcotte, 1983], and cohesive Coulomb wedge theory [Davis et al., 1983;Dahlen, 1984;Dahlen et al., 1984;Dahlen, 1990;Wang and Hu, 2006] have all been proposed to explain the geometry of the landward trench slope. Viscous-flow models show downwarping of the overriding plate as a direct consequence of dynamic coupling to the negative buoyancy force of the subducting slab through the wedge Gurnis, 1992, 1994]. This slab induced downwarping, however, is broad (>400 km), centered landward of the trench approximately above the slab center of mass [Hager, 1984], and is quite different from observed topography. The inclusion of a low-viscosity wedge in dynamic flow models reduces coupling through the wedge Gurnis, 2001, 2003; and inclusion of the plate boundary as a fault in viscous flow models is necessary to reproduce the short-wavelength, high-amplitude near-trench topography [Zhong and Gurnis, 1992;Zhong et al., 1998]. The gravity effect of the subducting slab beneath the fore arc is small, either because the density contrast is small or because the slab is regionally compensated at depth [Watts and Talwani, 1975;Davies, 1981].
Geochemistry, Geophysics, Geosystems 10.1002/2014GC005684 BASSETT AND WATTS CRUSTAL STRUCTURE AND SEISMICITY: 1. SUBDUCTING RELIEF Landward of the fore arc, the ensemble profile shows a broad free-air gravity high that decreases by 50 mGal across the back arc ( Figure 3b). This shallow gradient is flanked in some, predominantly oceanicoceanic, margins by a narrower (100-200 km), large amplitude (>100 mGal) high coincident with the island arc. These gravity anomalies reflect the combined effects of the topography and compensation of the island arc, and the dynamics of the back arc.   show the average Airy-Heiskanen isostatic gravity anomaly. For each data set, average profiles of each margin are normalized 300-600 km seaward of the trench axis and stacked at the bottom of each plot. The global average of these profiles is shown in black. FA 5 Fore arc. Red triangle with bars marks the location of most volcanic arcs [Syracuse and Abers, 2006].

10.1002/2014GC005684
The close correspondence between the ensemble average profiles calculated from isostatic and free-air gravity anomalies (Figure 3b) suggests that the residual gravity-anomaly does indeed represent the non-isostatic contribution to the trench-normal gravity signal. At short-wavelengths, the residual probably reflects the fore-arc, arc, and back-arc structure while at long-wavelengths it probably reflects dynamic processes such as those associated with sinking slabs and flow in the overlying mantle wedge.

Interpretations of Subducting Relief
The seafloor is far from uniform and comprises a rich mosaic of volcanic and tectonic structures including seamounts [Hillier and Watts, 2007], large igneous provinces [Coffin and Eldholm, 1994], fracture zones, and aseismic, active and extinct ocean ridges [Matthews et al., 2011]. When carried into subduction zones, these structures have been shown to modulate seismogenic processes [Kelleher and McCann, 1976;Bilek et al., 2003;Robinson et al., 2006;Mochizuki et al., 2008;Das and Watts, 2009;Sparkes et al., 2010] and the morphology, subsidence, and uplift history of the fore arc [Lallemand and Le Pichon, 1987;Lallemand et al., 1989;Dominguez et al., 1998;Dominguez et al., 2000;Oakley et al., 2008]. The morphological imprint on the fore arc is well expressed in grids of residual bathymetry and residual free-air gravity anomaly. In this section, these observations are presented alongside the relations shown with seismicity.

The Louisville Ridge Seamount Chain: Tonga Trench Collision Zone
At 26 S, the Tonga trench is intersected by the Louisville Ridge seamount chain [Hayes and Ewing, 1971;Lonsdale, 1986;Watts et al., 1988;Ballance et al., 1989;Koppers et al., 2004Koppers et al., , 2012. West of Osbourn, which is the oldest unsubducted seamount, the residual bathymetric grid shown in Figure 4a reveals a 70-80 km wide belt of positive bathymetric anomalies >1 km in amplitude extending 60 km across the fore arc.
Occupying the trench, the largest anomaly is associated with Mo'unga seamount, which was previously recognized as a conical, 1.5 km high, and 10 km wide uplift from swath bathymetry data [Lonsdale, 1986;Pontoise et al., 1986;Ballance et al., 1989;Stratford et al., 2015]. The residual bathymetric grid suggests a NE-SW elliptical geometry for Mo'unga with amplitude >2.5 km and semi-major and semiminor axes of 50 and 25 km, respectively. A second smaller residual bathymetric high (amplitude >1.5 km) is observed NW of Mo'unga. Both highs are located up to 50 km south of the projected track of the unsubducted Louisville Ridge, along which bathymetric anomalies are typically 1 km in amplitude.
The Louisville collision zone has been the subject of two recent marine-geophysical research cruises [Grevemeyer and Flueh, 2008;Peirce and Watts, 2011] and seismic reflection and refraction data calibrate the association of residual bathymetric anomalies with subducting relief. Figure 4c shows two preliminary forward P-wave velocity (V P ) models from Bassett [2014].
Section A-A' is dip parallel, traverses both Osbourn and Mo'unga, and shows both seamounts to be associated with a 2-3 km increase in thickness of a layer overlying the subducting Pacific plate crust with V P 4.2-5.3 km s 21 . A third increase in layer thickness is modeled at greater depth, beneath the smaller residual bathymetric high, but this model region is not well constrained.
Section G-G' is strike-parallel and traverses this smaller high, which is correlated with peaks in free-air gravity (160 mGal) and magnetic (1280 nT) anomalies. This profile reveals a sharp and local reduction in forearc velocity gradient, with seismic wave speeds between 2 and 12 km depth (b.s.f.) >1.0 km s 21 slower than those observed 20 km north and south of the bathymetric high at equivalent depth. A low-velocity body is modeled above the Pacific plate crust from wide-angle reflections, but wave speeds near the slabfore-arc-Moho intersection are not well constrained.
Despite alignment with the unsubducted Louisville Ridge, preliminary analysis of refraction data along Profile C resolves no subducting seamounts beneath the fore arc [Bassett, 2014]. This velocity-structure is consistent with the absence of both magnetic and residual bathymetric anomalies. A low-velocity Geochemistry, Geophysics, Geosystems   [Hayes et al., 2012]. The trench axis is dashed in gray. Black dash north of the contemporary collision zone marks the possible extent of seamount subduction in the last 1 Myr [Ballance et al., 1989]. Note the absence of large bathymetric anomalies within this section of the fore arc. Gray profiles show active source MCS/OBS profiles with the sections displayed in Figure 4c marked in black. Red dots along these profiles mark the subducting slab-fore-arc-Moho intersection. Historically active volcanoes are shown as red triangles [Simkin and Siebert, 2002]. Figure 4b EHB/ISC earthquake density and hypocenters for GCMT thrust earthquakes. Contours at 200 m increments show the geometry of the Louisville flexural moat assuming residual bathymetric anomalies in the outer fore arc and the geochemically anomalous arc volcano V [Timm et al., 2013] constrain the geometry of the subducting portion of the ridge [Bassett, 2014]. Note the correlation between the predicted boundaries of the flexural moat (black dash) with the northern and southern limits of the Louisville seismic gap (white dash). Figure 4c Preliminary forward P-wave (V P ) velocity models of crustal structure along profile segments A-A' and G-G' [Bassett, 2014]. The velocity model along Profile A-A' is similar to that published by Stratford et al. [2015] and differences are minimal in well-constrained model regions.
(V P 2.8-5.4 km s 21 ) layer overlying the Pacific Plate crust to 12 km depth is interpreted as subducting volcaniclastic and pelagic sediments, which accumulate and reach up to 1.5 km thickness in the Louisville Ridge's flexural moat [Contreras-Reyes et al., 2010].
These wide-angle seismic models and magnetic anomaly modeling in Bassett [2014] support the interpretation of residual bathymetric anomalies as reflecting the presence of subducting bathymetric relief. The subducting slab-fore-arc-Moho intersection is marked by red dots in Figure 4a, and this intersection may spatially limit the expression of subducting bathymetric features in fore-arc morphology.
Morphological characteristics of the Louisville collision zone are accompanied by a sharp 200 km wide reduction in the number of earthquakes known as the ''Louisville Gap'' [Habermann et al., 1986;Scholz and Small, 1997]. Earthquake density plotted in Figure 4b shows this gap centered on large residual bathymetric anomalies in the outer fore arc, landward which an isolated cluster of thrust earthquakes give the gap a two-lobbed appearance. This cluster may be analogous to earthquakes observed in-front of subducting seamounts in Japan [Mochizuki et al., 2008] and Cascadia [Tr ehu et al., 2012]. The absence of geodetic observations and the short-duration of instrumental earthquake catalogs make it hard to determine if this gap represents a creeping megathrust, or an interseismically locked seismogenic zone with each lobe capable of generating a future Mw 8.2 earthquake.
The distribution of magnetic and residual bathymetric anomalies may suggest an approximately East-West geometry for the subducting segment of the Louisville Ridge [Timm et al., 2013;Bassett, 2014]. This geometry is supported by Pb and Sr isotope ratios requiring up to 40% input of Louisville Ridge material into partial melts beneath the arc volcano ''V'' (Labeled in Figures 4a and 4b); and may also provide one possible physical explanation for the Louisville seismic gap.
The contours in Figure 4b show the geometry of the Louisville flexural moat calculated in 3-D by extending the mean cross-sectional structure of the ridge along the westerly geometry implied by the distribution of residual bathymetric anomalies and using an effective elastic thickness (T e ) of 10 km [Contreras-Reyes et al., 2010]. Correlating each lobe of the seismic gap with each flank of the flexural moat, Bassett [2014] suggest the locally higher flux of sediment and entrained fluid may reduce the shear strength of the subduction interface. But thicker sediment volumes have also been linked with large earthquakes [Ruff, 1989;Scholl et al., 2011] and although these scenarios remain hard to distinguish between, seismic-reflection, and morphological observations showing a significantly greater degree of active and pervasive upper plate extensional deformation within the seismic gap lead us to prefer interpretations where the underlying megathrust is weak.

Mariana Arc
The Mariana arc is commonly regarded as the type example of an erosive, extensional, and aseismic subduction margin [Uyeda and Kanamori, 1979;Hussong and Uyeda, 1982]. Residual bathymetry is shown in Figure 5b. Three large (diameter >50 km, amplitude 2 km) bathymetric anomalies occupy the central Mariana fore arc and continuity with bathymetric anomalies seaward of the trench axis strongly suggests an association with subducting seamounts. Along the southern Izu-Bonin trench, the subducting Ogasawara Plateau is expressed as a 2-3 km uplift of the outer fore arc, the amplitude and along-strike width (180 km) of which is well correlated with plateau dimensions seaward of the trench axis. All observations of subducting relief appear restricted to slab depths 20 km [Hayes et al., 2012], which is similar to 15-20 km crustal thickness estimates from seismic refraction data across the southern Mariana fore arc [Takahashi et al., 2007]. Figure 5c shows that earthquake density is largest along the southern Mariana megathrust and a 600 km long segment has ruptured in >70 intermediate magnitude (5Mw<7) thrust events. No positive residual bathymetric anomalies are observed within this seismic band, the northern termination of which is coincident with a large anomaly likely associated with a subducting seamount. Megathrust seismicity is sparse along the central and northern Mariana arc and all large residual bathymetric anomalies are characterized by low-earthquake density. The amplitudes, form, and relations with seismicity of these anomalies support earlier suggestions that along the Mariana arc, bathymetric relief onboard the Pacific plate may be subducting intact [Fryer and Smoot, 1985] and aseismically [Wang and Bilek, 2014].

10.1002/2014GC005684
One unique characteristic of the Mariana arc is the presence of fore-arc seamounts associated with serpentinite mud diapirism Fryer and Smoot, 1985;Fryer and Fryer, 1987;Fryer et al., 1999;Fryer, 2012]. These seamounts typically occur 30-100 km west of the Mariana trench axis [Fryer and Fryer, 1987], have little or no magnetic signature [Hussong, 1982], and are associated with the upwelling of serpentinite muds formed in the mantle wedge as it is hydrated by fluids from the subducting plate [Fryer, 2012]. Locations where serpentinite seamounts have been sampled since 1981 are shown as red pluses in Figure 5c (inset) and we note that almost all are near large positive residual bathymetric anomalies. This correlation supports emplacement mechanisms whereby vertical tectonic movements and fracturing of the overthrusting plate required to accommodate subducting oceanic seamounts provides faults and planes of weakness along which diapirs of serpentinized ultramafics can rise Fryer and Smoot, 1985].

Java
Offshore Java (Figure 6), the major bathymetric relief is associated with the Roo Rise oceanic plateau. The plateau has crustal thickness of 12-18 km [Shulgin et al., 2011], a sharp eastern margin, and a tiered western margin. The inner boundary defines the western margin of the shallower central segment of the plateau, which can be tracked across trench axis and associated with a 2 km reduction in fore-arc water depth. The location of highs and lows within this elevated segment are contiguous with plateau physiography seaward of the trench axis. Northwest of the Roo Rise, a second residual bathymetric anomaly with amplitude >2 km in the outer fore arc is likely associated with an isolated subducting seamount.
The bathymetric anomalies we associate with subducting relief are limited by the 20 km depth contour of the subducting slab [Hayes et al., 2012]. This contour is correlated with the location and geometry of the outer arc high, and the subducting slab-fore-arc-Moho intersection as constrained by wide-angle seismicrefraction profiles [Kopp et al., 2002;Planert et al., 2010;Shulgin et al., 2011] (Figure 6a). Thrust earthquakes in Java are sparse, of highly variable mechanism and typically low magnitude. Only two earthquakes have Mw7, but these events are significant for being both tsunami earthquakes and widely cited examples of interactions between earthquakes and subducting seamounts [Masson et al., 1990;Abercrombie et al., 2001;Bilek and Engdahl, 2007]. Epicenters and distributed slip models for the 1994 (Mw 7.8) and 2006 (Mw 7.7) tsunami earthquakes from Bilek and Engdahl [2007] are shown in Figure 6c.
The 2006 earthquake propagated east with slip distributed between the trench and mantle-wedge corner. Residual bathymetry within the rupture area is typically low (<1 km), although a positive anomaly near the trench is correlated with a small increase in coseismic-slip ( Figure 6c). The close proximity of the eastern rupture extent to a large residual bathymetric anomaly in the outer fore arc lead us to agree that the 2006 earthquake did not rupture or initiate near a subducting seamount, but was quite possibly arrested by one [Bilek and Engdahl, 2007].
The 1994 earthquake initiated on the southwestern flank of the largest residual bathymetric anomaly associated with the Roo Rise ( Figure 6c). This earthquake propagated 120 km east and west of the epicenter, extending between the trench axis and 40 km depth. A detailed local correlation between coseismic slip and residual bathymetry is not observed, but we do interpret a regional correlation between the area of large coseismic slip and the western high-standing block of the plateau. Coseismic slip diminishes between the tiered western plateau boundaries and we interpret a second correlation between the western extent of coseismic slip and the outer boundary of the Roo Rise. We agree with existing links between the 1994 earthquake and subducting relief [Abercrombie et al., 2001], but make this link with oceanic plateau rather than seamount subduction. This distinction may be significant given their contrasting modes of isostatic compensation and the lower height to length/width ratios of oceanic plateaus, which may influence megathrust roughness.

Japan
In 2011 offshore Honshu, Japan, the Mw 9.0 Tohoku-Oki earthquake ruptured 400 km of strike [Ide et al., 2011;Koketsu et al., 2011;Lay et al., 2011;Simons et al., 2011]. The Pacific plate seafloor adjacent to this rupture area is smooth [Wang and Bilek, 2014], but further south and of significance to this study, subducting seamounts of the Joban seamount chain have been interpreted near the trench axis [Lallemand et al., 1989;Nishizawa et al., 2009] and >5 km beneath the fore arc [Mochizuki et al., 2008]. Residual bathymetry for the Japan trench is shown in Figure 7. In the north, Erimo seamount is observed seaward of the trench axis and a second residual bathymetric high (1.5 km) of similar dimensions has been associated with a subducting seamount [Lallemand and Le Pichon, 1987]. In the southern Japan trench, residual bathymetric highs with 2 km amplitude are aligned with unsubducted Joban seamounts and can be tracked >100 km landward of the trench axis. We tentatively interpret a third subducting seamount beneath the Ty oshi spur (TS), part of a regional positive anomaly extending from Inubo Cape (IC). The shortwavelength character of these anomalies is limited by the 30 km slab contour [Hayes et al., 2012] and near fore-arc crustal thickness estimates (25 km) [Fujie et al., 2013].
Most thrust earthquakes along the Japan megathrust are distributed between 10 and 50 km depth. In the EHB/ISC catalog, regions of low recent earthquake density are located: near the trench in northern Honshu where the outer fore arc is characterized by positive residual bathymetry; north of this region within a 70 km wide trench-normal bathymetric trough bounded by the Hokkaido continental slope break; and in regions of active seamount subduction. The seamount we interpret offshore Hokkaido is located in a region where the outer fore arc has regionally low-earthquake density. The subducting Joban seamounts have similarly low-earthquake density, but appear so adjacent to, and define the southward boundary of, a recently active segment of the megathrust.
In Figure 7c, we show the distribution of seismic coupling on the Pacific plate megathrust prior to the Tohoku-Oki earthquake, coseismic slip contours from inversion of horizontal displacements [Loveless and Meade, 2010], and the epicenters of aftershocks occurring within 2 days (gray dots). The region of Joban seamount subduction was characterized by low coupling before the earthquake, experienced no coseismic slip, and almost all aftershocks occurred on the megathrust to the north. These observations strongly support the suggestion of Wang and Bilek [2014] that Joban seamounts subduct via aseismic creep and may have acted as a barrier to the southward propagation of the Tohoku-Oki earthquake. The absence of large amplitude residual bathymetric anomalies in the center or northern region of the rupture zone make it unlikely subducting seamounts were associated with either the nucleation or northern termination of earthquake rupture [Zhao et al., 2011].

Kamchatka
Residual bathymetry for Kamchatka and the western Aleutians is shown in Figure 8. This grid is generated by calculating the average trench-normal profile from Kamchatka, but extending the removal of this profile further east and west along the western Aleutian and Kuril trenches, respectively. The subducting Emperor seamount chain can be identified as a 1-2 km high bathymetric anomaly 200 km landward of the trench axis and beneath the Kronotsky peninsula (KP in Figure 8a). The landward extent of this anomaly is coincident with the 40 km slab contour (Figure 8a), and is close to crustal thickness estimates (30-40 km) along the east coast of Kamchatka from gravity measurements, active source seismic [Bogdanov and Khain, 2000] and receiver function modeling [Levin et al., 2002].
In the 20th century, the Kamchatka megathrust has ruptured in six large (Mw 7), three great (Mw 8), and one giant (1952 Mw 9.0) earthquake (Figure 8c) [B€ urgmann et al., 2005]. The region adjacent to the subducting Emperor seamount chain ruptured most recently in the 1997 Mw 7.8 Kronotsky earthquake, which propagated entirely southwest from an epicenter near the subducting ridge flank. Aftershock distributions for older events are uncertain, but these also suggest rupture extents restricted to one side of the subducting ridge, which appears a barrier to rupture propagation. Inversion of horizontal GPS velocities show the region of the subducting Emperor seamount chain to be strongly coupled [B€ urgmann et al., 2005].
The black line in Figure 8b shows the small circle best fitting the locations of Kamchatka arc volcanoes [England et al., 2004]. Where the Emperor/Meiji seamounts subduct, arc volcanoes are located 50 km further landward and have residual bathymetry (4 km) at least 1 km higher than adjacent arc segments. Geochemical data [Keller et al., 2000;Cottrell and Tarduno, 2003], the change in physiography of the seamount chain, and the reduction in width of the flexural moat suggest the Hawaiian hotspot may have been close to a spreading ridge during the formation of Detroit seamount. Subduction of the Emperor seamount chain may thus be more analogous to aseismic ridge subduction, other examples of which are discussed below.
Geochemistry, Geophysics, Geosystems 10.1002/2014GC005684 6. Aseismic Ridges 6.1. South-America The Nazca plate subducts beneath South America at rates increasing from 54 mm yr 21 beneath Ecuador and Columbia to 80 mm yr 21 beneath southern Chile [DeMets et al., 2010]. The active Chile Rise intersects the margin near 45 S and to the north, many of the fracture zones, seamounts, and aseismic ridges on the Nazca plate have been related to persistent seismogenic segment boundaries [Bilek, 2010;Sparkes et al., 2010] or the dynamics of individual earthquakes [Barrientos and Ward, 1990;Robinson et al., 2006]. The two largest bathymetric reliefs are the Nazca and Carnegie aseismic ridges (Figure 9).
Seaward of the trench, the Carnegie Ridge has residual bathymetry up to 3 km in height, crustal thickness of 14 km , and a width increasing toward the trench axis from 100 km to 300 km. Over an identical width, the adjacent fore arc has residual bathymetry of similar amplitude (2-3 km). The surficial expression of the subducting ridge is strongest across the outer fore arc and diminishes entirely within 90 km of the trench axis. There is no residual bathymetric expression of the Carnegie Ridge across the inner fore arc, but further landward and analogous to observations made in Kamchatka (Figure 8b), we note that both the width (100 km) and elevation (residual bathymetry >2km) of the volcanic arc increases by a factor >2 relative to adjacent regions.
The Grijalva fracture zone defines the southern boundary of the shallower northern segment of the Nazca plate, corresponding to a more seismically active segment of the megathrust. Where the Carnegie ridge subducts, thrust earthquakes with Mw 6 are located landward of residual bathymetric anomalies and have complex source-time functions [Bilek, 2010]. The southward propagation of large historical earthquakes in 1906 (Mw 8.8) and 1942 (Mw 7.8) was limited by the northern ridge flank [Kelleher, 1972;Collot et al., 2004].
EHB earthquake density is regionally low along the Peru, Ecuador, and Columbia segments of the convergent plate boundary. The ISC catalog, in contrast, shows higher densities where the Carnegie ridge is subducted, particularly in the outer fore arc near subducting ridge flanks, and near the volcanic arc (Figure 9b). The different magnitude ranges over which the EHB (mb 3) and ISC (mb 2) bulletins are constructed may suggest this increase in density is from small earthquakes.
The Nazca aseismic ridge is 200 km wide, trench-normal, and has a residual bathymetry high of up to 1.5 km amplitude. Across the fore arc, the surficial expression of the subducting ridge has amplitude up to 2 km, is strongly limited to the outer fore arc and has a width well correlated with ridge physiography seaward of the trench. The volcanic arc terminates east of this region and although the Peruvian Andes have residual bathymetry up to 1 km higher than the adjacent region further north, it is not clear if this additional elevation is associated with ridge subduction.  Figure  9f). Low slip is modeled southeast of this zone, is not required to fit InSAR and GPS observations [Chlieh et al., 2011] and at depths <50 km the north-western Nazca ridge flank acted as a barrier to rupture propagation.
The southern flank of the Nazca ridge ruptured in 1942 (Mw 8.2) and 1996 (Mw 7.7). The 1942 event is suggested to be complex, with southward rupture propagation and significant moment release away from the Nazca ridge [Swenson and Beck, 1996]. The 1996 earthquake had an epicenter on the southern ridge flank and also propagated south rebreaking the deeper portion of the 1942 rupture area [Chlieh et al., 2011]. Neither earthquake propagated north of the ridge crest and the south-eastern ridge flank may also act to limit rupture propagation.
Subduction of the Nazca ridge is similarly correlated with an increase in ISC earthquake density (Figure 9e). The highest densities are near the flanks of the subducting ridge and are separated by a region of lower density near the ridge-crest.
Bilek [2010] associate the 2007 earthquake with a gap left since the 1687 earthquake, but note that the 2 m of coseismic slip is an order of magnitude less than the 20 m of relative plate motion accrued. From interseismic GPS observations, Chlieh et al. [2011] show the megathrust immediately overlying the subducting Nazca ridge is weakly coupled (creeping) and up to 60% of relative plate motion may be accommodated aseismically (Figure 9g).
Intermediate magnitude thrust events in both of these collision zones are more complex, have longer durations [Bilek, 2010] and are consistent with observations of rupture complexity in large events [Swenson and Beck, 1996;Bilek, 2010]. These observations are interpreted as reflecting the structural complexity and lower rigidity of the megathrust, and correlation with elevated earthquake densities and geodetic observations of low coupling may suggest that upper plate deformation and small earthquakes accommodate a significant proportion of relative plate motion where aseismic ridges are subducted.

Central-America
The subducting plate offshore Costa Rica is extremely rough hosting the Cocos Ridge, Panama fracture zone, and >10 ridge flanking seamounts. Residual bathymetry in this region is shown in Figure 10.
Seaward of the trench axis, the Cocos Ridge is 200 km wide, 2 km high, and has thicknesses up to 20 km [Walther, 2003]. Residual bathymetric anomalies associated with this ridge, and a margin normal ridge and trough associated with the Panama fracture zone, are strongly limited to the outer fore arc within 40-60 km of the trench axis.
Landward of the Cocos Ridge, the Cordillera de Talamanca (CDT) has residual bathymetry >2 km, is a factor >2 wider than the adjacent Cordillera de Guanacste (CDG) and Serrania de Tabasara (SDT), and is the highest mountain range in Central America (Figure 10a). This observation is analogous to the wider and higher volcanic arcs landward of the Emperor/Meiji and Carnegie ridges.

10.1002/2014GC005684
Most thrust earthquakes are distributed across the inner fore arc and landward of large residual bathymetric anomalies (Figure 10b). Exceptions are a cluster of events aligned with the Panama fracture zone; and a second smaller cluster located west of a conical 2 km residual bathymetric high (circled in Figure 10) that is likely associated with a subducting seamount. A Mw 7.0 event occurred <50 km downdip of this high in 1990 [Bilek et al., 2003] and it is possible this seamount, or a deeper subducting geometrical irregularity, contributed to the termination of rupture in the 1950 Mw 7.7 earthquake near the eastern Nicoya Peninsula ( Figure 10c) [Protti et al., 1995].
Historical earthquakes landward of the Cocos Ridge occurred in 1941 (M 7.7) and east of the Panama fracture zone in 1934 (M 7.7) [Kelleher et al., 1973]. Source time functions for recent events in 1983 (Mw 7.4), 1990 (Mw 7.0), and 1999 (Mw 6.9) are complex, consisting of multiple subevents, often at different depths and/or with different focal mechanisms [Bilek et al., 2003]. This complexity is in stark contrast to the simple source time function associated with the 1992 (Mw 7.7) tsunami earthquake offshore Nicaragua [Kanamori and Kikuchi, 1993], which Wang and Bilek [2011] interpret as reflecting the roughness of the subducting plate and the structural complexity of the megathrust. The distribution of residual bathymetric anomalies supports this interpretation.

Fracture Zones: Southern Chile
Fracture zones are aseismic, off-axis, traces of the transform faults that segment mid-ocean ridge systems. The sediment cover of ocean basins is often sufficient to bathymetrically obscure fracture zones; however,  [Robb and Kane, 1975;Fox et al., 1976;Prince and Forsyth, 1988;Matthews et al., 2011].
In Southern Chile, the active Chile Rise is subducted near 45 S and a dense network of ENE striking fracture zones in the adjacent Nazca and Antarctic plates are clearly expressed in grids of residual free-air gravity anomalies ( Figure 11). The largest fracture zones (e.g., Mocha, Valdivia, Guafo, Darwin) are correlated with up to 50 km steps in the location of the trench axis and can be interpreted up to 300 km further landward (Figure 11b).
The obliquity of most fracture zones with respect to margin-strike decreases across the trench axis in response to the dip of the subducting slab. Across the fore arc, they are correlated with boundaries between structural segments characterized by positive or negative/near-zero free-air gravity anomalies and cusps or transitions in the physiography of the coastline. Figure 11c shows the distribution of slip in the 1960 Mw 9.5 Chilean earthquake [Barrientos and Ward, 1990]. This giant event ruptured the megathrust between the Mocha and Guamblin fracture zones and the northern and southern rupture limits are coincident with shallower fore-arc segments, the latter reflecting subduction of the Chile Rise (SCR). The distributed slip model shows three main asperities and Barrientos and Ward [1990] noticed that the boundaries between these asperities, and the northern and southern rupture limits, were close to the projected locations of fracture zones. Our constraints on fracture zone geometries and observations of margin segmentation support the suggestion of Barrientos and Ward [1990] that fracture zones may be an important source of margin segmentation along the southern Chile fore arc.
The influence of fracture zones on the segmentation of subduction boundaries is further demonstrated by observations in the Aleutian and southern Kuril Islands presented in part 2 of this study.

Discussion
It has long been suggested that the tectonic structure of the subducting plate exerts a strong influence on interseismic and coseismic deformation in subduction zones [Kelleher and McCann, 1976;Spence, 1977], but the physical and mechanical processes by which tectonic features subduct and the seismological consequences of this interaction remain a subject of debate and contrasting mechanical models have been proposed. One model suggests seamounts are accreted, either by decapitation at some level or by shearing along their base [Cloos, 1992;Cloos and Shreve, 1996], and act as seismic asperities generating large earthquakes within the seismogenic zone [Cloos, 1992;Scholz and Small, 1997]. An alternate model, proposed by Wang andBilek [2011, 2014], posits that seamounts subduct more or less intact, with the development of a dense fracture network near and above the subduction interface resulting in fault conditions favorable for many small earthquakes and creep, and unfavorable for the nucleation of large earthquakes. It is not possible in our grids to discriminate between actively subducting and accreted relief, but the morphology and spatial distribution of bathymetric anomalies we interpret within fore arcs, and the relationships they show with recent seismicity, may provide constraints discriminating between these models.

Distribution and Character of Subducting Seamounts
From residual bathymetric grids for all subduction zones on Earth, >200 bathymetric anomalies are identified within fore arcs that may be associated with subducting relief. Isolated and irregular bathymetric anomalies are interpreted as subducting seamounts, most of which are similar in wavelength, amplitude, and morphology to seamounts located seaward of the trench axis.
The area and height (h) of each bathymetric anomaly is estimated by determining the deepest closed residual bathymetric contour circumscribing the anomaly, before calculating the area and depth difference between the basal contour and anomaly peak. Most anomalies have height <1 km (Figure 12b) and their spatial distribution is colored for height and shown in Figure 12a. Anomalies with h1 km and area 500 km 2 are shown as red dots (N536) and, coupled with aseismic ridges (yellow pluses N>10), represent our most confident interpretations of subducting relief. Anomalies with 0.5h>1 km (N596) and h<0.5 (N572) are shown as blue and green dots, respectively and, although more speculative than larger anomalies, in many regions (e.g., Ryukyu, Mexico) the lower amplitude is consistent with seamount heights seaward of the trench axis. A strong correlation is observed between locations of residual bathymetric anomalies and the bathymetric character of the adjacent subducting plate seafloor, which supports Geochemistry, Geophysics, Geosystems 10.1002/2014GC005684 assumptions made in existing regional studies of the influence of bathymetric relief on seismogenic behavior. A catalogue containing the location, height, area, and volume of each bathymetric anomaly is provided in supporting information.
The histograms in Figure 12c show the distribution of residual bathymetric anomaly peaks with respect to slab-depth and trench distance. Most anomaly peaks (95%) are located within 70 km of the trench axis and above regions where the subducting SLAB1.0 depth is less than 17 km [Hayes et al., 2012]. The underlying density plot reflects the full area occupied by all residual bathymetric anomalies (excluding aseismic ridges) and reveals a similar restriction of most anomalies to the shallow and near trench region of subduction zones. In Japan and Kamchatka, bathymetric anomalies in line with the Joban and Emperor seamount chains are identified >100 km landward of the trench axis, associated with SLAB1.0 depths 30 km. These large lateral extents may reflect the shallow trajectory of the subducting Pacific Plate and/or fore-arc crustal thickness.
Density plots showing the distribution of slab-depths and trench distances associated with residual bathymetric anomalies for each active margin segment are shown in supporting information. Where known, the depth of the fore-arc Moho is shown as a dashed red line, which in Tonga, Java, and Mariana appears to bound the landward limit of most residual bathymetric anomalies. These observations are consistent with the lack of surface uplift observed above a possible deep (30-40 km) subducting seamount imaged below the fore-arc Moho in Sumatra [Singh et al., 2011]. The isolated and irregular nature of subducting seamounts makes it difficult to reliably constrain the dip-parallel extent of subducting relief. This difficulty is significantly reduced for the case of aseismic ridge subduction.

Character of Aseismic Ridge Subduction
Subducting aseismic ridges appear subject to similar observation limits as seamounts and large residual bathymetric anomalies are strongly limited to the outer fore arc. In contrast to seamounts, however, where aseismic ridges subduct it is simpler to assume the location, width, and height of the subducting Nankai -SW Japan Peru Figure 13. Aseismic ridge profiles. Profiles showing topography (gray), ensemble average topography (black), and residual bathymetry (red 5 gray-black) across five regions of aseismic ridge subduction. The subducting SLAB1.0 geometry is shown in black (vertical scale compressed) [Hayes et al., 2012]. Vertical black and gray (with arrow) lines mark the interpreted position of the trench-slope break for the ensemble average and residual bathymetric profiles, respectively, and illustrate the reduction in width of the outer-fore arc.
Geochemistry, Geophysics, Geosystems 10.1002/2014GC005684 bathymetric relief from ridge physiography seaward of the trench axis; and in most cases, the assumption that subducting ridges extend beyond the outer fore arc is likely to be similarly valid. The residual bathymetric expression of aseismic ridges can thus be interpreted as constraining the spatial limits over which subducting bathymetric anomalies directly influence fore-arc morphology, enabling the physical controls on these limits to be considered. Figure 13 shows profiles of bathymetry (grey solid line) and residual bathymetry (filled red line) traversing five subducting aseismic ridges. The ensemble-average profile, calculated across the 250 km of trench either side of the ridge, is shown in black and reflects the mean structure of the trench slope in the absence of subducting relief. The subducting SLAB1.0 geometry is also shown in black (note different depth scale) [Hayes et al., 2012]. These profiles reveal the contrasting architecture of the fore arc in regions of aseismic ridge subduction.
In general, the amplitude of residual bathymetric anomalies is maintained in the trench axis. The geometry of the seaward trench slope is parallel to the ensemble-average profile and the outer fore arc is similarly elevated. The trench-slope break, however, is located 45-70 km (factor >2) nearer the trench. Relative to the ensemble-average profile, the outer fore arc is steeper and the inner fore arc is shallower/near-horizontal resulting in a gradual dip-parallel reduction in the bathymetric expression of the ridge. We thus observe a strong correlation between the location of the trench-slope break, which occurs 10-50 km landward of the trench axis and above SLAB1.0 depths of 5-15 km, and the landward limit of large residual bathymetric anomalies. The amplitude reduction in ridge expression across the inner fore arc occurs over 50-100 km and residual bathymetric anomalies converge with ensemble average profiles, indicating there is no expression of the subducting ridge, 80-120 km from the trench axis.
In Costa-Rica, Ecuador, and Kamchatka, the active volcanic-arc landward of subducting ridges is elevated relative to adjacent regions by 1.5-2.0 km. The top of the slab beneath the arc/highlands is typically >60 km deep and near-arc uplift is juxtaposed against the notable absence of residual bathymetric anomalies associated with subducting relief across the inner fore arc.

Implications for Mechanical Models of Bathymetric Relief Subduction 8.3.1. Interpretation of Aseismic Ridge Profiles
The along-strike extent of outer fore-arc and near-arc uplift indicates both are clearly related to aseismic ridge subduction, but what might the margin normal transitions illustrated in Figure 13 reflect physically with regard to the mechanics of bathymetric relief subduction and the properties of the overthrusting wedge and/or megathrust? One simple interpretation comes from considering the influence of margin normal variations in rigidity on where and how the upper plate deforms, both to accommodate subducting relief and in response to stresses transmitted across the plate boundary.
Assuming steady slip of the plate interface against constant resistance, the classical critical wedge theory predicts a greater surface slope for a stronger basal fault and/or weaker wedge material [Davis et al., 1983;Dahlen et al., 1984]. The subduction interface beneath the outer wedge is unlikely to be stronger than beneath the inner wedge and we thus follow Byrne et al. [1988] in interpreting the structural contrast across the slope break as reflecting the relative strength of materials comprising the inner and outer fore arc. This interpretation is supported by inferences of rigidity contrasts from normalized thrust earthquake source durations [Bilek, 2010] and seismic velocities in north Chile [Sallarès and Ranero, 2005], the latter of which suggest declining pore-fluid pressure and less fracturing of the overthrusting plate across the inner fore arc. The gradient of the reduction in residual bathymetric amplitudes landward of the slope break ( Figure 13) may reflect the length scale over which the rigidity of the upper plate increases, thereby increasing resistance to buoyancy-driven upward flexure above subducting relief, or may alternatively reflect the increasing dip of the subduction interface.
In Ecuador (Figure 9b), Peru (Figure 9e), and Costa-Rica (Figure 10b), thrust earthquakes with Mw6 are located entirely landward of the trench-slope break and the large residual bathymetric anomalies of the outer fore arc. This distribution is consistent with dynamic Coulomb wedge theory, which associates the trench-slope break with the updip transition from aseismic velocity strengthening to seismogenic velocity weakening behavior on the subduction interface [Wang and Hu, 2006]. Geochemistry, Geophysics, Geosystems

10.1002/2014GC005684
The subducting slab beneath volcanic arcs is typically 80-120 km deep [England et al., 2004;Syracuse and Abers, 2006] and in regions of aseismic ridge subduction always exceeds the thickness of overthrusting plates. In contrast to uplift of the outer fore arc, the presence of a mantle wedge makes it unlikely that the additional elevation of arcs is directly associated with the excess topography of subducting ridges and direct influences of this nature would be expected to also extend across the inner fore arc, which is not observed. An alternate explanation might be that the added thickness of subducting crust increases the volatile flux to the mantle wedge enabling a greater degree of partial melting, resulting in more arc volcanism and a larger arc. Landward of the subducting Emperor seamount chain, Klyuchevskoy is the most active arc volcano on Earth (Figure 8c), but this also is not our preferred interpretation. The forces that drive deformation at subduction zones are transmitted across the plate interface and we suggest near arc uplift may be related to the influence of bathymetric relief subduction on the magnitude of these forces.
In the ''breaking through'' model of Wang andBilek [2011, 2014], the geometrical incompatibility of subducting relief is overcome by permanently deforming the surrounding rocks of both the subducting and overthrusting plates. At low temperature, the deformation is in the form of pervasive fracturing (Figure 14c), which is predicted by analog modeling [Dominguez et al., 1998;Dominguez et al., 2000], and likely reflected in the structural complexity of exhumed subduction zones [Vannucchi et al., 2006] and seismic images in regions of recent/ongoing seamount subduction [Kodaira et al., 2000;Ranero and von Huene, 2000;Bangs et al., 2006;Bell et al., 2010]. Any change in fault geometry in the slip-direction acts to resist fault slip and one significant aspect of the breaking through model relates to the strength of the megathrust. Although creeping faults are typically considered to be weak, and locked faults considered strong, along a rough fault the breaking and wearing of geometrical irregularities within a broad and internally complex fracture network likely results in an integrated geometrical resistance against creep that is likely to be stronger than the locking stress of a smooth fault [Wang andBilek, 2011, 2014]. Fault strength is linked to surface heat-flow by frictional heating and from heat-flow observations across 10 subduction zones, Gao and Wang [2014] show that rough creeping megathrusts dissipate more heat and therefore tend to be stronger than the smooth megathrusts that generate large earthquakes.
The stresses transmitted through the fore arc may thus be greater in regions where aseismic ridges subduct, but why is deformation associated with these stresses focused near the volcanic arc? In Kamchatka [Wang, 1996;Kulinich et al., 2007], Sumatra [Sieh and Natawidjaja, 2000], and Nicaragua [Turner et al., 2007], the partitioning of oblique convergence onto strike-slip faults coincident with the arc may reflect the higher heatfluxes and mechanical weakening of the lithosphere. We suggest thermal weakening may also result in compressional deformation being focused here, resulting in uplift and broadening of the volcanic arc. This hypothesis is consistent with higher earthquake densities and the presence of shallow (<25 km) thrust earthquakes in the overthrusting plate near the Ecuadorian volcanic arc (Figure 9b), and near the Cordillera de Talamanca in Costa Rica (Figure 10b). It is also consistent with force balance calculations relating the elevation of volcanic-arc/highlands to the horizontal push transmitted through the wedge, and hence the shear strength of the megathrust [Smith, 1981;Lamb and Davis, 2003;Lamb, 2006].
The inner fore arc separates regions of lower rigidity in the outer fore arc and near the volcanic arc, providing a stable and relatively strong environment for the formation of fore-arc basins, fore-arc sliver plates, and accumulation of the elastic stresses recovered in earthquakes.

Interpretation of Seamount Morphology and Distributions
The residual bathymetric anomalies we interpret as subducting seamounts are similar in wavelength, amplitude, and morphology to seamounts observed seaward of the trench axis. These observations are inconsistent with models where the primary mechanical effect of seamount subduction is to increase normal stress on the subduction interface [Scholz and Small, 1997]. High normal stresses are generated by the flexural rigidity of the fore arc and such resistance to upward flexure would quickly increase the wavelength and reduce the amplitude of the bathymetric anomaly observed, making any morphological consistency with unsubducted seamounts unlikely.
In Figure 4c, we show that residual bathymetric anomalies associated with the subducting Louisville Ridge seamount chain are correlated with a local reduction (V P up to 1 km s 21 ) in fore-arc wave speeds. In the Mariana fore arc, serpentinite mud volcanoes cluster above large residual bathymetric anomalies associated with Geochemistry, Geophysics, Geosystems 10.1002/2014GC005684 subducting seamounts (Figure 5b). Both observations support the model of Wang andBilek [2011, 2014], in which the subduction of bathymetric relief is accommodated by elevated degrees of upper-plate fracturing.
The form and distribution (Figure 12a) of residual bathymetric anomalies are also considered to be inconsistent with mechanical models proposing wholesale decapitation of subducting seamounts [Cloos, 1992;Scholz and Small, 1997]. The Mariana trench, for example, is located west of a portion of Pacific Plate characterized by a dense distribution of large (>2 km high) seamounts ( Figure 5). Only 3-5, however, are presently interpreted beneath the fore arc. This density is consistent with the frequency with which seamounts will enter the trench over the next 5 Myr and a considerably higher density would be expected if a significant proportion were being accreted.
The ambiguity in identifying regions where seamounts are likely to have recently subducted drops considerably where linear seamount chains are identified. Few seamount chains are parallel with relative plate motion vectors, causing most collision zones to migrate with time. In Tonga, the rate of southward migration may be as high as 180 mm yr 21 [Lonsdale, 1986;Ballance et al., 1989], suggesting that the fore arc north and within 180 km of the contemporary collision zone may have been a site of seamount subduction within the last million years. This region is dashed in Figure 4a and bathymetric anomalies of similar amplitude to those observed in the contemporary collision zone are not observed. The region to the north may contain small fragments of broken up seamounts, but identification of isotopic and trace element signatures in arc lavas requiring up to 40% of a slab-derived Louisville signature in arc volcanoes suggests the bulk of the Louisville Ridge is subducted [Timm et al., 2013].
In Japan, the collision zone between the Joban seamount chain and the Japan trench is migrating north, but large residual bathymetric anomalies are not observed south of the contemporary collision zone. This observation may be interpreted as suggesting seamounts are not accreted, or alternatively, that accreted seamounts lose all residual bathymetric expression. The later scenario is possible as seamounts lose their compensation, or are sufficiently fragmented that they lose their competence, but neither are preferred interpretations.
More generally, Wessel et al. [2010] used satellite-derived bathymetry to predict there may be as many 12,000 large seamounts with height, h, above the surrounding seafloor >1.5 km. Hillier and Watts [2007] used shipboard bathymetry data to predict there may be as many as 40,000 large (h > 1 km) seamounts, 60% of which remain to be discovered, and >200,000 probable seamounts of varying amplitude (0.1 < h < 6.7 km). Given these populations, if large-scale accretion were a process ubiquitous with seamount collision then the composition of both contemporary fore arcs and ancient subduction zones should be largely basaltic. Fragments of accreted seamounts in Japan are much smaller than the average dimensions of modern seamounts, which Isozaki et al. [1990] interpret as suggesting the bulk of colliding seamounts are subducted. The metamorphic grade of these fragments suggests deep lower crustal underplating of seamount material, as opposed to offscraping within an accretionary wedge, as the dominant mode of accretion with less than 1% of samples in the zeolite facies metamorphic grade. These observations are consistent with models where the majority of seamounts are subducted and while fractured seamounts may lose small pieces during subduction, as suggested by geochemical observations of rocks with OIB isotope and trace element signatures within fore arcs [Geldmacher et al., 2008], wholesale accretion is likely to be rare.

Influence on Seismicity
We have shown that regions of contemporary seamount subduction are typically associated with reduced levels of megathrust seismicity. Subduction of the Louisville ridge, for example, is associated with a 200 km wide gap in shallow (<65 km depth) seismicity, the width, and geometry of which suggest an association with the volcaniclastic filled flexural moat (Figure 4b) [Bassett, 2014]. Seamounts subducting beneath the central Mariana arc are similarly characterized by low-earthquake density and appear so against the more seismically active southern segment of the Mariana trench (Figure 5b). Tsunami earthquakes in Java were influenced by subducting relief, but rather than nucleating from or rupturing across a subducting seamount, the 2006 Mw 7.7 earthquake appears to have ruptured a comparatively smooth segment of the plate interface with eastward rupture propagation arrested by a subducting seamount ( Figure  6c). A similar situation is observed in Japan and subduction of the Joban seamount chain is characterized by low-earthquake density, weak interseismic coupling, and may have limited southward rupture propagation in the 2011 Mw 9.0 Tohoko-Oki earthquake (Figures 7b and 7c). The strength of conclusions that can be drawn from these observations is limited by the short duration of instrumental earthquake catalogs, but we do suggest that holistically, our observations of subducting relief are more consistent with associations with small earthquakes and creep than with large megathrust events. This translates into the question: how do bathymetric anomalies subduct aseismically?
The subducting seamounts interpreted above are 1000-3000 km 2 in area, have residual bathymetric elevations 2 km and with respect to wavelength and amplitude, represent the end-member of subducting seafloor roughness (Figure 15a). At seismogenic depths, the ''breaking through'' model suggests that the geometrical incompatibility of this roughness may be accommodated by the development of an upper plate fracture network, the aseismic evolution of which is described by Wang and Bilek [2014]. They suggest that cataclasis and pressure solution are the dominant mechanisms facilitating aseismic creep. Figure 15. Cartoons summarizing the possible influence of subducting seamounts and aseismic ridges on structure and seismicity in subduction zones. Blue and filled red dash represents regions of moderate and high-earthquake density, respectively. Dotted black line marks the fore-arc slope break. (a) Subducting seamounts represent seafloor roughness over a narrow (typically <100 km) wavelength and fore-arc uplift and megathrust complexity are expected over a similar length-scale. The flexural compensation of seamounts may increase the flux of sediments, potentially reducing the strength of the megathrust and resulting in weak creeping zones either side of the seamount. (b) Aseismic ridges are locally compensated, resulting in a strong correlation between ridge physiography its influence on seismicity and plate-interface coupling. Megathrust complexity is greatest on the ridge flanks, explaining high densities of small earthquakes and rupture termination of large earthquakes. The megathrust may be less rough beneath the ridge crest, but remains sufficiently complex to facilitate strong creep and drive the additional deformation of the volcanic arc/highlands.

10.1002/2014GC005684
Within the fracture network, the alternating failure of various parts of the system causes the network to deform almost continuously. Structural and stress heterogeneities limit the failure dimensions of individual fractures, resulting in large numbers of small magnitude earthquakes, as indicated by B-values where seamounts subduct offshore Costa Rica as large as 2 [Ghosh et al., 2008]. Extrapolation of high B-values to low magnitudes implies a large number of brittle events below the resolution of seismometers and brecciation and comminution associated with the fracturing process leads to catalastic flow [Wang and Bilek, 2014]. Highly fractured rocks also allow a wide distribution of fluids, which may facilitate pressure solution creep occurring on the many block and grain boundaries within the fracture network.
The reduced levels of seismicity we observe may reflect the magnitude completeness of the earthquake catalogs we analyze with the numerous small earthquakes and microseismicity predicted by this model likely falling below the resolution (Mw3) of global earthquake catalogs [Obana et al., 2009;Tr ehu et al., 2012]. Where large events are associated with subducting seamounts, these are likely, and have been shown, to be complex [Bilek et al., 2003;Wang and Bilek, 2011].
Subducting aseismic ridges or other large and laterally contiguous bathymetric features represent seafloor roughness over a larger scale than subducting seamounts, and the widths of ridges approach the rupture lengths of large (Mw 7) megathrust earthquakes (Figure 15b). The steepest gradients, and hence the maximum seafloor and megathrust roughness, are associated with the ridge flanks and contrary to subducting seamounts, the lateral continuity of aseismic ridges means they do not have a steep frontal slope during subduction, except when they first collide with a margin. As a consequence, the plate interface across the central section of subducting ridges may be comparatively smooth and sufficiently wide that complex large earthquakes can rupture within this zone. The maximum roughness is associated with subducting ridge flanks, which may explain the high densities of small magnitude earthquakes and the tendency for earthquakes propagating from epicenters above the subducting ridge/plateau to stay within its bounds as is observed for earthquakes in Java (1994) and Costa-Rica (1941 and1950). For earthquakes nucleating toward aseismic ridges, there are many examples of along-strike rupture propagation being arrested near the crest or flanks of these features [Bilek, 2010;Sparkes et al., 2010]. In South America, Bilek [2010] shows that the complexity and normalized source duration of moderate magnitude earthquakes (5.5> Mw >7.1) is greatest on the southern flanks of the subducting Nazca and Carnegie ridges, consistent with our suggestion of maximum roughness and complexity.
It thus seems likely that the mechanical processes by which oceanic seamounts, aseismic ridges, and oceanic plateau subduct are similar, with seafloor roughness promoting upper plate deformation and megathrust complexity enabling creep [Wang andBilek, 2011, 2014]. In addition to where and the wavelength over which this roughness occurs (Figure 15), one additional and potentially important distinction between these features may be the mechanism of isostatic compensation . Ridges and oceanic plateau are locally compensated because they formed on or near a mid-ocean ridge crest, which explains the close spatial correspondence between the residual bathymetric expression of these features and the extent of their influence on megathrust seismicity, geodetic locking (Dashed lines in Figures 9g and 15b) [Chlieh et al., 2011], and near arc physiography. Oceanic seamounts, in contrast, are often flexurally compensated providing local depressions where thick packets of volcaniclastic sediment may accumulate ( Figure  15a). In Tonga, the width of the Louisville seismic gap is a factor >3 wider than the dimensions of seamounts entering the margin, which Bassett [2014] relate to the width of the flexural moat and the influence of higher sediment fluxes on the shear strength of the subduction interface ( Figure 15a). The influence of subducting seamounts may thus extend beyond fracturing and deformation of the overthrusting plate, and what may be as important is their role in facilitating a larger flux of sediment into the seismogenic zone. This may occur either through subduction of a flexural moat [Bassett, 2014], or by altering the geometry of the megathrust [Bell et al., 2010].

Conclusions
An ensemble averaging technique is developed to isolate and remove the long-wavelength, large-amplitude (7.5 km and 200 mGal), trench-normal topography and gravity fields associated with subduction zones. This procedure lowers the dynamic range, removes the steep gradients associated with the seaward and landward trench-slopes and significantly increases the resolvability of short-wavelength/lower-Geochemistry, Geophysics, Geosystems 10.1002/2014GC005684 amplitude structure. Application of this technique to all subduction zones on Earth has resulted in the identification of >200 residual bathymetric anomalies within subduction fore arcs. Isolated anomalies with height 1 km and area 500 km 2 (N 5 36), and aseismic ridges (N > 10), represent confident interpretations of subducting relief. Subducting fracture zones are resolved by residual free-air gravity anomalies up to >200 km landward of the trench-axis.
Isolated residual bathymetric anomalies interpreted as subducting seamounts are located near the trench (within 70 km), at shallow subducting slab depth (SLAB1.0 depth < 17 km) and are similar in wavelength, amplitude and morphology to unsubducted seamounts seaward of the trench-axis. The morphology and spatial distribution of subducting seamounts is not consistent with mechanical models proposing wholesale decapitation and accretion, and the majority of seamounts are likely subducted intact. Subducting seamounts are associated with reduced levels of megathrust seismicity and may have acted as barriers to rupture propagation in large earthquakes in Java and Japan.
Subducting aseismic ridges are associated with uplift and steepening of the outer-fore arc trench-slope, a gradual reduction in residual bathymetric expression across the inner fore arc, and a local increase in the width and elevation of the volcanic-arc/orogen. These contrasting expressions may reflect the influence of margin normal variations in rigidity on where and how the upper plate deforms, both to accommodate subducting relief and in response to stresses transmitted across the plate boundary. The outerfore arc and near-arc regions may have lower rigidity due to fracturing and thermal weakening, respectively. Subducting aseismic ridges are associated with complex earthquakes and fault creep, and the, 'breaking through' model proposed to accommodate the geometrical incompatibility of subducting seamounts may also be applicable to the case of aseismic ridge subduction. This model proposes that bathymetric irregularities subduct by severely deforming and fracturing both the subducting anomaly and the overthrusting plate. The structural complexity and associated geometrical resistance against slip within this fracture network results in a strong megathrust and higher stress transmission across the subduction interface. Near-arc uplift landward of subducting ridges may reflect the predicted local increase in the magnitude of these stresses, which drive deformation that is focused near the arc due to the higher heat flux. Seafloor and megathrust roughness is greatest on the flanks of aseismic ridges, which explains the tendency for earthquakes propagating from epicenters within or away from regions of subducting aseismic ridges to terminate in these regions. The bathymetric expression of subducting ridges places strong constraints on the lateral extent over which subducting bathymetric anomalies influence the surface morphology of the fore arc trench-slope. Understanding the physical controls on these limits is an important subject worthy of future research. Morgan are thanked for helpful comments on an early version of this manuscript. David Sandwell is thanked further for providing an updated version (V23.1) of the global marine gravity grid. Cedric Twardzik and Johnny Hunter are thanked for helpful discussion. Figure construction and grid processing was conducted using GMT [Wessel and Smith, 1991]. This study is global in nature and the grids of residual free-air gravity anomaly and residual topography we generate are freely available as an electronic supplement and online. This study was supported in part by a University of Oxford Clarendon Scholarship and UK Natural Environment Research Council grants NE/F005318/1 and NE/1026839/1. Geochemistry, Geophysics, Geosystems 10.1002/2014GC005684