Earth's deepest earthquake swarms track fluid ascent beneath nascent arc volcanoes

Most of the world's explosive volcanoes are located in volcanic arcs, formed by fluid-fluxed melting of upper mantle rocks. The fluids that facilitate melting are released from subducted tectonic plates as they sink into the mantle. Yet, we have sparse knowledge of the migration pathways of melts through the upper mantle (i.e., between the surface of the subducted plate and arc volcanoes). We are also uncertain of the time required for this migration to occur. Here, we show evidence of two earthquake swarms that occur in the upper mantle beneath the Mariana and Izu-Bonin arc systems. The best-resolved swarm occurs beneath the Mariana arc, where the earthquakes define a sub-vertical pipe-like structure with a diameter of ∼50 km and occurs between depths of ∼10-250 km. To test the robustness of depth locations, we used a fully non-linear grid search algorithm, double-difference relocation, as well as an analysis of pP-P arrival times and depth sensitive phases. In addition, we calculated centroid moment tensor solutions using a 3D Earth model to understand the mechanism of failure within the swarms. These data demonstrate that the sub-vertical earthquake swarm occurs within the upper mantle between the subducted slab and the overriding volcanic arc, with seismicity concentrated within discrete dayto month-long swarms of activity over a single two-year period. The Izu-Bonin example shares a similar subvertical pipe-like geometry with seismic activity bracketed within a two-year period. We infer that these rare earthquake swarms record the ascent of hydrous melt and/or fluid, from dehydration of the subducting plate. This implies that hydrous minerals within subducted slabs continue to dehydrate to depths of at least 200-250 km. Also, the short duration of earthquake swarms implies that fluids/melts can be rapidly transported through the sub-arc mantle at rates in the order of km/hr. This is consistent with rapid melt ascent rates inferred from geospeedometry and experimental petrology and is reminiscent of patterns seen during episodic tremor and slip events. Publication Details White, L. T., Rawlinson, N., Lister, G. S., Waldhauser, F., Hejrani, B., Thompson, D. A., Tanner, D., Macpherson, C. G., Tkalcic, H. & Morgan, J. P. (2019). Earth's deepest earthquake swarms track fluid ascent beneath nascent arc volcanoes. Earth and Planetary Science Letters, 521 25-36. Authors Lloyd T. White, Nicholas Rawlinson, Gordon S. Lister, Felix Waldhauser, Babak Hejrani, David Thompson, Dominique Tanner, Colin G. Macpherson, Hrvoje Tkalcic, and Jason Morgan This journal article is available at Research Online: https://ro.uow.edu.au/smhpapers1/757 Earth and Planetary Science Letters 521 (2019) 25–36


Introduction
Arc volcanism is generated by flux melting of the sub-arc mantle. This is driven by progressive dehydration of a subducting oceanic plate as it experiences greater pressures and temperatures * Corresponding author at: GeoQuEST Research Centre, School of Earth, Atmospheric and Life Sciences, University of Wollongong, NSW, 2522, Australia.
during its descent (Coats, 1962;Schmidt and Poli, 1998;Rüpke et al., 2002Rüpke et al., , 2004Hacker et al., 2003;Tatsumi and Kogiso, 2003). However, we know little about the pathways and mechanisms that transport fluids and melts from the subducted plate through the upper mantle to the overriding volcanic arc. Tomographic studies have shown vertical zones with relatively high attenuation in the upper mantle beneath several island and continental arc volcanoes (Haberland and Rietbrock, 2001;Rychert et al., 2008;Syracuse et al., 2008;McGary et al., 2014). These highly attenu- ating zones have been interpreted to reflect transport pathways of fluids from the dehydrating slab beneath continental arc volcanoes (Haberland and Rietbrock, 2001;Husen and Kissling, 2001;Rychert et al., 2008;Syracuse et al., 2008;McGary et al., 2014). Such interpretations are intuitive based on our current understanding of fluid-fluxed melting of the sub-arc mantle (Rüpke et al., 2004). The magma ascent pathways inferred from seismic tomography (Haberland and Rietbrock, 2001;Rychert et al., 2008;Syracuse et al., 2008;McGary et al., 2014) show similar geometries to shallower (10-60 km) seismicity beneath continental arc systems that have also been inferred to document fluid/melt migration through sub-arc mantle (e.g., Špičák et al., 2004Chang et al., 2017Chang et al., , 2019Hotovec-Ellis et al., 2018;Halpaap et al., 2019). The idea of earthquakes occurring in the sub-arc mantle is counter-intuitive as most earthquakes occur in regions of the Earth known to have mechanical strength (e.g., in the lithosphere above and adjacent to subduction zones, and along fracture zones in oceanic crust). Deeper earthquakes are common within subducted slabs, either in Wadati-Benioff Zones, in paired seismic zones related to dehydration of phases such as serpentine (Green and Houston, 1995;Rüpke et al., 2002Rüpke et al., , 2004Hacker et al., 2003;Tatsumi and Kogiso, 2003;van Keken et al., 2011;Abers et al., 2013), or due to 'slow slip' or 'episodic tremor and slip' (ETS) events located near the interface of the surface of the subducting slab and its overriding plate (Katsumata and Kamaya, 2003;Rogers and Dragert, 2003;Obara et al., 2004;Brudzinski and Allen, 2007;Ghosh et al., 2010;Hyndman et al., 2015). Apart from Wadati-Benioff Zones, there are also reports of enigmatic earthquake occurrences beneath hotspot volcanoes and intra-plate rifts (Wolfe et al., 2003;Keir et al., 2009;Blanchette et al., 2018), or within the mantle at depths of up to ∼100 km (Chen and Molnar, 1983). Most of these intermediate to deep earthquakes occur within cold lithospheric mantle or are associated with weakening of hotter mantle or phase changes/dehydration. However, there is a growing body of evidence that indicates a new type of earthquake might exist in the upper mantle between subducted oceanic lithosphere and an overriding volcanic arc (e.g., Špičák et al., 2004Chang et al., 2017Chang et al., , 2019 and that these earthquakes may delineate fluid migration pathways between subducted, dehydrating lithosphere and overriding volcanic arcs. Such phenomena potentially provide the evidence required to support the theoretical work of Davies (1999) who proposed that rare pulses of fluids could drive fracturing and/or melt channels in the hottest part of the sub-arc mantle, rapidly transporting large volumes of melt towards the surface. These sub-arc earthquakes might also provide insight into the rate of fluid/melt transfer through the mantle, but before we investigate these issues further, we must first ensure that the earthquakes have been accurately located. We therefore report a robust and comprehensive analysis of the 1985-1986 earthquake swarm beneath the Izu-Bonin arc that was first reported by Špičák et al. (2009). We also report the discovery of a much clearer example of a sub-arc mantle earthquake swarm that occurred beneath the Mariana volcanic arc.

Initial analysis of earthquake hypocenter data
We systematically searched two global relocated hypocenter databases to identify geographically distinct zones of seismic activity in the sub-arc mantle. One of these databases was that of Engdahl et al. (1998), which we herein refer to as "EHB". The other database was sourced from the International Seismological Center (2013), which we herein refer to as "ISC". We used the program 'eQuakes' (Lister et al., 2014) to visualize the hypocenter data in time and space, and to make inferences about the geometry of the Wadati-Benioff Zone. Further details are provided in the Supplementary Text.
Our analysis of the EHB and ISC hypocenter data shows the distinct vertically elongated zone of seismicity beneath the Izu-Bonin arc that was reported by Špičák et al. (2009), as well as a similar, but more extensive example beneath the Mariana arc ( Fig. 1; Supplementary 3D Model; Interactive Map File). The two examples are defined by pulses of seismicity that delineate transient, sub-vertical pipes of earthquake activity between the slab and near-surface. The best-resolved example is beneath the Mariana arc/back-arc. Here, the sub-vertical earthquake pipe extends from ∼10 km beneath the surface, through the sub-arc mantle, intersecting the upper surface of the subducted Pacific Plate at a depth of ∼200-250 km ( Fig. 1; Supplementary 3D Model). These earthquakes delineate a sub-vertical, pipe structure ∼50 km in diameter ( Fig. 1; Supplementary 3D Model), i.e. a narrow column within the resolution of the method for earthquake location used in generating the hypocenter databases. The earthquakes defining this structure range in magnitude from 2.9 to 5.2 M w (ISC, 2013), with the majority occurring as a series of seismic events between January 2006 and December 2007 ( Fig. 1; Supp. Fig. 1; Supp. Data 1 and Supp. Video 1). Likewise, the events that define the earthquake pipe beneath the Izu-Bonin Arc are distributed from the near surface to depths of ∼100-150 km and occur within daily to monthly pulses of activity between April 1985 and October 1986 ( Fig. 1; Supp. Data 2).
Hypocenters of the Mariana and Izu-Bonin sub-arc earthquake pipes both occur within two-year time periods, with most events concentrated in pulses occurring over several days to a few months. This is best demonstrated with animations that show the spatial distribution of the earthquakes (Supp. Videos 1 and 2). The animations demonstrate that the sub-arc earthquakes are clearly distinguished from background seismic activity in the down-going slab and the overriding lithosphere. Another animation shows a spatio-temporal analysis of the sub-arc mantle earthquakes, demonstrating that these occur within discrete time periods and are of limited spatial extent (Supp. Video 3).

Hypocenter relocation to test interpretation of sub-arc mantle earthquakes
We hypothesise that the sub-arc mantle earthquakes identified in both the EHB and ISC hypocenter datasets ( Fig. 1) are associated with the ascent of sub-arc fluids/melts. 1 Yet, to test this hypothesis, it is crucial that the earthquake swarms are well-located and truly occur in the upper mantle between the upper part of the subducting slab and the near surface. For instance, the EHB and ISC hypocenter datasets do not include a robust estimate of location uncertainty for individual earthquakes that we consider in this study. As such, it is essential that we recompute the depth and geographic locations of all events to assess their accuracy, and ensure that the observed sub-vertical zone of seismicity are required to explain the data. To this end, we use an absolute location method based on a fully non-linear grid search as well as a relative relocation method based on double differences. As additional tests, we also (1) directly invert differential pP-P arrival times for hypocenter depth -this is the least well constrained parameter in event location, but essential to defining the vertical extent of the earthquake pipes; and (2) identify a number of depth sensitive phases in data recorded by the dense USArray Transportable Array, which 1 A key unknown is the physical and chemical properties of fluids ascending through the sub-arc mantle. We have referred to this material as 'fluid/melt' as it could be several phases over the depth-range of each swarm (e.g. supercritical brines, fluids with a high component of dissolved silicates and/or hydrous silicate melts).  we use to supply independent constraints on source depth (see Supplementary Information for further details).

Earthquake location using a fully non-linear grid search algorithm
We apply a fully non-linear earthquake location algorithm to available arrival times for the earthquakes identified in the ISC (2013) catalogue which define the location and geometry of both the Mariana and Izu-Bonin sub-arc earthquake pipes. The arrival time predictions were calculated using the program ttimes in the presence of the ak135 global reference model (Kennett et al., 1995). The inverse problem of reconciling observed and predicted arrival times was solved using a nested grid search, similar to the technique described by Pilia et al. (2013). Further details of the methodology are provided in the Supplementary Text.
We find that for both the Mariana and Izu-Bonin sub-vertical earthquake pipe swarms, the new event locations ( Fig. 2a and 2b) were sufficiently close to the ISC locations ( Fig. 1d and 1e) that the delineation of the pipes by the earthquakes was essentially unchanged. Furthermore, the average horizontal and vertical uncertainties for events in the Mariana pipe are 11 km and 17 km respectively, while for Izu-Bonin, the equivalent figures are 16 km and 29 km respectively. Despite a similar geographic distribution of picks (see Supp. Fig. 2), the earthquakes which characterize the Mariana pipe appear -based on the above findings -to be more accurately located than those which characterize the Izu-Bonin example. This is likely due, at least in part, to better station coverage and improvements in data quality (and hence more accurate picks) as a result of the modernization of seismic recording equipment over the 20-year period between the Izu-Bonin pipe (1985)(1986)(1987) and the Mariana pipe (2005)(2006). A comparison of the new results with the ISC locations (Supp. Fig. 1) for the Mariana example indicates that there is very little difference in the sub-vertical pipelike geometry, with the mean location difference for all events in latitude, longitude and depth being 6 km, 12 km and 20 km respectively, which is approximately within the 95% confidence region of the non-linear locations. For Izu-Bonin, the equivalent figures are 14 km, 17 km and 42 km, with the latter result reflecting the increased depth distribution of events in the non-linear relocation; the distribution of ISC hypocenters in this case reflect the fact that depths are often fixed at a shallow value (e.g., 33 km).
Based on these results, we conclude that the inferred, narrow, sub-vertical pipes that connect the surface to the subducting slab, as delineated by the earthquake hypocenters, are unlikely to be artifacts caused by relocation error. For additional details, including the potential effects of using a 3D reference model, a comparison of locations obtained using an L1 norm, excluding S phases and excluding sP and pP phases, and plots which only show the most accurately located events, refer to the Supplementary Text, Supp. Figs. 1-5 and Supp. Table 1.

Analysis of depth-sensitive phases
Differential pP-P arrival times were scrutinized to assess the robustness of the depth estimates produced by the catalogue and non-linear locations. For each individual station that records the P and pP phase of these events, an accurate estimate of source depth can be obtained that is relatively insensitive to large-scale lateral variations in wavespeed between the source and receiver.
To perform our analysis, we target earthquakes in the subvertical pipe that have at least one station with P and pP arrival time picks and directly invert the differential pP-P times for depth using the ak135 global reference model (Kennett et al., 1995). This work was also undertaken using the computer program 'ttimes' and correcting for ellipticity. P and pP traveltimes were calculated and their difference taken over a range of depths, which produced the predicted differential arrival time T m (d) as a single valued function of depth d. This was compared to the observed pP-P differential time T o , with the final depth corresponding to the value of d for Assuming an error of the order of 1 s on the observed pP-P times, the uncertainty in the depth estimate is approximately 5 km for an event at 100 km depth. When multiple stations constrain the depth of a single event, the final result is taken as the average, with the standard deviation of the ensemble taken as an estimate of uncertainty.
A comparison between the depths obtained using this method and the non-linear grid search is shown in Supp. Fig. 4. Clearly, the results for both the Mariana and Izu-Bonin examples are in good agreement, and support the presence of sub-vertical earthquake pipes within hot, viscously deforming, regions of the sub-arc mantle. However, the later arriving phases are difficult to pick, and the accuracy of the arrival times of such phases is not quantified in most earthquake catalogues including those provided by the ISC, so some form of independent verification is warranted in order to ensure that we have left no stone unturned.
To this end, we were also able to identify depth-sensitive phases from a sub-set of 11 Mariana pipe earthquakes in USArray data from 2006-2007. Due to the large number of stations available at a distance range of 75 • -93 • , it was possible to enhance the later arriving pP phase by applying adaptive stacking (Rawlinson and Kennett, 2004). By using the same differential pP-P approach for constraining earthquake depth as above, we found that ISC (2013) depths agree with our estimates to within 20 km for 10 out of 11 events (for additional details see Supplementary Text as well as Supplementary Data 1 and 3).

Double-difference earthquake relocation
As an additional measure to better delineate the two subvertical earthquake pipes, we improve the relative locations of the hypocenters using a multi-phase, teleseismic double-difference algorithm (Waldhauser and Schaff, 2007). The double-difference method is a weighted, iterative, least-squares procedure that takes advantage of the similarity in ray paths between a common station and neighbouring earthquakes. It relates the residual between the observed and predicted phase travel time difference for pairs of earthquakes to changes in the vector connecting their hypocenters through the partial derivatives of the travel times for each event with respect to the unknown model parameters.
First and later arriving regional and teleseismic phase arrival time picks from earthquakes listed in the bulletin of the ISC (2013) for the years 1964-2012 near the Mariana and Izu-Bonin pipes were used to form millions of pick delay times between nearby events observed at common stations (Supp. Fig. 6e-f). For each event pair we select only the best delay times (highest pick quality) within bins of 3 • × 3 • to avoid strong spatial clustering of partial derivatives in areas with dense seismic networks (e.g., arrays or local networks reporting to the International Seismological Center).
Delay times in our analysis are predicted using the global travel time model ak135 (Kennett et al., 1995). The use of later arrivals greatly improves the spatial sampling of the partial derivatives, leading to significant improvement in location precision. This is especially true for relative depths as they are constrained by the steep take-off angles of teleseismically recorded mantle and core phases.
During the final iterations of the double-difference procedure the events are linked with at least 5 delay times at inter-event distances less than 80 km. The final teleseismic double-difference locations have RMS residuals of 0.410 s (initial RMS = 1.3 s) in the case of the Mariana pipe, and 1.3 s (initial RMS = 2.0 s) in the case of the Izu-Bonin pipe (Supp. Fig. 6a-b). Relative location uncertainties are derived from a bootstrap analysis of the final residual vector including 200 samples for each event (see Waldhauser and Ellsworth, 2000 for further details). The medians of the major axis of the horizontal and vertical projection of the 90% confidence ellipsoids are 2.8 km and 1.4 km for the Mariana events, and 3.5 km and 3.7 km for the Izu-Bonin events (Supp. Fig. 6c-d). As for the non-linear relocation method, the higher residuals and lower resolution of Izu-Bonin pipe events are due to the lower data quality and fewer stations during the years 1985-1986.
The results of the double-difference analysis provide further support of the two sub-arc mantle earthquake pipes and yields a more refined image of the seismicity within each example ( Fig. 2c and d). The median relative location errors at the 90% confidence level are 2.1 km and 3.6 km for the Mariana and Izu-Bonin examples, respectively (Supp. Fig. 6).

Failure mechanisms and the location of centroid moment tensor solutions
To understand the mechanisms of the sub-arc mantle earthquakes further, we examined events in the Global Centroid Moment Tensor (GCMT) database (Ekström et al., 2012), and selected only those that occurred within either of the two sub-vertical earthquake pipes. Only two events that are recorded in the GCMT catalogue occur within the Izu-Bonin example between 1986 and 1988, so we focused on the Mariana example, where there are 20 pipe events with centroid moment tensor (CMT) solutions between 2006 and 2008. These all report extensional faulting mechanisms and were found to occur between 12-18 km depth (Supp. Table 2). However, there is a major discrepancy in the depth location of the GCMT solutions compared to the equivalent earthquake hypocenters (Engdahl et al., 1998;ISC, 2013), as well as the depths obtained from the relocations performed as part of this study (Supp. Table 2). We considered that these discrepancies may be due to the nature of how GCMT solutions are calculated as well as the earthquakes occurring in a structurally complex region. The primary issue is that GCMT (Ekström et al., 2012) uses the spherically symmetric Earth model PREM (Dziewonski and Anderson, 1981) to calculate synthetic seismograms which are then used to model the observed waveforms in the time windows around body and surface wave arrivals as follows: a lowpass filter with a cut-off period of 125 s was applied, whereas body and surface waves are filtered in the intervals of 120-40 s and 120-50 s (more details are given at: http://www.globalcmt . org /CMTfiles .html). Others have since adopted 3D velocity models in the processing of CMT solutions and in doing so, have shown that there may be tens of kilometer discrepancies in horizon-tal and vertical location, as well as variations in moment-tensor components compared to what is listed in the GCMT catalogue (Hjörleifsdóttir and Ekström, 2010;Valentine and Trampert, 2012;Hejrani et al., 2017). We therefore recalculated the CMT solutions for the 20 earthquakes in the GCMT catalogue that occur within the pipe-like structure. These calculations use the Earth mantle model S362ANI (Kustowski et al., 2008) combined with the Earth crustal model CRUST 2.0 (http://igppweb .ucsd .edu /~gabi /rem .html) and follow the technique used to calculate Continental Scale CMT using a 3D Earth model (CSCMT-3D) (Hejrani et al., 2017).
We adopted the methodology used by Hejrani et al. (2017) to perform a CMT inversion. According to this approach, six elementary ground-displacement seismograms are generated at trial depths according to six basic mechanisms, prior to the inversion: (1) where M m are the basics mechanisms from Equation (1) and a m are six corresponding coefficients. Therefore, the ground motion due to an arbitrary moment tensor can be written as follows: where D i (x) is the observed data (seismograms), G ij are the elastodynamic Green's functions between the source (x 0 ) and receiver (x). Equation (3) can be simplified as: where D is a matrix of the observed displacement seismograms, E is a matrix of the synthetic data due to six basic mechanisms: M 1 , . . . , M 6 , and A is the vector of six coefficients. We estimate the vector A by the least squares method and grid-search the depth and time of the events (Kikuchi and Kanamori, 1991;Hejrani et al., 2017). In this study, we calculate the synthetics, E in Equation (4), on 16 trial points with intervals of 10 km, starting from the depth of 15 km and ending at the depth of 165 km. We use the spectral element code, SPECFEM3D GLOBE 7.0.0 (Komatitsch andTromp, 2002a, 2002b;Komatitsch et al., 2015) published under the GPL 2 license, to calculate synthetics for 15 broadband stations in distance ranges of 500-2800 km surrounding the focal area. The data for these stations are accessible via the IRIS data management center (https://ds .iris .edu /wilber3 /find _event).
The quality of all waveforms were checked manually and poor quality waveforms were removed prior to the inversion. We then performed a moment tensor inversion for 20 events reported in the GCMT catalogue (Ekström et al., 2012) in several period bands. We tested period bands of 100-33 s, 100-28 s and 100-25 s. By systematically increasing the frequency content to capture shorter periods, we improved the accuracy of the depth parameter. An example of the grid search for periods 100-28 s for the event which occurred at 07:27:55 on 2007/03/26, is provided in Fig. 3. This event was reported to have occurred at the depth of 12 km (fixed during the inversion) in the GCMT catalogue. Our analysis shows the depth of 95 km to be optimal ( Fig. 3a-b). Our CSCMT-3D grid search shows a local minimum at the depth of 15 km as well, but depths between 46-125 km show higher overall correlation ( Fig. 3a-b). A good fit was obtained for all waveforms in the inversion (Fig. 3c-d).
This analysis resulted in two potential depth and tensor solutions for most of the events (Fig. 3; Supp. Table 2; Supp. Data 4). The fact that two solutions were obtained for many of the events is due to the limited data coverage in this remote part of the Pacific, which precludes a more definitive result. However, this does not mean that the two solutions should be weighted equally. For example, we obtained 'deep' solutions (∼>30 km) for all of the analysed earthquakes, but 'shallow' solutions (∼<20 km) were only possible for 12 of the 20 earthquakes (Supp. Table 2). Of these, there is good agreement between 16 'deep' solutions when compared to the depths reported for the equivalent hypocenters listed in the EHB and ISC catalogues as well as depths obtained from the non-linear and double-difference methods discussed above (Supp. Table 2). Three of the 'shallow' CSCMT-3D solutions record depths that are similar to those reported in the GCMT database (Ekström et al., 2012), and these also correspond with depths previously reported for the equivalent hypocenters (Engdahl et al., 1998;ISC, 2013), as well as those obtained from the non-linear and double-difference methods discussed above (Supp. Table 2). The fact that multiple location methods yield similar depths provides strong support for the existence of the pipe-like structure within the sub-arc mantle.
We also note that the waveform inversion procedure, largely as a result of the trade-off between depth and focal parameters, modified the source mechanism of each event identified in the GCMT database (Ekström et al., 2012). Instead of all events recording an extensional mechanism as per GCMT, the 'deep' solution indicates that these events show vertical compensated-linear-vector-dipole (CLVD) focal mechanisms with dominant vertical tension axes. The double-couple percentage of these events is typically lower than 90%, with most events yielding the double-couple values of 10-50% (Fig. 4a). Such solutions have been found to be associated with melt/fluid/gas migration before volcanic unrest near the Earth's surface (e.g. Shuler et al., 2013), but other explanations such as rupture on a ring fault or superposition of rotating reverse faults have been suggested (e.g., Tkalčić et al., 2009). There is no apparent correlation between the double-couple value and the depth of the events within the sub-vertical pipe (Fig. 4b), yet events that yield a preferred depth solution of 0-20 km yield moment tensor solutions with a vertical pressure axis ( Fig. 4c; Supp. Data 4), while deeper events yield vertical tension axes (Fig. 4d) or indicate strike-slip failure (Supp. Data 4).

Do the earthquakes delineate the path of fluid/melt transport?
The earthquake swarms reported here (and in Špičák et al., 2004Chang et al., 2017Chang et al., , 2019Hotovec-Ellis et al., 2018;Halpaap et al., 2019) suggest that ruptures can occur within the upper mantle above a subducted slab. This differs from conventional understanding that the relatively hot upper mantle is capable of deforming by viscous creep (e.g. Karato and Wu, 1993;Burgmann and Dresen, 2008), at rates that would preclude the buildup of significant levels of deviatoric stress. The mantle earthquake pipes must therefore have been driven by a perturbation that was occurring at a time scale faster than the time required for viscous relaxation to take place, likely driven by the addition of fluid into the sub-arc mantle that are ultimately introduced from dehydration reactions of hydrous minerals such as chlo-  (Ekström et al., 2012) solution is also shown for a comparison. (B) The color of the 'beach balls' marks the double couple (DC) percentage of the moment tensor. The optimum CMT solution is reflected by higher correlation values, yet we report the range of depths with correlation values that fall within 5% of the optimum value as is indicated in (A) and Supp. Table 2. (C) The location of broadband seismic stations used in to calculate the CSCMT-3D moment tensor; (D) The fit between the observed (thick gray) and the CSCMT-3D solution's synthetic waveform (black), as well as the GCMT solutions synthetic waveform calculated using the 3D Earth model S36ANI (red). The numbers above and below each waveform indicate the correlation between the observed and synthetic values for the CSCMT-3D solution (black) and the GCMT solution (red). This demonstrates that the GSCMT-3D moment tensor solution provides a better fit to the observed waveforms, compared with the GCMT solutions. (For interpretation of the colors in the figure, the reader is referred to the web version of this article.) rite, antigorite, phengite, lawsonite and zoisite (e.g. Peacock, 1990;Ulmer and Trommsdorff, 1995;van Keken et al., 2011), or phases that are stable to even greater depths (Gemmi et al., 2011).
We therefore propose that the earthquake swarms record the near-vertical pathways of slab-derived fluid induced melts, during their upward migration towards volcanic centres in the overrid- ing plate (Fig. 5). The fluid/melt migrated rapidly upwards through the upper mantle and beneath the back-arc spreading system (Fig. 5), intruding within a short-enough time to induce deformation of the mantle (i.e. shorter than the viscoelastic relaxation time τ R = μ E ≈ 10 18 Pa-s/70-100 GPa = 10 7 s or 0.3 y for a representative viscosity of 10 18 Pa-s and Young's modulus 100 GPa).
Release of slab fluids into the upper mantle will ultimately flux melting. Some have theorized that fluid/melt transport could be driven by vug-waves (Morgan and Holtzmann, 2005) or a process akin to 'hydro-fracking' (Shaw, 1980;Nicolas and Jackson, 1982;Davies, 1999). We consider that the sub-vertical earthquake pipe seismicity records these processes, where failure is driven by the pressure difference between the wall rock and the fluid/melt within the overall sub-vertical fracture network. The pressure differences between the dyke-wall and the fluid/melt in a connected crack system will increase rapidly upwards due to the density differences between the wall rock and the pressurised fluid/melt. Therefore melt-fracking should be expected to commence once the vertical component of the connected length of the fracture network exceeds a critical magnitude. As the majority of mechanisms obtained from the CSCMT-3D data have vertical CLVD focal mechanisms with dominant vertical tension axes ( Fig. 4d; Supp. Data 1 and 4), we propose these represent fluid/melt migration along what would be equivalent to a growing dyke network nearer the surface or that the earthquake locations could otherwise represent failure and collapse of wallrock into spaces created once the fluid/melt phases have migrated upward.
To test our hypothesis further, we examined major and trace element geochemical data obtained from igneous rock samples collected from the seafloor above the Mariana and Izu-Bonin earthquake pipes (downloaded via the EarthChem portal (http://ecp . iedata .org, accessed 2 August, 2016)). These data were assessed to determine what lithologies existed above and in the vicinity of each example, as well as to determine if these rocks show any influence from slab-derived fluids. We examined 219 samples in the region above the Mariana pipe and 123 samples in the region above the Izu-Bonin pipe. We reclassified the lithology of each sample within the database according to the lithology field defined on a total alkali vs. silica (TAS) plot (Supp. Fig. 7). The reclassified results were imported into ArcGIS to examine their spatial distribution. We then plotted the trace element data for samples that were collected within the total radial extent of the Mariana and Izu-Bonin pipes. These data were plotted on 'spider plots' to examine relative abundances with respect to N-MORB (Sun and McDonough, 1989) (Supp. Fig. 7). The geochemical data show that a range of lithologies exist above each earthquake pipe; these include, basalt, basaltic-andesite, andesite and dacite (Supp. Fig. 7). The majority of samples above the two earthquake pipes are basalt and basaltic-andesite, yet these display similar trace element abundances to the more fractionated andesite and dacite samples. The spider plots demonstrate that these samples are depleted in high field strength elements and are enriched in heavy rare earth elements, as is typical for subduction derived melts. This implies an arc-magma influence on magma compositions in this region, although, more complex interactions between recent magmas and former arc-crust remain a potential alternative origin for this signal.
The rocks found in the crust above each sub-arc mantle earthquake pipe appear to show a recent influence of arc volcanism. These areas are active zones of arc-volcanism and hydrothermal activity and there are volcanic edifices at the surface within the vicinity of the earthquake pipes (Supp. Fig. 8). However, there is no direct record of volcanic eruptions in the ocean crust or volcanoes above the identified earthquake pipes at the time of their seismic activity (http://volcano .si .edu -accessed 16 March, 2016). The volcanoes in these regions are typically submarine, so observing volcanism linked to swarm hypocenters would require monitoring or searching for recent melt emplacement on the seafloora prospective target for future seafloor surveys. However, it is also possible that ascending fluids or melts evidenced by sub-arc mantle seismicity could pond within the overriding oceanic lithosphere for some time before eruption (Macpherson et al., 2010). We therefore propose that the earthquake pipes may be recording fluid/melt migration beneath established arc volcanoes (e.g. the Izu-Bonin example) or beneath newly forming arc volcanoes with no record of eruptions or no topographic expression on the surface of the ocean-floor (e.g. the Mariana example) (Supp Fig. 8).
It is also worth considering why the sub-arc mantle earthquakes are apparently rare features since they have only been recognized as discrete pipes beneath the Aleutian Islands, Equador/ Columbia and the Izu-Bonin-Mariana arc (Špičák et al., 2004Chang et al., 2017Chang et al., , 2019Hotovec-Ellis et al., 2018, this paper). Their rarity either relates to the relative sparsity of data in the available ∼50-year catalogue of hypocenter locations, or because particular conditions must be met to cause failure. We are not aware of distinct features or slab geometries that clearly explain the reported sub-arc mantle earthquake pipes. However, one could speculate that the stress regime in the crust above the earthquake pipes may influence their location (e.g., decompression associated with back-arc spreading) or that a sudden release of fluids from the down-going slab may drive such features. One interesting point is that the Izu-Bonin and Mariana earthquake pipes consist of earthquakes with magnitudes of ∼5 Mw; we speculate that similar sub-arc mantle earthquakes are more common, but rupture slowly or have smaller source dimensions and therefore are not recorded in global hypocenter catalogues. We cannot resolve such issues here, but expect that additional examples will emerge in future, and that together, these will shed further light on the physical processes that are responsible for sub-arc mantle seismicity.

Fluid/melt ascent rate
If we assume that a burst of seismic activity within the pipe represents transport of fluid/melt across the full depth of the Mariana conduit (a distance of ∼200 km) over a two-year period we can obtain an estimate for the maximum velocity of fluid transport of ∼0.01 km/hr. This estimate is slower than the 0.04-36 km/hr range of velocities reported for propagation of mantle-derived dykes estimated from (a) computed settling rates of entrained xenoliths, and (b) laboratory-derived rates of depressurization reactions (Dalton and Wood, 1993), but is faster than the approximate (∼0.0001 km/hr) minimum rates estimated from 226 Ra-230 Th isotopic data collected from arc volcanoes (Turner et al., 2001).
Such rates are similar to the ∼0.003 km/hr ascent rates proposed from geospeedometry of Ni zonation in olivine mantle xenocrysts (Ruprecht and Plank, 2013). However, our estimate assumes that fluids migrate from the base to the top of the pipe over the course of the observed lifetime of an individual swarm. It is also possible that fluid/melt could ascend at higher velocities in a series of pulses within the two-year periods of activity that we observe beneath the Izu-Bonin and Mariana arcs.
To investigate this further, we examined the spatial and temporal distribution of the double-difference relocated hypocenters, paying particular attention to periods where multiple ruptures occurred over a range of depths and over one to several days (Supp. Fig. 9, 10, 11, Supp. Data 5). We show that "bursts" of activity exist within the two-year periods of seismic activity in both the Izu-Bonin and Mariana pipes -similar to shallower earthquake swarms beneath Mammoth Mountain, California (Hotovec-Ellis et al., 2018). The "bursts" occur along thin sub-pipes and show a general trend of shallowing depths over the course of days, particularly if the shallowest events (depths of <50 km) are excluded from the Mariana example (Supp. Fig. 10). We determined the range of possible rates of shallowing within each of these "bursts" by calculating the velocity between multiple hypocenters, where a velocity was determined for each hypocenter pair provided that an event was younger and shallower than another (see Supp. Data 5).
What was unclear from the catalogue data was whether these bursts of activity were travelling from deeper to shallower levels through time, i.e. in the pattern we would anticipate if the earthquakes reflect upward fluid/melt transport from the subducted slab to the near surface. To test for patterns associated with ascending fluid movement, we plotted the time of each earthquake in the relocated seismicity as a function of hypocenter depth, as well as in map and section views (Supp. Fig. 9-11). These spatio-temporal plots show that activity within the pipes occurs in bursts, some lasting only a few days, with the events occupying the entire depth range of the pipes. We identified several such bursts and investigated the corresponding location of the events within the pipes. The burst events locate within narrow conduits of approximately 10 km in diameter that extends from the near slab interface to the surface (Supp. Fig. 9-10), with events near the surface (<30 km depth) showing a broader lateral distribution.
We find that the data within each "burst" also contains noise, and recognize that there may be some inherent bias in selecting the connectivity between different earthquakes through time. We therefore devised an unbiased routine to determine the velocity of fluid/melt transport. Our method for calculating the ascent rate for each earthquake pipe is underpinned by the assumption that sequential earthquakes represent fluids/melts ascending through the Earth's mantle that originated from the dehydrating subducting plate. This assumes that fluid/melt migration occurs in a net upward trajectory. This is consistent with observations made from high-resolution studies of seismicity of upper crustal volcanic systems, even where there is some evidence of minor downward fluid migration (e.g., Shelly et al., 2013).  6. Box plots showing a summary of the ascent rates that were calculated for day-to-several-daily bursts of seismic activity identified within the Izu-Bonin and Mariana pipes. The median ascent rates for each "burst" (identified in Supp. Fig. 9 and 10) are shown by a bold horizontal line within each box. The height of each box represents 50% of the data; the height of the thin vertical line above and below each box represents 95% of the data; the colored dots indicate values that fall outside the 95% confidence limit. We worked on the basis that fluid/melt ascent could occur between any two earthquake hypocenters, provided that the later event was shallower than the earlier event. This means that in the majority of cases there are multiple solutions for any given hypocenter. We calculated the velocity (change in vertical height/ time) for each of the hypocenter pairs. The results of this procedure are summarised in Supp. Table 3 and are also plotted as box-plots to compare the median velocity values and range of values obtained from both the Izu-Bonin and Mariana pipes (Fig. 6). These calculations do not formally consider the effects of uncertainties in space or time; however, these uncertainties should have minimal impact on the calculated ascent rates because the timing of each event is well defined and because the medians of the major axis of the vertical projection of the 90% confidence ellipsoids is 1.4 km for the Mariana events and 3.7 km for the Izu-Bonin events (Supp. Fig. 6 and 11). We also conducted additional tests to determine if the results of these calculations might be generated using a series of randomized date/time and depth values (see Supplementary Information, Supp. Figs. 12-14 and Supp. File 5 for further details and discussion of the suitability of the randomly generated data).
The results of potential ascent rates that were obtained for each pipe are summarised in Fig. 6 and Supp. Fig. 14. This image shows that similar ascent rates were obtained from the Mariana and Izu-Bonin earthquake pipes. The median velocity for the four Izu-Bonin "bursts" is 1.95 km/hr, similar to a median velocity of 1.27 km/hr for the four Mariana "bursts" (equivalent to ∼0.3 m/s). These velocity estimates assume that the seismic events are connected and record upward fluid migration on the basis of the vertical-tension axes obtained for the revised moment tensor solutions for the majority of events ( Fig. 4d; Supp. Data 4). It is also possible that these represent cascading failure events that sweep through a subarc mantle that is already prone to failure, perhaps due to minor perturbations in differential stress in sub-arc mantle that already contains slowly migrating fluids. However, the calculated ascent rates merit consideration, if only as an estimate of the fastest possible rates of fluid transport in a sub-vertical pipe swarm. It is also worth considering how these rates compare to transport rates proposed from independent observations. For instance, our calculated median rate of ∼1 km/hr is in agreement with the proposed monthly to yearly time scale required for mixing of mantle melts, ascent through the crust, hybridisation in the upper crust and andesitic eruptions such as Irazu in Costa Rica (Ruprecht and Plank, 2013), and approximately double the inferred minimum rates of ∼0.7 km/hr inferred from dyke transport modelling required to preserve the observed Ni zonation in olivine (Ruprecht and Plank, 2013). To add further perspective, the median ascent rate velocities that we infer from the earthquake pipe events are an order of magnitude slower than the several meters per second (∼>10 km/hr) rates proposed for kimberlite ascent (Sparks, 2013). Therefore, the rates that we propose for melt transport on the basis of seismicity are certainly viable. However, these rates are much faster than those that have been applied to numerical modelling of melt migration within sub-arc mantle (e.g., 0.8-8.0 cm/yr: Dahm, 2000).

Conclusion
Earthquake hypocenter data were examined in hundreds of cross-sections of modern-day island arc subduction zones. This work shows that there are two clear examples of sub-vertical swarms of earthquakes within sub-arc mantle (i.e., mantle wedge). One would assume that the physical properties of this part of the Earth would prevent earthquakes from occurring. As such, we used a fully non-linear grid search algorithm, double-difference relocation, as well as an analysis of pP-P arrival times and depth sensitive phases to assess the veracity of the location of the earthquakes that define these pipe-like sub-arc mantle earthquake structures. The results of these tests provide further support for earthquake activity within the upper mantle. We propose that these rare examples of mantle seismicity are related to the dehydration of the subducting plate which introduces fluids into the upper mantle, modifying the physical properties of these rocks and causing them to fracture in a brittle, rather than ductile manner. This assertion is supported by the available geochemical data obtained from previous dredge samples collected from the seafloor above the earthquake pipes. In addition, as the seismic activity that defines each pipe occurred within a two-year timeframe, we can tentatively use these observations to infer the potential rates of fluid/melt ascent between the subducted slab and the overriding lithosphere. These calculations can be used to independently assess existing inferences of fluid/melt transport that have been based on different isotopic data (and the assumptions made in the interpretation of these data).