TY - JOUR
AU - Peyrat, Sophie
AB - Abstract We use seismic tomography to investigate the state of the supraslab mantle beneath northern Chile, a part of the Nazca-South America Plate boundary known for frequent megathrust earthquakes and active volcanism. We performed a joint inversion of arrival times from earthquake generated body waves and phase delay times from ambient noise generated surface waves recorded by a combined 360 seismic stations deployed in northern Chile at various times over several decades. Our preferred model shows an increase in Vp/Vs by as much as 3 per cent from the subducting slab into the supraslab mantle throughout northern Chile. Combined with low values of both Vp and Vs at depths between 40 and 80 km, we attribute this increase in Vp/Vs to the serpentinization of the supraslab mantle in this depth range. The region of high Vp/Vs extends to 80–120 km depth within the supraslab mantle, but Vp and Vs both increase to normal to high values. This combination, along with the greater abundance of ambient seismicity and higher temperatures at these depths, suggest that conversion from basalt to eclogite in the slab accelerates and that the fluids expelled into the supraslab mantle contribute to partial melt. The corresponding maximum melt fraction is estimated to be about 1 per cent. Both the volume of the region affected by hydration and size of the wave speed contrasts are significantly larger north of ∼21°S. This latitude also delimits large coastal scarps and the eruption of ignimbrites in the north. Ambient seismicity is more abundant north of 21°S, and the seismic zone south of this latitude is offset to the east. The high Vp/Vs region in the north may extend along the slab interface to depths as shallow as 20 km, where it corresponds to a region of reduced seismic coupling and overlaps the rupture zone of the recent 2014 M8.2 Pisagua earthquake. A potential cause of these contrasts is enhanced hydration of the subducting oceanic lithosphere related to a string of seamounts located on the Iquique Ridge of the Nazca Plate. Interferometry, Seismic tomography, Subduction zone processes, Continental margins: convergent, South America 1 INTRODUCTION It is generally agreed that subducting oceanic lithosphere transports an abundance of aqueous fluid and hydrated minerals into the Earth (e.g. Ulmer & Trommsdorff 1995; van Keken et al.2011). While much of this fluid is released at shallow depths, it is thought to be responsible for hydration of the mantle wedge beneath the overriding plate (Peacock 1993; Peacock & Hyndman 1999; Iwamori 2000; Kamiya & Kobayashi 2000; Bostock et al.2002). Depending on the age and speed of the subducting lithosphere, a significant amount of water may be carried to depths in excess of 200 km (van Keken et al.2011). The effects of the transportation and release of water in subduction zones are profound. For example, water can lower the melting temperature of the mantle by some 400 °C, resulting in the generation of magma that eventually finds its way to the surface in the form of a volcanic arc (Gill 1981; Moran et al.1992; Peacock 1993; Plank & Langmuir 1993; Yogodzinski et al.1995; Stern & Kilian 1996; Kelemen et al.1998; Kerrick & Connolly 2001; Grove et al.2006, 2009). Dehydration reactions and phase transitions within the subducted slab are believed to be responsible for the intermediate-depth seismicity within slabs (e.g. Kirby et al.1996; Peacock 2001; Omori et al.2002, 2004; Hacker et al.2003; Abers et al.2013). Water released in subduction zones by phase transition from basalt to eclogite (between 1.2 and 3.3 per cent, Hacker 1996) appears to be important in forming serpentine in the mantle wedge above the slab (Fyfe & McBirney 1975; Bostock et al.2002; Hacker et al.2003; Hacker 2007). If the supraslab mantle is serpentinized in this corner of the uppermost mantle, its low shear strength may control the down-dip rupture limit of great megathrust earthquakes (e.g. Hyndman et al.1997; Bostock et al.2002). Because elastic wave speeds are sensitive to both serpentinization and partial melt (e.g. Hammond & Humphreys 2000; Kerrick 2002; Carlson & Miller 2003; Hacker et al.2003; Christensen 2004; Watanabe et al.2007; Karato 2012; Ji et al.2013), empirical evidence for the role of water in subduction has largely come from seismic imaging and earthquake locations. For example, seismic velocity models such as those from Cascadia (Bostock et al.2002), Costa Rica (DeShon & Schwartz 2004), New Zealand (Eberhart-Phillips & Bannister 2010) and Japan (Kamiya & Kobayashi 2000; Zhang et al.2010), suggest that serpentinization in the mantle wedge is on the order of 10–20 wt%. While these studies are enlightening, neither the extent to which the supraslab mantle is hydrated nor the amount of the partial melt generated within that part of the mantle is well known. In this paper, we investigate the pervasiveness of hydration and partial melt in the supraslab mantle beneath the Nazca-South America Plate boundary in northern Chile by generating elastic wave speed images through joint inversion of observations of local earthquake body waves and ambient noise surface waves. Our motivation for choosing this particular area is partly the availability of an extensive data set directly above a mantle wedge, but also because the frequent occurrence of megathrust earthquakes and the prevalent arc volcanism along this plate boundary allows an opportunity to investigate the consequences of variations in supraslab wave speed anomalies for Andean margin seismicity and volcanism. 1.1 Overview of the northern Chile Plate boundary The central Andes evolved over the past ∼200 Ma as a series of magmatic arc systems related to the subduction of oceanic lithosphere beneath the South American continental margin (e.g. Isacks 1988; Scheuber et al.1994; Allmendinger & Gubbels 1996; Allmendinger et al.1997). Large volumes of calc-alkaline lavas and related plutons have been emplaced along trench-parallel provinces since the Jurassic, showing a general trend of eastward migration (Coira et al.1982). The current magmatic arc commenced in Middle Miocene times (∼27 Ma) (Baker & Francis 1978), and deposited significant volumes of volcanic material in northern Chile during Neogene times. The arc remains very active along its ∼2000 km length. The forearc region of the central Andes typically is subdivided into four morphotectonic units (e.g. Reutter et al.1988; Fig. 1). From west to east, these include: (1) the Coastal Cordillera, which is the Jurassic age magmatic arc, (2) the Central Depression, a Jurassic back arc basin and mid-Cretaceous magmatic arc covered by younger sedimentary rocks, (3) the Precordillera, a late-Cretaceous magmatic arc with Mesozoic sediments lying unconformably upon metamorphic basement and plutonic rocks and (4) the Western Cordillera, the currently active volcanic arc with peaks in excess of 6000 m. The missing Jurassic forearc is presumed to have been tectonically eroded; the eastward migration of the volcanic front from the Coastal Cordillera to the Western Cordillera suggests that an EW section of more than 200 km has been consumed (Rutland 1971; Ziegler et al.1981; Scheuber et al.1994). Figure 1. Open in new tabDownload slide Simplified geologic map of northern Chile. The arrow to the left of the trench indicates the direction of the Nazca plate relative to South America. Topography in the forearc is highly variable. There is no well-developed accretionary wedge near the trench (Von Huene et al.1999), and an extensional forearc is linked to the coast by a ∼1000 m escarpment. The subdued relief of the Coastal Cordillera and Central Depression give way in the east to the steep western flank of the Precordillera and high altitudes of the Western Cordillera. The Nazca-South American Plate boundary currently has one of the fastest convergence rates on Earth at about 6–7 cm yr−1 (DeMets et al.1990, 1994; Altamimi 2002; Sella et al.2002). As a result, this boundary has produced on average one earthquake with magnitude Mw ≥ 7.9 every 10 yr over the past 100 yr. The segment of the boundary in northern Chile presently is considered to be a seismic gap, having accumulated more than 9 m of slip deficit during the last 130 yr of interseismic coupling, and believed capable of generating an Mw 9.0 earthquake (Fig. 2). The recent 2014 Mw 8.2 Pisagua earthquake in northern Chile is considered by some to be a potential foreshock of that larger event. Figure 2. Open in new tabDownload slide Seismotectonic setting of northern Chile and southern Peru. Left-hand panel: estimated ruptures of historic seismicity are shown at the far left; the bold purple line corresponds to the last rupture of the northern Chile seismic gap in 1877. Recent ruptures are indicted on the map by the year in which they occurred. Map colours correspond to bathymetry of the Nazca plate as shown in the palette at the lower right. Ages of the subducting Nazca plate are shown next to dashed lines. Thin rectangle locates the map on the right. Right-hand panel: interseismic plate interface backslip rate, coseismic slip (white contours) for 2014 (Mw 8.2) Pisagua main shock and largest (Mw 7.7) aftershock, post-seismic slip (blue contours) induced by Pisagua main shock event, all determined from GPS observations by Ortega-Culaciati et al. (2015). The epicentre (blue star) and focal mechanism are shown for the 2014 Pisagua main shock event. Coseismic slip (white contours) is also shown for the other recent large earthquakes in the study area: the Antofagasta (Mw 8.1) 1995 (Pritchard & Simons 2006) and Tocopilla (Mw 7.7) 2007 (Béjar-Pizarro et al.2010) events. Northern Chile and southern Peru together have experienced great earthquakes in the past (e.g. the 1868 Mw 8.7 southern Peru and 1877 Mw 8.9 northern Chile events (Comte & Pardo 1991)) (Fig. 2). By contrast, few shallow earthquakes related to the interplate coupled region have occurred within the forearc. Among rare examples is the 2007 Mw 5.7 Pisagua earthquake, located in the outer forearc region during the interseismic stage of the subduction seismic cycle (Allmendinger & González 2010). Modelling of GPS and InSAR data suggests that during the interseismic phase of the earthquake cycle, the plate contact in northern Chile is locked to 35 km depth with a transition zone (brittle–ductile) between 35 and 55 km depth, but with variations in those limits along strike (Norabuena et al.1998; Bevis et al.2001; Khazaradze & Klotz 2003; Chlieh et al.2004; Pritchard & Simons 2006; Béjar-Pizarro et al.2010; Métois et al.2013). Recently, Ortega-Culaciati et al. (2015) re-analysed decades of geodetic data and obtained a detailed map of interplate coupling along the seismic gap in northern Chile (Fig. 2). To first order, they observe that the part of the boundary between Pisagua (∼21°S) and Antofagasta (∼23°S) is highly coupled, while that north of ∼21°S is heterogeneous but generally less coupled. They found an isolated, highly coupled patch offshore of Pisagua surrounded by low coupled areas, in agreement with the findings of Métois et al. (2013). The coseismic slip of the 2014 Mw 8.2 Pisagua earthquake occurred in this isolated patch while the postseismic slip occurred in the neighbouring low coupled regions. Despite this variability, interseismic coupling models suggest that the northern Chile seismic gap remains capable of generating a large megathrust earthquake. 2 DATA AND METHODOLOGY The data analysed in this study are a combination of arrival times from earthquake generated body waves and phase delay times from ambient noise generated surface waves recorded by seismic stations deployed in northern Chile for different periods of time over the past three decades (Fig. 3). Figure 3. Open in new tabDownload slide Locations of seismic stations used in this study, with symbols indicating provenance as shown in the insert at the upper right. The location of the northern Chile seismic gap is shown by the dashed purple line. 2.1 Body waves from local earthquakes The abundance of local earthquake data collected from this region (Fig. 3 and Table 1) is largely in response to the frequent megathrust events and the anticipated end of the current seismic cycle in northern Chile. For example, a rapid deployment of 10 broad-band stations before, and 16 broad-band stations after, the 2014 Mw 8.2 Pisagua main shock by the Oficina Nacional de Emergencia (ONEMI), the Departamento de Geofísica of the Universidad de Chile (DGF) and the Chilean Centro Sismológico Nacional (CSN), together with the Iquique Local Network (ILN-IPOC) and GeoForschungsZentrum (GFZ) deployments, resulted in the highest density of stations ever deployed in this region. Table 1. Summary of the deployments that recorded the data used in this investigation. Operator and network abbreviations refer to the Universidad de Chile (UCH), Universidad de Antofagasta (UA), Institut de Recherche pour le Developpement (IRD) the Institut de Physique du Globe de Strasbourg (IPGS), GeoForschungsZentrum (GFZ), Freie Universität Berlin (FUB), the Integrated Plate Boundary Observatory Chile (IPOC), the Chilean Centro Sismológico Nacional (CSN) seismic networks, the Oficina Nacional de Emergencia (ONEMI), the Departamento de Geofísica of the Universidad de Chile (DGF), and the Iquique Local Network (ILN-IPOC). Instrument types are short period (sp) and broad band (BB). Start and end times of deployment are in year and Julian day (JD). Data from IPOC (doi:10.14470/PK615318) is available online via the GFZ website. All other data available on request from DGF. Seismic network . Start . End . . Operator . Name . Year . JD . Year . JD . Number and type of sensor . IRD-UCH-UA ANTOFAGASTA PERMANENT 1990 313 1995 209 10 sp DGF IQUIQUE 1991 1991 194 1991 227 21 sp DGF-GFZ-FUB CALAMA 1994 1994 88 1994 114 34 sp + 5 BB IRD-IPGS-DGF ARICA 1996 1996 173 1996 225 46 sp UTA ARICA PERMANENT 1997 91 2002 273 16 sp DGF TMP 1 2005 349 2006 38 10 sp DGF TMP 2 2005 182 2005 212 16 sp DGF TMP 3 2006 127 2006 265 15 sp DGF TMP 4 2006 274 2007 3 14 sp DGF TMP 5 2007 244 2007 320 12 sp DGF TMP 6 2007 324 2008 5 12 sp DGF TMP 7 2007 146 2007 241 17 sp IPOC-CSN-GFZ IPOC + CSN 2007 318 2014 91 16 BB + 4 BB DGF TMP 8 2009 183 2009 365 23 sp DGF TMP 9 2011 207 2012 211 15 sp DGF TMP 10 2011 319 2011 348 15 sp DGF TMP 11 2012 183 2012 270 10 sp DGF PISAGUA 2013 2013 143 2013 280 26 sp IPOC-GFZ-CSN-ONEMI IPOC + CSN + CSNTemporary 2014 92 2014 240 46 BB + 12 BB GFZ-IPOC ILN 2012 91 2013 120 10 BB Seismic network . Start . End . . Operator . Name . Year . JD . Year . JD . Number and type of sensor . IRD-UCH-UA ANTOFAGASTA PERMANENT 1990 313 1995 209 10 sp DGF IQUIQUE 1991 1991 194 1991 227 21 sp DGF-GFZ-FUB CALAMA 1994 1994 88 1994 114 34 sp + 5 BB IRD-IPGS-DGF ARICA 1996 1996 173 1996 225 46 sp UTA ARICA PERMANENT 1997 91 2002 273 16 sp DGF TMP 1 2005 349 2006 38 10 sp DGF TMP 2 2005 182 2005 212 16 sp DGF TMP 3 2006 127 2006 265 15 sp DGF TMP 4 2006 274 2007 3 14 sp DGF TMP 5 2007 244 2007 320 12 sp DGF TMP 6 2007 324 2008 5 12 sp DGF TMP 7 2007 146 2007 241 17 sp IPOC-CSN-GFZ IPOC + CSN 2007 318 2014 91 16 BB + 4 BB DGF TMP 8 2009 183 2009 365 23 sp DGF TMP 9 2011 207 2012 211 15 sp DGF TMP 10 2011 319 2011 348 15 sp DGF TMP 11 2012 183 2012 270 10 sp DGF PISAGUA 2013 2013 143 2013 280 26 sp IPOC-GFZ-CSN-ONEMI IPOC + CSN + CSNTemporary 2014 92 2014 240 46 BB + 12 BB GFZ-IPOC ILN 2012 91 2013 120 10 BB Open in new tab Table 1. Summary of the deployments that recorded the data used in this investigation. Operator and network abbreviations refer to the Universidad de Chile (UCH), Universidad de Antofagasta (UA), Institut de Recherche pour le Developpement (IRD) the Institut de Physique du Globe de Strasbourg (IPGS), GeoForschungsZentrum (GFZ), Freie Universität Berlin (FUB), the Integrated Plate Boundary Observatory Chile (IPOC), the Chilean Centro Sismológico Nacional (CSN) seismic networks, the Oficina Nacional de Emergencia (ONEMI), the Departamento de Geofísica of the Universidad de Chile (DGF), and the Iquique Local Network (ILN-IPOC). Instrument types are short period (sp) and broad band (BB). Start and end times of deployment are in year and Julian day (JD). Data from IPOC (doi:10.14470/PK615318) is available online via the GFZ website. All other data available on request from DGF. Seismic network . Start . End . . Operator . Name . Year . JD . Year . JD . Number and type of sensor . IRD-UCH-UA ANTOFAGASTA PERMANENT 1990 313 1995 209 10 sp DGF IQUIQUE 1991 1991 194 1991 227 21 sp DGF-GFZ-FUB CALAMA 1994 1994 88 1994 114 34 sp + 5 BB IRD-IPGS-DGF ARICA 1996 1996 173 1996 225 46 sp UTA ARICA PERMANENT 1997 91 2002 273 16 sp DGF TMP 1 2005 349 2006 38 10 sp DGF TMP 2 2005 182 2005 212 16 sp DGF TMP 3 2006 127 2006 265 15 sp DGF TMP 4 2006 274 2007 3 14 sp DGF TMP 5 2007 244 2007 320 12 sp DGF TMP 6 2007 324 2008 5 12 sp DGF TMP 7 2007 146 2007 241 17 sp IPOC-CSN-GFZ IPOC + CSN 2007 318 2014 91 16 BB + 4 BB DGF TMP 8 2009 183 2009 365 23 sp DGF TMP 9 2011 207 2012 211 15 sp DGF TMP 10 2011 319 2011 348 15 sp DGF TMP 11 2012 183 2012 270 10 sp DGF PISAGUA 2013 2013 143 2013 280 26 sp IPOC-GFZ-CSN-ONEMI IPOC + CSN + CSNTemporary 2014 92 2014 240 46 BB + 12 BB GFZ-IPOC ILN 2012 91 2013 120 10 BB Seismic network . Start . End . . Operator . Name . Year . JD . Year . JD . Number and type of sensor . IRD-UCH-UA ANTOFAGASTA PERMANENT 1990 313 1995 209 10 sp DGF IQUIQUE 1991 1991 194 1991 227 21 sp DGF-GFZ-FUB CALAMA 1994 1994 88 1994 114 34 sp + 5 BB IRD-IPGS-DGF ARICA 1996 1996 173 1996 225 46 sp UTA ARICA PERMANENT 1997 91 2002 273 16 sp DGF TMP 1 2005 349 2006 38 10 sp DGF TMP 2 2005 182 2005 212 16 sp DGF TMP 3 2006 127 2006 265 15 sp DGF TMP 4 2006 274 2007 3 14 sp DGF TMP 5 2007 244 2007 320 12 sp DGF TMP 6 2007 324 2008 5 12 sp DGF TMP 7 2007 146 2007 241 17 sp IPOC-CSN-GFZ IPOC + CSN 2007 318 2014 91 16 BB + 4 BB DGF TMP 8 2009 183 2009 365 23 sp DGF TMP 9 2011 207 2012 211 15 sp DGF TMP 10 2011 319 2011 348 15 sp DGF TMP 11 2012 183 2012 270 10 sp DGF PISAGUA 2013 2013 143 2013 280 26 sp IPOC-GFZ-CSN-ONEMI IPOC + CSN + CSNTemporary 2014 92 2014 240 46 BB + 12 BB GFZ-IPOC ILN 2012 91 2013 120 10 BB Open in new tab The initial body wave data set derived from these deployments consists of P and S arrival times from 33 351 events recorded over a period of 25 yr by a variety of networks comprising 360 stations (Table 1; Figs 3 and 4). From this combined data set, we selected earthquakes that were recorded by at least 10 stations and with predicted arrival times within an initial outlier residual threshold of the larger of 2 s and 10 per cent of the total traveltime (essentially this means we presume the initial model is, on average, within about ±10 per cent of the actual structure). Application of these criteria resulted in a reference data set of 11 874 events with 110 640 P and 106 680 S wave arrival times. These arrival times and their associated uncertainties are estimated manually. Nominal uncertainties for P wave arrival times are between 0.1 and 0.5 s while those for S waves typically are about twice that amount. Figure 4. Open in new tabDownload slide Locations of earthquakes used in this study, indicated by circles filled with colours corresponding to depth ranges shown in the insert at the lower left. White squares joined by dashed lines and labelled P1–P16 locate the ends and midpoint origins (X = 0) in the cross sections shown in Figs 10–12. 2.2 Ambient noise data and preprocessing Most of the ambient noise analyzed in this study was recorded by 18 broad-band stations operated by the Integrated Plate Boundary Observatory Chile (IPOC), CSN and ONEMI over a period of 3 yr (2012 January 1–2014 December 31). This primary dataset is supplemented by recordings from 29 broad-band stations deployed by the DGF, CSN and ONEMI near Pisagua for 4 months (2014 March 29–2014 July 31) and 12 broad-band stations of the Iquique network operated by GFZ Potsdam for 7 months (2012 January 10–2013 September 5). The 59 total stations in these combined networks (Fig. 3) provide 1422 contemporaneous station pairs from which Green's functions may be estimated. For the most part we followed the pre-processing steps described in Bensen et al. (2007) to generate estimated Green's functions (EGFs) from the vertical component of cross-correlated noise, using a version of the CU-Boulder ANCC software,1 modified to permit sample rates and frequency bands appropriate for our data. Running-mean normalization was used to reduce contamination by coherent signals (mostly from earthquakes), and we pre-processed a band between periods of 2 and 150 s. Cross-correlations for each contemporary station-pair were generated for each day of record, and then stacked over the duration of co-recording. In an attempt to improve the signal-to-noise ratio (SNR) of the final stack we iteratively removed daily correlations that deviated by more than one standard deviation from the overall stack. However, we found that generally neither the shape of the EGF nor the SNR were significantly different from those obtained from a simple stack of the entire data set. We measured phase velocity dispersion using a technique based on that described in Yao et al. (2006). The EGFs are narrow-band filtered at 1 s intervals between 4 and 30 s in a pass-band of ±0.2 s of the central period and windowed within 3 cycles of an estimated arrival time. Noise is sampled in a part of the seismogram before the expected arrival, and seismograms with an SNR less than a given threshold (usually 10) are excluded. Those remaining are transformed to a velocity versus period space and a phase velocity dispersion curve is determined by tracking amplitude maxima over available periods. In practice, simply following the maxima does not guarantee immunity from cycle skips, and hence as a pre-processing step we generated a pilot dispersion curve from a least squares fit to all of the original curves that appear to fall within the correct cycle. In this case the pilot curve for phase velocity c (in km s−1) as a function of period T (in seconds) was represented by the polynomial \begin{eqnarray} c\left( T \right) &=& 3.142226\ + \ 0.03087T\nonumber\\ && - \ 0.0002410516\ {T^2} - 0.00000065{T^3}. \end{eqnarray} (1) For each EGF we choose the cycles that most closely follow this curve. Points without a clear maximum near the pilot curve were visually re-examined and discarded if ambiguous. Dispersion curves were generated for a band between 5 and 20 s period. Examination of dispersion curves from several station pairs suggests that cycle skips could be reliably identified at periods longer than 5 s; hence we adopt that as a minimum period. The choice of 20 s maximum period comes from abiding by the three-wavelength rule of Bensen et al. (2007) which in this case excludes nearly all but the N–S paths at longer periods. Because of the proximity of the stations to the Pacific Ocean (Fig. 3), and the strong dependence of path length on azimuth (most N–S paths are long, most E–W paths are short), one may be concerned that an azimuthally inhomogeneous distribution of noise could significantly bias the phase velocities determined from the EGFs. This concern was reinforced by both the asymmetry of some of the EGFs and discrepancies in dispersion curves generated by the causal and acausal sides of some of the cross correlations. To estimate and correct for this bias we employed a technique based on that described by Yao & van der Hilst (2009). To summarize their approach, we first determine the azimuthal distribution of noise E(θ) by considering that the observed cross correlation Cn(ω, t, dn) for station pair n as a function of frequency ω, time t, and interstation distance dn can be constructed as a weighted sum of monochromatic plane waves. Taking the Fourier transform (FT) of this construction gives \begin{eqnarray} &&{FT\left[ {{C_n}\left( {\omega ,t,{d_n}} \right){W_n}\left( {t,{d_n}} \right)} \right] = \sum\nolimits_{m = 1}^M E\left( {{\theta _m}} \right)}\nonumber\\ &&\times FT\{ {\cos[ {\omega ( {t - \delta {t_{nm}}} )} ]{H_{nm}}( {t,\delta {t_{nm}}} ){W_n}( {t,{d_n}} )} \} , \end{eqnarray} (2) where θm is the mth azimuth and δtnm is the phase delay time, calculated by integrating the inverse of the phase velocity along the path between the two stations. Wn and Hnm are time windows for the surface wave train in the EGF and the obliquely arriving monochromatic waves, respectively, assuming a range of wave speeds from 2 to 5 km s−1. The distribution of noise E(θ) is then determined by minimizing the energy in the difference between the observed and calculated real and imaginary parts of the Fourier transform of the cross correlations for all station pairs. We slightly modify this approach by noting that, for any pair of stations, the cross correlation at each station depends only on that part of E within the 180° of azimuth that travels towards that station. Hence, the sum in eq. (2) is a concatenation of two independent equations relevant to the causal and acausal parts of the cross correlation. Separating the two is advantageous for many station pairs where one side of the cross correlation has little energy above the noise level as the other side can still be used. Hence, we split the single eq. (2) into two sums over the relevant azimuths and add a row to the normal equations for each side of the cross correlation. Having estimated E(θ), we use (1) to calculate a phase shift δϕn between an observed cross-correlation and one which would have resulted from an isotropic distribution of noise (i.e. E(θ) is a constant) and define a bias μ as \begin{equation} \mu = - \frac{{\delta {\phi _n}}}{{\omega {t_n}}} = \frac{{t_n^{iso} - {t_n}}}{{{t_n}}}, \end{equation} (3) where tn is an estimate of the phase delay time between station-pair n derived from a 2-D phase velocity map using dispersion curves generated from averaged cross correlations (causal and acausal). The estimates of the ‘true’ phase delay time, tiso, are then used to update the phase velocity maps. This process may be iterated until the estimates of tiso stabilize. For the data set used here, stability was achieved after a single iteration. As may be expected, the noise distribution we determined for northern Chile (Fig. 5) is dominated by sources from the west, although the direction of propagation for the maximum is slightly north of east (about 70° azimuth). The biases resulting from this noise distribution are variable but generally small (<0.6 per cent) (Fig. 6). These estimated biases were used to adjust the surface wave phase delay times and new phase velocity maps (Fig. 7) were computed. The phase delay times derived from these maps are then used in the joint inversion with the body wave arrival times. Figure 5. Open in new tabDownload slide Normalized noise distribution E as a function of propagation azimuth for a selection of periods (5–9 s). Note that most of the noise propagates to the ENE, as may be expected from the proximity of the recording stations to the Pacific coast of South America. Figure 6. Open in new tabDownload slide Phase velocity bias determined from the noise distributions shown in Fig. 5 as a function of period. Despite the heterogeneity in the noise distribution, the expected bias is small (<0.6 per cent) for all station pairs. Figure 7. Open in new tabDownload slide Phase velocity maps at selected periods (indicated at the lower right of each panel) determined from ambient noise generated Rayleigh waves. The coastline of northern Chile is plotted as a reference. Values are percent differences from an average 1-D background, and the contour interval is 2 per cent. 2.3 Joint inversion methodology Like most joint inversions, the simultaneous fitting of body wave arrival times and surface wave phase delay times is motivated by complementary sensitivities of these observations to different parts of the model. In this case, surface waves are included mostly to reduce potential trade-offs in crustal and mantle structure caused by the lack of local earthquake activity in the crust. Surface waves also provide additional constraints on shear wave speeds that compensate for the lower number of, and higher uncertainties in, shear wave arrival times. The joint inversion methodology is based on an approach described in Nunn et al. (2014), modified for the local event case using procedures described in Roecker et al. (2004, 2006). The subsurface is parameterized by specifying P and S wave speeds on a 3-D grid of nodes spaced at 5 km intervals in depth and latitude and about the same distance in longitude (increments in longitude are everywhere 0.04717°). Intragrid wave speeds are determined by trilinear interpolation. Body wave traveltimes within the medium are calculated using a 3-D eikonal equation solver in a spherical (Earth-centred) coordinate system (Li et al.2009; Zhang et al.2012). Surface wave phase delay times at a given frequency ω in a given 3-D Vs model are determined by first assuming that such a model can be constructed by combining 1-D models at each areal grid point (Montagner 1986). With this assumption, we calculate phase velocities c and partial derivatives ∂c/∂Vs for the 1-D model at each areal point using the locked-mode method of Gomberg & Masters (1988). The phase delay times between any two points are then calculated by integrating the reciprocal of the phase velocity along the great circle path between them. The inverse problem is a standard linearized expansion of the forward problem, and partial derivatives for each observation are calculated along the ray paths. Following a procedure used in Roecker et al. (2004) we solve for perturbations to P slowness (Up = 1/Vp) and the ratio r = Vp/Vs = Us/Up rather than to shear wave slowness itself by setting ΔUs = Δ(rUp) = UpΔr + rΔUp and including a shear wave observation using \begin{eqnarray} \ {T_{so}} - {T_{sc}} = \sum\nolimits_{k = 1}^4 {\frac{{\partial {T_{sc}}}}{{\partial {h_k}}}\Delta {h_k}} + \sum\nolimits_{i = 1}^m {\frac{{\partial {T_{sc}}}}{{\partial {U_{si}}}}\left[ {{U_{pi}}\Delta {r_i} + {r_i}\Delta {U_{pi}}} \right]} ,\!\!\!\!\nonumber\\ \end{eqnarray} (4) where Tso and Tsc are the observed and calculated shear wave traveltimes, respectively. The partial derivatives are taken with respect to the 4 hypocentre parameters hk and the m values of Usi that exhibit sensitivity to the shear wave time (i.e. those that bound the elementary volumes through which the ray passes). A revised shear wave speed model is then determined from Us = rUp after modifying r and Up with Δr and ΔUp. We take this approach partly for technical reasons—it allows the S wave model to take advantage of the greater resolving power of P waves and provides a means to couple the sensitivity to Us of the surface waves to perturbations in Up—but also because we ultimately seek to interpret variations in r as proxies for hydration and melt fraction. Determining r by dividing individually determined Vp and Vs is problematic because of differences in the level of resolution for these quantities. Phase delay times of surface waves at any frequency ω are included in a similar way by relating derivatives of shear wave slowness Us with respect to phase velocity c(ω) by the chain rule: \begin{eqnarray} {T_{so}}\left( \omega \right) - {T_{sc}}\left( \omega \right) &=& \sum\nolimits_{i,j} {\frac{{\delta T\left( \omega \right)}}{{\delta {c_{ij}}\left( \omega \right)}}} \sum\nolimits_k \frac{{\delta {c_{ij}}\left( \omega \right)}}{{\delta {U_{sijk}}}}\nonumber\\ &&\times \left[ {{U_{pijk}}\Delta {r_{ijk}} + {r_{ijk}}\Delta {U_{pijk}}} \right], \end{eqnarray} (5) where the indices i, j and k refer to the latitude, longitude and depth coordinates, respectively, of each node in the model. Each row of the system of linear equations is weighted to reflect the uncertainties in the observations. Body wave arrivals are weighted by the inverse of the square of the associated time uncertainty. Surface wave phase delay times are similarly weighted by estimating an uncertainty based on the SNR of the EGF used in generating the dispersion curve, multiplied by a factor related to the period (longer periods generally have less well defined maxima). Uncertainty estimates thus obtained range from a few tenths of a second to several seconds but generally are about 3 per cent of the phase delay time. Outliers in the surface wave data set are identified and removed during the creation of phase velocity maps as part of the noise bias analysis. For body waves, data quality criteria are enforced at each iteration in order to eliminate potential outliers and to disqualify less well constrained hypocentres. Thresholds for outlier identification are applied at each iteration, and observations that exceed these thresholds are disqualified for that iteration (but may be reintroduced in subsequent iterations if a wave speed adjustment brings the residual under the threshold). These thresholds are generous at the start but gradually reduced over several iterations to the larger of 0.5 seconds or 5 per cent of the traveltime. Hypocentres with fewer than 10 acceptable P and S wave observations after the criteria are applied are excluded. For most runs this culling process retained about 95 per cent of the original observations over the course of inversion. We also applied different strategies for relative weighting of body and surface wave observations in a joint inversion, but in the end found that the best result was obtained by first inverting the surface wave phase delay times alone for several iterations prior to jointly inverting both the body and surface waves with even weighting. This appears to be due to both the greater sensitivity of surface waves to shallow structure and the much larger number of body wave observations available. The resulting system of linear equations is solved iteratively using the LSQR algorithm of Paige & Saunders (1982). In addition to regularization by damping, perturbations to the P and S models are smoothed after each iteration by averaging neighbouring nodes with a moving window. After trial-and-error tests of trade-off between variance reduction and model complexity we decided to use a window with 5 nodes in the E–W direction (the centre node and two nodes either side), 11 nodes in the NS direction (±5 nodes from the centre) and 3 nodes (±1 node from the centre) in the Z direction. The larger window in the N–S direction reflects an anticipated lesser degree of heterogeneity parallel to the strike of the Andes, but also compensates for the longer ray paths in that direction. Note that at each iteration we are smoothing the perturbations to the model as opposed to the model itself; hence smaller scale variations in wave speed still can appear over several iterations. Hypocentre coordinates are included as variables in the inversion (eq. 3) but are also relocated in updated models prior to any iteration. The models shown here generally will have gone through about 10 iterations with only surface wave times followed by an additional 10 iterations with combined surface and body wave times. Generally, no significant change in variance nor in values of wave speeds are realized beyond that point. 2.4 Preliminary inversions and starting models As mentioned above, we determined that the best way to incorporate surface wave observations was to allow them to adjust the model prior to the joint inversion. In developing a starting model for the surface wave only (SWO) inversion, we began with a body wave only (BWO) inversion that started with an adaptation of the 1-D model of Husen et al. (1999) for the Antofagasta region (Fig. 8). The resulting 3-D BWO model for Vp and Vp/Vs, obtained after 13 iterations, was averaged laterally to obtain a 1-D estimate of Vs. This step ensures some compatibility with the body waves, and also provides a means for determining appropriate values for regularization and thresholds for data misfit. Because surface wave inversions tend to retain biases introduced by interfaces like the Moho, we smoothed the transition from crustal to mantle velocities over a range of depths (40–70 km) suggested by previous estimates of crustal thickness in this area (e.g. Beck et al.1996; Yuan et al.2000, 2002; McGlashan et al.2008). This 1-D Vs model (Fig. 8) was then used as a starting model for the SWO inversion. Figure 8. Open in new tabDownload slide 1-D starting models used in this study. Blue lines show Vp and Vs from Husen et al. (1999) used in the BWO inversion. Red lines show a model deduced from a 1-D average of the BWO inversion. Vs from this model is used as a starting model for the SWO inversion. After 10 iterations of the SWO inversion, the wave speeds in the resulting 3-D Vs model were then multiplied by a constant Vp/Vs ratio to construct a starting 3D Vp and Vs model for the joint inversion. Based on both the results from the BWO inversion and a linear fit to Tp versus Ts–Tp for the entire body wave data set (Fig. 9), we take Vp/Vs = 1.76 as the preferred value. As discussed below, while the absolute values for the wave speeds and their ratios in the final models are sensitive to the choice of initial value of Vp/Vs, there are metrics that argue in favour of 1.76. Figure 9. Open in new tabDownload slide A Wadati plot [P traveltime (Tp) versus S-P traveltime (Ts–Tp)] for all arrival times in the dataset. The red line is a least-squares fit corresponding to a Vp/Vs of 1.7622 ± 0.0002. 3 DESCRIPTION OF THE PREFERRED MODEL Several trial inversions, discussed below, were performed to evaluate the effects of a variety of assumptions and choices of parameters on the final model. Those chosen for our preferred model (Figs 10–12) manage to fit the observations well while suppressing short wavelength artefacts. Relative to the 3-D starting model generated by the SWO inversion described above, the preferred model reduces the variance in the residuals from 0.063 to 0.044 s2, or by about 30 per cent. The final variance is larger than the anticipated noise level variance of 0.024 s2, suggesting that the observations are not overfit. The body wave residual variance is slightly less than that generated by the original (and less regularized) BWO inversion (0.046 s2) suggesting that the inclusion of the surface wave observations reduces the space of models that fit the body wave arrival times equally well. At the same time, as discussed below, these models are generally similar. Figure 10. Open in new tabDownload slide (a) Cross-sections P1–P8 of Vp for the preferred model. Locations of sections are shown in Fig. 4, and correspond to the labels in the lower left-hand corner in each section. Depth is relative to sea level, and distance is with respect to the midpoints of each section (plotted as white squares in Fig. 4). Values are indicated by the colour palette in the upper left part of the figure, and are shown for thick contours in between the sections. Thin contours are at 0.2 km s–1 intervals. The plots are shaded to reflect density of sampling, with the brightest regions corresponding to regions samples by at least 10 ray paths. Cross-sections of hypocentres within ±15 km perpendicular to a vertical plane defined by the section are shown as small open circles. (b) Cross-sections P9–P16 of Vp for the preferred model. Meanings of symbols and colours are the same as in (a). Figure 10 Open in new tabDownload slide (Continued.) Figure 11. Open in new tabDownload slide (a) Cross-sections P1–P8 of Vs for the preferred model. Meanings of symbols and colours are the same as in Fig. 10(a), except that the thin contour interval is now 0.5 km s–1. (b) Cross-sections P9–P16 of Vs for the preferred model. Meanings of symbols and colours are the same as in Fig. 10(a), except that the thin contour interval is now 0.5 km s–1. Figure 11 Open in new tabDownload slide (Continued.) Figure 12. Open in new tabDownload slide (a) Cross-sections P1–P8 of Vp/Vs for the preferred model. Meanings of symbols and colours are the same as in Fig. 10(a), except that the thin contour interval is now 0.005 and the thick contour interval is 0.02. Thick contour values start at 1.72. (b) Cross-sections P9–P16 of Vp/Vs for the preferred model. Meanings of symbols and colours are the same as in Fig. 10(a), except that the thin contour interval is now 0.005 and the thick contour interval is 0.02. Thick contour values start at 1.72. Figure 12 Open in new tabDownload slide (Continued.) Because most of the trial runs we performed were designed to test the robustness of certain features of interest in the mantle, we mostly discuss wave speeds in the subducted slab and the supraslab mantle (a discussion of finer scale crustal structure will appear in a companion paper). We also focus on larger scale features that our trial inversions show to be robust, in the sense that they are not overly sensitive to subjective choices concerning how the inversion is carried out. One of the more conspicuous features of the preferred model is an ubiquitous positive gradient in Vp/Vs ratio, increasing from about 1.75 within the slab to 1.78–1.82 in the supraslab mantle (Fig. 12). This gradient is in most cases perpendicular to the trend defined by the seismic zone and, while it varies in magnitude along strike, extends along much of the length of the slab within the mantle. Specifically, the contrast in the northern sections (P1–P8 in Fig. 12) is about 3 per cent while that in the south (P9–P16) is closer to 2 per cent. The high Vp/Vs region also appears to extend over a significantly larger volume in the sections north of about 21°S. Vp/Vs within the slab itself is higher in the north (∼1.76) than in the south (∼1.74). In most sections, the supraslab region of high Vp/Vs extends from about 40 km depth to the end of the seismic zone at about 120 km depth, and may extend to depths as shallow as 20 km in some of the northern sections (sections P1-P3 in Fig. 12a). Note that while Vp/Vs is high in the supraslab region, both Vp and Vs are relatively low in the 40–80 km depth range (Figs 10 and 11). Vp and Vs both increase within the slab at depths of 80–100 km, and then remain high to the end of the seismic zone. Sections with a well-defined double seismic zone at depths greater than 90 km (P1–P4 in Fig. 12) are characterized by a relatively low value of Vp/Vs in the area between the two zones. Vp and Vs at these depths are both lower in the lower part of the double seismic zone than in the upper part (e.g. sections P1–P3 in Figs 10 and 11). 4 TESTS OF RESOLUTION AND ROBUSTNESS As no single type of trail inversion provides conclusive evidence about the quality of an image, a variety of tests with both real and synthetic data were performed to estimate the level of robustness of certain features in the preferred wave speed model. In order to account for the effects of observational uncertainty, each synthetic data set is modified by random noise using a normal distribution with standard deviations that reflect the expected uncertainties in the actual data. Our main concerns are the extent to which the data can resolve the features we wish to interpret, and the dependence of the model on subjective choices made while carrying out the inversion (e.g. degree of regularization, weighting of observations, or starting model assumed). 4.1 Recovery of prismatic anomalies A standard approach to estimating resolution in arrival time inversion is the ‘checkerboard’ test that attempts to recover alternating positive and negative perturbations in abutting prisms. We conducted versions of this test in which ±5 per cent perturbations relative to a background 1-D model are applied to individual rectangular prisms that are 20 × 20 km in area and 10 km in depth. The perturbed elements are separated from one another by at least one prism with no perturbation so that the potential for smearing may be observed. Results of these tests show that, in general, regions sampled with at least 10 ray paths are reasonably well recovered (Fig. 13). Moreover, despite recovered perturbations being smoothed at every iteration, the amount of smearing outside the anomaly is small over most of the region. In the test shown here (Fig. 5), Vp and Vs are both perturbed by ±5 per cent, with a zero Vp/Vs anomaly. The near identical recovery of Vp and Vs in this example shows that Vs can be resolved when solving for Vp/Vs rather than Vs directly, even when Vp/Vs itself is associated with a null signal. Figure 13. Open in new tabDownload slide Map view of results of the prism reconstruction test at 50 km depth for Vp (left-hand panel) and Vs (right-hand panel). Values are indicated by the colour palette in the upper left-hand part of the figure. The plots are shaded to reflect density of sampling, with the brightest regions corresponding to regions samples by at least 10 ray paths. The rectangles locate the bounds of the prisms that are perturbed from the background 1-D model. Note that the reconstructed images are nearly identical for Vp and Vs. 4.2 Reconstruction of the model While recovery of an isolated perturbation can be a fair representation of resolution of that variable, reconstruction of a prismatic pattern includes effects of coupling between perturbations that would be appropriate only if the Earth were indeed composed of alternatingly perturbed prisms. For example, such a representation requires 3-D ray paths to alternately focus and defocus around each perturbation. An alternative way to demonstrate the ability of a data set to recover interpretable features at the scale and amplitude as those that appear in the model is to reconstruct the same distribution of imaged anomalies (e.g. Prevot et al.1991). Hence, we generated a synthetic dataset using the same source–receiver distribution as the actual dataset and followed the same iterative procedure as with real data to attempt to recover the preferred model from a 1-D starting model. The results (Fig. 14) show that the recovered model generally looks very much like the true model both in pattern and amplitude, but there are some exceptions. First, while the recovered model manages to reproduce the absolute wave speeds rather well, it tends to slightly underestimate (by as much as a percent) the anomalies in the southern sections and overestimate them by a similar amount in the north. This is likely a result of the greater density of both earthquakes and stations in the north. Second, there is some lateral smearing of Vp in the northernmost sections at depths greater than about 90 km that prevents the appearance of the isolated low Vp/Vs in the true model. The reason for this lateral smearing is not obvious; indeed, one might expect that the relative abundance of seismicity in this part of the seismic zone would lead to better resolved features at depth. This smearing suggests that real earth structure in the north and at depths greater than 90 km is likely to be laterally smeared in the modelling of real data as well, and hence that the Vp and Vs anomalies are in reality smaller in lateral extent than they appear in the image. Figure 14. Open in new tabDownload slide Results of the reconstruction tests for Vp (top two panels), Vs (middle two panels) and Vp/Vs (bottom two panels). In each case the ‘True’ model, taken from the preferred model, is shown at the top, and the reconstructed model directly beneath. These sections correspond to P12 in Figs 4 and 10–12. 4.3 Thin slab test One of the main features we seek to interpret is the contrast between wave speeds in the slab and the supraslab mantle, and to test the robustness of this gradient we constructed hypothetical models by placing a thin (10 km thick) body with a relatively low Vp/Vs (1.67 versus 1.75) ratio representing that part of the slab just beneath the mantle seismicity. The purpose of this test is to investigate the extent to which (1) localized anomalies within the slab spread across steep gradients into surrounding regions and (2) intraslab anomalies may be imaged by body waves. Several bodies with different lateral and vertical dimensions where tested; the salient results can be illustrated by models reconstructed from a Vp/Vs anomaly located in a thin (15 km thick), curved zone located along the underside of the (upper) seismic zone from about 50 km depth to the base of the model (Figs 15 and 16). We examine the resolution of lateral variations in the slab by restricting the lateral extent of one of the slab models to about 1° of latitude. As with all tests, we mimic the body wave arrival time data set used in the real inversion, but as only structure below 60 km is perturbed we do not include surface waves in a joint inversion. In this case, Vp remains at its 1-D value while Vs is perturbed to give an appropriate Vp/Vs contrast. In addition to Vp/Vs, this allows us to examine how well Vs can be recovered from estimates of Vp and Vp/Vs. Figure 15. Open in new tabDownload slide Map views of Vp/Vs for the short (top panel) and long (bottom panel) hypothetical slab models shown at 100 km depth. The true model is at the right, and the reconstructed model is at the left. Figure 16. Open in new tabDownload slide Cross-section views reconstruction of Vp/Vs for the long (top panel) and short (bottom panel) slab models shown at 20°S. The true model is shown in the middle. Bright areas are those sampled by 10 or more rays. This section corresponds to P5 in Fig. 4. The Vp/Vs and Vs anomalies in the reconstruction (Figs 15 and 16) are both well recovered in both models at depths shallower than the deepest earthquakes (about 120 km depth) but, due to a lack of ray paths, no deeper. The thickness and location of the zone are well recovered, with very little contamination of the surrounding areas by the lower Vp/Vs in the slab, or of higher Vp/Vs from outside the slab. The upper termination of the slab at 50 km depth is also well constrained. Finally, the lateral variation in the extent of the slab is well reconstructed. At the same time, while the variations in Vp/Vs are evident in the reconstruction, the actual values are underestimated in much of the slab and show intraslab heterogeneity on the order of 0.026 along its strike and dip. 4.4 Smoothing, weighting and starting models Among the subjective assumptions made as part of the joint inversion are the extent to which the model is smoothed, the relative influence of different date types, and the starting model chosen. Several tests were performed varying each of these assumptions individually, and the general conclusions can be summarized in a few comparative tests. Specifically, the BWO inversions discussed above use a different set of data (no surface waves), a different starting model (the 1-D Antofagasta model of Husen et al.1999), and/or a shorter smoothing window (±1 node about the centre) than that used to generate the preferred model. A comparison of BWO models with different moving average window lengths (top and bottom panels in Fig. 17) shows that, as may be expected, the effects of smoothing are akin to passing the rough model through a low pass filter. All of the first order features that we see in the smoothed model can be identified in the rough one, but the fit to the data using the rough model is only marginally better (by about 5 per cent in variance). The difference between the BWO and joint models arises partly from the difference in starting model; in particular, the Moho transition is more gradual in the joint model than the BWO models. Moreover, the inclusion of surface waves helps constrain the wave speeds in the upper few tens of km. These choices account for differences in the supraslab mantle wave speeds determined by the BWO inversion. For example, the lower wave speeds in the upper mantle in the BWO model are likely compensating for crustal thickening towards the central Andes. Indeed, the added constraint on crust and uppermost mantle structure provided by surface waves is a principle motivation for a joint inversion. Note that while the rough BWO model manages to fit the body waves better than the joint model by about 2 per cent, the joint model does a slightly better job (by 3 per cent) than the equally smoothed BWO model. Of course, neither of the BWO models manages to fit the surface wave data nearly as well as the joint model (surface wave variance in the joint model is reduced by about 70 per cent compared to that in either of the BWO models). Figure 17. Open in new tabDownload slide Comparison of models generated with body waves only (BWO) with minimal smoothing (top row), the preferred model (middle row), and BWO with the same smoothing as the preferred model (bottom) row. Shown are Vp (left-hand column), Vs (middle column) and Vp/Vs (right-hand column). Bright areas correspond to regions sampled by at least 10 rays. Each cross section corresponds to P3 in Fig. 4. Creating a starting model for a joint inversion from the results of an SWO inversion requires a choice of Vp/Vs ratio. As discussed above, a Vp/Vs of 1.76 is consistent with the arrival times in the body wave data set. At the same time, tests with a range of values from 1.73 to 1.78 shows that while the patterns and gradients in Vp/Vs are virtually invariant, the average value of Vp/Vs remains close to its starting value. The fit to the entire data set is better with a Vp/Vs of 1.76 than with either 1.73 or 1.78 (by about 3 per cent), and is nearly the same as 1.75. Moreover, there are about 1100 earthquakes (or roughly 10 per cent of the entire dataset) distributed throughout the study area, for which the fit to the data improves on an event by event basis with a Vp/Vs of 1.76 (Fig. 18). Since events are relocated prior to each iteration, the minimization of residuals on an event by event basis is a valid metric for judging the significance of a given model, and on this basis, as well as the improvement in overall variance, we take the model generated with a starting value of 1.76 as the preferred model. Figure 18. Open in new tabDownload slide Histogram of the number of earthquakes with standard deviations within bins of 0.05 s for different choices of Vp/Vs in the starting model. Note that while results from Vp/Vs of 1.75 (blue) and 1.76 (red) are nearly the same, these choices reduce the standard deviation for more than 1100 events when compared to low (1.73) and high (1.78) values. The starting model refers to the 1-D SWO model shown in Fig. 8. 4.5 Inferences from trial results As a result of the trials described above, we conclude that the first order features in the preferred model are robust. Specifically, the increase in Vp/Vs, from the subducted slab into the supraslab mantle, as well as the size and location of the high Vp/Vs regions, appear well constrained. These features are reasonably independent of the starting model, but the actual values of Vp/Vs tend to stay close to the starting value. Nevertheless, we can appeal to a preference of the data for a certain range of Vp/Vs to argue for an optimal value. While large gradients in Vp, Vs and Vp/Vs are well constrained, localized regions with small changes (on the order of 1 per cent in velocity and 0.017 in Vp/Vs) are probably not significant. We are less confident in features in the model at depths greater than about 110 km. Although tests that recover specific features such as a low Vp/Vs in the region between the deeper parts of the double seismic zone were successful, other tests suggest lateral smearing may blur the image. 5 DISCUSSION 5.1 Hydration of the supraslab mantle The most conspicuous first order features in our preferred model are (1) the gradient from low to high Vp/Vs from the slab into the supraslab region, mostly at depths greater than about 40 km, along the entirety of the northern Chile margin, and (2) the extensive region, particularly north of about 21°S, of Vp/Vs > 1.8 located mostly above the interface between the slab and the supraslab mantle as represented by the Wadati-Benioff Zone (WBZ) earthquakes. This high Vp/Vs region is also characterized by relatively low values of Vp and Vs to depths of about 80 km. Because the introduction of fluids into the mantle can decrease both Vp and Vs and increase Vp/Vs by either partial melting (e.g. Hammond & Humphreys 2000; Karato 2012) or serpentinization (e.g. Christensen 1996, 2004; Ji et al.2013), a viable explanation of these features is the dehydration of the subducted slab and the release of volatiles into the supraslab mantle. Variations on this scenario have been suggested in previous studies of subduction zones (e.g. Kamiya & Kobayashi 2000; Bostock et al.2002; DeShon & Schwartz 2004; Eberhart-Phillips & Bannister 2010; Zhang et al.2010) including northern Chile (e.g. ANCORP Working Group 2003). Our results suggest that the way in which the mantle is hydrated can depend not only on the type of convergent margin, but on variations in the subduction environment within a particular margin type. Locating the supraslab mantle along the Andean margin requires an estimation of both the continental Moho and the top of the subducting Nazca Plate. The latter is assumed to coincide with the upper boundary of the WBZ earthquakes. The depth of the continental Moho is known to some extent from previous studies in this area (Beck et al.1996; Yuan et al.2000, 2002; McGlashan et al.2008) that suggest that the thickness of the continental crust increases from about 40 km beneath the Coastal Cordillera and Central Depression to about 60 km under the Western Cordillera. Hence, the supraslab mantle ‘wedge’ at depths of 40–60 km is more like a channel with a width of perhaps a few 10 s of km. Most sections of our model (Figs 10 and 11) show a prominent east-dipping gradient in Vp and Vs between 40 and 60 km depth about 25 km above the top of the WBZ that we interpret as the crust to mantle transition. Vp and Vs are both low (∼7.2–7.6 km s−1 and 4.1–4.3 km s−1, respectively), and Vp/Vs is generally high (1.78–1.80), within the supramantle at depths between about 40 and 80 km depth, with the more extreme values found in the northern sections. These trends are consistent with serpentinization of this part of the mantle. Thermomechanical modelling of the subducting Nazca Plate (Dorbath et al.2008) suggest that temperatures are cold enough at these depths (<700 °C) for serpentinite to be stable (Ulmer & Trommsdorff 1995). In most cross sections, a gradient in Vp, Vs and Vp/Vs collocates with the upper part of the WBZ, and correlates quite well with the ‘Nazca Reflector’ seen in both receiver function (Yuan et al.2000) and active source (ANCORP Working Group 2003) imaging. Note that while there is a positive gradient in Vp and Vs from the supraslab mantle into the interior of the slab, Vp and Vs within the slab are low. In the northern sections, Vp/Vs within this part of the slab is high, suggesting the presence of serpentinite within this part of the slab. These contrasts are greater north of 21°S, suggesting more significant hydration of the subducted slab in the north. While estimates derived from the existing literature disagree on the extent to which this part of the mantle is serpentinized, the reduced Vp and Vs combined with increased Vp/Vs is consistent with higher levels of serpentinization in the north. For example, the high P-T model of Ji et al. (2013) would predict ∼45 per cent serpentinite in the south and ∼70 per cent in the north, while predictions from the dunite-antigorite model of Christensen (2004) are closer to 30 and 42 per cent (Fig. 19). Figure 19. Open in new tabDownload slide Values of Vp, Vs and Vp/Vs taken from the mantle wedge region of the preferred model for regions north (red symbols) and south (blue symbols) of 21°S, plotted on the Dunite-Antigorite model of Christensen (2004) (squares) and the high PT measurements of Ji et al. (2013) (circles). Note that while the estimates of serpentine per cent vary significantly, both correspond to a higher percentage north of 21°S. Vp and Vs both increase in the WBZ between 80 and 100 km and reach 8.2 and 4.6 km s–1, respectively. Vp/Vs remains high below 80 km, and in fact is higher in the 80–120 km depth range than anywhere else in the model. For the most part the highest Vp/Vs is in a narrow zone that runs more or less parallel to the WBZ. As with the shallower parts of the model, Vp/Vs is in general much higher in the north (1.82) than in the south (1.78). Combined with the normal to high Vp and Vs in the supraslab mantle below 80 km depth, plus the instability of serpentinite at higher temperatures in this part of the mantle, we infer that this part of the supraslab mantle is partially melted. Based on the model of Hammond & Humphreys (2000), the observed increase in Vp/Vs could be achieved with a melt fraction of about 1 per cent in the north and less than half of that in the south. We interpret the large, downdip positive gradient in Vp and Vs within the slab at depths between about 80 and 100 km to be a consequence of accelerated phase transitions from basalt to eclogite in the crust and from serpentinite to peridotite in the mantle of the slab. The region of high Vp/Vs extends nearly the length of seismic zone but returns to average values at about 120 km depth, suggesting that phase transitions within the slab have ceased at this depth. The increased level of seismic activity between 100 and 130 km depth is indicative of higher rates of phase transition. Therefore, as similarly suggested by Yuan et al. (2000), the lessening of wave speed contrasts between the WBZ and the supraslab mantle at these depths could signify the completion of phase transitions in the slab. 5.2 The north–south contrast and surface geology As discussed above, while the influence of water in the slab and supraslab mantle can be seen all along the northern Chile subduction zone, there are some significant differences between the northern and southern parts of the model, mostly having to do with the size of the anomalies and the volumes over which they occur. These differences are gradational, but the largest gradients occur at about 21°S. The most significant of these is the extent of the region of high Vp/Vs, which, as discussed above, we associate with hydration resulting either in serpentinization of peridotite or the presence of partial melt. The difference in the size of this region (Fig. 21) north and south of 21°S suggests a more profound influence of water in the north. We note that the oldest (>50 Ma) oceanic lithosphere subducted along the Nazca-South American plates boundary is located along this part of northern Chile (Fig. 2). Because the effective elastic thickness of oceanic lithosphere and the depth to the brittle‐ductile transition increases with age, older oceanic lithosphere has a higher water reservoir capacity (Watts 2001; Contreras-Reyes & Osses 2010; Contreras-Reyes et al.2011). Hence, the effects of dehydration are likely to be more apparent here than anywhere else along the Andean margin. The northern region is much more seismically active in the 70–145 km depth range (prior to the 2104 Pisagua earthquake, the rate of activity was about twice as high in the north), consistent with earthquake activity being related to rates of dehydration reactions (Kirby et al.1996; Omori et al.2002, 2004). Moreover, there is a shift in the locations of the deeper hypocentres to the east south of 21°S (Fig. 4), meaning that the dip of the slab is less in the south. An explanation consistent with dehydration is that the greater loss of volatiles in the north leads to increased production of eclogite within the subducted slab (e.g. Hacker et al.2003; Green et al.2010), which in turn increases its negative buoyancy. There is also a north–south difference in patterns of volcanic activity. In particular, the northern part is associated with an anomalous volume of Neogene ignimbrite deposits (Wörner et al.1994; Wörner et al.2002), which is consistent with a larger percentage of volatiles in the ascending magma. Interestingly, the southern boundary of the large hydrated region has an extended region of low Vp/Vs near the surface that corresponds to the location of the Pica gap in volcanism (Wörner et al.1994; Figs 1 and 2). The northern part also has prominent linear topographic scarps that deform Neogene to Pleistocene deposits along the forearc, including portions that overlie the interplate zone and the Precordillera piedmont. While these are dramatic topographic features, their activation mechanism remains poorly understood (i.e. Armijo & Thiele 1990; González et al.2003; Allmendinger et al.2005; Farías et al.2005; Loveless et al.2005; González et al.2008). The inferred regions of enhanced hydration correlate with the occurrence of megathrust earthquakes along this part of the margin. The high Vp/Vs at depths below 40 km is consistent with the downdip extent of large megathrust events being limited by sepertinization (e.g. Hyndman & Peacock 2003). High Vp/Vs also extends up to 20 km depth between 19° and 19.8°S (sections P1–P3 in Fig. 12), suggesting the presence of significant amounts of water at shallower depths along the plate boundary. The importance of fluids for interplate coupling and hence for rupture mechanics is well known (e.g. Tsuji et al.2008; Nakajima et al.2009; Uchida et al.2009; Abers et al.2013), and appears here in a seismic gap capable of generating a Mw ∼9.0 megathrust earthquake. Recent geodetic studies (Métois et al.2013; Ortega-Culaciati et al.2015, Fig. 2) show a heterogeneous distribution of interseismic coupling within this gap, with the north (18°–20.5°S) being significantly less coupled than the south (20.5°–23°S). This less coupled region correlates well with the shallow distribution of high Vp/Vs (Fig. 20), suggesting that hydration of the interface along the interplate contact controls the level of coupling. Figure 20. Open in new tabDownload slide Tectonic context of the region of high Vp/Vs (>1.78) as shown by a blue surface that envelopes the volume. Most of this region is located in the supraslab mantle north of 21°S, and may extend close to the surface in the region where there 2014 Pisagua earthquake occurred. Seismicity is plotted as small circles filled with colours representing depth range as shown in Fig. 4. Active volcanos are shown by red circles. The dashed red line is the Pica gap. The volume within the blue surface represents a region of enhanced hydration leading to supraslab hydration and partial melting, particularly north of 21°S. These contrasts in the effects of hydration should be related to variations in the level of hydration of the Nazca plate. While it is difficult to know the initial conditions of the subducted plate, there presently is a succession of seamounts located on a larger swell on the Nazca plate, referred to as the Iquique ridge, that intersects the trench at about 21°S (Contreras-Reyes & Carrizo 2011; Figs 1 and 2). An extrapolation along the azimuth of this ridge into the subducted slab corresponds to the large region of high Vp/Vs that we see in the north. Because seamounts should be accompanied by enhanced hydrothermal circulation in the lithosphere, the subduction of the Iquique ridge could account for what we interpret as the elevated hydration of the supraslab mantle in this region. Note that this in turn suggests an additional role for subducted seamounts in the generation of megathrust earthquakes beyond the expected roughening of the plate boundary (e.g. Cloos & Shreve 1996; Scholz & Small 1997; Contreras-Reyes & Carrizo 2011; Geersen et al.2015). 5.3 Comparison with other studies in northern Chile A number of local seismic investigations of the crust and upper mantle beneath this part of the Andean margin have been published (e.g. Comte et al.1994; Schmitz & Kley 1997; Graeber & Asch 1999; Patzwahl et al.1999; Schmitz et al.1999; Yuan et al.2000, 2002; ANCORP Working Group 2003; Schurr et al.2003; Rietbrock & Waldhauser 2004; Koulakov et al.2006; Dorbath et al.2008; Bloch et al.2014). We note that the tomographic study of Graeber & Asch (1999) of the PISCO network deployed between 21°S and 25°S shows (1) a similar contrast of Vp/Vs from about 1.70 in the slab to 1.80 in the supraslab mantle, (2) extensive regions of low Vp/Vs (<1.70) in the crust, consistent with a largely felsic (i.e. granitic) lithology (Christensen 1999) and (3) a large increase in Vp from about 7.5 to >8.5 km s−1 within the slab in the 80–100 km depth range. A similar change in Vp within the slab was determined as well by Comte et al. (1994). As discussed above, the Nazca Reflector (NR) imaged by the Andean Continental Research Project (ANCORP) EW seismic reflection profile at 21°S (ANCORP Working Group 2003) aligns with the gradient in Vp and Vs found in our model. The high amplitude reflections of the NR terminate abruptly at 80 km depth, corresponding in our model to the increase in Vp and Vs to levels found in the surrounding mantle. Bloch et al. (2014) interpret the NR as an indicator of enhanced hydration, and suggest that a region of increased reflectivity above the NR to about 20 km depth corresponds to pathways for fluids within the crust and supraslab mantle, similar to those conjectured in previous studies (Springer 1999; Schurr et al.2003; Koulakov et al.2006). We note that this shallower region of enhanced reflectivity corresponds to a region of high Vp/Vs seen in sections P7 and P8 (Fig. 12a) that are closest to the ANCORP line. The tomographic study of Dorbath et al. (2008) overlaps the northern section of our region (∼18°–20°S). Their estimates of Vp and Vs in the upper WBZ show the same increases from low to high values at about 80 km depth as our model. In fact, Vp in their section P3 looks quite similar to section P3 in Fig. 10(a), including the decrease in Vp in the lower seismic zone at depths >90 km. At the same time, their estimates of Vs in the lower seismic zone, generally >5 km s−1, are significantly higher than our estimate. (∼4.4 km s−1). As a result, they deduce (by division) a low Vp/Vs (∼1.68) in this part of the double seismic zone, in contrast to our value of 1.76. Our model does show a region of low Vp/Vs, but it appears in the region between the two zones. Dorbath et al. (2008) also do not find any coherent increase in Vp/Vs in the supraslab region, which is a dominant feature in our model. While the causes are not easily discerned, these differences in results are most likely due to differences in data and methodology. 6 CONCLUSIONS Based on the results of this investigation and a comparison with results from previous studies in this and other subduction zones, we infer that the main features we observe are a consequence of dehydration of the subducting slab, phase transitions from basalt to eclogite and serpentinite to peridotite within the slab, and the hydration of the supraslab mantle. As summarized in Fig. 21, fluids from the subducting Nazca plate serpentinize the supraslab mantle between 40 and 80 km depth. At depths between 80 and 100 km, the phase transformations within the slab accelerate, resulting in higher wave speeds and more intense seismic activity. Water released into the supraslab mantle at these depths induces partial melt which makes its way up to the active arc in the Western Cordillera. At the end of the seismic zone at 140 km depth, the slab is dehydrated and the phase transformations are complete. Figure 21. Open in new tabDownload slide Summary of the principle inferences of this study as a simplified cartoon (top panel) and plotted on section P3 for Vp and Vp/Vs (bottom panel). Low Vp and high Vp/Vs within the supraslab mantle at depths between about 40 and 80 km correspond to serpentinization of the mantle wedge. The gradient to high Vp with high Vp/Vs between about 80 and 120 km transitions to partial melt that ascends (symbolized by the large arrow) to the volcanic arc in the Western Cordillera. Within the subducting Nazca Plate, a basaltic crust and serpentinized mantle at shallower depths transitions to an ecolgitic crust and peridotitic mantle at greater depths. We note that this scenario is similar in many respects to that proposed by Oncken et al. (2003) for the ANCORP results at 21°S from an entirely different set of observations. The current study reveals a significant contrast in the effects of hydration north and south of this latitude. We attribute the more intense effects in the north to a gradient in the level of hydration in the subducted Nazca plate that eventually results in a release of larger volumes of water into the supraslab mantle. These larger quantities of water also appear to influence the nature of volcanism in the magmatic arc and the coupling of the plate boundary at depths of 20–40 km. We suggest that the cause of this greater influx of water is the subduction of the Iquique ridge on the Nazca Plate. 1 Available at: http://ciei.colorado.edu/Products/. The ambient noise processing was greatly facilitated by the availability of the CU Boulder software and we especially thank Mike Ritzwoller for helpful advice at various stages of the processing. We also acknowledge GFZ/IPOC (doi:10.14470/PK615318) for use of their high quality data and as well the efforts of the many people who collected the data from the many other networks used in this study. We thank Frederik Tilmann and Ingo Grevemeyer for thoughtful reviews of this manuscript. This study was funded by FONDECYT project 1130071. The geodetic portion of this project was done in collaboration with FONDECYT project 11140904. REFERENCES Abers G.A. Nakajima J. van Keken P.E. Kita S. Hacker B.R. 2013 Thermal-petrological controls on the location of earthquakes within subducting plates Earth planet. Sci. Lett. 369 178 187 Google Scholar Crossref Search ADS WorldCat Allmendinger R.W. González G. 2010 Invited review paper: Neogene to Quaternary tectonics of the coastal Cordillera, northern Chile Tectonophysics 495 1–2 93 110 Google Scholar Crossref Search ADS WorldCat Allmendinger R.W. Gubbels T. 1996 Pure and simple shear plateau uplift, Altiplano-Puna, Argentina and Bolivia Tectonophysics 259 1 13 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Allmendinger R.W. Jordan T.E. Kay S. Isacks B.L. 1997 The evolution of the Altiplano-Puna Plateau of the central Andes Annu. Rev. Earth planet. Sci. 25 139 174 Google Scholar Crossref Search ADS WorldCat Allmendinger R.W. González G. Yu J. Hoke G. Isacks B. 2005 Trench-parallel shortening in the northern Chilean Forearc: tectonic and climatic implications Geol. Soc. Am. Bull. 117 89 104 Google Scholar Crossref Search ADS WorldCat Altamimi Z. 2002 ITRF2000: a new release of the International Terrestrial Reference Frame for earth science applications J. geophys. Res. 107 B10 2214 doi:10.1029/2001JB000561 Google Scholar Crossref Search ADS WorldCat ANCORP Working Group 2003 Seismic imaging of a convergent continental margin and plateau in the central Andes (Andean Continental Research Project 1996 (ANCORP’96)) J. geophys. Res. 108 B7 2328 doi:10.1029/2002JB001771 OpenURL Placeholder Text WorldCat Armijo R. Thiele R. 1990 Active faulting in northern Chile: ramp stacking and lateral decoupling along a subduction plate boundary? Earth planet. Sci. Lett 98 40 61 Google Scholar Crossref Search ADS WorldCat Baker M.C.W. Francis P.W. 1978 Upper Cenozoic volcanism in the central Andes: ages and volumes Earth planet. Sci. Lett. 41 175 187 Google Scholar Crossref Search ADS WorldCat Beck S. Silver P. Drake L. Zandt G. Myers S. Wallace T. 1996 Crustal-thickness variations in the central andes Geology 24 407 410 Google Scholar Crossref Search ADS WorldCat Béjar-Pizarro M. et al. 2010 Asperities and barriers on the seismogenic zone in North Chile: state-of-the-art after the 2007 Mw 7.7 Tocopilla earthquake inferred by GPS and InSAR data Geophys. J. Int. 183 390 406 Google Scholar Crossref Search ADS WorldCat Bensen G.M. Ritzwoller M. Barmin M. Levshin A. Lin F. Moschetti M. Shapiro N. Yang Y. 2007 Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements Geophys. J. Int. 169 1239 1260 Google Scholar Crossref Search ADS WorldCat Bevis M.E. Kendrick R. Smalley B. Brooks R. Allmendinger R. Isacks B. 2001 On the strength of interplate coupling and the rate of backarc convergence in the Central Andes: an analysis of the interseismic velocity field Geochem. Geophys. Geosyst. 2 doi:10.1029/2001GC000198 Google Scholar OpenURL Placeholder Text WorldCat Bloch W. Kummerow J. Salazar P. Wigger P. Shapiro S.A. 2014 High-resolution image of the North Chilean subduction zone: seismicity, reflectivity and fluids Geophys. J. Int. 197 3 1744 1749 Google Scholar Crossref Search ADS WorldCat Bostock M.G. Hyndman R.D. Rondenay S. Peacock S.M. 2002 An inverted continental Moho and serpentinization of the forearc mantle Nature 417 536 538 Google Scholar Crossref Search ADS PubMed WorldCat Carlson R.L. Miller D.J. 2003 Mantle wedge water contents estimated from seismic velocities in partially serpentinized peridotites Geophys. Res. Lett. 30 5 1250 doi:10.1029/2002GL016600 Google Scholar Crossref Search ADS WorldCat Chlieh M. de Chabalier J.B. Ruegg J.C. Armijo R. Dmowska R. Campos J. Feigl K.L. 2004 Crustal deformation and fault slip during the seismic cycle in the northern Chile subduction zone, from GPS and InSAR observations Geophys. J. Int. 158 695 711 Google Scholar Crossref Search ADS WorldCat Christensen N.I. 1996 Poisson's ratio and crustal seismology J. geophys. Res. 101 3139 3156 Google Scholar Crossref Search ADS WorldCat Christensen N.I. 1999 Physical and chemical properties of South Island, New Zealand rocks, Report Madison Univ. of Wis. 60 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Christensen N.I. 2004 Serpentinites, peridotites, and seismology Int. Geol. Rev. 46 9 795 816 Google Scholar Crossref Search ADS WorldCat Cloos M. Shreve R.L. 1996 Shear-zone thickness and the seismicity of Chilean-and Marianas-type subduction zones Geology 24 2 107 110 Google Scholar Crossref Search ADS WorldCat Coira B. Davidson J. Mpodozis C. Ramos V. 1982 Tectonic and magmatic evolution of the Andes of northern Argentina and Chile Earth Sci. Rev. 18 303 332 Google Scholar Crossref Search ADS WorldCat Comte D. Pardo M. 1991 Reappraisal of great historical earthquakes in northern Chile and southern Peru seismic gaps Nat. Hazards 4 23 44 Google Scholar Crossref Search ADS WorldCat Comte D. Roecker S. Suárez G. 1994 Velocity structure in northern Chile: evidence of subducted oceanic crust in the Nazca Plate Geophys. J. Int. 117 625 639 Google Scholar Crossref Search ADS WorldCat Contreras-Reyes E. Osses A. 2010 Lithospheric flexure modeling seaward of the Chile trench: implications for oceanic plate weakening in the Trench Outer Rise region Geophys. J. Int. 182 1 97 112 Google Scholar OpenURL Placeholder Text WorldCat Contreras-Reyes E. Carrizo D. 2011 Control of high oceanic features and subduction channel on earthquake ruptures along the Chile–Peru subduction zone Phys. Earth planet. Inter. 186 1 49 58 Google Scholar Crossref Search ADS WorldCat Contreras-Reyes E. Grevemeyer I. Watts A. Flueh E. Peirce C. Moeller S. Papenberg C. 2011 Deep seismic structure of the Tonga subduction zone: implications for mantle hydration, tectonic erosion, and arc magmatism J. geophys. Res. 116 doi:10.1029/2011JB008434 Google Scholar OpenURL Placeholder Text WorldCat DeMets C. Gordon R.G. Argus D.F. Stein S. 1990 Current plate motions Geophy. J. Int. 101 425 478 Google Scholar Crossref Search ADS WorldCat DeMets C. Gordon R.G. Argus D.F. Stein S. 1994 Effect of recent revisions to the geomagnetic reversal timescale on estimates of current plate motions Geophys. Res. Lett. 21 2191 2194 Google Scholar Crossref Search ADS WorldCat DeShon H.R. Schwartz S.Y. 2004 Evidence for serpentinization of the forearc mantle wedge along the Nicoya Peninsula, Costa Rica Geophys. Res. Lett. 31 L21611 doi:10.1029/2004GL021179 Google Scholar Crossref Search ADS WorldCat Dorbath C. Gerbault M. Carlier G. Guiraud M. 2008 Double seismic zone of the Nazca plate in northern Chile: high-resolution velocity structure, petrological implications, and thermomechanical modeling Geochem. Geophys. Geosyst. 9 Q07006 doi:10.1029/2008GC002020 Google Scholar Crossref Search ADS WorldCat Eberhart-Phillips D. Bannister S. 2010 3-D imaging of Marlborough, New Zealand, subducted plate and strike-slip fault systems Geophys. J. Int. 182 73 96 Google Scholar OpenURL Placeholder Text WorldCat Farías M. Charrier R. Comte D. Martinod J. Hérail G. 2005 Late Cenozoic deformation and uplift of the western flank of the Altiplano: evidence from the depositional, tectonic, and geomorphologic evolution and shallow seismic activity (northern Chile at 19°30′ S) Tectonics 24 TC4001 doi:10.1029/2004TC001667 Google Scholar Crossref Search ADS WorldCat Fyfe W.S. McBirney A.R. 1975 Subduction and the structure of andesitic volcanic belts Am. J. Sci. A275 285 297 Google Scholar OpenURL Placeholder Text WorldCat Geersen J. Ranero C. Barckhausen U. Reichert C. 2015 Subducting seamounts control interplate coupling and seismic rupture in the 2014 Iquique earthquake area Nat. Commun. 6 8267 doi:10.1038/ncomms9267 Google Scholar Crossref Search ADS PubMed WorldCat Gill J. 1981 Orogenic Andesites and Plate Tectonics Springer Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gomberg J.S. Masters T.G. 1988 Waveform modelling using locked-mode synthetic and differential seismograms: application to determination of the structure of Mexico Geophys. J. Int. 94 193 218 Google Scholar Crossref Search ADS WorldCat González G. Cembrano J. Carrizo D. Macci A. Schneider H. 2003 The link between forearc tectonics and Pliocene-Quaternary deformation of the Coastal Cordillera, northern Chile J. South Am. Earth Sci. 16 321 342 Google Scholar Crossref Search ADS WorldCat González G. Gerbault M. Martinod J. Cembrano J. Carrizo D. Allmendinger R. Espina J. 2008 Crack formation on top of propagating reverse faults of the Chuculay Fault System northern Chile: insights from field data and numerical modelling J. Struct. Geol. 30 6 791 808 Google Scholar Crossref Search ADS WorldCat Graeber F. Asch G. 1999 Three dimensional models of P wave velocity and P-to-S velocity ratio in the southern central Andes by simultaneous inversion of local earthquake data J. geophys. Res. 104 20 237–20 256 Google Scholar Crossref Search ADS WorldCat Green H.W. Chen W.P. Brudzinski M.R. 2010 Seismic evidence of negligible water carried below 400-km depth in subducting lithosphere Nature 467 828 831 Google Scholar Crossref Search ADS PubMed WorldCat Grove T.L. Chatterjee N. Parman S.W. Médard E. 2006 The influence of H2O on mantle wedge melting Earth planet. Sci. Lett. 249 74 89 Google Scholar Crossref Search ADS WorldCat Grove T.L. Till C. Lev E. Chatterjee N. Medard E. 2009 Kinematic variables and water transport control the formation and location of arc volcanoes Nature 459 694 697 erratum 460, 1044 Google Scholar Crossref Search ADS PubMed WorldCat Hacker B.R. 1996 Eclogite formation and the rheology, buoyancy, seismicity, and H2O content of oceanic crust Dynamics of Subduction 337 246 Bebout G.E. et al. AGU Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hacker B.R. 2007 Ascent of the ultrahigh-pressure Western Gneiss Region, Norway Convergent Margin Terranes and Associated Regions: A Tribute to W.G. Ernst Special Paper 419 171 184 Cloos M. Carlson W.D. Gilbert M.C. Liou J.G. Sorensen S.S. Geological Society of America Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hacker B.R. Peacock S.M. Abers G.A. Holloway S.D. 2003 Subduction factory: 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions ?, J. geophys. Res. 108 B1 2030 doi:10.1029/2001JB001129 Google Scholar OpenURL Placeholder Text WorldCat Hammond W.C. Humphreys E.D. 2000 Upper mantle seismic wave velocity: effects of realistic partial melt geometries J. geophys. Res. 105 B5 10 975–10 986 Google Scholar OpenURL Placeholder Text WorldCat Husen S. Kissling E. Flueh E. Asch G. 1999 Accurate hypocenter determination in the seismogenic zone of the subducting Nazca plate in north Chile using a combined on-/offshore network Geophys. J. Int. 138 687 701 Google Scholar Crossref Search ADS WorldCat Hyndman R.D. Peacock S.M. 2003 Serpentinization of the forearc mantle Earth planet. Sci. Lett. 212 417 432 Google Scholar Crossref Search ADS WorldCat Hyndman R.D. Yamano M. Oleskevich D.A. 1997 The seismogenic zone of subduction thrust faults Island Arc 6 244 260 Google Scholar Crossref Search ADS WorldCat Isacks B.-L. 1988 Uplift of the central Andean Plateau and bending of the Bolivian Orocline J. geophys. Res. 93 3211 3231 Google Scholar Crossref Search ADS WorldCat Iwamori H. 2000 Deep subduction of H2O and deflection of volcanic chain towards backarc near triple junction due to lower temperature Earth planet. Sci. Lett. 181 41 46 Google Scholar Crossref Search ADS WorldCat Ji S. Li A. Wang Q. Long Ch. Wang H. Marcotte D. Salisbuty M. 2013 Seismic velocities, anisotropy, and shear-wave splitting of antigorite serpentinites and tectonic implications for subduction zones J. geophys. Res. 118 1015 1037 Google Scholar Crossref Search ADS WorldCat Kamiya S. Kobayashi Y. 2000 Seismological evidence for the existence of serpentinized wedge mantle Geophys. Res. Lett. 27 819 822 Google Scholar Crossref Search ADS WorldCat Karato S. 2012 On the origin of the asthenosphere Earth planet. Sci. Lett. 321–322 95 103 Google Scholar Crossref Search ADS WorldCat Kelemen P.B. Hart S.R. Bernstein S. 1998 Silica enrichment in the continental upper mantle via melt/rock reaction Earth planet. Sci. Lett. 164 387 406 Google Scholar Crossref Search ADS WorldCat Kerrick D. 2002 Serpentinite seduction Science 298 1344 1345 Google Scholar Crossref Search ADS PubMed WorldCat Kerrick D.M. Connolly J.A.D. 2001 Metamorphic devolatilization of subducted marine sediments and the transport of volatiles into the Earth's mantle Nature 411 293 296 Google Scholar Crossref Search ADS PubMed WorldCat Khazaradze G. Klotz J. 2003 Short- and long-term effects of GPS measured crustal deformation rates along the south central Andes J. geophys. Res. 108 2289 doi:10.1029/2002JB001879.361 Google Scholar Crossref Search ADS WorldCat Kirby S. Engdahl E.R. Denlinger R. 1996 Intermediate-depth intraslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs Subduction: Top to Bottom: American Geophysical Union Geophysical Monography 96 195 214 Bebout G.E. et al. AGU Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Koulakov I. Sobolev S. Asch G. 2006 P- and S-velocity images of the lithosphere-asthenosphere system in the central Andes from local-source tomographic inversion Geophys. J. Int. 167 106 126 Google Scholar Crossref Search ADS WorldCat Li Z. Roecker S. Wei B. Wang H. Schelochkov G. Bragin V. 2009 Tomographic image of the crust and upper mantle beneath the western Tien Shan from the MANAS broadband deployment: possible evidence for lithospheric delamination Tectonophysics 477 1–2 49 57 Google Scholar OpenURL Placeholder Text WorldCat Loveless J.P. Hoke G.D. Allmendinger R.W. González G. Isacks B.L. Carrizo D.A. 2005 Pervasive cracking of the northern Chilean Coastal Cordillera: new evidence for forearc extension Geology 33 973 976 Google Scholar Crossref Search ADS WorldCat McGlashan N. Brown L. Kay S. 2008 Crustal thickness in the central Andes from teleseismically recorded depth phase precursors Geophys. J. Int. 175 1013 1022 Google Scholar Crossref Search ADS WorldCat Métois M. et al. 2013 Revisiting the North Chile seismic gap 255 segmentation using GPS-derived interseismic coupling Geophys. J. Int. 194 1283 1294 256 Google Scholar Crossref Search ADS WorldCat Montagner J.-P. 1986 Regional three-dimensional structures using long-period surface waves Ann. Geophys. 4 283 291 Google Scholar OpenURL Placeholder Text WorldCat Moran A.E. Sisson V.B. Leeman W.P. 1992 Boron depletion during progressive metamorphism—implications for subduction processes Earth planet. Sci. Lett. 111 331 349 Google Scholar Crossref Search ADS WorldCat Nakajima J. Tsuji Y. Hasegawa A. 2009 Seismic evidence for thermally-controlled dehydration in subducting oceanic crust Geophys. Res. Lett. 36 L03303 doi:10.1029/2008GL036865 Google Scholar Crossref Search ADS WorldCat Norabuena E. et al. 1998 Space geodetic observations of Nazca-South America convergence across the Central Andes Science 279 358 362 Google Scholar Crossref Search ADS PubMed WorldCat Nunn C. Roecker S. Priestly K. Liang X. Gilligan A. 2014 Joint inversion of surface waves and teleseismic bodywaves across the Tibetan collision zone: the fate of subducted Indian lithosphere Geophys. J. Int. 198 1526 1542 Google Scholar Crossref Search ADS WorldCat Omori S. Kamiya S. Maruyama S. Zhao D. 2002 Morphology of the intraslab seismic zone and devolatilization phase equilibria of the subducting slab peridotite Bull. Earth. Res. Inst. Univ. Tokyo 76 455 478 Google Scholar OpenURL Placeholder Text WorldCat Omori S. Komabayashi T. Maruyama S. 2004 Dehydration and earthquakes in the subducting slab: empirical link in intermediate and deep seismic zones Phys. Earth planet. Inter. 146 297 311 Google Scholar Crossref Search ADS WorldCat Oncken O. et al. 2003 Seismic imaging of a convergent continental margin and plateau in the central Andes (Andean Continental Research Project 1996 (ANCORP’ 96)) J. geophys. Res. 108 2328 doi:10.1029/2002JB001771 Google Scholar Crossref Search ADS WorldCat Ortega-Culaciati F. et al. 2015 Imaging the seismic cycle in the central Andean subduction zone from geodetic observations Proceedings of the AGU Fall Meeting T43C-3024 San Francisco, CA, USA Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Paige C.C. Saunders M.A. 1982 LSQR: an algorithm for sparse linear equations and sparse least squares ACM Trans. Math. Softw. 8 1 43 71 Google Scholar Crossref Search ADS WorldCat Patzwahl R. Mechie J. Schulze A. Giese P. 1999 Two-dimensional velocity models of the Nazca plate subduction zone between 19.5°S and 25°S from wide-angle seismic measurements during the CINCA95 project J. geophys. Res. 104 7293–7317 Google Scholar Crossref Search ADS WorldCat Peacock S. Hyndman R.D. 1999 Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes Geophys. Res. Lett. 26 2517 2520 Google Scholar Crossref Search ADS WorldCat Peacock S.M. 1993 The importance of blueschist-eclogite dehydration in subducting oceanic crust Geol. Soc. Am. Bull. 105 684 694 Google Scholar Crossref Search ADS WorldCat Peacock S.M. 2001 Are the lower planes of double seismic zones caused by serpentine dehydration in subducting oceanic mantle? Geology 29 4 299 302 Google Scholar Crossref Search ADS WorldCat Plank T. Langmuir C.H. 1993 Tracing trace-elements from sediment input to volcanic output at subduction zones Nature 362 739 743 Google Scholar Crossref Search ADS WorldCat Prevot R. Roecker S.W. Isacks B.L. Chatelain J.L. 1991 Mapping of low P wave velocity structures in the subducting plate of the central New Hebrides, southwest Pacific J. geophys. Res. 96 19 825–19 842 Google Scholar Crossref Search ADS WorldCat Pritchard M.E. Simons M. 2006 An aseismic slip pulse in northern Chile and along-strike variations in seismogenic behavior J. geophys. Res. 111 B8 1 14 Google Scholar Crossref Search ADS PubMed WorldCat Reutter K.-J. Giese P. Goetze H.-J. Scheuber E. Schwab K. Schwarz G. Wigger P. 1988 Structures and crustal development of the Central Andes between 21° and 25°S The Southern Central Andes 17 231 261 Bahlburg H. Breitkreuz C. Giese P. Lecture Notes in Earth Sciences Springer-Verlag Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rietbrock A. Waldhauser F. 2004 A narrowly spaced double seismic zone in the subducting Nazca plate Geophys. Res. Lett. 31 L10608 doi:10.1029/2004GL019610 Google Scholar Crossref Search ADS WorldCat Roecker S. Thurber C. McPhee D. 2004 Joint inversion of gravity and arrival time data from Parkfield: new constraints on structure and hypocenter locations near the SAFOD drill site Geophys. Res. Lett. 31 L12S04 doi:10.1029/2003GL019396 Google Scholar Crossref Search ADS WorldCat Roecker S. Thurber C. Roberts K. 2006 Refining the image of the San Andreas Fault near Parkfield, California using a finite difference travel time computation technique Tectonophysics 426 189 205 Google Scholar Crossref Search ADS WorldCat Rutland R.W.R. 1971 Andean orogeny and ocean floor spreading Nature 233 252 255 Google Scholar Crossref Search ADS PubMed WorldCat Scheuber E. Bogdanic T. Jensen A. Reutter K.-J. 1994 Tectonic development of the north Chilean Andes in relation to plate convergence and magmatism since Jurassic Tectonics of the Southern Central Andes 121 139 Reutter K.J. Scheuber E. Wigger P.J. Springer-Verlag Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Schmitz M. Kley J. 1997 The geometry of the central Andean backarc crust: joint interpretation of cross section balancing and seismic refraction data J. South Am. Earth Sci. 10 99 110 Google Scholar Crossref Search ADS WorldCat Schmitz M. et al. 1999 The crustal structure beneath the Central Andean forearc and magmatic arc as derived from seismic studies; the PISCO 94 experiment in northern Chile (21°–23°S) J. South Am. Earth Sci. 12 237–260 Google Scholar Crossref Search ADS WorldCat Scholz C.H. Small C. 1997 The effect of seamount subduction on seismic coupling Geology 25 6 487 490 Google Scholar Crossref Search ADS WorldCat Schurr B. Asch G. Rietbrock A. Trumbull R. Haberland C. 2003 Complex patterns of fluid and melt transport in the central Andean subduction zone revealed by attenuation tomography Earth planet. Sci. Lett. 215 105 119 Google Scholar Crossref Search ADS WorldCat Sella G.F. Dixon T.H. Mao A. 2002 REVEL: a model for recent plate velocities from space geodesy J. geophys. Res. 107 B4 2081 doi:10.1029/2000JB000033 Google Scholar Crossref Search ADS WorldCat Springer M. 1999 Interpretation of heat-flow density in the Central Andes Tectonophysics 306 377 395 Google Scholar Crossref Search ADS WorldCat Stern C.R. Kilian R. 1996 Role of the subducted slab, mantle wedge and continental crust in the generation of adakites from the Andean Austral volcanic zone Contribut. Mineral. Petrol. 123 3 263 281 Google Scholar Crossref Search ADS WorldCat Tsuji Y. Nakajima J. Hasegawa A. 2008 Tomographic evidence for hydrated oceanic crust of the Pacific slab beneath northeastern Japan: implications for water transportation in subduction zones Geophys. Res. Lett. 35 L14308 doi:10.1029/2008GL034461 Google Scholar Crossref Search ADS WorldCat Uchida N. Nakajima J. Hasegawa A. Matsuzawa T. 2009 What controls interplate coupling?: evidence for abrupt change in coupling across a border between two overlying plates in the NE Japan subduction zone Earth planet. Sci. Lett. 283 1–4 111 121 Google Scholar Crossref Search ADS WorldCat Ulmer P. Trommsdorff V. 1995 Serpentine stability to mantle depths and subduction-related magmatism Science 268 858 861 Google Scholar Crossref Search ADS PubMed WorldCat van Keken P.E. Hacker B.R. Syracuse E.M. Abers G.A. 2011 Subduction factory: 4. Depth‐dependent flux of H2O from subducting slabs worldwide J. geophys Res. 116 B01401 doi:10.1029/2010JB007922 Google Scholar Crossref Search ADS WorldCat Von Huene R. Weinrebe W. Heeren F. 1999 Subduction erosion along the north Chile margin J. Geodyn. 27 345 358 Google Scholar Crossref Search ADS WorldCat Watanabe T. Kasami H. Ohshima S. 2007 Compressional and shear wave velocities of serpentinized peridotites up to 200 MPa Earth planet. Space 59 233 244 Google Scholar Crossref Search ADS WorldCat Watts A.B. 2001 Isostasy and Flexure of the Lithosphere Cambridge Univ. Press 458 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wörner G. Moorbath S. Horn S. Entenmann J. Harmon R.S. Davidson J.P. Lopez-Escobar L. 1994 Large- and fine-scale geochemical variations along the Andean arc of northern Chile (17.5–22S) Tectonics of the Southern Central Andes 69 76 Reutter K.-J. Scheuber E. Wigger P. Springer-Verlag Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wörner G. Uhlig D. Kohler I. Seyfried H. 2002 Evolution of the West Andean Escarpment at 18°S (N. Chile) during the last 25 Ma: uplift, erosion and collapse through time Tectonophysics 345 183 198 Google Scholar Crossref Search ADS WorldCat Yao H. Van der Hilst R.D. 2009 Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to SE Tibet Geophys. J. Int. 179 2 1113 1132 Google Scholar Crossref Search ADS WorldCat Yao H. Van der Hilst R.D. De Hoop M.V. 2006 Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis. I. Phase velocity maps Geophys. J. Int. 166 732 744 Google Scholar Crossref Search ADS WorldCat Yogodzinski G.M. Kay R.W. Volynets O.N. Koloskov A.V. Kay S.M. 1995 Magnesian andesite in the western Aleutian Komandorsky region: implications for slab melting and processes in the mantle wedge Geol. Soc. Am. Bull. 107 505 519 Google Scholar Crossref Search ADS WorldCat Yuan X. et al. 2000 New constraints on subduction and collision processes in the central Andes from P-to-S converted seismic phases Nature 408 958 961 Google Scholar Crossref Search ADS PubMed WorldCat Yuan X. Sovolev S. Kind R. 2002 Moho topography in the Central Andes and its geodynamic implications Earth planet. Sci. Lett. 199 389 402 Google Scholar Crossref Search ADS WorldCat Zhang H. Roecker S. Thurber C.H. Wang W. 2012 Seismic imaging of microblocks and weak zones in the crust beneath the southeastern margin of the Tibetan plateau Earth Sciences Dar I.A. InTech Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Zhang N. Zhong S. Leng W. Li Z.-X. 2010 A model for the evolution of the Earth's mantle structure since the Early Paleozoic J. geophys. Res. 115 B06401 doi:10.1029/2009JB006896 Google Scholar OpenURL Placeholder Text WorldCat Ziegler A.M. Barrett S.F. Scotese C.R. 1981 Palaeoclimate, sedimentation and continental accretion Phil. Trans. R. Soc. Lond., Ser. A 301 253 264 Google Scholar Crossref Search ADS WorldCat © The Authors 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © The Authors 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Three-dimensional elastic wave speeds in the northern Chile subduction zone: variations in hydration in the supraslab mantle
JO - Geophysical Journal International
DO - 10.1093/gji/ggw318
DA - 2016-11-01
UR - https://www.deepdyve.com/lp/oxford-university-press/three-dimensional-elastic-wave-speeds-in-the-northern-chile-subduction-gUAAOBq5eO
SP - 1080
EP - 1105
VL - 207
IS - 2
DP - DeepDyve
ER -