A subduction and mantle plume origin for Samoan volcanism

The origin of Samoan volcanism in the southwest Pacific remains enigmatic. Whether mantle melting is solely caused by a mantle plume is questionable because some volcanism, here referred to as non-hotspot volcanism, defies the plume model and its linear age-progression trend. Indeed, non-hotspot volcanism occurred as far as 740 km west of the predicted Samoan hotspot after 5 Ma. Here we use fully-dynamic laboratory subduction models and a tectonic reconstruction to show that the nearby Tonga-Kermadec-Hikurangi (TKH) subduction zone induces a broad mantle upwelling around the northern slab edge that coincides with the non-hotspot volcanic activity after 5 Ma. Using published potential mantle temperatures for the ambient mantle and Samoan mantle plume, we find that two geodynamic processes can explain mantle melting responsible for intraplate volcanism in the Samoan region. We propose that before 5 Ma, the volcanism is consistent with the plume model, whereas afterwards non-hotspot volcanism resulted from interaction between the Subduction-Induced Mantle Upwelling (SIMU) and Samoan mantle plume material that propagated west from the hotspot due to the toroidal component of slab rollback-induced mantle flow. In this geodynamic scenario, the SIMU drives decompression melting in the westward-swept plume material, thus producing the non-hotpot volcanism.

plume-induced volcanoes such as in Hawaii are unlikely to explain the rejuvenated volcanism west of Savai'i because in Hawaii it ceases at most ∼5 Myrs after the shield volcano is born 10 . Hence, an alternative mechanism to the plume model is required to explain the abnormally young volcanism in the western Samoan region, and    Fig. 1a,b were created using the GEBCO_2014 Grid (version 20150318, www.gebco.net) with 0.5 arc-minute resolution and all maps were plotted using the GMT software developed by Wessel and Smith (1991) 60 (version 4.5.14, http://gmt.soest.hawaii.edu). The drawings in tectonic reconstruction maps of Fig. 1b,c,d were done using the GPlates software (version 2.0, http://www. gplates.org). A data bundle provided with the software allowed us to retrace the evolution of the TKH trench, the Pacific subducting plate, the Samoan volcanoes, the Hikurangi Plateau, and the Chatham rise [52][53][54][55] , as well as the Samoan hotspot 1,56-58 , using the global rotation model of Matthews et al. 52 and the Indo-Atlantic moving hotspot reference frame as a reference 59 . possibly the voluminous rejuvenated volcanism on Savai'i and Upolu. To complicate matters further regarding the volcano age distribution, Malulu and Rose, two volcanoes located east of the predicted present-day hotspot, are thought to be older 8 whereas the plume model predicts that these volcanoes should not exist yet but only be built in the future or recently if one considers the relative proximity with the present-day hotspot. Another striking feature of the Samoan region relates to the spatial distribution of volcanoes of which some tend to break the linear trend (Fig. 1a). For example, Uo Mamae and Papatua are clearly offset 120 km south of Upolu and 60 km south of Tutuila, respectively. Interestingly, Uo Mamae experienced volcanic activity at 0.97 Ma 11 , similar to the late rejuvenated volcanism on Upolu. The western Samoan region also seems less linear than the well aligned en-échelon volcanic trails of the eastern Samoan region 8 . Moreover, a particularity of the Samoan volcanism is that in general the volcanic rocks of rejuvenated origin have a distinctive isotopic geochemistry when compared with the shield lavas 9,11 . This was an additional source of motivation that led to the proposal of alternative mechanisms to explain the abnormal volcanic rejuvenation throughout the chain. Suggested mechanisms include multiple consecutive mantle plume sources 7 , and lithospheric flexure and cracking due to the nearby TKH subduction zone in combination with shear melting at the base of the lithosphere 8,9,11 . Proximity of the Samoan region to the northern termination of the TKH subduction zone (Fig. 1b) indeed suggests that subduction related processes could provide an alternative origin for the Samoan volcanism. A first examination of our tectonic reconstruction (see Methods) shows a phase of slab rollback starting at 10 Ma and accelerating at 5 Ma, in remarkable concomitance with the rejuvenation of volcanism throughout the volcanic chain after 5 Ma (Fig. 1c,d).
Quasi-toroidal mantle return flow driven by slab rollback is known to occur around the lateral slab edges of subduction zones [12][13][14][15] , including around the northern edge of the Tonga slab 16,17 . This return flow has a component of upwelling that is produced away from the lateral slab edges in the trench parallel direction, thus in an intraplate setting [18][19][20] . The TKH subduction zone is unique because it is characterised by high trench retreat velocities in the north, reaching 8.2 cm/yr averaged over the last three Myrs as indicated by our kinematic reconstruction, that progressively decrease to the south and reverse to trench advance (Fig. 1b). This asymmetry in trench retreat reflects an along-trench gradient in slab sinking and rollback, which coincides with faster trenchward subducting plate motion, leading to more slab material being consumed in the north as seen on tomographic images 21,22 . Mantle flow that develops in such a subduction setting may differ from the more common symmetrical trench retreat cases. We thus developed fully-dynamic, laboratory-based, subduction models of the TKH subduction zone to simulate asymmetric subduction and to quantify the subduction-induced mantle flow and associated upwelling (SIMU) around the northern lateral slab edge (see Methods and Fig. 2). In our models the kinematic asymmetry is caused by presence of the positively buoyant Hikurangi plateau and Chatham rise born by the Pacific plate in the south of the TKH subduction zone. Because they are positively buoyant, the Hikurangi plateau and Chatham rise resist slab sinking and rollback whereas further north the Tonga-Kermadec slab segment quickly migrates eastward since there is no other resistance to slab motion than viscous mantle drag. We moreover performed a tectonic reconstruction (see Methods) to compare with our modelling results, allowing us to test whether Samoan volcanism can possibly be related to the SIMU occurring north of the Tonga slab.

Results
Geodynamic modelling of Subduction-Induced Mantle Upwelling (SIMU). Out of 10 subduction experiments that were performed, we report here the results of one model (Figs 3, 4 and Table S1) that best reproduces the kinematics of trench motion as determined by our tectonic reconstruction, and best approximates the present-day slab geometry in the upper mantle. The experiment evolves in three subduction phases, as commonly observed in fully-dynamic subduction models, as a result of the interaction between the slab and upper-lower mantle transition zone simulated by the rigid bottom of the tank 13,15 (Fig. 3). The three subduction phases are as follows: (1) An initial slab free sinking phase with progressive acceleration of the subduction process due to increasing volume of negatively buoyant slab material diving in the relatively less dense upper mantle; (2) A slab folding phase due to interaction with the rigid bottom of the tank simulating the 660 km mantle discontinuity; (3) A steady-state slab rollback phase during which velocities become approximately constant over time. During the steady-state rollback phase of our model (phase 3 on Figs 3d and 4a,b), the trench retreats quickly in the north, with scaled velocities reaching 8 cm/yr, and it progressively slows down to the south to reverse to trench advance in the middle of the Hikurangi plateau (Fig. 4c). This trench kinematics is indeed similar to the natural case (Figs 1c,d and 4a). A simplification of our model is the absence of a lower mantle. Despite this simplification, the similarity in trench kinematics indicates a good reproduction of slab sinking and lateral migration in the upper mantle.
In our model, the fast slab rollback in the north produces a quasi-toroidal upper mantle flow around the northern slab edge, with a component of upwelling (Fig. 3). The SIMU is broad, initiates soon after subduction inception and is always present during the subduction process (Figs 3 and 4b). Variations in upwelling rate correspond to the 3 subduction phases (Fig. 4b). We note that the SIMU occurs both outboard of the sub-slab and mantle wedge domain (Fig. 3c,d) whereas in models with uniform trench retreat it is predominantly focused outboard of the sub-slab domain 23 . The asymmetric slab rollback of the TKH subduction zone thus promotes the development of a broad SIMU extending perpendicularly to the trench, with scaled upward velocities averaging 1.2 cm/yr during the steady-state slab rollback phase.
Comparison with tectonic reconstruction. We compare our model with the last 14 Myr of the natural prototype kinematics. We use the slab folding phase (phase 2) of our model as an approximation of the subduction dynamics that occurred between 14-10 Ma because it replicates the dominance of down-dip slab motion over minimal slab rollback (Figs 1d and 5c). We then use the steady-state phase (phase 3) of our model as it We find that the broad SIMU of our subduction model overlaps with all intraplate volcanoes in the Samoa region that are active at 0-1 Ma (Fig. 5a). In particular, areas of high SIMU velocities match very well with volcanoes that were active recently despite their location away from the predicted hotspot, such as Savai'i at 0 Ma and Wallis at ∼1 Ma. The correlation starts to be less clear at 5 Ma (Fig. 5b) and the SIMU is of low magnitude for older times due to very slow slab rollback (Fig. 1d), which thus compares with the slow trench retreat observed during phase 2 of our model (Fig. 5c). These results strongly suggest that the SIMU around the northern slab edge of the TKH subduction zone contributes to causing the intraplate volcanism in the Samoan region after 5 Ma, in particular the rejuvenated volcanism occurring west of the predicted hotspot. The location of the Samoan hotspot after 5 Ma is remarkably distant from the active non-hotspot volcanoes, but it can explain older Samoan volcanism (Fig. 5c).

Discussion
The cause for mantle melting that ultimately leads to formation of Samoan volcanism can be discussed in light of previous research on potential temperature T P for the ambient mantle and mantle plumes. T P is the extrapolation at the Earth's surface of adiabatic temperature-depth profiles, and is determined from the primary melt calculated using forward thermodynamic modelling of adiabatic decompression melting 24,25 . Using the intersection between the peridotite solidus curve and the adiabatic temperature profile at a given T P , one can assess the depth for the onset of mantle melting in the Samoan region. We compute this mantle melting onset depth that corresponds to the primary melt depth and we use it as a proxy to evaluate the efficiency of three different geodynamic  scenarios in producing mantle melting (Fig. 6). We assume a peridotite solidus curve and an adiabatic temperature gradient as determined in previous research 24 , and we use estimated lower and upper bounds of 75 km and 100 km 26 for the depth of the lithosphere-asthenosphere boundary (LAB) below the Samoan region. We use a T p of 1340 ± 60 °C for the ambient mantle (T P (MOR) in Fig. 6a ) as constrained from estimates at mid-oceanic ridges 27 since decompression melting there results from passive upwelling with no initial temperature anomaly. In contrast, T P estimates for the Samoan mantle plume (T P (PLU) in Fig. 6b,c) range between 1395 °C and 1524 °C 28 .
As can be expected, the Samoan mantle plume is hotter than the surrounding ambient mantle, and the excess  temperature is in agreement with T P estimates for other plumes 29 . We assume an adiabatic gradient of 0.55 °C/km, which is an average of the values ranging between 0.45-0.7 °C/km that are commonly used in thermodynamic studies 24,30 . Calculation of the dry peridotite solidus is done using the following equation, which was determined experimentally 24  with P(T ) S the pressure in GPa at a given temperature T S in °C, and T S 0 the extrapolation of temperature at the surface, which is 1,100 °C. We consider a dry peridotite solidus because the Samoan volcanism is located away from the arc volcanism, thus hydration processes should be minimum. In contrast, we infer that wet peridotite solidi should be considered to study the adakite-like volcanism that occurs as a result of interaction between mantle upwelling and slab possibly carrying sediments. Such volcanism does not occur in the Samoan region but does occur in the mantle wedge close to the northern slab edge of the TKH subduction zone 31 . The process of SIMU (Fig. 6a) occurs in the ambient mantle, hence with the lowest T P without bringing any thermal anomaly. The onset of mantle decompression melting in the asthenosphere is possible via this process only if we consider minimum values for the LAB and maximum values for T P (orange arrow and star in Fig. 6a). In sharp contrast, the Samoan mantle plume (Fig. 6b) brings a heat excess that allows the onset of melting in the asthenosphere below the LAB, and thus consistently explains the origin of intraplate volcanism occurring near the predicted hotspot. We propose that before 5 Ma, the Samoan mantle plume is the only process responsible for the intraplate volcanism because our tectonic reconstruction shows that volcanic activity always occurred near the predicted hotspot location and far from the northern Tonga subduction segment (Figs 1d and 5c). The plume might also continue to produce Samoan volcanism near the predicted hotspot location after 5 Ma.
To explain the non-hotspot volcanism occurring west of the hotspot after 5 Ma (Fig. 1c), we suggest that the SIMU interacts with Samoan mantle plume material to promote decompression melting (Fig. 6c). Noting that active non-hotspot volcanism became progressively more distant west of the predicted Samoan hotspot between 1-5 Ma (Figs 1c,d and 5a,b), we propose that Samoan mantle plume material has been progressively dragged west of the hotspot after 5 Ma by the toroidal component of subduction-induced mantle flow. Toroidal return flow has already been suggested around the northern slab edge of the Tonga slab by azimuthal anisotropy 17 , consistently with trench-parallel anisotropic fast directions in the sub-slab domain 32 . Such toroidal mantle flow is known to strongly deform the mantle plume tail and head 33 , and it has been suggested for the Samoan plume using geochemistry and modelling 16,31,34 . Variations of the isotropic shear velocity at 250 km depth north of the Tonga slab are also in agreement with the suggestion that hot material from the Samoan plume is swept westward 35,36 . Considering the onset time of the interaction between the subduction-induced toroidal component of mantle flow and the Samoan plume at around 5 Ma, the propagated plume material may have originated from the Samoan plume head that may have been present at that time. Alternatively, the dragged plume material could come from plume tail material ponded underneath the base of the lithosphere. Our geodynamic subduction model demonstrates that the rollback-induced return flow around the northern edge of the Tonga slab has a strong upwelling component over a broad area. We propose that this SIMU is the driver of decompression melting in the hot westward-swept Samoan plume material after 5 Ma.
Our two-phase geodynamic scenario reconciles a range of observations and geochemical data. It satisfactorily explains the occurrence of intraplate volcanism due to the Samoan hotspot at all times and also of non-hotspot volcanism west of the hotspot after 5 Ma. It is moreover consistent with the geochemistry of shield and rejuvenated lavas from the regions that have a common Samoan mantle plume source 5 despite the great distance of the rejuvenated lavas from the hotspot. It also accounts for the different isotopic signatures of the shield (hotspot) and rejuvenated (non-hotspot) lavas that could thus be due to the two different mechanisms triggering decompression melting in the upper mantle. The triggering of decompression melting in the westward-swept Samoan plume material by SIMU provides a new alternative origin to intraplate volcanism that can be explored in other comparable tectonic settings. An interesting parallel can be done with the triggering of decompression melting in the Hawaiian plume head by upwellings resulting from small-scale lithospheric convection that might produce volcanism that is too young to be attributed to the Hawaiian hotspot 37 . Another type of interaction that has been proposed is between mantle plume, mid-oceanic ridges and background mantle flow to explain volcanism occurring far away from the Réunion 38 and Tristan Da Cunha 39 hotspots. In this case, decompression melting occurs near the mid-oceanic ridges where the oceanic lithosphere is thinner.

Methods
Fully-dynamic laboratory subduction models. The laboratory modelling builds on previous studies that use a fully-dynamic approach to study the subduction process 13,20,23,40,41 . The models are self-consistent and driven by buoyancy forces only, except that we initiate subduction by creating a small slab perturbation of about 3 cm, scaling to 150 km, at the beginning of the experiments. The driving force is the slab negative buoyancy and the resistive forces are assumed to be dominated by viscous resistance and drag in the ambient mantle, slab, and lithospheric plate. The fully-dynamic (buoyancy-driven) modelling approach allows us to investigate the mantle flow that results solely from the dynamics of slab sinking and rollback. In this approach, the rheological response of the lithosphere and asthenosphere is assumed to be dominantly viscous over geological timescales. Our scaling procedure follows previous subduction modelling studies that use the Stokes' settling law to describe the velocity of the sinking slab 40,42,43 , defined as follows: where v is the slab sinking velocity, C is a constant, ρ ∆ is the density contrast between the slab and the ambient mantle, l is a characteristic length, g is the gravitational acceleration, and η is the dynamic shear viscosity of the sublithospheric upper mantle. To achieve dynamic similarity our scaling uses identical force balance between model (subscript m) and natural prototype (subscript p), such that:  41 , with a density of 3,200 kg · m −3 for the upper mantle in nature, a density of 1,420 kg · m −3 for the glucose syrup in the models, and density contrasts between the slab and mantle of 80 kg · m −3 in nature and 97 kg · m −3 in the models (see details in paragraph below). With these numbers, we calculate a sublithospheric upper mantle viscosity p η of ~2.19 × 10 20 Pa·s, which is in the range of natural estimates of 10 19 -10 21 44 . Another criterion that allows us to achieve dynamic similarity is the use of a small Reynolds number, estimated between 7.03 × 10 −6 and 2.53 × 10 −4 . Such a small Reynolds number ensures that the models are in the laminar symmetrical flow regime and guarantee the dominance of viscous resistive forces over negligible inertial forces 43 .
The models involve four linear-viscous (Newtonian) layers contained in a rectangular tank that is 140 cm long and 140 cm wide (see Fig. 2). One layer is made of transparent glucose syrup and fills the tank up to 13.2 cm to simulate the sub-lithospheric upper mantle, the rigid bottom of the tank approximating the 660-km discontinuity. The glucose syrup has a density of 1,420 kg·m −3 and a viscosity of 202 ± 8 Pa·s at 20 °C. Three layers allow to simulate the Pacific subducting lithosphere that carries the Hikurangi plateau and the Chatham rise to the south of the TKH subduction zone. A 100-km-thick subducting oceanic plate is simulated using a 2-cm-thick layer of silicone (Wacker silicone, from Dow Corning company) mixed with fine iron powder. A 100-km-thick plate is a good first order approximation of the Pacific subducting plate that shows variation in age along the TKH trench, with ages of ∼85 Ma near the Osbourne Trough that increase in both directions to reach ∼125 Ma at the southern termination and ∼110 Ma at the northern termination of the subduction zone 45 . At these ages, 100 km is an upper bound for the depth of the lithosphere-asthenosphere boundary (LAB) because of seafloor flattening at ages greater than 70 Ma 26,46 . Moreover, for oceanic lithosphere older than 80 Ma, geophysical estimates of the Gutenberg discontinuity seem to indicate a LAB at 75 km on average 26 . The subducting plate layer is 72 cm (scaling to 3,600 km) in the trench-parallel direction and 60 cm (scaling to 3,000 km) in the trench-normal direction, leaving a distance of 34 cm between either lateral edge of the slab and the sidewalls to minimise sidewall effects 23 . It has a density of 1517 kg·m −3 and a viscosity of 6.32 ± 0.1 × 10 4 Pa·s, as measured using an Anton Paar Physica MCR301 rheometer. The subducting plate-upper mantle density contrast of 97 kg·m −3 reflects natural conditions of a 100-km-thick oceanic lithosphere 47 , although it is slightly higher in the models (17 kg·m −3 ) to negate surface tension effects that are negligible in nature, and it makes the subducting plate negatively buoyant 23 . Both the trailing edge and the lateral edges of the subducting plate are free, representing a mid-oceanic ridge and strike-slip faults, respectively, that offer negligible resistance to subducting oceanic plate motion.
The Hikurangi plateau, an oceanic plateau formed of mafic igneous crust 48 , is simulated using a mixture of Wacker silicone and Dow Corning 3176 Dilatant compound (referred to as pink putty) in weight proportion of 10.66% and 89.34%, respectively. The Wacker silicone is Newtonian with a viscosity of 4.8 × 10 4 Pa·s. The pink putty has a more complex power-law behaviour but it approximates a Newtonian rheology with a viscosity of 6.2 × 10 5 Pa·s at strain rates lower than 10 −3 s −1 , which are similar to those occurring in our models 49 . The Hikurangi plateau layer has a density of 1,120 kg·m −3 , which corresponds to a density contrast of 300 kg·m −3 in relation to the glucose syrup, and reflects an upper mantle-oceanic plateau density contrast as estimated in nature with an average oceanic crust density of 2,900 kg·m −3 50 . The Chatham rise, a continental fragment, is also simulated using a mixture of Wacker silicone and pink putty but in weight proportion of 72.25% and 27.75%, respectively. The Chatham rise layer has a density of 1,020 kg·m −3 , which corresponds to a density contrast of 400 kg·m −3 in relation to the glucose syrup, and reflects an upper mantle-Chatham rise density contrast as estimated in nature with an average continental crust density of 2,800 kg·m −3 51 .
To quantify the upper mantle flow produced by slab sinking, we use a stereoscopic Particle Image Velocimetry (sPIV) technique that allows to compute the 3 components (3 C) of the mantle flow velocity field in a plane. The glucose syrup layer simulating the upper mantle is randomly filled with 20-50 μm diameter polymethylmethacrylate-rhodamine B phosphorescent particles that are illuminated with a horizontal laser sheet at 5.5 cm depth (275 km). More details on the sPIV method applied to subduction models can be found in a previous publication 23 . Two high-resolution cameras (2,046 × 2,046 pixels) placed above the model allow to compute ScIentIfIc REPORTS | (2018) 8:10424 | DOI:10.1038/s41598-018-28267-3 the mantle flow velocity field in map view. The 3 C mantle flow velocity field was computed with a spatial resolution of 1.75 pixels/mm, a seeding density of ~40 particles/cm 2 , an interrogation window of 64 × 64 pixels with an overlap of 50%, and a time lapse between two images of 45 s. After computation and post-processing, a few artefacts remain on the vertical component of the 3 C velocity field because the sPIV system may detect motion of particles other than the phosphorescent markers, such as tiny bubbles anywhere in the glucose syrup. These artefacts cannot be removed as they have the same magnitude as the signal we want to image, but they are easy to distinguish from it as they are of smaller scale than, for instance, the broad scale upwellings (Fig. 3). These artefacts should not be taken into account when interpreting the mantle flow velocity field images. Two other cameras allow to take pictures from one side, allowing to follow the slab geometry in the north, and from the front (mantle wedge) side, allowing to see the slab geometry along the trench. An additional camera placed above the model allows to take pictures when the laser does not emit, to record pictures in map view where passive tracers placed on the plate are seen. The subducting plate velocity is measured using passive markers equidistantly placed both on top of and at the sides of the subducting plate. The trench retreat velocity is measured using the colour contrast between the glucose syrup and the subducting plate.
We performed a set of 10 experiments where we varied the thickness of the Hikurangi plateau and Chatham rise, and also compared results with models that do not include them. Some experiments were repeated to ensure reproducibility (see Table T1 in Supplementary Information). We then selected one model that best reproduces the kinematics of trench motion and slab geometry at depth of the Tonga-Kermadec-Hikurangi subduction zone, to discuss results focusing on the subduction-induced mantle flow around the northern lateral slab edge.
Tectonic reconstruction. The tectonic reconstruction maps (Figs 1b-d, 5 and Fig. S1 in Supplementary   Information) were performed using the GPlates software (version 2.0, http://www.gplates.org). A data bundle provided with the software allowed us to map the TKH trench, the Pacific subducting plate, the Samoan volcanoes, the Hikurangi Plateau, and the Chatham rise [52][53][54][55] . A data bundle was also used 1,56-58 , from which the location of the present-day Samoan hotspot was determined. The predicted position of the above geological features back in time was then reconstructed using the global rotation model of Matthews et al. 52 and the Indo-Atlantic moving hotspot reference frame as a reference 59 . In this reference frame, the Samoan hotspot moves slowly southward between 10-0 Ma (<0.65°, Fig. 1c) and in the east-west direction between 14-10 Ma (<0.4°).
Data Availability. The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.