Intraplate volcanism originating from upwelling hydrous mantle transition zone

Most magmatism occurring on Earth is conventionally attributed to passive mantle upwelling at mid-ocean ridges, to slab devolatilization at subduction zones, or to mantle plumes. However, the widespread Cenozoic intraplate volcanism in northeast China1–3 and the young petit-spot volcanoes4–7 offshore of the Japan Trench cannot readily be associated with any of these mechanisms. In addition, the mantle beneath these types of volcanism is characterized by zones of anomalously low seismic velocity above and below the transition zone8–12 (a mantle level located at depths between 410 and 660 kilometres). A comprehensive interpretation of these phenomena is lacking. Here we show that most (or possibly all) of the intraplate and petit-spot volcanism and low-velocity zones around the Japanese subduction zone can be explained by the Cenozoic interaction of the subducting Pacific slab with a hydrous mantle transition zone. Numerical modelling indicates that 0.2 to 0.3 weight per cent of water dissolved in mantle minerals that are driven out from the transition zone in response to subduction and retreat of a tectonic plate is sufficient to reproduce the observations. This suggests that a critical amount of water may have accumulated in the transition zone around this subduction zone, as well as in others of the Tethyan tectonic belt13 that are characterized by intraplate or petit-spot volcanism and low-velocity zones in the underlying mantle. The widespread intraplate volcanism in northeast China and the unusual ‘petit-spot’ volcanoes offshore Japan could have resulted from the interaction of the subducting Pacific slab with a hydrous mantle transition zone.

The Cenozoic intraplate volcanism in northeast China is located more than 1,000 km westward of the Japan Trench 1 , while the young alkaline basalts (0-6 Ma) known as petit-spots outcrop up to 600 km eastward of the trench 4 (Fig. 1). The formation mechanism of these types of onshore and offshore volcanism is still debated, as there is no geological and geophysical correlation with mantle plumes or arc volcanism 6,14 .
Seismic tomography models indicate that in this region the Pacific Plate is currently stagnant in the mantle transition zone (MTZ), extending continuously up to nearly 1,000 km to the inland of northeast China 3,8,10 . Thus, it has been proposed that the Cenozoic intraplate magmatism is related to the dehydration of the Pacific slab in the MTZ 2,15 .
The primary petit-spot magma has been determined to be volatilerich with extremely enriched mantle (EM1-like) isotopic compositions 6,7 . The lack of hotspot tracks in this region excludes a contribution from a mantle plume. It has been postulated that the petit-spot magma forms in the asthenosphere and migrates upward through the oceanic lithosphere by reactive porous flow in response to plate flexure 4,6,7 .
Based on electrical conductivity surveys, the MTZ probably holds about 0.1 wt% water 16 . The MTZ below northeast China and Japan is particularly wet, with at least 0.5-1 wt% water 17 . The MTZ is primarily composed of wadsleyite and ringwoodite minerals that can accommodate 1-3 wt% water, which is 1 to 2 orders of magnitude higher than the water (hydrogen) solubility in upper-and lower-mantle minerals. Given the large contrast in water solubility between the MTZ and upper/ lower mantle, it is reasonable to expect deep dehydration melting when subducting slabs excite vertical flow in the nearby wet MTZ 18 . Indeed, seismic low-velocity zones (LVZs) above 410 km and below 660 km have been observed not only in Japan [8][9][10][11][12]14,15 , but also around subduction zones in Europe 19 and the western United States 20,21 .
To test this hypothesis, we construct two-dimensional numerical experiments in which a self-sustained oceanic plate subduction is characterized by trench retreat and slab stagnation into a homogeneously or heterogeneously wet MTZ (see Methods). The subducting plate and entrained dry upper mantle push the adjacent wet MTZ downward to the lower mantle such that a partially molten layer forms between 700 km and 800 km depth (Fig. 2a, region labelled M2) (Supplementary Video 1). On the other hand, MTZ material uplifted to the upper mantle starts to partially melt above 410 km (Fig. 2a, M3). Slab stagnation and retreat is accompanied by sub-slab MTZ upwelling and new melting (M1). These partially molten regions above and below the MTZ cause large seismic LVZs (Fig. 2c). When melt percolation is active (see Methods), extraction to the surface occurs, forming intraplate and petit-spot volcanisms ahead of and behind the trench, respectively. Figure 3 shows the spatial and temporal trend of modelled volcanics for the reference model in Fig. 2. The first intraplate volcanism occurs about 500 km away from the trench, then spreads in two opposite directions. The mantle water content decreases after melt extraction, which precludes further deep (≥200 km) melting of the residual peridotite (Fig. 2b). As the slab rolls back, more distal wet MTZ is sucked into the upper mantle wedge, such that partial melting and volcanoes will form further away from the trench. The new generated volcanism is not homogeneously distributed as it is strongly influenced by mantle flow and trench movement. Furthermore, a heterogeneous distribution of water in the MTZ would prevent the formation of any temporal-spatial magmatic sequence as wetter portions would melt earlier than drier regions at the same pressure and temperature (P-T) conditions. It is noteworthy that intraplate volcanism also occurs a few hundred kilometres in front of the slab tip. After about 12 Myr of modelled subduction, petit-spot volcanoes appear behind the trench. They are located up to about 300 km seaward of the trench and exhibit a similar magmatic activity trend to the intraplate volcanism.
We further test the influence of initial water content in the transition zone and other parameters on the genesis of asthenospheric melting (Extended Data Fig. 5). Melting commences 40 km above the transition zone, and no petit-spot volcanoes are formed for 0.2 wt% initial water. The thickness of the partially molten layer could range from a few tens of kilometres to more than 100 km, depending on the melt extraction efficiency and water content. A petit-spot volcano might be located more than 600 km from the trench if the melt extraction process is efficient (Extended Data Fig. 5b). Given the assumed homogeneous distribution of water in the MTZ, these models provide upper-bound estimates on the volumes of volcanics and melt. However, the results also hold for a more realistic heterogeneous distribution of the water in this mantle level (Extended Data Fig. 5e).
When comparing the model results with seismic and geological observations, we note that around the Pacific slab, three remarkable seismic LVZs outside the transition zone are clearly imaged (Fig. 1b). These are well correlated with the locations of intraplate and petit-spot volcanoes, and the modelled partially molten zones (Figs. 1a and 2). Although seismic low-velocity anomalies are generally attributed to thermal effects 14 , or to the presence of water 22 , melt 11,20 and/or major element compositional heterogeneities 16,23 , it has recently been argued that some of these LVZs could be artefacts induced by seismic anisotropy 24 . Nevertheless, the authenticity of the sub-slabs LVZ1 and LVZ2 appearing in tomographic models has been confirmed by other independent studies using, respectively, an accurate scrutiny of the seismic ray paths 25 sampling the LVZ1, and receiver functions in the case of LVZ2 11 . The LVZ3 sits below the active Changbai volcano and appears to extend down to 410 km as revealed by multiple high-resolution tomography models 3,8,15,26 . A thermal anomaly from a non-hotspot upwelling, if it hypothetically exists, is difficult to reconcile with the large velocity drop of LVZ1. The hot material will rapidly cool when flowing upward, because of adiabatic decompression and the latent heat of the wadsleyite-to-olivine reaction. Laboratory experiments show that seismic wave speeds are insensitive to moderate (<1 wt%) water contents for olivine 27 and wadsleyite 22 ; thus, the LVZs are very likely to be caused by partial melting and/or compositional heterogeneities. The presence of basalts at the bottom of the upper mantle can be excluded as it would generate a positive seismic anomaly 28 . On the other hand, basalts accumulating at the base of the MTZ could be effectively dragged by the slab into the uppermost lower mantle and generate the LVZ2. However, receiver functions indicate that the lower-mantle LVZs are within the 750-780 km depth range 11 , which is likely to be below the post-garnet phase transition where basalts are seismically faster than mantle rocks.
The presence of melt in the deep mantle, which is mostly catalysed by the involvement of volatiles, decreases seismic velocities and provides a magmatic source for intraplate/petit-spot volcanism. Our numerical models thus suggest that a hydrous transition zone with at least 0.2-0.3 wt% water beneath northeast China and offshore Japan can

Article
comprehensively explain the LVZs and the intraplate and petit-spot volcanism. This model does not exclude the devolatilization of the stagnant Pacific slab as a mechanism to explain the LVZ3 region and the overlying intraplate volcanism 15 , which favours the upwelling of volatile-rich plumes from the MTZ 29 as envisaged by the Big Mantle Wedge model 15,30 . However, the same slab-derived volatiles cannot obviously be the cause of both the LVZ1 and LVZ2 and of the petit-spots, implying the presence of a metasomatized MTZ before the last subduction episode. The accumulation of water in the MTZ could be caused by, for example, delamination of volatile-rich lithospheric roots 31 , or by previous slab dehydration episodes in the MTZ and subsequent absorption of the water by wadsleyite and ringwoodite. Alongside with water, reduced (by redox-freezing) carbonated sources and restitic K-hollandite-bearing sediments are required to explain the volatile-rich, alkaline and EM1-type petrological and geochemical signature of the basalts 32,33 . This is not surprising, as the MTZ, a graveyard for stagnating slabs, is the most likely candidate to host volatiles and subducted sediments, and long-term isolation of these MTZ domains would be consistent with the ancient metasomatizing episodes estimated for intraplate basalts [32][33][34] . Subsequent subduction events would mobilize the wet and (carbon + alkali)-bearing MTZ rocks, promoting the formation of silica-undersaturated magmas in the upper mantle. It is important to note that the addition of these components is not critical to our results, because the location and amounts of partial melting above and below the MTZ will still be dictated by the distribution of wet MTZ domains, while reduced carbonated sources are expected to experience redox melting at shallower depths (<250 km) 35 .  The process proposed here could potentially explain also the Cenozoic anorogenic volcanism in the Mediterranean 13 and intraplate volcanism in the Turkish-Iranian Plateau 36 regions characterized by the long-term subduction of the Tethys Ocean. Together with surface intraplate/petit-spot volcanism, constraints on deep seismic low velocities and/or high electrical conductivity may thus indicate a volatile-rich and/or partially molten mantle within and around the transition zone.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-020-2045-y.

Modelling approach
The 2D petrological-thermomechanical numerical code I2VIS used in this study is based on a finite difference method using a marker-in-cell technique on a staggered grid 37 . It solves mass, momentum and energy conservation equations (1)-(3) on the Eulerian grid and interpolates physical properties to the markers for advection accordingly. v x where v i is velocity, x i coordinate, P dynamic pressure, ρ density, g gravity acceleration, c p heat capacity, T temperature, k thermal conductivity, H r radioactive heating,

Model configuration
The initial model set-up (6,000 × 1,000 km discretized with 1,501 × 501 nodes) is composed of a 3,500-km subducting plate and a 2,500-km overriding plate. The model imposes free-slip mechanical boundary condition at the top with 30-km-thick and viscosity of 10 18 Pa s 'stickyair' to mimic free surface; the bottom boundary is no slip, and side boundaries are periodic. The bottom no-slip condition is needed to define an initial horizontal velocity from which finite differences can be computed for this variable. Comparison with results from a model with a bottom free-slip condition and closed vertical walls indicate that the bottom no-slip boundary condition does not affect the subduction dynamics at all, as it is confined above the lower mantle. The initial thermal structure is defined by the half-space cooling age for the plates (50 Myr old) and an adiabatic thermal gradient of 0.5 K km −1 for the underlying mantle. The thermal boundary conditions are isothermal on the top and bottom, while side boundaries are periodic, consistent with the mechanical boundary conditions.
To initiate subduction, the subducting slab extends down to about 200 km in the upper mantle together with a rheologically weak zone on top of it which lubricate the initial contact between the plates. The high numerical resolution (4 km × 2 km) used here is needed to ensure plate contact lubrication at shallow depths and localized, bendingrelated hydration at the trench outer-rise. Tests at a lower resolution (4 km × 4 km) result in less-localized slab mantle hydration, whereas with a resolution of 8 km × 4 km, self-sustained subduction and slab rollback do not appear spontaneously.

Viscous-plastic rheological model
The rock mechanical behaviour is represented by the effective viscosity, which combines ductile (dislocation, diffusion and Peierls creep) and brittle (Drucker-Prager) deformation. The effective ductile viscosity is given by the harmonic average of the combined rheologies (parameters and physical meaning are defined in Extended Data Table 1): where the dislocation and diffusion creep are given by 38 : For hydrated (wet) mantle, viscosity is reduced by and C w , C w0 are water content and reference water content (100 ppm, which is the water content for the dry upper mantle), respectively. The Peierls creep η Peierls is given by 39 :

Petrological modelling
Petrological solid-solid phase changes are included through the density and enthalpy look-up tables for basalt and pyrolite obtained from PERPLE_X 40 . Therefore, phase transition boundaries at 410 km and 660 km have been considered. The solidus (T s = f(P, T, H 2 O)) and liquidus (T l = f(P, T)) temperatures for the upper mantle and MTZ are taken from high-pressure experiments 41 (Extended Data Fig. 1). At lower mantle conditions, T s and T l vary considerably among different experiments. Here we adopt the dry solidus and liquidus of chondritic mantle 42 , as these are more compatible with the results of KLB-1 peridotite 43 , while the wet solidus 44 was measured on samples with an estimated water content of 400 ppm wt.
A conservative estimate of the melt fraction in the wet upper mantle 18 is applied: where W a , W ol and W m (10 wt%) are water mass fraction of the ambient mantle, olivine and melt, respectively. Note that the water solubility in olivine increases with pressure, so the melt fraction will decrease with depth if W a remains constant (Extended Data Fig. 3). The silicate melt density (Extended Data Fig. 2) is taken from highpressure sink-float experiments 45 , which show that the melt becomes denser than the surrounding mantle at around 400 km (refs. 45,46 ) owing to the increased compressibility. However, the presence of water generally reduces the melt density such that it becomes buoyant relative to solid mantle (see Extended Data Fig. 2), rendering melt extraction at this depth possible. We also test the melt density from molecular dynamics simulations at high-pressure conditions 47 (Extended Data Figs. 4f, 5d).

Melt extraction timescale
The distance over which the compaction rate decreases by a factor of e is the characteristic length scale of the compaction process and is known as the compaction length, δ c : c 4 3 f where ζ and η are the effective bulk and shear viscosities, respectively, of the partially molten rock; η f is the fluid viscosity; K is the permeability given by the empirical equation: where ϕ is the porosity (melt fraction); K 0 (10 −12 m 2 ) is the permeability at the reference porosity ϕ (0.01) 0 ; and n = 3. The relative migration velocity between the melt and the solid matrix is w: Indeed, a small migration time implies a large ∆ρ (that is, high buoyancy force) and/or high melt fraction (that is, high permeability) and/or weak solid matrix which can be easily deformed during compaction/decompaction processes. In this study, we use a reference value t ref = 6 kyr.
Melt extraction depends not only on these three parameters (and fluid viscosity), but also on the dihedral angle (that is, melt interconnectivity). Previous experiments showed that the dihedral angle decreases systematically with increasing pressure such that it probably allows for complete wetting at about 400 km (ref. 48 ) where the dihedral angle is <5° (ref. 49 ). At these depths, melt interconnectivity is high even for low amounts of melt (<1%), which makes melt extraction possible provided that the extraction timescale is sufficiently low (that is, the melt migration process is efficient and not hindered by the solid matrix). As such, when the extraction timescale t is smaller than t ref , the material is extracted at the surface forming plutonic intrusions or volcanics 50 . Note that when the density of the melt is higher than that of the solid surrounding mantle, as occurs between 11.5 GPa and 13.5 GPa when the water content is low 45 (Extended Data Fig. 2b), there will be no melt extraction. In these conditions, the denser melt should percolate downward and accumulate over the 410-km discontinuity 18 . However, dry melting generally does not occur at ambient mantle conditions, except when there is an abnormal heat source associated with a mantle plume. As hydrogen partitions preferentially into the melt, the water content in the melt would be quite high, decreasing its density 45,46,49 . As a result, hydrous melt should be less dense than the solid matrix throughout the upper mantle 45 .
The melt migration process is illustrated here with more realistic models accounting for visco-elastoplastic deformation in a two-phase flow regime. These models demonstrate that melt migration from the deep upper mantle to the surface should occur through several mechanisms: viscous diapirism, viscoplastic decompaction channels and elastoplastic dyking 51,52 (Extended Data Fig. 6). For weak host rocks where viscous deformation dominates, such as the asthenosphere, magma migrates by diapirism. When the magma moves through the lithosphere-asthenosphere boundary (or the lower crust in continents) where both ductile and brittle deformation occur, the fluid compaction pressure might reach the tensile strength, and magma could migrate by channelling. If the host rock is completely elastoplastic, such as the core of lithospheric mantle and upper crust, magma migrates by dyking.

Water budget
The phase diagram reporting the maximum water content that can be hosted in hydrated or wet mantle rocks (that is, absorbed by nominally anhydrous minerals, NAMs) is built upon the compilation from refs. 41,53 (Extended Data Fig. 3). It is often assumed that heterogeneously serpentinized mantle rocks below the oceanic Moho can contain up to 2 wt% H 2 O (refs. 54,55 ). As such, the maximum water content in hydrated rocks 53 is scaled accordingly.
As the oceanic crust completely dehydrates at about 300 km depth 41 , generating fluids that fuel arc-volcanism only, we assume a dry crust for the sake of simplicity. On the other hand, dehydration of the underlying mantle within the transition zone is thought to cause intracontinental magmatism 15 . Consequently, we allow for serpentinization by bending-related deformation when the strain of mantle rocks is greater than 0.1 (ref. 55 ).
When the rock water content exceeds the saturation limit, decomposition of hydrous minerals or water exsolution in NAMs occurs, and fluid markers are generated and migrate according to Darcy's law 55,56 until they are absorbed by dry markers: ρ s , ρ f are the densities of solid and fluid, respectively; V 0 is a constant percolation velocity; g is the gravity acceleration vector as defined in equation (2), and g y is its vertical component. Upon partial melting and extraction, the water is partitioned into the extracted melt and water in the residual peridotite as: w res w w melt where D = 0.01 is the hydrogen partition coefficient for olivine polymorphs.
Falling block tests. The validity of the petrological model used here can be easily tested with a simple model in which a falling block (simulating the subducting slab) sinks into the wet MTZ, exciting wet upwellings to the upper mantle and squeezing water into the lower mantle (Extended Data Fig. 4). These tests indicate that the melt layer gets thicker (>100 km) when melt extraction is not efficient, owing to very small amounts of melt/water and/or denser melt phase. This might explain the thick low-velocity layers above 410 km in many regions 9 . After melt extraction, less water remains above the transition zone, causing higher viscosity and less melt fraction, which yields a larger extraction timescale: that is, the melt preferentially ponds above 410 km depth.
Seismic velocity anomalies. The seismic velocity perturbation in Fig. 2c have been computed as:

Article
The change of seismic wave velocities caused by the existence of a fluid phase is given by: with A ss , A sl the area of solid-solid contact and solid-liquid contact, respectively 58 .
Extended Data Fig. 7 shows K b /k and N/μ for the equilibrium geometry model at various dihedral angles.

Data availability
The dataset generated during the current study is available at https:// figshare.com/articles/Yang_Faccenda_Nature2019/9933056.