Article Highlights

  • We present imaging of the Tonga–Samoa region from intermediate and deep-focus earthquakes

  • Distinct anomalies in the Vp/Vs ratio span the whole upper mantle to ~ 600 km depths close to the Samoa plume

  • The anomalies trace deep seated magmatic processes evoked by the interaction of the Samoa plume with the Tonga slab

1 Introduction

Mantle plumes, first suggested by Wilson (1963) and Morgan (1971), represent relatively hot, low-density mantle material that ascends because of its buoyancy. They were introduced to explain the intraplate oceanic island chains ageing progressively along the chain in the direction of plate motions. The surface manifestation of the plumes forms large hotspots relatively fixed to each other which represent focused zones of melting characterized by high heat flow. Their thermal origin together with geochemical imprints such as high 3He/4He ratios observed in hotspot basalts (Kellogg and Wasserburg, 1990) indicate that the plumes are rooted deep in the mantle and transport primordial mantle material upwards (e.g., Condie, 2001). This is one of the fundaments that discriminate convective flows of plate tectonic movements primarily driven by sinking of cold plates in favour of deep-seated plumes driven by heat exchange (Morgan, 1971; Foulger and Natland, 2003; Foulger, 2010). Mantle plumes are supposed to consist of two parts, a large bulbous head at the top and relatively narrow tail connecting the plume with its source. They are distributed irregularly in the Earth, mostly inside tectonic plates; however, some of them occur near lithospheric plate margins reflecting the large diversity of interactions (Foulger, 2010).

The plume behaviour near subducting plates differs significantly from a typical mantle-plume hotspot model in the intraplate environment and is still poorly understood. Samoa in the south-west Pacific (Fig. 1) represents such a hotspot in a proximity to a subducting plate margin. It sits in a place where the Pacific plate tears with one part moving over the hotspot of upwelling magma forming a chain of volcanic islands, while the other one subducting beneath the Australian Plate. Early studies disputed timing of the Samoan volcanic chain and advanced a non-plume origin such as lithospheric bending (Natland, 1980; Foulger and Natland, 2003). Recent Samoa studies favoured its plume origin; however, complex interactions with the nearby Tonga subduction zone discriminated it from a classical hotspot model (Hart et al., 2004; Koppers et al., 2008; Chang et al., 2016). Strak and Schellart (2018) debated two subsequent geodynamic processes responsible for the intraplate volcanism in the Samoan region: the plume model in the early stages, followed by a non-hotspot volcanism resulting from the interaction between subduction and the plume in the later stages. Laboratory experiments confirm that the plume-driven flow can be altered in subduction zones, where the slab motion plays a dominant role in the mantle (Druken et al., 2014). In all cases, the crucial question remains, how deep can the fertile source reach.

Fig. 1
figure 1

The Tonga–Samoa region with superimposed deep seismicity. a The epicentres of earthquakes in a depth range of 150–670 km (according to the Global CMT catalog) for the period of 1992–2018 with the magnitude mb > 4.7 colour-coded by depth (from orange to pink). Seismic stations MSVF (Fiji Islands) and RAR (Cook Islands) are marked by yellow triangles. Two white crosses mark the location of events shown in Fig. 4. The inset shows the position of the Tonga-Fiji-Kermadec subduction zone (red box). b The epicentres of all detected earthquakes in a depth range of 70–670 km (according to the USGS catalogue), colour-coded by depth (modified after Martin, 2014). FP, Fiji; S, Samoa. Note a distinct inflexed shape of the earthquake foci in the northern Tonga region

Full size image

Due to the proximity to the subducting plate, the Tonga–Samoa region (Fig. 1) is enigmatic and displays numerous unique features. The Tonga subduction zone is the fastest convergent plate boundary in the world at a rate of 24 cm/year. It exhibits the most intense mantle seismicity being 10 times larger than any other subduction zone in the world (Bevis et al., 1995; Gurnis et al., 2000). Its extensional Lau Basin is the fastest opening back-arc at a rate of 17 cm/year (Bird, 2003). Anomalous ocean-island-basalt records in Samoa and north-western part of the Lau Basin show extremely high 3He/4He ratios and isotopic Sr–Nd–Pb element enrichment which appear to originate from the Samoa plume (Lupton et al., 2009; Price et al., 2014). On the other hand, debates on the Samoan volcanic lineament as a plume-driven or non-plume-driven hotspot chain involves limitations due to predicted age of some volcanoes (Natland, 1980).

The origin of anomalies connected with melt can essentially be traced through imaging of compressional (Vp) and shear (Vs) seismic waves and their ratio (Vp/Vs) by seismic tomography. However, detecting narrow vertical structures that may represent plume tails over the entire mantle is a challenging task. In the Samoan region, the plume-related mechanisms were indicated by low‐velocity anomalies traced deep to the lower mantle imaged by recent global teleseismic tomography studies (Bijwaard and Spakman, 2000; Hall and Spakman, 2002; Montelli et al., 2004; French and Romanowicz, 2015; Maguire et al., 2018). However, global whole-mantle tomography studies suffer from limited resolution and details in seismic structure cannot be resolved with the required accuracy. Thus, the upper mantle dynamics and interactions between the Samoa plume and the Tonga slab remain enigmatic.

Beside mantle models from the global tomography imaging, the intense seismic activity in the Tonga subduction region offers the opportunity for regional tomography studies with regional-scale details utilizing intermediate and deep-focus earthquakes in the slab. Between 1993 and 1995, elaborate land experiment SPASE and passive seismic ocean bottom experiment LABBATS were designed across the northern Tonga and Fiji (Wiens et al., 1995) providing information on the structure in the mantle wedge and the back-arc Lau Basin (Zhao et al., 1997; Xu and Wiens, 1997; Koper et al, 1999; Conder and Wiens, 2006). These experiments improved our knowledge of the region and emphasized the importance of regional studies; however, they retained the following limitations.

  • Zhao et al. (1997) analysed the arrival times of P-waves of earthquakes that occurred in the Tonga-Fiji region during the experiment. They imaged the mantle wedge with low-velocity anomalies of −2 to −4% at depths of 100–400 km and very low velocity anomalies of −5 to −7% above (at depths < 100 km) referenced against the IASP91 model (Kennett and Engdahl, 1991). The slab was imaged by velocities of 4–6% higher than in the surrounding mantle. The lowest velocities were found westwards of the Lau Basin and beneath the Tonga volcanic arc in agreement with the attenuation studies, which showed extremely high attenuation beneath the Lau back arc (Roth et al., 1999). However, Zhao et al. (1997) did not include the S-waves in their interpretation and, therefore, no Vp/Vs constraints on a possible distribution of melt and volatiles were provided.

  • Koper et al. (1999) used the P- and S-wave traveltimes from selected earthquakes to investigate a fractional change in the S-wave velocities relative to the P-wave velocities (ratio dlnVs/dlnVp). They found that this ratio was 1.2–1.3 for the whole mantle wedge region, consistently with the experimentally predicted values mapping the effect of temperature on the P and S velocities. However, they only solved for the bulk mantle wedge and provided no information on the lateral Vp/Vs variations within the wedge.

  • Conder and Wiens (2006) expanded the P- and S-wave data from regional studies by data from the 2001–2002 SAFT experiment (Tibi and Wiens, 2005) and jointly inverted for the Vp and Vp/Vs structure. Their distribution of Vp anomalies down to depths of 700 km were similar to Zhao et al. (1997), while the Vs velocity variations exhibited a fast slab (0.2–0.4 km/s) and a slow back arc and spreading centre (0.3–0.5 km/s) in respect to the best fitting model. Their Vp/Vs structure highlighted positive anomalies associated with the arc (~ 0.05–0.1) and the spreading centre (~ 0.1), and the large negative anomalies associated with the cold down-going slab (~ 0.04–0.1). They provided detailed imaging down to depths of 700 km; however, they confined their results to 2-D profile along the main E–W trending line (Zhao et al., 1997) across Fiji.

  • Thus, the differences between the northern and southern Tonga segments were not addressed and potential information on melt generation connected with the Samoan plume, its interactions with the slab and delineation of areas containing thermal anomalies are still missing.

A less costly but still very powerful alternative to regional studies might be employing data recorded at permanent stations of the Global Seismographic Network (GSN) operated jointly by IRIS and US Geological Survey. This network consists of 152 broadband stations, globally distributed and providing open access data (https://www.iris.edu/hq/programs/gsn). Although the number of suitable nearby stations is limited, it can be overcome by the advantage of numerous earthquakes observed during a long period of time and thus covering densely the whole active slab. Moreover, focusing on earthquakes reported in the Global CMT catalog (Dziewonski, et al., 1981; Ekström, et al., 2012) offers an ideal source of accurately known locations and origin times. Together with an advanced approach in data processing, the time residuals extracted from theoretical traveltimes predicted by global velocity models can be recovered with a high accuracy. Such residua converted to velocities can afford for imaging of the upper mantle surrounding the slab.

The aim of this study is not to surpass regional tomography models but to provide additional information on seismic velocities and especially their Vp/Vs ratio around the slab. Because of the vast amount of high-quality data, a new perspective is incorporated by comparing the upper mantle structure of the northern and southern Tonga segments, which was not considered in previous studies. When confronted with geochemical results on ocean island basalts, the distribution of seismicity, deformation and extension in the area, the constraints on thermal/compositional nature of heterogeneities can be obtained.

2 Tectonic and Structural Setting

2.1 Tectonics

Tonga belongs to the 2,550-km-long Tonga-Kermadec subduction system, which is the outcome of the Pacific plate submerging beneath the Australian plate (Fig. 1). Tonga subducts steeply, dipping ~ 60° into the mantle (Fig. 2) and is the fastest converging and the most seismically active system in the world. The very high subduction rates arise from a combination of the Australian-Pacific plate convergence in the east and the back-arc spreading with the Lau Basin opening in the west (Bevis et al., 1995). The convergence rate is assessed to ~ 15 cm/year, but due to the trapezoidal Lau Basin opening from N to S, the northern Tonga Trench indicates the convergence rate of 24 cm/year which decreases to 6 cm/year southwards along the Tonga–Kermadec arc (Martinez et al., 2006; Smith and Price, 2006). In the north, the Tonga Trench ends south of Samoa where the system bends perpendicularly to the west, intersects the Fiji Fracture Zone and continues to the inactive Vitiaz Trench. Further to the west, another subduction system of Vanuatu is located. To the south, the Louisville Ridge, a striking seamount chain of submerged volcanoes, intersects the Tonga–Kermadec Trench at 25.6°S and separates the Tonga and Kermadec subduction segments (Fig. 1).

Fig. 2
figure 2

The 3-D cubic-spline interpolation of 1687 earthquake foci (10-km grid) in a depth range of 150–670 km and mb > 4.7 (black dots). The locations are from the Global CMT catalog for the period of 1976–2017. Depth is colour-coded; azimuth/elevation of the view angle is 295°/12°. Transparent horizontal planes indicate the 410-km (green) and 660-km (violet) discontinuities bordering the upper mantle transition zone. Note the westward upwelling of foci in this zone. The grey dots mark the foci bottom projections; note their inflexed shape in the north. The grey plane represents the surface with the location of station MSVF (red triangle), and the grey triangle represents its bottom projection. Station RAR (indicated by the arrow and red triangle) is aside of the imaged area. Note that the imaged surface does not represent the actual slab geometry

Full size image

The extension-dominated, non-accretionary convergent margin of the Tonga-Kermadec system led to a formation of new oceanic microplates. Their formation started in Miocene ~ 4–5 Ma, when the region went through crustal extension and formation of numerous different spreading centres leading to a separation of the Australian and Pacific plates (Bird, 2003; Smith and Price, 2006). Prior to that, the 100–140 Ma old and relatively cold Pacific plate was subducting beneath the Australian plate generating the Lau–Colville Ridge. The crustal extension in Miocene initiated opening of the back-arc Lau Basin and its propagation southwards. Nowadays, due to the opening, new crust is produced in front of the Tonga system, while the old crust is absorbed behind in the Tonga Trench due to subduction. The subduction is attributed to the excess weight of the old Pacific plate entering the mantle, while by contrast, the extension in the Lau Basin is driven by clockwise rollback of the Tonga plate (Schellart and Spakman, 2012).

The Samoa hotspot lies ~ 130 km from the northern termination of the Tonga Trench, at a place where the Pacific plate pulls apart and its southern portion subducts beneath the Australian plate, while its northern portion continues westwards. The deep mantle plume concept of the Samoa hotspot is supported by the Samoan volcanic lineament indicating the mantle plume volcanic chain gradually ageing to the west. Seismic low-velocity structure in the mantle imaged in global tomography models could represent a plume originating at the core-mantle boundary (Gurnis et al., 2000; Courtillot et al., 2003; Montelli et al., 2004; French and Romanowicz, 2015; Chang et al., 2016; Maguire et al., 2018). In addition, numerous geochemical studies are in favour of the plume model (Hart et al. 2004; Koppers et al., 2008; Jackson et al., 2010; Price et al., 2014). However, the plume concept also involves considerable limitations, with the most striking occurrence of rejuvenated volcanism in the western Samoa where several volcanoes show a volcanic activity too young to fit the plume predicted age (Natland, 1980).

2.2 Structural Setting

The Tonga wedge and the associated back-arc Lau Basin are imaged by low-velocity anomalies (Zhao et al., 1997; Koper et al., 1999; Conder and Wiens, 2006) in agreement with the results from attenuation studies detecting an extensive zone of high attenuation in the uppermost mantle (Roth et al., 1999; Wei and Wiens, 2018). This low-velocity zone extends westwards afar off the active Lau spreading centre and imply high temperatures in the upper mantle related to the magma production generating the fast-spreading back-arc Lau Basin (Zhao et al., 1997; Martinez and Taylor, 2002). Thermal modelling suggested a counter-flow originating in the subduction mantle wedge and its upwelling in the back-arc (Koper et al., 1999). On the other hand, the proximity of the Samoa hotspot influences the mantle-wedge setting. The plume-related material was identified in the northern Lau Basin southwards of Samoa (Pearce et al., 2007; Chang et al., 2016) with inflows suggested by the Tonga trench retreat (Strak and Schellart, 2018). North–south flows promoted by the Samoa presence were also identified from shear‐wave splitting (Smith et al., 2001). Such hotter material increases the temperature in the northern Tonga back arc and can also change the topography of the thermally controlled 410-km discontinuity as suggested by Li et al. (2019).

At the base of the upper mantle, global seismic tomography models (Hall and Spakman, 2002; Schellart and Spakman, 2012) revealed increased seismic velocities interpreted as a large amount of slab-like material retained in the mantle transition zone. Combined spatial distribution of the P- and S-wave velocities with deep seismicity and focal mechanisms allocated high velocity anomalies of about 1–1.5% in the transition zone to the juxtaposition of a steeply dipping Tonga slab and a sub-horizontal remnant of a slab not any more connected with the actively subducting lithosphere (Chen and Brudzinski, 2001; Hall and Spakman, 2002; Burkett and Billen, 2010; Fukao and Obayashi, 2013). Gurnis et al. (2000), Richards et al. (2011), and Chang et al. (2016) discussed a concave folded shape of this detached slab beneath the North Fiji and Lau Basins. They suggested plume-related processes to explain it; however, the spatial extent of the Samoa plume interactions at the contact with the Tonga slab remains undisclosed.

3 Method and Data Processing

3.1 Data

We focused on earthquakes located in the Tonga subduction zone (latitude 14.9–26°S, longitude 174°W-178°E) and reported in the Global CMT catalog (Dziewonski et al., 1981; Ekström et al., 2012). The data set consists of 1687 earthquakes in a depth range of 150–670 km with magnitudes mb > 4.7 that occurred during the period of 1976–2017 (Fig. 2). Together they form steeply dipping Wadati–Benioff zone with a saddle shape and marked decrease in dip at ~ 410 km depth followed by a dip increase at 550–600 km depth (e.g., Giardini and Woodhouse, 1984; Myhill, 2013). The seismic activity is very high in shallow parts and generally decreases from depths of ~ 100–150 km. At intermediate depths, the seismic energy is not released uniformly but in zones forming clusters and double seismic zones (Chen et al, 2004; Brudzinski et al. 2007). The most intense seismicity is disclosed in the mantle transition zone in a range of 500–670 km (Fig. 3) with 3-D morphology and complex seismicity and deformation patterns interpreted as a slab upwelling, buckling, tearing or a remnant slab (Chen and Brudzinski, 2001; Warren et al., 2007; Bonnardot et al., 2009; Fukao and Obayashi, 2013; Myhill, 2013; Hayes et al., 2018; Billen et al., 2003; Billen 2020). The northernmost part of the Tonga seismicity is influenced by the inactive Vitiaz lineament and by tearing of the Pacific plate near Samoa (Fig. 1), (Giardini and Woodhouse, 1984); to the south, the seismicity in the Kermadec slab is less intense and terminates at ∼480 km depth.

Fig. 3
figure 3

taken from the Global CMT catalog for the period of 1976–2017 with the selected depth range of 150–670 km and mb > 4.7. Note the world most intense seismicity in the Tonga upper mantle transition zone at depths 490–670 km. b Distribution of the earthquake magnitudes (mb) in Tonga

a Distribution of foci of the Tonga earthquakes with depth (red bars) superimposed on global seismicity where the Tonga earthquakes are excluded (blue bars). The data are

Full size image

We processed waveforms recorded at the IRIS-USGS Global Seismographic Network stations (seismic network codes II and IU) which we retrieved from the IRIS Data Management Center. We selected two closest stations, MSVF and RAR (Figs. 1a, 2), operating since 1994 and 1992, respectively, having long-time recordings and being in a favourable position to the Tonga subduction. The MSVF station was located at Fiji Islands (latitude 17.7448°S; longitude 178.0528°E), westwards of the Tonga slab (mean epicentral distance of 5.95°, in the range of 1.7–16.7°); the RAR station was located in the Pacific Ocean, at Rarotonga, Cook Islands (latitude -21.212°S, longitude -59.773°W), eastwards of the Tonga slab (with the mean epicentral distance of 17.5° in the range of 14.2–20.7°). Both stations were equipped with the Geotech broadband borehole seismometer KS-54000 designed for low-noise monitoring in frequencies of 0.003–5 Hz and buried at 100 m depths.

To avoid complexities in the uppermost mantle structure contaminated by lithosphere and shallow tectonic processes, we eliminated the earthquakes down to the depths of 150 km and focused on a depth range of 150–670 km. In these depths, the mb magnitude range was of 4.8–7.0; however, the prevailing magnitudes were between 5.0 and 6.0 (Fig. 3). The selection resulted in the extraction of 538 events for station MSVF recorded during the period of 1994–2017 and in 868 events for station RAR recorded during 1999–2017.

3.2 Waveform Pre-Processing

The aim was to retrieve high-precision velocity perturbations and their ratio in the upper mantle around the Tonga slab. This required extracting accurate traveltime residua between observed and calculated direct P- and S-wave arrivals for individual events at a given station. We analysed vertical displacement recordings sampled with frequency of 20 Hz. The waveforms were filtered with a zero-phase bandpass filter, which eliminated time delays between the filtered waveform and the unfiltered signal. To optimally suppress noise and spurious high-frequency phases, we tested a series of Butterworth filters in several frequency bands. The low corner frequency was alternatively 20, 30 and 40 s, and the high corner frequency was 4, 6, 8 and 10 s. The optimum balance between waveform smoothing needed for increasing a similarity of processed signals and preserving maximum information resulted in the bandpass filter of 8–40 s (Fig. 4) similar to filtering of 0.02–0.2 Hz applied by Conder and Wiens (2006) to the regional Tonga data set. Too narrow filters produced rather oscillating signals with no clear arrival onsets and too broad filters kept high frequencies responsible for dissimilarity of the waveforms.

Fig. 4
figure 4

Example of vertical displacement recordings for the P- and S-waves for two earthquakes with different magnitudes at station MSVF. Unfiltered waveform (blue lines) and signal filtered with a zero-phase Butterworth filter with the frequency band of 8–40 s (red lines). a, b Event 33 with mb 5.7 at a depth of 611 km. c, d Event 63 with mb 5.3 at a depth of 601 km. Locations of events are depicted in Fig. 1. Note that high seismic noise prevents accurate manual picking of the first arrival and that the smoothed filtered signal preserves a consistent position of the maximum amplitude in the unfiltered signal

Full size image

3.3 Waveform Cross-Correlation and Principal Component Analysis

Since the traveltime residua produced by the velocity perturbations are rather small (up to several seconds), we cannot determine them using manual picking of the arrival times. The errors of manual picking are due to noise in data up to ± 2 s (see Fig. 4); hence, they would completely mask or seriously deteriorate the effects of the velocity perturbations. Therefore, the traveltime residua were retrieved from the cross-correlation of waveforms further refined by correlating with the common wavelet obtained by the principal component analysis. The cross-correlation of waveforms for calculating accurate differential traveltimes became popular after developing the double-difference location algorithm by Waldhauser and Ellsworth (2000) and the double-difference tomography method by Zhang and Thurber (2006). The waveform cross-correlation can provide arrival times with the accuracy higher by more than one order compared to the manual picking, and found numerous applications in calculating precise earthquake locations (Schaff and Waldhauser, 2005; Hauksson et al., 2012; Trugman and Shearer, 2017), as well as in tomography studies on all scales (Hauksson and Unruh, 2007; Houser et al., 2007; Li et al., 2008; Buehler and Shearer, 2012; Obayashi et al., 2013). The basic condition for a successful application of the waveform cross-correlation method is a similar shape of the analysed waveforms (Schaff and Waldhauser, 2005). Therefore, this method is particularly suitable for analysing intermediate and deep-focus earthquakes, which produce rather simple pulse-dominated waveforms.

The differential traveltimes are usually calculated by the cross-correlation of waveforms with a reference noise-free waveform of a typical shape, for which the arrival time is well determined and which yields high correlations with other waveforms. However, selection of the reference waveform is somewhat arbitrary. In our study, we applied a mathematically more correct way and we extracted the common wavelet of all analysed waveforms (i.e. the most typical waveform) by applying the principal component analysis (PCA). The PCA is a statistical tool for computing principal components of a large data set, which reduces the number of variables and finds the first few principal components preserving as much independent information as possible (Jackson, 1991). In this study, the first principal component represents a ‘stacked’ wavelet that is common to all processed events (Vavryčuk et al., 2017, 2021). In fact, it represents a centre of a cluster formed by waveforms in a multidimensional space, in which the distance between the individual waveforms is measured by their cross-correlation. Higher principal components comprise deviations of the individual waveforms from the common wavelet including noise. The PCA method is fast, robust and works even for waveforms with mixed polarities. The retrieved common wavelet was a result of the PCA applied to the observed waveforms, which were cross-correlated, aligned and windowed independently for both direct P- and S-waves. Consequently, the cross-correlation of waveforms with the common wavelet maximizes the accuracy of the differential traveltimes.

The waveform processing consisted of the following steps:

  1. (1)

    We aligned the waveforms according to the theoretical P- and S-wave arrivals (Fig. 5). For calculation of theoretical traveltimes, we used hypocentral locations and origin times taken from the Global CMT catalog. The arrival times of individual phases were computed by the TauP method (Crotwell et al., 1999). As the reference Earth model, we used the ak135 velocity model (Kennett et al., 1995).

  2. (2)

    We selected the time window around the theoretical traveltimes and independently extracted the direct P- and S-waves. In this way, we suppressed the differences between individual traces produced mainly by scattered waves following the direct waves.

  3. (3)

    Since the locations, magnitudes and focal mechanisms of the analysed earthquakes varied, the waveforms were not identical and thus not fully correlated. In order to remove anomalous events and to suppress errors due to dissimilarity of waveforms, we excluded wavelets, for which the absolute value of the correlation coefficient was less than 0.8. Some authors consider a sufficient correlation of similar waveforms higher than 0.6 (Hauksson et al., 2012), but we applied stricter criteria to ensure that only highly correlated waveforms were involved and only highly accurate differential traveltimes were obtained.

  4. (4)

    The cross-correlated and windowed waveforms were aligned independently for both direct P- and S-waves.

  5. (5)

    We applied the PCA method to the aligned waveforms to retrieve the best common wavelet (Fig. 5c,d). The time differences between the common wavelet and the cross-correlated wavelets of individual events served for the refinement of the direct P- or S-wave time residuals. Since we focused on the cross-correlation of the P- and S-wavelets irrespective of their amplitudes and polarities, the procedure was rather insensitive to the radiation pattern of individual earthquakes.

  6. (6)

    We finally calculated residuals between the theoretical and observed arrival times for both P- and S-waves.

Fig. 5
figure 5

Example of waveform processing of the Tonga-Kermadec earthquakes at seismic station MSVF. a, b The vertical displacement recordings aligned according to theoretical P- and S-wave arrivals (vertical red lines), calculated by the TauP code (Crottwell et al., 1999) in the 1-D reference ak135 velocity model (Kennett et al., 1995). Black dashed lines indicate the window for the cross-correlation and calculation of the PCA common wavelet. c, d Waveforms of the P- and S-waves, respectively, aligned according to the cross-correlations and the retrieved common wavelet. The red curve represents the common wavelet determined by the PCA method. The waveforms were filtered by the bandpass in the frequency range of 8–40 s and comply with the correlation threshold of 0.8. Note that the plots do not show the entire set of the analysed waveforms

Full size image

The waveform processing at the MSVF station resulted in 308 and 358 wavelets for the P- and S-waves, respectively (Fig. 5), ready for assessing the time residuals Δt. The time residuals of the P-waves ranged mostly from −3 to 3 s; the residuals of the S-waves ranged mostly between −7 and 3 s. Similarly, processing at the RAR station resulted in 325 and 404 wavelets for the P- and S-waves, respectively. The time residuals of the P-waves at the RAR station ranged mostly from −3 to 1 s; the residuals of the S-waves ranged mostly between −4 and 1 s. The time residuals Δt associated with the individual earthquake foci were recalculated to the velocity perturbations Δv and normalized to the average velocity v along the whole ray paths. The normalized velocity perturbations Δv/v were projected onto the foci of earthquakes along the latitude from N to S (Fig. 6). In such a way, the presented velocity perturbations are related to the structure sampled by the rays along the path between the earthquake foci and the station, while supressing local properties of the slab.

Fig. 6
figure 6

The colour-coded velocity perturbations Δv/v in the upper mantle around the Tonga slab bordered by the transform fault between Tonga-Kermadec at latitude of 26°S. a, b The velocity perturbations of the P- and S-waves west of the Tonga slab (station MSVF); c, d velocity perturbations of the P- and S-waves east of the Tonga slab (station RAR). The perturbations Δv are normalized to the average velocity v along the whole ray paths and projected onto the foci. Face view in 2-D along the latitude from N to S. Note that the velocity perturbations Δv/v reflect the structure along the whole ray paths between the earthquake foci and the station

Full size image

3.4 Errors in Traveltime Residua and Velocity Perturbations

The retrieved time residuals reflect not only velocity anomalies in the target area, but also errors that might be of random and/or systematic origin. The random errors comprise uncertainties in picking of the arrival times due to noise and dissimilarities in waveforms, and uncertainties in the earthquake locations. Because of strict selection criteria (cross-correlation coefficient higher than 0.8), the errors in picking by the waveform cross-correlation should be less than 0.2 s. In the case of manual picking, the errors would be of one order higher and thus they would mask the effects of the velocity perturbations. The estimates of the location errors reported in the Global CMT catalog (Ekström et al., 2012) are less than 30 km with the majority between 0 and 20 km (Fig. 7) and less than 0.3 s for the earthquake origin times. Since these errors differ for individual earthquakes, they affect the time residuals and velocity perturbations randomly and produce their scatter. Such scatter is visible in Fig. 6 for both the MSVF and RAR stations (e.g., P-wave residuals for foci deeper than 500 km at station MSVF in Fig. 6a). Nevertheless, the random noise in the residuals is not so high to deteriorate significantly the overall pattern, which shows clear systematic trends. This confirms that the random location and picking errors do not seriously affect the time residuals and the velocity perturbations determined by our analysis. This applies to both stations MSVF and RAR (Fig. 6).

Fig. 7
figure 7

Histograms of location errors of the earthquakes used for determination of the P- (a) and S- (b) wave velocity perturbations at station MSVF in the Tonga back-arc region

Full size image

Similarly, the systematic location errors cannot be responsible for the observed trends in the residuals in Fig. 6. The reasons are as follows: (1) we observe quite high negative residuals for the P- and S-waves for the MSVF station at depth less than 300 km. If such anomaly was produced by systematic location errors, the location errors should be about 30 km or more for the whole slab (latitudes from −14° to −26°). Such systematic location errors in such a large range of foci are fully unrealistic. (2) Fig. 6 indicates that the patterns in residuals for the MSVF and RAR stations are uncorrelated. We observe high negative velocity perturbations for the P- and S-waves for the MSVF station, but no anomaly is observed for the similar foci at the RAR station. If the foci at these depths were mislocated, appreciable residuals (either positive or negative) should be present at both stations. Similar conclusions were reported by Conder and Wiens (2006) who tested different location models producing different initial locations and found little impact on final results.

In addition, the Vp/Vs ratio as the main outcome of this study is even less sensitive to systematic location errors. This is confirmed by numerical tests (provided in Sect. 4. Results). These tests reveal that the uncertainties in the Vp/Vs ratio produced by the potential systematic location errors are even more suppressed compared to the uncertainties in the time residuals and velocity perturbations of the P- and S-waves.

3.5 Method Precis

Unlike standard first-arrival or full-waveform inversions, our concept is based on several factors resulting in a new insight into distribution of the P- and S-wave velocities and the Vp/Vs ratios. (1) The intermediate and deep-focus earthquakes are characterized by rather simple pulse-dominated waveforms, which can be effectively cross-correlated. (2) The bandpass filtering of waveforms with the zero-phase filter eliminates the time shifts between filtered and unfiltered signals. The final bandpass filter of 8–40 s represents the optimum waveform smoothing needed for the effective cross-correlation in the waveform processing while preserving maximum information (Fig. 4). (3) The time residuals determined by the cross-correlation of waveforms are more accurate and less influenced by noise and other effects than those determined by the standard picking of the first arrival times. (4) The accuracy of the time residuals is further improved by cross-correlating windowed P or S wavelets of individual earthquakes with the common wavelet constructed by the principal component analysis (Fig. 5). (5) The inexpensive data acquisition, but yet, due to long-time observations reported in the Global CMT catalog, provides a large amount of high-quality data ready to process and invert.

4 Results

4.1 Velocity Anomalies Around the Tonga Slab

The P- and S-wave velocity anomalies in the mantle surrounding the Tonga region were traced by traveltime residuals observed at two nearby stations MSVF and RAR. The rays to the MSVF station sampled the upper mantle westwards of the Tonga slab beneath the Lau Basin (Figs. 1, 8a); the rays to the RAR station sampled the upper mantle of the Pacific plate eastwards of the Tonga region (Figs. 1, 8b). The time residuals were recalculated to the normalized velocity perturbations Δv/v imaged along the latitude (Fig. 6). Further, they were projected onto the surface obtained by a cubic-spline interpolation of the earthquake foci in a 25-km grid (Fig. 9). Note that the presented velocity perturbations are not related to local properties of the subducting slab, but they reflect properties of the structure along the whole ray paths between the earthquake foci and the station.

Fig. 8
figure 8

The ray coverage for station MSVF (a) and station RAR (b). The red dots mark the epicentres of foci; the green lines mark the horizontal projections of rays. The stations are marked by blue triangles

Full size image
Fig. 9
figure 9

The colour-coded velocity perturbations Δv/v in the upper mantle around the Tonga slab bordered by the transform fault between Tonga-Kermadec at latitude of 26°S. a, b The velocity perturbations of the P- and S-waves west of Tonga (station MSVF); c, d the velocity perturbations of the P- and S-waves east of Tonga (station RAR). The perturbations Δv are normalized to the average velocity v along the whole ray paths (similarly to Fig. 6) and projected onto the surface imaged by the interpolated foci in a 25-km grid. Note that the velocity perturbations Δv/v are not related to local properties of the slab but to the structure along the whole ray paths between the earthquake foci and the station. The azimuth/elevation of the view angle is 280°/5°

Full size image

Beneath the back-arc Lau Basin to the west (station MSVF), the velocity perturbations Δv/v ranged from −3 to + 3% for the P-waves (Fig. 9a), and from −5 to + 1% for the S-waves (Fig. 9b). Patterns of the P and S perturbations were roughly similar and with consistent trends. Strong negative velocity anomalies were confined to depths of 150–250 km; positive anomalies were mostly located at depths of 500–670 km. This indicates a highly inhomogeneous structure in the mantle wedge and transition zone due to processes associated with the slab subduction. Besides, distinct lateral variations with negative velocity anomalies were recovered in the northern part of the Tonga region, which points to additional processes connected with the nearby Samoa plume at the northern rim of the Tonga back-arc region.

By contrast, the P- and S-wave velocity perturbations Δv/v in the upper mantle eastwards of the slab (station RAR) were much smaller than in the back-arc region. They attained values between −1 and 2% for the P-waves (Fig. 9c) and between −2 and 1% for the S-waves (Fig. 9d). They also displayed rather positive values at depths of 150–250 km. This indicates that the Earth’s mantle eastwards of the Tonga slab does not deviate significantly from the ak135 reference velocity model. In addition, the smooth and rather small RAR anomalies compared to large systematic MSVF variations prove that the uncertainties due to location errors are insignificant and the methodology based on the waveform cross-correlation is highly accurate for both P- and S-waves.

The notable lateral velocity variations in the Tonga subduction area (station MSVF) point to the presence of two essentially different segments: the northern and southern segments delimited by the latitude of 20°S. The difference between segments is clearly visible in Figs. 10a, 11a, SI Movie S1 showing the Vp/Vs ratios for station MSVF averaged along the whole ray paths and projected along the latitude from N to S. The high Vp/Vs ratios are concentrated mostly in the northern segment, which is characterized by a complicated slab geometry. Interestingly, the high Vp/Vs values spread across the whole upper mantle down to depths of ~ 500–600 km. The southern segment is more regular except for a localized low Vp/Vs anomaly between 350 and 550 km depths at around 22°S latitude.

Fig. 10
figure 10

The colour-coded Vp/Vs ratio averaged along the whole ray paths at station MSVF and projected onto the foci. a The observed Vp/Vs values. Note a pronounced high Vp/Vs ratio anomaly (red colour) in the northern segment spreading to ~ 600 km depths. b The reference Vp/Vs values of the ak135 velocity model. The face view in 2-D along the latitude from N to S

Full size image
Fig. 11
figure 11

The colour-coded Vp/Vs ratio averaged along the whole ray paths at stations MSVF and RAR and projected onto the surface imaged by the interpolated foci in a 25-km grid. a, b The observed Vp/Vs values. Notice a pronounced high Vp/Vs ratio anomaly (red colour) in the northern segment spreading to ~ 600 km depths for the MSVF station (panel a). c, d The reference Vp/Vs values of the ak135 velocity model. The imaged Vp/Vs values are not local values at the surface, but they represent averaged values along the whole ray paths from foci to the station. The images of the theoretical Vp/Vs values slightly differ for the MSVF and RAR stations, because the ray geometry is different for both stations. The azimuth/elevation of the view angle is 280°/5°

Full size image

Compared to the subduction area (Fig. 11a), the Vp/Vs ratio in the upper mantle eastwards of the slab (station RAR, Fig. 11b) behaves quite differently and it does not display any pronounced anomaly. This is consistent with the observed low-velocity perturbations shown in Fig. 9c,d. However, the comparison of variations of the Vp/Vs ratios at stations MSVF and RAR is not so straightforward, because both stations have different epicentral distance and the Vp/Vs ratio for longer rays (station RAR) might be more smoothed than for shorter rays (station MSVF) due to averaging along the ray path. In order to demonstrate that the absence of visible anomalies at the RAR station is not produced by this smoothing effect, we show histograms of residuals between observed and theoretical Vp/Vs ratios at both stations (Fig. 12). The residuals for the RAR station (Fig. 12b) are artificially amplified by a factor of 2.5 to eliminate a potential smoothing effect produced by roughly twice larger epicentral distance at RAR than at MSVF. The histograms indicate that even the amplified Vp/Vs residuals at RAR are about 3–4 times lower than at MSVF. Hence, the different Vp/Vs patterns at the MSVF and RAR stations in Fig. 11a,b are not the artefacts of data processing but they really reflect remarkably different properties of the upper mantle westwards and eastwards of the slab.

Fig. 12
figure 12

The histograms of the residuals between the observed and theoretical Vp/Vs ratios shown in Fig. 11. The Vp/Vs residuals at station RAR (lower plot) are amplified by a factor of 2.5 to eliminate a potential smoothing effect produced by calculating the Vp/Vs ratios for about twice longer rays than those at station MSVF (upper plot)

Full size image

4.2 Additional Numerical Tests

In order to make sure that the observed high Vp/Vs anomaly in the northern segment is not an artefact of systematic location errors, we performed additional numerical tests. We assumed that the earthquake foci in the latitude range between -20°N and -14°N were systematically mislocated due to a complex geological structure in this area. We shifted the foci by 30 km to the east and recalculated the P- and S-wave velocity perturbations, and subsequently the Vp/Vs ratio based on these new locations. The effect of the systematic shift in the foci locations is demonstrated in Fig. 13. Compared to Fig. 9a,b, the P- and S-wave velocity perturbations for shifted foci (in the north) are different and the negative velocity anomalies disappeared for both types of waves (see Fig. 13a,b). However, if we compare the Vp/Vs ratios calculated for the original and shifted foci (Fig. 13c,d), we find that their differences are almost negligible. The tests were performed not only for systematic foci shifts in the E–W direction but also for shifts in the N–S direction and in depth. In all cases, the Vp/Vs ratio was almost unaffected by the shift, similarly as reported by Conder and Wiens (2006). This points to the reliability and robustness of the retrieved variations in the Vp/Vs ratio.

Fig. 13
figure 13

Numerical tests for location errors. a, b The colour-coded P and S velocity perturbations Δv/v for foci with the latitude > -20°N shifted eastwards by 30 km. c, d The colour-coded Vp/Vs ratio for shifted and original foci. The velocity perturbations and the Vp/Vs ratios are calculated from time residuals observed at the MSVF station. For other details, see caption of Fig. 9. Note an almost identical variation of the Vp/Vs ratio for the original and shifted foci (see (c) and (d))

Full size image

4.3 Depth-Dependent Velocities in the Tonga Subduction Region

To illustrate essential differences in the Vp and Vs velocities and the Vp/Vs ratio between the northern and southern Tonga segments, we calculated the representative P- and S-wave velocities in the subduction region west of the Tonga slab. The traveltime residua were inverted for 1-D depth-dependent velocities with the inversion performed in iterations. In the first iteration, the rays were computed by the TauP code (Buland and Chapman, 1983; Crotwell et al., 1999) in the reference ak135 model. The observed residua were inverted for smooth velocity perturbations in three different sectors bounded by the velocity change in the ak135 model at 210, 410 and 660 km depths. The reference model was updated after each iteration; the iterative process converged fast after 4 iterations with no further improvement.

The inverted representative depth-dependent velocities and the Vp/Vs ratios calculated from traveltime residua (Fig. 14) do not surpass the P- and S-wave velocity anomalies displayed in Figs. 9, 11. However, they illustrate the overall differences in the Tonga subduction region and especially the anomalous behaviour of the northern segment. These velocities, referenced to the ak135 model, show the most remarkable deviations for depths less than 210 km (Fig. 14a). Such deviations attain values of −3% and −5% for the P- and S-wave velocities, respectively, and are also associated with an extremely high Vp/Vs ratio calculated from retrieved P- and S-wave velocities (Fig. 14b). By contrast, the deviations at greater depths are positive, being about 2.5% and 1.5% for the P- and S-wave velocities at depths of 210–410 km, and about 3% for both P- and S-wave velocities at depths of 410–660 km. The deviation of the observed Vp/Vs ratio at depths greater than 210 km is quite small being less than 1%; similarly as the Vp/Vs results for depths greater than 200 km achieved by Conder and Wiens (2006, Fig. 4).

Fig. 14
figure 14

The velocity models in the Tonga subduction region calculated from the time residuals at station MSVF. a, b The whole Tonga back-arc region; (c,d) the northern Tonga segment delimited by the latitude 20°S. a, c The P- and S-wave relative velocity differences from the ak135 reference as a function of depth. Blue line: P-waves; red line: S-waves; dashed black line: P- and S-waves in the ak135 model. b, d The Vp/Vs ratio as a function of depth compared to the reference ak135 model. Magenta line: the retrieved model; dashed black line: the ak135 model. Note a higher Vp/Vs ratio down to 600 km depth in the northern segment (magenta line in (d)). The grey area indicates the depths of eliminated events

Full size image

Lateral variations are observable in Fig. 14c,d showing the velocity differences and the Vp/Vs ratio averaged along the whole ray paths for the northern Tonga segment. Greater velocity differences result in a higher Vp/Vs ratio and, interestingly, high values of the Vp/Vs ratio spread down to depths of ~ 600 km.

5 Discussion

Distinct patterns of velocity anomalies (Figs. 9, 11, SI Movie S1) point to the strong differences between structures affected by the subduction (west to the slab) and the essentially unaffected Pacific plate (east to the slab). The differences are consistent for both P- and S-waves. The subduction region west to the slab (station MSVF) is particularly interesting and is characterized by two opposing zones: (1) a low-velocity zone at depths less than ~ 250 km and (2) a high-velocity zone at greater depths. Since this pattern is not observed on the other side of the slab (station RAR), it must be related to the slab subduction processes. Moreover, this region is laterally inhomogeneous (Figs. 9, 11) with notable anomalies in the north. In our interpretation, we concentrated on isotropic interpretations; however, anisotropy persists in such a complex environment (e.g., Smith et al., 2001; Vavryčuk, 2004, 2006) and can partly influence our results.

5.1 Low Velocities in the Tonga Wedge and Back-Arc

First, pronounced negative velocity anomalies in the mantle wedge and the Lau back-arc region at depths less than 250 km (Figs. 9a,b, 11a) exhibit a broad low-velocity zone extending westwards from the mantle wedge beyond the active Lau Basin. This observation was also confirmed by other authors (Zhao et al., 1997; Koper et al., 1999; Martinez and Taylor, 2002; Stern, 2002; Wei and Wiens, 2018). This zone is explained by the increased temperatures in the upper mantle related to the magma production generating the oceanic crust at the fast-opening Lau back-arc centre.

Second, pronounced negative P- and S-wave anomalies accompanied by changes in the Vp/Vs ratio reaching depth of ~ 500 km are detected around 22°S latitude and coincide with the gap in the seismicity (Fig. 1). Benefiting from the favourable sampling of the structure between station MSVF and the old Louisville Ridge (Fig. 1), these anomalies show the imprints of the old Louisville subduction. Since the area coincides with the northern termination of the double seismic zones at greater depths (Fig. 1b), they can be perceived as a whole upper-mantle contrasting feature.

5.2 High Velocities in the Transition Zone

The observed positive velocity anomalies for the P- and S-waves of ~ 3% at depths of 410–670 km reflect the velocity rise in the back-arc transition zone (Fig. 9a). Similar results were obtained by Bijwaard and Spakman (2000) who used the P-wave inversion and revealed anomalies of + 1.5% relative to the reference ak135 model limited not only to the slab but extending west of the Tonga subduction. The mantle tomography (Hall and Spakman, 2002) revealed high-velocity anomalies of about + 1% stretched westwards of Tonga. However, compared to our results, the global tomographic models are well resolved for larger spatial scales and smeared by the presence of the slab. Our imaging focuses on regions aside of the slab and results in more distinct positive anomalies of ~ 3% with lateral variations.

The transition zone, extending between 410 and 660 km depths, represents the boundary between the upper and lower mantle (Kárason and van der Hilst, 2000). Here, rocks go through a pressure-induced isochemical phase transition from olivine to its high-pressure polymorphs (to wadsleyite at 410 km and then to ringwoodite at 520 km) (Ringwood, 1975). The rocks also become denser due to the transformation of pyroxene to majoritic garnet. However, at conditions within a cold subducting slab, this reaction is extremely slow and the slab remains buoyant relative to the surrounding mantle and stagnates (e.g., Brudzinski and Chen, 2003). On the other hand, a stagnating slab will never float above the 410-km discontinuity since warmer, less dense olivine of the ambient mantle controls it. The eventual descent of slab into the lower mantle is explained by another phase transformation of pyroxene to a high density akimotoite at the base of the transition zone (van Mierlo et al., 2013).

This explains why the velocity anomalies in the transition zone are of different origin than those above 410 km and that they mainly reflect variations in petrology and temperature. The low-temperature and high-velocity patterns indicate a stagnant slab (Chen and Brudzinsky, 2001). The oceanic lithosphere subducting deeper to the mantle flattens, bends and stagnates in the transition zone (Fukao and Obayashi, 2013) or even is detached (Bonnardot et al., 2009) prior to final descending through the 660-km discontinuity into the lower mantle. In addition, modelling studies show that the slab stagnation is often accompanied by the retreat of trench (van der Hilst and Seno, 1993) which also promotes slab flattening in the transition zone (Griffiths et al., 1995; Čížková et al., 2002; Čížková and Bina, 2013).

The presence of a plume, however, can change the topography of the transition zone. The upper discontinuity at 410 km depth represents a phase change from olivine to wadsleyite with a thermally controlled equilibrium. In the Tonga–Samoa region, the 410-km discontinuity was found to be ~ 30 km shallower inside the Tonga slab and ~ 20 km deeper further to the northwest under Fiji. Such a difference might be of thermal origin of the plume material entering the northern Lau Basin (Li et al., 2019). However, deflected topography (Li et al., 2019, their Fig. 3b) is mostly aside the positions of our rays sampling the 410-km discontinuity (Fig. 8) and thus, it has only little impact on our results.

5.3 Contact of the Samoa Plume with the Tonga Slab

The seismic velocities and the Vp/Vs ratio are strongly laterally heterogeneous and differ for the southern and northern segments of the Tonga subduction region (Figs. 9, 11). In the southern segment, the low-velocity anomaly is confined to depths less than 250 km. In the northern segment, the low-velocity anomalies extend to depths of ~ 600 km and are associated with the high Vp/Vs ratio (Fig. 8a, SI Movie S1). The high Vp/Vs found in the low S-wave velocity zone can be associated either with compositional or thermal variations. Since large compositional anomalies cannot be segregated from thermal anomalies (Afonso et al., 2010), this points to the affection by the nearby Samoa plume with a hot material advected from the lower mantle or even the core/mantle boundary. Recent global seismic tomography studies depicted a low‐velocity conduit as a plume beneath the Samoa hotspot (Montelli et al., 2004; French and Romanowicz, 2015; Chang et al., 2016; Maguire et al., 2018). Our results, benefiting from a favourable sampling of the whole upper mantle in the target area (Fig. 1), delimit the lateral and depth extent of the plume-related processes not imaged before.

The obtained results are consistent with other findings. (1) A Tonga slab deformation in the north, demonstrated by distribution of seismicity (Fig. 1b) and focal mechanisms (Giardini and Woodhouse 1984; Vavryčuk, 2006), point to the Benioff zone sharply bending near Samoa. (2) The extensional back-arc Lau Basin reflects intense and complex magmatic processes; its trapezoidal shape induced by extremely fast spreading in its northernmost part (Martinez and Taylor, 2002) evince progressive opening from the north to the south. (3) Extremely high 3He/4He ratios and isotopic Sr–Nd-Pb elements, detected in the ocean island basalt in Samoa and north-western part of the Lau Basin, appear to originate in the Samoa plume (Lupton et al., 2009; Price et al., 2014, 2017; Jackson et al., 2010). The highest geochemical enrichment values (Price et al., 2014, their Figs. 5–7) were detected along the termination of the earthquakes in the northern Lau Basin (Fig. 1b). Price et al. (2014) suggested that the underplated Samoan material was advected into the Lau Basin as the Tonga trench moved to the east via rollback during the past 4 Ma and enabled the hot plume material to entrain. As the trench progressed eastwards, the plume signs in the Lau centre became intense. They were manifested geochemically and by a strong slab deformation induced by plume-related processes near Samoa acting as a backstop barrier.

The high Vp/Vs ratio detected in the northern segment reaches depths of ~ 500–600 km and underlines the whole upper-mantle magma production within the extensional environment enhanced by deeper-seated plume-related processes. These processes were indicated also by Chang et al. (2016) who pointed to the upper- and middle-mantle interactions between the Samoa plume and the Tonga slab from the global anisotropic tomography and 3-D petrological-thermo-mechanical numerical simulations. In their model, the Samoa plume collided with the Tonga slab at the bottom of the mantle transition zone at ~ 4–10 Ma, being both coincident at that time. The collision promoted the fastest slab retreat (roll-back) accompanied by dragging of plume materials into the mantle wedge from the northern rim of the Tonga slab, resulting nowadays in the high Vp/Vs ratio. In a similar way, a high‐resolution mapping of the 410‐km topography under the north Lau Basin was used as a proxy for temperature changes resulting from thermally controlled phase transition from olivine to wadsleyite (Li et al., 2019). A downward deflection of the 410-km discontinuity supports the hypothesis of a hot material mantle entering the place where the descending Tonga slab interacts with the ascending Samoa plume. Since the high Vp/Vs ratio, connected with the low S velocities, images thermal anomalies, it provides the missing evidence of the whole upper-mantle processes at the contact of Samoa with Tonga suggested by other authors (Chang et al., 2016; Li et al., 2019; Smith et al., 2001) as a mantle-plume inflow.

6 Conclusions

Based on processing of a large data set of earthquakes with foci in the Tonga subduction slab recorded at two nearby GSN stations, we provide information on depth and lateral distributions of the Vp, Vs, and Vp/Vs ratio in the upper mantle in the Tonga–Samoa region (Figs. 9a,b, 11a). Even though our study does not provide 3-D tomography images, the obtained results are new and unique, because the previous regional studies (Zhao et al., 1997; Xu and Wiens, 1997; Koper et al, 1999; Conder and Wiens, 2006) did not include the S-wave data in their imaging and/or the resulting structure was mainly confined along the 2-D W–E trending profile across Fiji. Thus, the differences between the northern and southern Tonga segments were not addressed.

Our results image anomalies in the Vp, Vs and Vp/Vs ratio that are crucial for tracing thermal nature of heterogeneities and delineate the areas containing partial melt. We confront these images with geochemical findings on ocean island basalts, the distribution of seismicity, the deformation and extension in the area, and provide the depth extent of the Samoa–Tonga contact.

  1. 1.

    Large low-velocity anomalies of -3% for the P- and -5% for the S-waves, are detected beneath the Lau back-arc region west of the Tonga slab at depths less than 250 km. They represent increased temperatures in the upper mantle related to the magma production generating the origin of the oceanic crust at the fast-opening back-arc centre. The low-velocity anomalies trace a pressure-induced melt extracted in the uppermost mantle wedge reversed and re-introduced in the Lau back-arc centre by the subduction-generated corner flow and indicate a relation of the back-arc opening with the convective mantle wedge circulation.

  2. 2.

    Positive velocity anomalies of ~ 3% for the P- and ~ 1% for the S-waves in the mantle transition zone at depths of 410–670 km point to the retention of high amount of the subducted material. It can be explained by folding of a coherent slab, slab retention or detached slab fragment retained at the mantle transition level (Chen and Brudzinski, 2001; Hall and Spakman, 2002) connected with the fast retreat (roll-back) of the Tonga slab. This is in agreement with modelling studies, where high trench retreat rates advance the slab stagnation and flattening in the transition zone (Griffiths et al., 1995; Čížková et al., 2002; Čížková and Bina, 2013).

  3. 3.

    The most important result, however, is achieved through the detailed Vp/Vs imaging of the subduction region, which points to strong differences between the northern and southern Tonga segments. The high Vp/Vs ratio anomalies in the north extend to depths of 500–600 km. They are associated with the thermal variations and accentuate the whole upper-mantle magma production. Together with a progressive opening of the Lau Basin demonstrated by its south-wedging trapezoidal shape, the Tonga slab deformation in the north demonstrated by the distribution of seismicity and focal mechanisms (Giardini and Woodhouse, 1984), geochemical imprints with extremely high 3He/4He ratios and isotopic Sr–Nd–Pb element enrichments (Price et al., 2014), they trace deeper-seated magmatic processes evoked by the interaction of the Tonga slab with the Samoa plume.

To conclude, our imaging of the upper mantle surrounding the Tonga subduction reveals a prominent Vp/Vs anomaly down to ~ 600 km depth in the northern segment of the region, which is close to the Samoa plume. The occurrence of this anomaly in the vicinity of a plume points to its thermal nature. Such a result provides missing evidence of the whole upper-mantle plume-related processes not imaged before and constrains the active contact of Samoa with Tonga.