## INTRODUCTION

Measurements of seismic anisotropy have opened up ways of dynamically imaging Earth’s mantle, showing how it flows as part of deep mantle convection, plume activity, and the motion of the tectonic plates (*1*). A key assumption here is that there is a relation between this anisotropy and the direction of mantle flow, mainly determined by the lattice-preferred orientation (LPO) in olivine aggregates (*1*, *2*). The seismically fast direction often defines the shear direction in melt-free olivine that has undergone simple shear (*1*), although this is not always the case, depending on the stress state, temperature, pressure, and water content of the olivine (*2*). Five distinct LPO fabrics (A to E types) for an olivine mineralogy are commonly associated with seismic observations (*2*), observed in laboratory experiments and samples of mantle rocks such as the ultramafic parts of ophiolite complexes, and xenoliths in volcanic eruption or those dredged from the ocean floor (*3*, *4*). However, samples of mantle rocks have also yielded an additional LPO fabric, referred to as the AG fabric (fig. S1), which has not previously been identified from direct seismic observations of the upper mantle (*4*).

The AG fabric in forsterite-rich olivine (Fo_{90}) aggregates is characterized by a girdle of high P- and S-wave speeds (8.7- to 8.8-km/s P-wave and 5.1-km/s S-wave speeds), with a slower speed orthogonal to this (7.8-km/s P wave and 4.7-km/s S wave) (fig. S1). These high P-wave speeds are notably different from the global average for (continental) P_{n} of 8.09 ± 0.2 km/s (*5*) and are higher than could be accounted for by measurement uncertainties (*5*). On the basis of the locations of where samples with this fabric have been found, it has been suggested that it may form locally in trench and mantle wedge environments and near the base of the lithosphere, possibly where nonsegregated melt phases are present (*4*, *6*). Here, using seismic refraction techniques, we present the first in situ evidence for this fabric in the upper mantle. We show that it occurs on a regional (~1000 km) scale in the shallow mantle beneath a large igneous province (LIP) and can be simply explained in terms of the anticipated late-stage strain field in a large mantle plume undergoing gravitational collapse. This way, a regional-scale AG fabric may be a diagnostic fingerprint of ancient, and now dispersed, mantle superplumes.

## RESULTS

Our new seismic data were acquired during active- and passive-source seismic refraction experiments across two fragments of LIP in the southwest (SW) Pacific, comprising the Manihiki (MP) and Hikurangi (HP) oceanic plateaux (Fig. 1). Plate tectonic reconstructions indicate that these, together with the Ontong-Java Plateau (OJP), originally formed an even larger oceanic plateau, referred to as the Ontong-Java-Manihiki-Hikurangi Plateau (OJMHP), created over about 2.5 million years (Ma) in a brief period of intense and prolific igneous activity at about 120 Ma ago (*7*, *8*). Melting in a hot mantle plume is widely considered the most viable explanation for the distinct chemistry of the lavas (*9*, *10*), but the anomalous subsidence history of the Ontong-Java fragment requires additional factors that remain controversial (*11*).

The western edge of the HP is now subducting along the Hikurangi Margin beneath eastern North Island, New Zealand (*12*–*14*). Previous active- and passive-source seismic experiments along the length of the Hikurangi Margin (*15*–*17*), for northeast-SW (NE-SW) ray paths, identified unusually high mantle P-wave speeds in the range of 8.7 to 9.0 km/s at depths of 30 to 70 km within the HP (fig. S2). Velocity models derived from seismic tomography, using dominantly northwest-southeast (NW-SE)–oriented ray paths, also yield P-wave speeds as high as 9 km/s (*13*), suggesting that these high speeds are not just confined to ray paths that traverse a single azimuth. New refraction data, described below and acquired as part of the Seismic Array Hikurangi Experiment (SAHKE) (*12*) (Figs. 1B and 2), show how these high wave speeds are part of a more regional and consistent seismic anisotropy.

SAHKE 04 was an onshore ~85-km-long NW-SE line in the SAHKE experiment (*12*), deploying 872 vertical seismographs spaced at ~100 m and 291 three-component instruments spaced at about 300 m. Twelve large explosions (located in drill holes) were used as seismic sources to image the crustal structure of the accretionary wedge where the underlying Pacific plate dips on average ~9°NW (Figs. 1B and 2A) (*14*, *18*). The subducted oceanic crust has a thickness of 10 ± 1 km, forming the western edge of the subducting HP (*12*, *14*, *19*). The SAHKE 04 line was too short to allow P_{n} wave speed determinations using active-source techniques, but several offshore earthquakes recorded by the array during the 2 weeks that it was operating made it possible to make simultaneous measurements of both P- and S-wave speeds from refracted phases beneath the subducted crust (Figs. 1B and 2).

GeoNet earthquake #3511915, which occurred ~42 km offshore and approximately on the line of the profile at a depth of 50 ± 3 km, was recorded by three-component seismometers (Fig. 2, A and D). A linear fit to first arrivals at stations on basement rock for P waves that have traveled through the uppermost mantle, after correction for the 9° ± 1°NW dip of the Moho (figs. S3 and S4), yields a speed to 8.8 ± 0.2 km/s. We confirmed this result with a full ray-tracing model, which includes the crustal structure, constraining the uppermost mantle P-wave speeds to be 8.7 to 8.8 km/s for rays traveling across the Hikurangi Margin (Figs. 1B and 2A), identical within error to those traveling at right angles along the margin (fig. S2) (*15*, *17*).

High uppermost mantle P-wave velocities are also indicated by a two-dimensional (2D) interpretation of two ocean bottom seismometer (OBS) gathers (OBS 14 and 15) for an offshore extension of SAHKE 04, referred to as SAHKE 01, comprising an array of OBS and airgun seismic sources (Figs. 1B and 2B) (*12*). Here, a ray-tracing model indicates a speed of 8.8 ± 0.2 km/s for P_{n} phases (Fig. 2B and figs. S5 and S6) in the gently dipping upper mantle of the HP both within and to the SE of the subduction zone (Fig. 2B). This is within the range of 8.1 to 8.6 km/s reported in an analysis of the SAHKE 01 line (*20*), but we note also that mantle velocities are likely to be variable here because they may become perturbed in cracks where the slab bends into the trench because of infiltration of water (*21*).

S waves with unusually high mantle speeds were recorded on all three seismograph components across the SAHKE 04 array, but they are most clearly seen in the horizontal components (Fig. 2, E and F, and fig. S7). We identify two mantle S-wave arrivals (Fig. 2, E and F); one is weak, yet persistent, over the first 30 km of the array, whereas the second is stronger and can be traced right across the array. Their true speeds are estimated to be 5.1 and 4.7 km/s (±0.15 km/s), respectively, from ray tracing constrained by the crustal structure model based on previous work (*14*) (Fig. 2D). We interpret these to result from a split of the S_{n} phase due to about 10% radial anisotropy (fig. S7): The slower S wave is well defined in both the vertical and radial horizontal seismograph components (fig. S4), indicating that it is the vertically polarized phase (S_{v}). The faster S wave has a lower amplitude, particularly for rays that have come through the 2 to 3 km of stratified and wet Late Cenozoic sediments east of kilometer 35 (fig. S4), as would be anticipated for the higher attenuation of the horizontally polarized phase (S_{h}) when the ray path steepens toward the surface (*22*). There are also slower S-wave arrivals right across the SAHKE 04 array with a dip-corrected speed of 3.6 ± 0.2 km/s, which are interpreted as the direct S waves through the oceanic crust.

We also analyzed uppermost mantle P-wave speeds beneath the MP, which forms another fragment of the OJMHP (figs. S8 and S9). Here, OBS receiver gathers from an active-source seismic experiment (*23*), on two lines at roughly 60° to each other (lines 0100 and 0200, fig. S8), indicate apparent upper mantle horizontal P-wave speeds of 8.79 ± 0.1 km/s for line 0100 and 8.76 ± 0.07 km/s for line 0200 (table S1). These estimates are based on averaging results for 33 OBS gathers that show clear P_{n} phases (fig. S8 and table S1). The effect of local Moho dip on the upper mantle P-wave velocities can be corrected for by sorting these gathers into pairs where reversed seismic profiling has sampled the same part of the upper mantle, yielding within error the same horizontal P-wave speeds of 8.86 ± 0.1 and 8.77 ± 0.03 km/s for lines 0100 and 0200, respectively (fig. S8 and table S1). This result is confirmed by ray-tracing models for four OBS gathers along line 0200, again requiring upper mantle P-wave speeds for subhorizontal rays of 8.9 ± 0.1 km/s (fig. S9).

Seismic velocities beneath the OJP are not as well resolved as those in our analysis of the MP and HP but are consistent with our results from these other plateaux. For example, a Rayleigh-wave dispersion study for the center of the OJP shows high sub-Moho S-wave speeds in the upper mantle of 4.75 km/s (*24*), similar to S_{v} speeds observed along the SAHKE 04 array, as described above. An active-source survey of the OJP from the 1970s (Fig. 1) reported unusually high P_{n} speeds from beneath the plateau around 8.6 km/s (*25*), but the uncertainty is at least ±0.2 km/s given the sparse spacing of receivers. We note that mantle xenoliths that may be derived from deep beneath the OJP (*3*) are sampling a deeper region than our seismic observations, because we are using P_{n} and S_{n} phases that are refracted in the top 10 to 20 km of the mantle.

## DISCUSSION

When all the seismic observations of uppermost mantle P-wave speeds are plotted on the reconstructed OJHMP (*8*), there appears to be no obvious correlation with azimuth, given a minimal ±0.1 km/s uncertainty in each wave speed determination (Fig. 3, A and D). We conclude that there is negligible azimuthal anisotropy (<1%) in the shallow mantle beneath the OJMHP over a horizontal distance greater than 1000 km. The Observed P- and S-wave speeds and polarizations, do however, match those for aggregates of olivine (Fo_{90}) with an AG fabric under uppermost mantle conditions (fig. S1) (*4*). We rule out eclogite (*26*), because the observed P-wave speeds and their anisotropy are too high (*27*), and geochemical and petrographic studies show that it cannot be the source of the lavas (*9*).

Numerical modeling of aggregates of forsterite-rich olivine subjected to uniaxial shortening under upper mantle conditions, using experimental data on the activity of all the possible slip systems, has robustly reproduced the strong symmetry and seismic velocities of the AG fabric, as would be anticipated from the symmetry of the imposed strain field (*4*, *28*). This indicates that the regional AG fabric that we have documented beneath the OJMHP could be explained in terms of a 3D strain field, on the 1000-km scale, that is dominated by vertical compaction and isotropic horizontal extension. If the OJMHP was the product of a superplume, then our numerical and analog models show that exactly this 3D strain regime would dominate in the shallow part of the plume head during the later stages, in cases with or without significant interaction with the overlying lithosphere (Fig. 3, B and C, and figs. S10 and S11).

In the early stage of plume head development, when the plume head is fed by steady-state vertical flow through the plume tail, previous studies have predicted that an axially symmetric Poiseuille-like flow develops midway through the plume head (Fig. 3B) (*29*, *30*). We reproduce this flow fabric with a 3D finite-element model (figs. S10 to S12) for the evolution of a finite-volume, buoyant mantle upwelling. A key property of the model is that as the material flux in the plume tail declines (Fig. 3C and figs. S10 to S12), the flow field in the vicinity of the top surface of the plume switches to one of vertical compaction, with equal radial and tangential extension, giving rise to the AG fabric (fig. S1). Subsequently, as the rocks cool, the late-stage AG fabric is frozen in and preserved over geological time.

Experiments show that in the presence of melt, the AG fabric may develop in olivine even if there is a component of simple shear, because under these conditions, multiple slip systems can be activated (*6*), potentially extending the regional AG fabric to a wider part of the plume head over a longer period of its development. We note, however, that although simple shear is likely to be a widespread feature of mantle flow, the AG fabric has not been observed before, supporting our conclusion that it is mainly formed under the special conditions of uniaxial flattening in a spreading plume head. Moreover, melt extraction in the plume head will also drive compaction of the residue mantle, intensifying the AG fabric at the observed horizontal length scale. It is likely, however, that the fabric is only preserved in the largest of plumes, the so-called superplumes (*7*), because the high rate and volume of upwelling and melting of mantle is required to dwarf the effect of shearing on the plume head due to the absolute plate motion (*31*). In this way, the flow field remains axially symmetric (Fig. 3D) and preserves the distinct AG fabric on length scales of 100 to 1000 km. Last, the coherency of our observed seismic fabric indicates that large LIPs can be the product of a single superplume, rather than a cluster of smaller plume heads (*32*), with a more complicated spatial distribution of strain fabrics.

We conclude that the AG fabric with its characteristic high P- and S-wave velocities may be an important feature of Earth’s uppermost mantle, created by late-stage collapse and spreading of large mantle plume heads. The AG fabric stands in contrast to the more common simple shear–related fabrics that arise from plate tectonics (*28*). Identifying this fabric elsewhere in the world is therefore important as it may help to quantify the distribution and extent of fossil superplumes, providing better constraints on the role of these plumes in mantle convection and the evolution of Earth. For example, the high sub-Moho P-wave seismic speeds of 8.7 to 9 km/s reported for the lithospheric mantle beneath the Siberian traps (*27*), constituting the second largest known LIP (*7*, *33*), also suggest a superplume origin.

## MATERIALS AND METHODS

### Seismic velocity anisotropy of type AG fabric

Seismic velocity anisotropy that may be produced by type AG fabric for an aggregate of forsterite-rich olivine can be extracted from its elastic properties. We used the elastic tensor form of this AG fabric within Christoffel equations to solve for anisotropic P- and S-wave velocities in all seismic propagation directions (*34*). This AG tensor, as described below, was oriented in accordance with radial plume flow. In fig. S1, we illustrate the resulting *V*_{p}, plus the fast and slow *V*_{s} draped on spheres that represent all azimuth and inclination directions, calculated for a pressure of 20 kbars and a temperature of 1100°C (*35*). We also show the percent anisotropy for *V*_{P} and *V*_{s} based on this expression where *V*_{average} = (*V*_{max} + *V*_{min})/2, where percentage anisotropy = (*V*_{max} − *V*_{min})/*V*_{average}.

The *V*_{p} sphere illustrates the radial symmetry of the P-wave anisotropy, with *V*_{p} slower in the vertical direction (7.81 km/s) in contrast to the horizontal direction (8.65 to 8.77 km/s). In the horizontal direction, the faster *V*_{s} is in the range of 4.99 to 5.11 km/s, and its wave particle motion direction is horizontal (red bars on *V*_{s}-fast sphere). In contrast, the slower *V*_{s} ranges between 4.65 and 4.68 km/s in the horizontal direction with vertical polarization motion. Correspondingly, the largest amounts of *V*_{s} % anisotropy is in the horizontal direction.

The AG fabric elastic tensor for an aggregate of olivine (Fo_{90}) under uppermost mantle conditions (*P* = 20 kbars and *T* = 1100°C) has the following tensor elements in Voigt index notation and is diagonal symmetric (*35*): *C*_{11} = 250.88, *C*_{22} = 258.19, *C*_{33} = 204.73, *C*_{44} = 73.46, *C*_{55} = 72.50, *C*_{66} = 83.82, *C*_{12} = 79.51, *C*_{13} = 74.71, *C*_{14} = −0.030, *C*_{15} = 0.088, *C*_{16} = −1.050, *C*_{23} = 74.34, *C*_{24} = 1.760, *C*_{25} = 0.130, *C*_{26} = 1.430, *C*_{34} = 0.0600, *C*_{35} = −0.0320, *C*_{36} = 0.130, *C*_{45} = −0.020, *C*_{46} = 0.052, and *C*_{56} = 0.043.

### Seismic experiments at the Hikurangi Margin

Estimates of upper mantle P-wave speeds for trench-parallel azimuths along the Hikurangi Margin (fig. S2) have been reported over the past 35 years (*15*–*20*). Most of these studies involve using earthquake waves as an energy source, although an active-source seismic study was conducted in 1991 as a joint project between Victoria University of Wellington, Leeds University (United Kingdom), and the Canadian Geological Survey (*15*). The latter used six large borehole explosions, recorded with 86 seismographs over a length of 300 km (fig. S2). The key finding of the active-source study was P-wave speeds of 8.7 ± 0.1 km/s in the upper mantle at depths greater than 32 km (fig. S2).

The SAHKE project was a joint New Zealand–Japan–USA experiment in 2010 to 2012 to image the subduction zone beneath the southern North Island of New Zealand (*12*). SAHKE consisted of both onshore and offshore seismic profiling across the subducting margin using active- and passive-source techniques. For the 90-km-long onshore segment of SAHKE 04 (Fig. 1B), across the lower North Island, a dense array of 874 vertical and 291 three-component seismographs was used, spaced at 50 to 100 m. These recorded 12 borehole explosions to image the deep mantle structure of the Pacific plate (*19*).

### Subduction zone structure at Hikurangi Margin

A detailed analysis of the southern Hikurangi subduction zone interface (megathrust) beneath SAHKE 04, using seismic energy from onshore explosions, together with onshore-offshore shooting, shows that the megathrust increases in dip along the profile, from about 5°NW at the SE end to 14° at the NW end (*14*). The dip of the megathrust in the midsection beneath the SAHKE 04 profile is 9° ± 1°NW. Wide-angle reflections and refractions image the Moho of the Pacific plate (i.e., the HP) beneath SAHKE 04, giving an oceanic crustal thickness of 10 ± 0.5 km (*14*).

The crustal structure in the overlying plate has been determined in previous analyses of seismic arrivals along the SAHKE array (*12*, *14*, *18*, *19*). This includes a number of shallow sedimentary basins and a low-velocity (5 km/s) 1- to 3-km-thick underplated sediment layer at the plate interface, with a large (7 to 8 km thick) hump of underplated sediments at the NW end of the line (*12*, *14*).

### Upper mantle P-wave speed beneath SAHKE 04

Several earthquakes were recorded during the 2-week period when the onshore SAHKE 04 array was operating (fig. S4). The most useful of these was the moment magnitude 2.5 event (GeoNet event #3511915) with a hypocentral depth of 50 ± 3 km, and located 42 km offshore to the NW, almost directly in line with the NW extension of the SAHKE profile (Fig. 1). This event was well located by 16 stations of the GeoNet network using 19 phases and had a standard error for location of 0.3 s (www.geonet.org.nz/earthquake/technical/3511915). P- and S-wave arrivals from the earthquake were recorded across the array (fig. S4).

For the first 33 km of the line, all stations are on Mesozoic greywacke basement rock, and the arrivals form roughly linear travel time plots. Beyond kilometer 35, the profile crosses the Wairarapa where thick Paleogene to Pliocene sedimentary basins up to 2.2 km exist (*36*). The effect of the low wave speeds in these basins is visible in the short wavelength (kilometer scale) features of the travel time plots, giving rise to marked changes in apparent wave speed, but an overall linear trajectory is evident across the whole array.

A least-squares fit to picks across the SAHKE 04 line yields apparent upper mantle P- and S-wave speeds of 10.15 ± 0.06 and 5.30 ± 0.08 km/s, respectively (table S1). These speeds are apparent rather than real for two structural reasons. First, the sedimentary sequences occur along the SE half of the line, delaying the arrivals, biasing downward the apparent speed. Counteracting this is the effect of the NW dip on the subduction zone interface, which biases upward the wave speeds. We can eliminate the effect of sedimentary basins by measuring the speeds at distances of 3 and 31 km along the array, where the stations are all on basement. This yields apparent P- and S-wave speeds of 11.1 and 5.7 (± 0.07) km/s, respectively.

Our first approach to resolving the upper mantle P-wave speed beneath the SAHKE array is a simple analytical approximation, where we assume that the wave paths approximate critically refracted waves, which travel up-dip along the interface of a dipping half-space (fig. S3A). This assumption is reasonable given the subhorizontal trajectory of the mantle ray paths, evident in the full ray tracing (Fig. 2A). Analytical solutions (*37*) exist for the simple approximation of a dipping half-space with speed *V*_{2} overlying a layer of lower speed *V*_{1,} and the dip of the layer *V*_{2} is δ (fig. S3A)

(1)where *V*_{2up} is the apparent up-dip wave speed and θ is the critical angle. The solution for *V*_{2} is then

(2)

Our first estimate for apparent up-dip P-wave speed (*V*_{2up} in Eq. 3) across the first 31 km of the array is 10.15 km/s (table S1). But the low wave speeds in the sedimentary basins at the far end of the profile will bias the apparent velocity downward. Accordingly, we make a delay-time correction for the sedimentary basins along the profile and recompute the apparent wave speed to be 11.1 ± 0.06 km/s (table S1). Confidence for this corrected apparent wave speed is gained because an identical wave speed is obtained (table S1) if we measure the apparent wave speed stations over the first 30 km of the line at the northwestern end, which are located on basement greywacke. The shape and velocity structure of the basins are taken from previous studies (fig. S2).

We solve Eq. 1 for *V*_{2up} = 11.1 km/s and a plausible range of values of *V*_{1} (5.5 to 5.9 km/s) and dip angle (9° ± 1°NW, fig. S3). The largest source of uncertainty is the dip of the interface (fig. S3). The associated predicted range of *V*_{2} is 8.6 to 9.2 km/s or an average of 8.85 km/s. The average speed is consistent with our previous determination of mantle P-wave speed from ray tracing and with the values measured on margin parallel (NE-SW; see above) profiles (*15*).

We next take a forward modeling approach by building a ray-tracing model (Fig. 2A) that uses GeoNet event #3511915 as a seismic source. Ray tracing is done interactively and is based on ray theory (*38*). This approach has allowed us to include the dip of the subduction zone and the previously determined crustal structure in the overlying plate (*12*, *14*, *18*). We searched for a best-fit P-wave speed for the upper mantle by considering four upper mantle velocity models that have uppermost mantle wave speeds of 8.1, 8.3, 8.7, and 8.9 km/s, respectively. A vertical gradient of 5 × 10^{−3} s^{−1} is imposed so that the waves turn and come back to the surface (that is, the speed increases downward by 0.3 km/s from top to bottom of ~80-km-thick plate). We showed that having mantle P-wave speed gradients in the range of 2 to 10 × 10^{−3} s^{−1} made no significant difference to travel-time fits.

The ray tracing (Fig. 2A) shows that the waves turn within a vertical zone 5 to 10 km thick so P-wave speeds no greater than the 0.05 km/s above the initial wave speed are encountered. Figure 2C shows observed P-wave arrival picks compared with calculated model values. The earthquake energy source is located at the subduction interface, and the radiating seismic energy forms nearly horizontal turning waves that mimic P_{n} and S_{n} critically refracted waves (Fig. 2A). As we are modeling the wave speed, and not absolute travel time, we arbitrarily fix the calculated travel times at the northwestern end of the profiles and compare their slope on the reduced travel-time plot. From visual inspection, the model with the upper mantle P-wave speed of 8.7 km/s at the top gives the best fit of the four trial models. We considered a range of models with different choices of the wave speed at the top of the mantle, ranging between 8.1 and 8.9 km/s. This indicates that the maximum uncertainty in the P-wave velocity determination for our models is ±0.2 km/s, consistent within the uncertainties in the analytical methods discussed above. Thus, our “best” estimate for the speed of upper mantle and subhorizontal P-wave rays in the HP, regardless of azimuth, is 8.85 ± 0.2 km/s.

### Upper mantle P-wave speeds beneath SAHKE 01 with OBS gathers

High upper mantle P-wave speeds are observed from gathers of seismic arrivals recorded by OBS for the offshore south-eastern part of SAHKE 01 line (Fig. 1B) (*12*). A multichannel seismic (MCS)reflection profile was shot along the SAHKE 01 line, and the interpretation of that (*39*) was used to constrain the crustal structure of our interpretation (Fig. 2B). The MCS data show the top of the Pacific plate at about 10 km beneath the sea floor but cannot resolve structure beneath that. OBS gathers 14 and 15 (fig. S5) were shot in opposed directions, which gives control on the dip of interfaces, but not necessarily on overlapping segments. OBS 14 and 15 are shown at reduction speeds of 6 and 8 km/s, respectively (fig. S5). Both gathers show a triplication in the travel time branches at an offset of ~60 km from the OBS, which is typical of a thin, medium-speed layer sandwiched between a slower and faster layer (*37*). A simple 2D ray-tracing solution (Fig. 2B), using the bathymetry and upper sedimentary structure (*39*), shows that the data are consistent with a regular Moho at about 24-km depth beneath a thin (3.5 km) layer of high-speed (7.4 to 7.6 km/s) lower crust (fig. S6A). However, we could also change the thin-layer P-wave speeds to be 8 to 8.1 km/s (i.e., regular upper mantle speed) and get the same fit. This latter interpretation, of regular mantle P-wave speeds overlying faster ones, would be in line with that of Chadwick (fig. S2). We tried fits with P_{n} speeds that varied from 8 to 9.2 km/s to the OBS 14 data (fig. S7, B to E), and the best fit is a P-wave speed of 8.8 ± 0.2 km/s, with a vertical speed gradient of 3.75 × 10^{−3} s^{−1}. The uncertainty in P_{n} is because the section of P_{n} picks only occurs on a short segment of about 40 km long.

### Upper mantle polarization of S waves at the Hikurangi Margin

Earthquake #3511915 generated clear S waves. We found the optimal display for these data is in the band-pass range of 0.1 to 4 Hz (fig. S4). We also applied static corrections for elevation and then plotted against trace number (Fig. 2E). We identify various phases. The slowest that we interpret are S waves that have traveled in the oceanic crust. Arrivals for a faster S wave indicate an apparent speed 5.7 km/s. We also identify an even faster, and much weaker, S-wave phase with an apparent speed of ~6.6 km/s with arrivals that can only be traced to kilometer 33 along the profile.

Ray tracing, using the same subduction zone structure as our P-wave model, but with velocities reduced to appropriate S-wave speeds, indicates that the faster phases have traveled in the upper mantle with true speeds of 5.1 and 4.65 km/s and uncertainties of ±0.15 km/s (Fig. 2D). We interpret these two phases as the horizontally polarized S_{h} (5.1 km/s) and the vertically polarized S_{v} (4.65 km/s) The wave speeds are consistent with those predicted for the AG fabric, as discussed in the main text (fig. S1). A feature of the S-wave seismogram is the dominance in amplitude of the S_{v} phase over that of the S_{h} phase at all offsets (fig. S6). This could represent a dominance of S_{v} in the earthquake source or that the dominate energy is created by a converted P to S_{v} at a regular interface near the source (*37*). S_{h} energy cannot be traced east of kilometer 33 where layered shallow sedimentary sequences would be predicted to selectively attenuate S_{h} energy that would be coming to the surface at a high angle (*22*).

We also checked our interpretation of S_{v} and S_{h} by rotating the seismogram data from the field orientations into the transverse and radial components with respect to the azimuth of the ray path to the seismic stations (fig. S6). The first arrivals on the transverse component become sharper in the offset range of 3 to 18 km, which is consistent with the first arrivals being S_{h}. But beyond offsets of 35 km, the radial component becomes sharper and stronger, and this is consistent with our interpretation that S_{v} is more dominant at these distances. In reality, we can expect both S_{v} and S_{h} on all three components as both the source and ray path through the crust will have complexities, such as crustal anisotropy, which will cause S_{v} and S_{h} to be seen on all three components.

A check on the validity of our S-wave anisotropy analysis is to measure the splitting times at a well-determined point along the ray. The best locality for this is at about kilometer 33, where there is a 1-s split (fig. S6 and Fig. 2, E and F). The relationship between this split time (Δ*t*) and the vertical (S_{v}) and horizontal (S_{h}) shear wave speeds (*V*_{sv} and *V*_{sh}), given the length of trave path in the mantle is (Δ*x*) is

(3)

For *V*_{sh} = 5.1 km/s (see above), then *V*_{sv} = 4.65 km/s, and Δ*t* = 1 s, and the predicted travel distance (Δ*x*) in the mantle is 52 km. From our ray-tracing model (fig. S8), the total travel path for arrivals at kilometer 33 is 63 km. This would imply that *V*_{sh} is 5.02 km/s to satisfy Eq. 3, well within the ±0.15-km/s uncertainty of our estimate of *V*_{sh}.

### P-wave speeds in the MP

The crustal structure of the MP was investigated in 2012, with active-source seismic profiling along two ~450-km-long lines (fig. S8A) (*23*, *39*). Although the published interpretations focus on the crustal, rather than mantle, structure (*23*), the OBS receiver gathers reveal high apparent P-wave speeds in the upper mantle of the MP (fig. S8). We examined 60 gathers (*40*), identifying P_{n} phases on 33 of them where the apparent wave speed is >8.2 km/s. These results are summarized in table S1, where we tabulate average apparent P_{n} speeds for each line. Although the averages for the two lines are roughly the same (8.86 and 8.77 km/s for lines 0100 and 0200, respectively), the SDs of the estimates vary between the two lines (0.42 km/s for line 100 and 0.26 km/s for line 0200). This difference is likely to be an effect of the markedly smoother ocean floor along line 0200 compared to that along line 0100. Irregularities in the ocean floor would be predicted to introduce noise due to scattering of seismic energy, giving rise to increased uncertainty in estimates of apparent wave speed.

For seismic refraction data collected on reversed lines, the true wave speed for any layer is close to the average of those speeds obtained by shooting from each direction (*37*). We use this property to better define the true wave speeds on each line by selecting pairs of P_{n} apparent wave speeds that are measured in opposing directions and are spaced at the appropriate distance such that the P_{n} waves, from different directions, are sampling the same part of the upper mantle. This way, we identify nine pairs of data for each line, and their results are summarized in table S1. The true P_{n} estimates (with standard error, assuming that true speed is the same for all pairs) are now 8.86 ± 0.1 and 8.77 ± 0.03 km/s for lines 0100 and 0200, respectively.

We made a simple ray-tracing model of the OBS gathers at stations 7, 15, 23, and 29 for line 0200 (fig. S9). P_{n} is seen as first breaks on all gathers, constraining the P_{n} speed to be 8.9 ± 0.1 km/s, for a crustal thickness that varies between 24 and 27 km. We also plotted the travel-time curves for P_{n} wave speeds of 8.3, 8.5, 8.9, and 9.2 km/s, calculated for the OBS #29 gather (fig. S9). Our visual best fit upper mantle P-wave speed for subhorizontal rays is 8.9 ± 0.1 km/s, identical within error to those determine along line 0100 and in the HP (Fig. 2A).

### Plume model

We solve Stokes equation in 3D for a spherical buoyancy anomaly at depth, which develops into an axisymmetric plume, rising and spreading out beneath the lithosphere under gravity. We nondimensionalize and simplify the equations under the assumption of cylindrical symmetry, valid when the plume ascent rate is greater than any surface plate motions. This reduces the number of dimensions in the computational model to 2, while producing a solution that is strictly equivalent to the 3D axisymmetric case. We consider both end member cases for the near-surface boundary condition, with either a melt layer at the top of the plume (i.e., free slip) or a rigid zero-displacement boundary condition. As this is a dynamic model, we do not need to assume any internal boundary conditions in the mantle, such as those required by kinematic models (*29*), which a priori enforce a Poiseuille flow regime. All quantities scale with the layer thickness, *L*, and the buoyancy flow time scale, τ, representing the e-fold time for a mantle instability to evolve (*41*)

(4)and τ = *t*/*t*′ where *t* is absolute time and *t*′ is dimensionless time. For the model of Fig. 1 (A and B), *t*′ = 80 and 300, respectively. Thus, for adopted values of η = 1 × 10^{20} Pa s, Δρ = 30 kg/m^{3}and *L* = 1 × 10^{6} m, the respective time steps are 1 and 4.4 Ma.

We assume a Newtonian rheology, as we are interested in the long-time behavior of the system. The finite-element Lagrangian “particle-in-cell” code Underworld (https://www.underworldcode.org/) is used to solve the Stokes equation, assuming incompressibility. Under the axisymmetric assumption, a Cartesian mesh (coordinates, radial distance from plume center and depth) is used with model domain length to height ratios of 4:1 and a uniformly spaced element resolution of 1024 × 256. Each element is populated with at least 32 particles, representing a maximum particle spacing of 2 km.

We tested plume head collapse in an analog model using “silly putty”; this is also effectively a test of the role of rheology because silly putty has a non-Newtonian viscosity, whereas our numerical simulations were for a Newtonian fluid. Photographs were taken using a Canon EOS 60D equipped with an EF100mm f/2.8 L Macro lens and Speedlight 580EX II, ISO 200 with

$\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$10$}\right.$-s exposure at f/16. Underlying graph has 1-mm grid spacing. A time-lapse video produced with time dilation factor of 60 (i.e., 1 s/min) can be seen here: www.dropbox.com/s/cgb8042xhty8v7i/PuttyTimelapse.mp4?dl=0.

The analog model replicated the late-stage flow patterns in the numerical models, indicating that gravitational collapse and spreading is a robust feature of plume head evolution, independent of rheology.