Evidence of Volatile‐Induced Melting in the Northeast Asian Upper Mantle

A seismic low velocity layer (LVL) above the mantle transition zone (MTZ), often thought to be caused by volatile‐induced melting, can significantly modulate planetary volatile cycles. In this work, we show that an LVL observed beneath northeast Asia is characterized by small, 0.8 ± 0.5 vol%, average degrees of partial melting. Seismically derived P‐T conditions of the LVL indicate that slab‐derived CO2 , possibly combined with small amounts of H2 O, is necessary to induce melting. Modeling the reactive infiltration instability of the melt in a stationary mantle above a stalled slab, we demonstrate that the volatile‐rich melt slowly rises above the stalled slab in the MTZ, with percolation velocities of 200–500 μ m/yr. The melt remains stable within the LVL for this geologically significant period of time, potentially transferring up to 52 Mt/yr of CO2 from the subducting slab to the mantle for an LVL similar in areal extent ( 3.4×106km2 ) to the northeast Asian LVL. Reaction between the melt channels and the LVL mantle precipitates up to 200 ppmw solid C in localized zones. Using the inferred small melt volume fraction to model trace element abundances and isotopic signatures, we show that interaction between this melt and the surrounding mantle can over the long‐term produce rocks bearing a HIMU like geochemical signature.

SUN ET AL.

10.1029/2021JB022167
2 of 17 wave triplications (Gao et al., 2006;Han et al., 2021;Song et al., 2004), and seismic tomography (Obayashi et al., 2006). Marked by a sharp interface with the overlying mantle and a 2%-6% reduction in shear wave velocity, the LVL is generally recognized as a chemical, rather than thermal, anomaly. A number of recent rock physics analyses suggest that the LVL is characterized by an average of 0.5-1 vol% partial melt embedded in a silicate matrix (Freitas et al., 2017;Hier-Majumder & Courtier, 2011;Hier-Majumder & Tauzin, 2017;Hier-Majumder et al., 2014;Xiao et al., 2020). While it is plausible that the LVL is a global phenomenon (Tauzin et al., 2010), the highest resolution of its internal structure is obtained from regional studies, such as in the western US (Eagar et al., 2010;Fee & Dueker, 2004;Hier-Majumder & Tauzin, 2017;Jasbinsek & Dueker, 2007;Jasbinsek et al., 2010;Schmandt & Humphreys, 2011;Song et al., 2004;Tauzin et al., 2013) and northeast Asia (Bagley et al., 2009;Han et al., 2021;Revenaugh & Sipkin, 1994;Liu et al., 2016). Especially useful in this respect are seismic results obtained from receiver function studies that, in addition to the reduction in shear wave velocities, provide information regarding the thickness of the MTZ (Chevrot et al., 1999;Lawrence & Shearer, 2006;Tauzin et al., 2008). Exploiting the Clapeyron slopes of phase transitions of olivine polymorphs, the MTZ thickness can be converted into mantle temperature, thus providing an additional set of constraints on the physical properties of mantle aggregates. These conversions indicate that the average temperature of the LVLs is far below the solidus of dry peridotite but is closer to the solidi of carbonated basalt and hydrated peridotite (Hier-Majumder & Tauzin, 2017; Petrological and geochemical evidences provided by sub-lithospheric diamonds (Harte, 2010;Shirey et al., 2021;Walter et al., 2008) and Cenozoic basalts from eastern Asia (Qian et al., 2020) record volatile subduction beyond the volcanic front to the transition zone and lower mantle (Shirey et al., 2021). Carbon and oxygen isotopic anomalies in the sublithospheric diamonds (Burnham et al., 2015;Ickert et al., 2015;Thomson et al., 2014), and Mg isotopic anomaly in basalts (S. Li et al., 2017) demonstrate that slab-derived, volatile-rich fluids interact with the surrounding mantle, producing diamonds and magmas bearing a slab signature. Clinopyroxenes from late Cenozoic Wudi lavas in eastern China reveal the presence of a HIMU (high 238 U/ 204 Pb ratio) component in the source of these lavas. Geochemical analyses also reveal that these basaltic melts are likely derived from partial melting of carbonated peridotites (Qian et al., 2020). While these petrological and geochemical observations require a reservoir of volatile-rich melts in the mantle, and the rock physics analyses suggest the LVL as a potential candidate, the efficiency of the LVL in sequestering volatiles still remains debated. Additionally, the seismic evidence of melting provides a snapshot in time, lacking constraints on the temporal evolution of volatile-rich melt and reactions between these melts and the mantle.
In this article, we carry out the following steps to quantify the LVL observed beneath northeast Asia. First, we conduct a new rock physics interpretation of published seismic data to yield an estimate of the melt volume fraction in the LVL in Section 2. Then, we compare seismic observations of the observed depth and inferred temperature from both the western and eastern Pacific margins with solidi of volatile-rich and dry mantle assemblages to determine the feasibility of volatile-induced melting in Section 3. Next, we use the seismically inferred melt volume fraction to model modes and timescales of reactive transport of carbonate-rich melts in an LVL in Section 4. Furthermore, based on the inferred melt volume fraction, we model the possible trace element abundances and isotopic signatures resulting from carbonate-rich melts in the LVL and compared them with the concentrations of these elements in natural lavas and xenoliths in Section 5. SUN ET AL.

10.1029/2021JB022167
3 of 17 Finally, we present the discussions and our conclusions in Section 6. Given the diversity of techniques employed in this study, we provide a short description of each step in the main article, and provide more detailed description of methods and additional results, when appropriate, in the Supporting Information S1.

Calculating the Melt Fraction
In this article, we analyze the seismic signals of P-to-S wave conversions above the northeast Asian MTZ. In the following two subsections, we present a brief description of the methods of seismic data collection and rock physics inversion to calculate the melt fraction of the LVL.

Seismic Data
The seismic data on the amplitude and depth of P-to-S (Ps) conversions are derived from Tauzin et al. (2017) for the western margin of the Pacific. We present a detailed description of the steps involved in the seismic data processing, picking, reconstruction of spatial variations, and detections, in the Supporting Information S1. This material also outlines the statistics of absolute depths and amplitudes over the region. We reconstruct the structure of the LVL from the analysis of Ps conversions recorded in the P-wave coda of teleseismic earthquakes (Langston, 1979;Vinnik, 1977). Far away from the source, P-waves arrive close to the surface approximately as a plane-wave with a steep incidence. These P-waves convert a few percent of their energy into shear waves at local discontinuities in elastic properties underneath the receivers. The amplitude of Ps conversions roughly relates to the shear wave velocity contrast across the discontinuity, whereas the time difference between the P-wave and Ps conversion gives information about the depth of conversion. Constraints on the thickness of the LVL and wave velocities within the layer come from the differential analysis of travel-times and amplitudes of Ps conversions from the top and the bottom of the LVL (the olivine-wadsleyite transition at 410 km depth). Local P-T conditions are also obtained from the analysis of the temperature-dependence of phase transitions in the MTZ.
The seismic processing workflow turns individual seismic records from earthquakes into a data cube containing amplitudes of Ps converted waves in the spatial domain. We reconstruct the structure over a grid ranging from 34°N to 49°N in latitude and 115°E-145°E in longitude, with a cell dimension of 0.5° × 0.5°. The spatial extent of this reconstructed grid is shown in the map of MTZ depth in Figure 1a. The seismic section shown in Figure 1b is a representative image of conversion amplitudes in this data cube, along the A-a profile at 37°N. Yellow colors mark positive amplitudes associated with sharp increases of seismic wave velocities such as at the 410-km discontinuity. Blue colors mark negative amplitudes associated with reductions of shear-wave velocities with depth. We use an automated method (see Supporting Information S1) to locate the negative LVL peaks within a depth range of 300-410 km (red dots in Figure 1c) and positive amplitudes indicating the top of the MTZ within a depth range of 360-450 km (blue dots in Figure 1c). Along the cross section in Figure 1c, the thickest part of the LVL beneath northeast Asia appears around the subducting slab, especially seaward on the underside of the subducting slab. This observation is spatially related with the deeper than average 410-km discontinuity at the exact same location, and is consistent with earlier interpretation of higher than normal temperatures from seismic tomography (Obayashi et al., 2006). Our seismic processing leads to a total of 1,490 data points on the grid with information on the amplitude of Ps conversions above the LVL, 410-km discontinuity, and 660-km discontinuity. This data set also contains the depths at which these peaks were observed, from which we calculate the MTZ thickness. In Section 2.2, we discuss how two independent seismically derived constraints, the normalized amplitude of Ps conversion above the LVL and the MTZ thickness, allow us to separate the thermal and chemical signatures of the anomaly.
We further filter the seismically processed data to retain only the most robust signals. To select these signals, we eliminate data points where the peaks of the LVL and the 410 km discontinuity are less than 10 km apart. We then discard points where the 410 km discontinuity was undetected, or detected with an amplitude below a threshold level of 0.5%. This typically occurs in areas with poor data coverage on the sides of the analyzed region. Next, we select data points where the absolute amplitude is higher than the standard error on amplitude for the top of the LVL. We also repeat this data filtering for the top of the MTZ. This 4 of 17 selection procedure reduces the initial grid containing 1,490 data points to a grid of 597 data points. In Table 1, we list the values of the key seismic parameters obtained from these 597 data points.

Rock Physics Analysis
In the rock physics analysis, we separate the seismic signature arising from temperature and bulk composition of the solid from the observed seismic signature. Any residual negative seismic anomaly is then considered to arise from the presence of melt. We calculate the melt volume fraction from this residual seismic anomaly over a large range of input  Note. The original, larger data set is presented by Tauzin et al. (2017). The data in this table is filtered from the larger data set as described in Section 2.1. an LVL consisting of 27% recycled slab eclogite (basalt fraction E X = 0.4), and a Clapeyron slope of 3 MPa/K. As shown in Figure 2a, the residual shear wave velocity anomaly, compensated for the effect of temperature and composition, is largely negative in this region, with percentage values ranging between −6.83 and 0.41. Despite the presence of a few highly negative signals, the median value of (%) is −2.13 with a standard deviation of 1.11. This largely negative shear velocity anomaly, not surprisingly, arises from the negative amplitudes of the Ps conversion. The normalized amplitude, norm E R is entirely negative in our filtered data set, ranging between dimensionless values of −1.83 and −0.18, with a median of −0.68 (Figure 2b).
The key outcome of our analysis is the detection of pervasive low-degree melting in the northeast Asian LVL. The histogram in Figure 2c outlines the frequency distribution of the melt volume fraction. The median value corrected from the effect of temperature is 0.76 vol% with a propagated uncertainty of 0.53 vol%. Within the LVL, the calculated melt volume fraction varies between 0 and 2.71 vol%. The influence of uncertainties in temperature, dihedral angle, Clapeyron slope, and bulk solid composition are considered in evaluating the uncertainty in the reported melt volume fraction. Here, we only discuss the robustness of our results with respect to temperature in the main text, while the trade-offs between the calculated melt volume fraction and each of these parameters are also shown in the Supporting Information S1 and Hier-Majumder et al. (2014).
Both temperature and melt reduce the effective shear wave velocity, thus rendering the separation of any related effects difficult. In other words, a warmer LVL, characterized by a warmer potential temperature, will register a lower melt volume fraction due to this trade-off. In addition, spatial variations in temperature, calculated from a fixed potential temperature and MTZ thickness, also contribute to this trade-off. For example, the LVL temperature in the data in Figure 2, corresponding to a mantle potential temperature of 1500 K (1223 °C), varied between 1189 and 1534 °C. We quantify the magnitude of this trade-off by the variation in the median melt volume fraction ( E  ) as a function of prescribed mantle potential temperature ( 0 E T ). Our error analysis reveals that this variation, . In other words, over or underestimation of the mantle potential temperature by 100 K will under or overestimate the melt by 0.3 vol%. Due to the lack of knowledge of this value, it has been traditionally difficult to quantify the melt volume fraction from seismic interpretations. Our technique of incorporating temperature in the reference shear wave velocity alleviates this problem.
The plots in Figure 3 demonstrate that our method successfully removes the effect of temperature. In Figure 3a, we plot the normalized amplitude of Ps conversion, norm E R , from the seismic observations as a function of calculated melt vol%. The data points are colored by the temperature at each grid point. For a given value of the amplitude, warmer data points register lower melt volume fractions, due to the trade-off discussed above. Figure 3b shows the residual (%) as a function of the calculated melt volume fraction. Notice that scatter in the data is absent in this plot, as the variations in temperature are already accounted for. Thus, using the temperature and composition compensated residual S E V  from Equation 1 to calculate the melt fraction allows us to isolate the effect of melting and calculate the melt volume fraction in a robust manner. . Geographical distribution of this data set is mapped in the Supporting Information S1.

of 17
In addition to the effect of temperature, we investigate which average melt fraction is supported by varying the seismic data, the dihedral angle and the bulk solid composition on the calculated melt volume fractions, as shown in Figures 4a and 4b. As noticed in the previous similar analyses (Hier-Majumder & Courtier, 2011;Hier-Majumder & Tauzin, 2017;Hier-Majumder et al., 2014;Xiao et al., 2020), for a given potential temperature and bulk solid composition, the predicted melt volume fraction increases with an increase in the solid-melt dihedral angle. Lower dihedral angle melts wet a larger fraction of intergranular contact, reducing the strength of the skeletal mineral framework (Hier-Majumder & Abbott, 2010). As a result, smaller volume fractions of melt are required to explain the reduction in matrix strength depicted by the reduction in shear wave velocity. In turn, the effect of a recycled slab component in the bulk composition, reflected by the bulk basalt fraction, is opposite. A solid consisting of a larger eclogitic slab component is stronger, hence more melt is needed to explain an observed seismic anomaly. In this work, we use a conservative estimate of 27% recycled slab component (40% basalt in the bulk) in the solid. In reality, the LVL above a stalled slab can contain a larger eclogite fraction and register a higher melt volume (Long et al., 2019;Motoki & Ballmer, 2015). 3.4 10 km E  , which is almost twice of the size of the western US LVL, 6 2 1.8 10 km E  reported by Hier-Majumder and . This result is possibly caused by differences in the age of the onset of the subduction and the size of the stagnant slab in the transition zone (Fukao et al., 2009). Long-term stable subduction provides a continuous supply of volatile components for the LVL, while the slab stalled in the transition zone contributes to devolatilization of the subducted slab and LVL generation (Hier-Majumder & Tauzin, 2017;.

Volatile Assisted Melting
Although the rock physics analysis predicts the presence of small degree of melting within the LVL, we are unable to constrain the composition of the melt from this analysis. The seismically observed depth and inferred temperature of the LVL, however, provides us with constraints on the nature of melting. Here, we combine the data for the northeast Asian LVL from this study with the previously analyzed data for the western US from Hier-Majumder and . As illustrated in Figure 5, converted temperatures and depths to the tops of the LVLs are compared to experimentally determined melting curves for carbonated basalt (Kiseeva et al., 2013;Zhang et al., 2020), hydrated peridotite (Ohtani et al., 2004), and dry peridotite (Hier-Majumder & Hirschmann, 2017). This comparison is built on two entirely independent sets of constraints, since the rock physics analysis in this work or Hier-Majumder and Tauzin (2017) did not use any constraints from these experimentally determined solidi to infer melt fractions.
Comparison between the seismic data and the experimentally determined solidi suggests that volatiles play a crucial role in melting the LVL mantle. The temperatures in both regions are considerably lower than the dry peridotite solidus and straddle the melting curves of volatile-rich solidi, especially that of carbonated basalts. This result supports the idea of small degree melting above the subducted slabs, aided by the presence of 2 CO E and possibly 2 H E O-rich fluids.
Supporting the role of carbonate-rich melting from the slab are Ti-rich 3 CaSiO E and majoritic garnet inclusions in sublithospheric diamonds. These inclusions are predominantly of an "eclogitic" or "pyroxenitic" paragenesis, generally marked by enriched trace element abundances indicative of equilibration with a low-degree melt, and are derived from pressures consistent with the low-temperature carbonated eclogite solidus, as shown in the histogram on Figure 5 ( Thomson et al., 2021;Walter et al., 2008). The overlap between the part of the high-pressure peak in the histogram and the LVL suggests that some of these garnets may have formed by reaction between slab-derived, carbonate-rich melts and the LVL mantle.
The reducing ambient mantle predicts that carbonate-rich melts in the LVLs are likely unstable in the peridotite mantle due to the solidi of carbonated peridotite (Brey et al., 2008;Dasgupta et al., 2013;Foley et al., 2009) and redox freezing (Rohrbach & Schmidt, 2011;Sun & Dasgupta, 2019). However, these studies do not consider the rate of solidification. To reconcile this apparent contradiction, our numerical model of reactive porous flow in  and this study considers the efficiency of redox freezing and the stability of carbonate-rich melts in the reduced mantle. As demonstrated by Sun, Hier- , carbonate-rich melts in the LVL can potentially survive redox freezing over long geological time scales (up to hundreds of Ma).

Melt Migration Modeling
We model the percolation rate of the buoyant, carbonate-rich melt through the mantle peridotite using a reactive porous flow model Sun, Payton, et al., 2020). Previous numerical simulations of the LVL (Leahy & Bercovici, 2007, 2010 treated the LVL as a homogeneous fluid surrounded by a high viscosity medium, with reactive infiltration driven by a phase change reaction. We now know that the partial melt in the LVL is rather small (0.8 E  0.5 vol%) and the LVL is a multiphase layer, so our model considers a process of reactive porous flow. To match the theory of "redox freezing" (Rohrbach & Schmidt, 2011), the reaction in this study is the redox reaction between the melt and the mantle matrix, rather than the phase change reaction. In this section, we primarily focus on the influence of the chemical  (Kiseeva et al., 2013;Zhang et al., 2020), wet peridotite (Ohtani et al., 2004), and dry peridotite (Hier-Majumder & Hirschmann, 2017), overlain by the slab geotherm (Turcotte & Schubert, 2002, Ch. 4). The scatter plot is from temperature and height of the top of the low velocity layers (LVLs) in northeast Asia and the western US. The histogram of the formation pressure of majoritic garnet inclusions in diamond, determined by the barometer of Thomson et al. (2021), is shown on the left panel. The solid gray curve marked "T16" is the solidus of carbonated average MORB with 2.5% 2 CO E from . The broken curve, marked "Z20," is the solidus of basalt with 2.5 wt% 2 CO E from Zhang et al. (2020). The broken black curve, marked "K13" is a solidus of an altered oceanic crust containing 10 wt% 3 CaCO E , which lost some siliceous component to subarc volcanism from Kiseeva et al. (2013). 9 of 17 reaction rate and of the slab-derived carbonate component input rate on the flux of the carbonate component in the melt and the concentration of solid C precipitation. A brief description of the model set-up and the key results are provided in the following two subsections, while details of the model including the dimensional parameters are presented in the Supporting Information S1.

Numerical Model
In our model, we focus on the LVL above a stalled slab as shown in Figure 6a. Since the carbonate-rich melt released by the subducting slab will be carried by the convective mantle, the mantle of the LVL with the melt will move as a whole horizontally, especially above the stalled slab in the MTZ. However, more important for this study is the vertical migration of the melt and the volatile within the melt, which is driven by porous flow. Under these considerations and assumptions, our model is simplified as slab-derived carbonate-rich melts percolating through a stationary mantle. During its ascent, the metastable carbonate-rich melt reacts with the reduced mantle, leading to the precipitation of solid carbon or carbonate. The freezing of carbonate-rich melts in the ambient mantle can involve two possible mechanisms: carbonation freezing and redox freezing. During carbonation freezing, carbonate-rich melts react with silicates in the mantle to produce carbonate minerals (Sun & Dasgupta, 2019), while redox freezing of carbonate-rich melts lead to the crystallization of solid diamond and metal oxides (Rohrbach & Schmidt, 2011). Carbonation freezing depends on slab-surface temperature-pressure conditions, local volatile influx, and the infiltrating melt composition. Redox freezing, in turn, is controlled by the reducing conditions in the ambient mantle. While both types of reaction are plausible, the numerical model of reactive infiltration in this article focuses on the redox freezing reaction, which is better characterized and evidenced more widely by superdeep diamonds. In future studies, the role of carbonation freezing can be modeled using the same theoretical construct presented in this section. While a number of plausible chemical reactions can be used to describe redox freezing, we consider the following reaction between the dissolved carbonate component in the melt and free Fe in the mantle (Dorfman et al., 2018;Rohrbach & Schmidt, 2011), This chemical reaction is strongly coupled to the reactive porous flow of the melt, described by a set of partial differential equations. We solve this system of five governing equations for five variables; the concentration of the carbonate component in the melt ( E c ), the concentrations of solid Fe and C in the mantle ( F E c and C E c , respectively), the percolation velocity of the melt ( E u ), and the fluid pressure ( E p ). The nondimensional form of governing equations is given by, The nondimensional numbers indicate the relative rates of the different processes involved in the reactive infiltration instability. The Péclet number, E e  , characterizes the ratio between the rates of dissolved carbonate transport by advection and diffusion. The Damköhler number ( E a  ) represents the ratio between the rates of chemical reaction and percolation. In this study, we select a narrower range, similar to that suggested by Spiegelman et al. (2001), and report results for E e  ranging between 10 and 100 and E a  ranging between 1 and 100, while the study of Sun, Payton, et al. (2020) (Sun, Payton, et al., 2020, also see Supporting Information S1). We carried out simulations for E  varying over a range between 0.1 and 6 in this study.
The model set-up and boundary conditions are displayed in the schematic diagram in Figure 6a. We solve the governing equations by adopting a set of Dirichlet, Neumann, and periodic boundary conditions on the primary unknowns, E u and E c . The complete set of boundary conditions are given by, (11)  . Finally, we set a zero diffusive flux Neumann condition on E c at the top and bottom boundaries, Since this is a time-dependent problem, the necessary initial conditions are prescribed as zero carbonate in the melt, 1 wt% free Fe and zero solid carbon in the mantle prior to the onset of melt percolation, namely, where 0 [Fe] E = 0.01 is the metallic iron mass fraction in the mantle matrix.
We discretize the governing PDEs using the finite element method. Then, based on the open source code MuPoPP, we generate the mesh, solve the discretized algebraic equations using Krylov sparse solvers, and output the data for further visualization and analysis.

Melt Channels, Residence Time, and Fluxes
The series of visualizations in Figure 6 display three snapshots in time, spanning over 80 Ma, demonstrating that carbonate-rich melts can survive redox freezing. Under these conditions ( 100, 10, 2 E e a       ), melt migration occurs through localized channels. Channelization also leads to the formation of reaction zones containing up to 200 parts per million by weight (ppmw) solid carbon (Figure 6, right), marking the sites of diamond formation. Mantle-derived melts arising from the LVLs can thus inherit geochemical components from the reaction zone and the nonreactive parts of the mantle, indicating heterogeneous geochemical signatures of their sources, as sometimes observed in basalts from northeast Asia (Qian et al., 2020;S. Li et al., 2017).
The flux of the carbonate component in the melt and the concentration of solid C, mapped over the parameter space in Figure 7, provide important information about carbon storage in the LVL. First, the flux of the carbonate component in the melt is strongly dependent on the carbonate input from subduction, as shown by the plot in Figure 7a. While this result is physically intuitive, it contradicts the conclusion that redox freezing completely consumes carbonate-rich melts (Rohrbach & Schmidt, 2011). To fully understand the stability of these melts, the rate of input from subduction must also be taken into account. We calculated the dimensional value of the flux using the calculated value of the northeast Asian LVL surface area, approximately . Second, this plot also demonstrates that the control of subduction input on the flux is substantially modified by the E e  number. While the dependence is nearly linear for 100 E e  , the flux remains much lower for E e  = 10 when 4 E   . Third, the influence of melt reactivity, characterized by the E a  number, is more pronounced on the concentration of average solid carbon concentration than the flux of melt carbonate, as revealed by a comparison between Figures 7a and 7b. With a tenfold increase in E a  , the C concentration in ppmw in the mantle also increases by nearly an order of magnitude. Over the range of parameters shown in the plot, the LVL sequesters between 0.1 and 80 ppmw solid C on average, while the LVL can locally rate to the LVL. For reference, we also show the total (molten and solid) carbon concentrations in the OIB and MORB sources (Hirschmann & Dasgupta, 2009). Error bars in this plot are smaller than the symbol sizes.
12 of 17 sequester up to 200 ppmw solid C (Figure 6b). Finally, the solid C concentration in the LVL increases with an increase in the subducted carbonate input. This result outlines the fact that the accumulation of solid C in the mantle might have varied over time, as the rate of subduction and the amount of slab-derived carbon varied.
An important outcome of the numerical models is the relatively sluggish velocity of melt percolation. The simulations in Figure 6 show that the melt channels reach the top of a 40 km thick LVL in 80 Ma, at a velocity of 0.5 km/Ma (500 E  m/yr). This velocity is even lower, E  Table 1). These estimates of sluggish melt flow bear two important consequences. First, the relatively immobile slab-derived melt can be entrained into faster moving upwelling convective flow and transported to shallower levels. Second, the several hundred million years of the melt residence time will allow extensive mantle metasomatism and in-growth of isotopic animalies.

Geochemistry
Inclusions in diamonds and petrological modeling provide evidence for the kinds of melts that may exist above a subducted slab in the transition zone, leading to metasomatism of mantle peridotite and generation of geochemically distinct sources. Inclusions of Ca(Ti,Si) 3 O E in diamond, interpreted as former Ca-perovskite, and many majoritic garnets have major and trace element compositions indicating derivation from subducted oceanic crust (Bulanova et al., 2010;Thomson et al., 2014;Walter et al., 2008). In particular, extreme enrichments in many incompatible elements (e.g., Th, U, Nb, Ta, REE) are inconsistent with a subsolidus origin from basaltic crust, but are best modeled by equilibration with a low-degree melt, likely of carbonatitic affinity, derived from subducted eclogitic crust in the transition zone (Bulanova et al., 2010;Walter et al., 2008). Further support for an origin in or close to the transition zone is a major peak in the crystallization pressures at E  14 GPa in the global database of majoritic garnet inclusions of eclogitic affinity (Thomson et al., 2021), in good agreement with estimates of the pressure at which slab top geotherms are expected to cross the solidus of carbonated eclogite shown in Figure 5 (Kiseeva et al., 2013;Walter et al., 2008;Zhang et al., 2020).
To illustrate the result of such interactions between a carbonate-rich melt and the mantle, Figure 8a shows the composition of a Ca(Ti,Si) 3 O E inclusion from Juina, Brazil (J1) (dashed green line) that has been identified as crystallizing from a primitive melt of oceanic crust (Bulanova et al., 2010;Walter et al., 2008). Also shown is the estimated trace element composition of the coexisting melt (solid green line) on the basis of experimental partitioning coefficients, where melt E C = C D min min melt / / and using partition coefficients from . The red dashed line represents a model of subducted, processed oceanic crust where a total of 10% hydrous component has been extracted at 2.5, 4.5, and 6 GPa, mimicking crust dehydration and using experimentally derived mineral/fluid partition coefficients (Kessel et al., 2005;Klimm et al., 2008;. Dehydration mobilizes elements like Sr and Pb relative to other trace elements. The solid red line represents a low-degree melt (0.1%) of processed MORB at transition zone conditions calculated using the methods, partition coefficients, and melting reactions described by . Also shown is a shaded field of calculated melts that could coexist with clinopyroxene xenocrysts from the Cenozoic Wudi basalt from eastern China (Qian et al., 2020). We note the striking similarity in the overall abundance patterns, and especially the Sr, Pb, Hf, and Zr depletions, as well as the predicted high U/Pb and Th/Pb ratios of all of these independently estimated melt compositions that potentially represent slab-derived melts present in the LVL above a subducted slab at transition zone depths. The black line shows the trace element abundances expected in primitive peridotitic mantle after the addition of 1% of the model low-degree melt derived from the subducted processed oceanic crust. are not meant to be dispositive, we note that HIMU-type lavas can be generated from the model melt-metasomatized peridotitic source in E  500 Ma. On the basis of isotopic and geochemical constraints, Qian et al. (2020) suggested that melts from dehydrated, carbonated oceanic crust have the chemical characteristics required for the HIMU component recorded in the Wudi xenocrysts, and indicated a mantle source metasomatized by carbonatitic melt from subducted crust no older than 700 Ma, consistent with our modeling. Likewise, high 206 Pb/ 204 Pb and low 207 Pb/ 204 Pb lavas from Bermuda have also been suggested to be caused by deep mantle metasomatism from a volatile-rich melt no earlier than E  500 Ma ago (Mazza et al., 2019), again consistent with our modeling of metasomatized mantle by carbonated melts derived from subducted oceanic crust. We note that calculated values at 500 Ma for 87 Sr/ 86 Sr (∼ 0.7030), 143 Nd/ 144 Nd (∼ 0.5130) and 176 Hf/ 177 Hf ( E  0.2829) using our model parameters are also generally consistent with a HIMU-like mantle source (Stracke et al., 2005). Overall, our modeling indicates that deep metasomatism of the mantle through the interaction with melts derived from oceanic crust is a viable mechanism to produce a HIMU-like component in eastern China Wudi lavas and likely anywhere that deep subduction results in melting of carbonated oceanic crust (Castillo, 2015;Weiss et al., 2016;Y. Xu et al., 2012;H. Li et al., 2016).

Discussions and Conclusions
Our rock physics analysis of the Ps conversion data indicates that the northeast Asian LVL is marked by an average of 0.8 vol% melt, with occasional patches exceeding 2 E  vol% melt (Figure 2). Comparison of the inferred temperature and depth of the seismic data from this study and a similar study of the LVL beneath the western US reveals, independently, that the melting is assisted by  . The thick black line shows a model two-stage growth curve using model parameters from Stacey and Kramers (1975). The thin black lines are growth curves for mantle metasomatized by addition of 1% of the low-degree melt produced from processed MORB, where age labels represent the time of metasomatism.  from the subducting slab to the upper mantle ( Figure 7a). As the two most representative LVLs associated with the subduction zones, the LVLs beneath northeast Asia and the western US are estimated to have a total area of 6 2 5.2 10 km E  and can transfer up to 78 Mt/yr 2 CO E in total. This value of carbon flux through the top of the LVLs is of similar magnitude to the global subducted carbon fluxes between 80-176 Mt/yr (Dasgupta & Hirschmann, 2010) and 300 E  50 Mt/yr (Plank & Manning, 2019) 2 CO E , suggesting that the LVLs related to subduction can be significant reservoirs of subducted carbon in the deep upper mantle, and an important component of the global carbon budget. At the same time, as the simulations predict in Figures 6b and 7b, the LVL can immobilize up to 80 ppmw solid C on average to the deep upper mantle by redox freezing, and even up to 200 ppmw solid C distributed in localized reaction zones. Compared to the carbon content of either the OIB source with 120-1,830 ppm 2 CO E or the MORB source with 25-95 ppm 2 CO E (Hirschmann & Dasgupta, 2009), this process can effectively increase the regional carbon content fixed in the mantle, demonstrating the potential of the LVL as a deep carbon reservoir. These results bear implications for the big mantle wedge (BMW) beneath northeast Asia. For example, the geochemical signature in the source of the Wudi basalts could have resulted from reaction between slab-derived melt and the ambient mantle during an early Cretaceous subduction event, which also led to the formation of the big mantle wedge beneath northeast Asia (Ma & Xu, 2021). The melt residence time in the LVL in our simulations is consistent with the timing of these two events in the northeast Asian upper mantle. Finally, using the estimated melt fractions, our geochemical model calculates the abundances of trace elements resulting from the melt-mantle interaction (Figure 8), demonstrating that such interaction can produce a HIMU like source component. When the HIMU like metasomatized mantle is carried by upwelling convective flow of the upper mantle to the base of the lithosphere, they will sustain the intraplate volcanism and form basalts, which inherit components or geochemical characters of the LVL, like the HIMU component in eastern China Wudi lavas.
While the results of this study, combining several different techniques, bears important implications for the global carbon cycle, a few factors need to be taken into account. First, future seismic investigations, employing a denser array of sensors will greatly improve the resolution of the seismic data and provide a more detailed picture of the internal structure of the LVL. Second, the calculated melt volume fraction depends on a number of assumed parameters. While we carried out a rigorous analysis spanning a large parameter space and reported the uncertainties associated with each of the parameters, future inversions using similar rigorous constraints on different seismic datasets will provide additional constraints on the melt fraction in the LVL. Third, our numerical simulation of melt migration is carried out in a two-dimensional domain, greatly simplifying the complex geometry of the LVL. In addition, these simulations do not account for possible temporal variations in the carbonate flux associated with subduction. Future 3D melt migration simulations, based on the tectonic history of the western Pacific subduction zone, will greatly enhance our understanding of the fluxes and residence times of carbonate-rich melts in the mantle beneath northeast Asia. Finally, such simulations will also aid in understanding the evolution of the geochemical signature of the mantle over time periods longer than those considered in this work.