The lateral variation of P- and S-wave velocities in the mantle is dominantly caused by temperature, water content, and chemical composition anomalies. Our mineral physics modeling builds a bridge between observed velocity anomalies and these physics and chemical factors. In addition, only considering anharmonic effect in the inversion of thermal structure would introduce a large degree of partial melt and a much thicker oceanic crust than observation (Karato, 2008); hence our mineral physics modeling will consider both the anharmonic effect and anelastic effect on velocity.

Using the improved mineral physics method, we predict the P and S velocity perturbations at different mantle conditions relative to a reference model (Figure 1). We choose a one-dimensional 1,400°C adiabatic geotherm as the reference temperature model (Supplementary Figure S1B); the pyrolite model (Carmichael et al., 1974; McDonough and Sun, 1995) as the reference mantle composition model, which consists of 82.5 mol% harzburgite and 17.5 mol% basalt (Supplementary Table S2). In addition, according to the research of hydrous wadsleyite and ringwoodite (Wang et al., 2019; Wang et al., 2021), we assume the reference water model is wadsleyite with 0.5 wt% water for the upper MTZ, ringwoodite with 0.25 wt% water for the lower MTZ, and other secondary minerals are dry. Furthermore, the reference Q model for anelastic effect is derived from PREM (Dziewonski and Anderson, 1981). The velocity anomalies caused by anharmonic effect (about 0.5%–3%) are all larger than those caused by anelastic effect (lower than 0.5%) since the energy dissipation in the MTZ is much smaller than that in the upper mantle (Figure 1). However, the anelastic effect is not negligible in the MTZ. For instance, the S wave anelastic anomaly caused by 200 K temperature anomaly is up to 0.5% in the MTZ (Figures 1E,F). Compared with the velocity anomalies caused by ∼1 wt% water content and 200 K temperature anomaly, the velocity anomaly caused by 7.5 mol% harzburgite anomaly is minor and lower than 1% in the MTZ (Figures 1C,D), which indicates the composition effect in the MTZ is minor. There are several possible composition anomalies in the MTZ, such as basalt content anomaly, Fe content anomaly, and Al content anomaly. Since the following inversion procedure does not have enough information to distinguish the differences between those anomalies, we use basalt content anomaly to represent mantle composition anomaly.

Due to the complicated phase transformations caused by temperature variations, the profile has some slight fluctuations. For instance, when temperature decreases 200 K, the phase boundary of β to γ olivine will elevate (Supplemenatry Figure S1A), causing the velocity anomalies around 450 km to increase largely (Figure 1E). Moreover, since the volume fraction of ilmenite increases with decreasing temperature around 550 km (Supplemenatry Figure S1A), the velocity anomalies fluctuate in the lower MTZ (Figure 1E).

Because of the significant water storage capacity of wadsleyite and the large water sensitivity of ringwoodite, the reduction of anharmonic velocity can reach 1%–3% in the water-saturated situation, which is comparable with those from 200 K temperature anomaly (Figure 1B). Unlike olivine (Mao et al., 2010), the anharmonic water effect of wadsleyite and ringwoodite is not largely suppressed by increasing pressure, so there are still remarkable velocity reductions in the deep MTZ. The velocity reduction caused by water in the lower MTZ is more significant than in the upper MTZ since ringwoodite is more sensitive to water than wadsleyite. In addition, P wave anomalies caused by water effects are slightly larger than S wave anomalies (Figures 1A,B), which is different from those caused by temperature and composition effects.

The 3-D sensitivity kernels of P- and S-wave velocity anomalies to water, temperature, and harzburgite abundance for the upper MTZ (462 km depth) and lower MTZ (587 km depth) are calculated (Figures 2, 3). At the upper MTZ, both P- and S-wave anomalies show small sensitivities to water, and they only change 1% in the range of 0–2 wt% water content when dT >0. However, they become more and more sensitive to water content when temperature decreases (Figures 2A,B) since the volume fraction of ringwoodite, which is more sensitive to water than wadsleyite, increases with decreasing temperature at 462 km (Supplementary Figure S1A). At the lower MTZ, water content significantly affects velocity anomalies because of the high sensitivity of ringwoodite to water content. For instance, shear wave velocity anomalies change from 6% to −2% with increasing water content when dT = −400 K (Figure 3B). When dT >400 K, velocity anomalies do not change with water content because when the temperature is greater than 2000 K, wadsleyite will directly decompose into garnet and magnesiowustite instead of transforming into ringwoodite (Figure 10.17 in Gasparik, 2003), which leads to the insensitivity of velocity anomalies to water content.

The temperature effect is comparable with the water effect in the MTZ. P and S wave velocity anomalies change about 4% and 8%, respectively, in the whole range of temperature variation (Figures 2C,D, 3C,D). Although temperature effect decreases with depths due to the pressure effects on thermal expansion and attenuation (Karato, 2008), it is prominent in the MTZ.

Compared with temperature and water effects, the chemical composition effect is minor. The variations of velocity anomalies are lower than 1% when harzburgite abundance changes from 70 mol% to 90 mol% (Figures 2E,F, 3E,F). Since the chemical effect is not suppressed by pressure and becomes more visible in the deeper mantle (Karato, 2008), we find that the velocity variations in 587 km caused by composition anomaly are slightly larger than those in 462 km (Figures 2E,F, 3E,F).

Based on the above mineral physics modeling, we can invert water content, temperature, and composition structures in the MTZ from a seismic tomography model. The grid search method is implemented to find the optimal temperature, composition, and water anomalies for observed P- and S-wave velocity anomalies. We define a misfit function χ ( m ) and regularization function R ( m ) as: χ ( m ) = | ( δ ⁡ log ⁡ V P ( m ) )   s y n − ( δ ⁡ log ⁡ V P )   o b s | + | ( δ ⁡ log ⁡ V S ( m ) )   s y n − ( δ ⁡ log ⁡ V S )   o b s | ( δ ⁡ log ⁡ V P ( S ) ( m ) )   s y n = ( δ ⁡ log ⁡ V P ( S ) ( m ) )   a h + ( δ ⁡ log ⁡ V P ( S ) ( m ) )   a n R ( m ) = | ( m − m r e f ) / m r e f | m = ( T , W , H ) where the subscripts s y n and o b s represent the synthetic velocity anomalies from our mineral physics modeling and the observed velocity anomalies from seismic model respectively. The synthetic velocity anomalies equal anharmonic velocity anomalies (subscript ah) plus anelastic velocity anomalies (subscript an). m is the unknown model, and the degree of freedom of m is 3 representing temperature T, water content W, and harzburgite mole fraction H, respectively. Since the composition effect is minor, we only introduce harzburgite anomalies when the observed velocity anomalies cannot be explained well by temperature and water anomalies. Moreover, if there are multiple optimal solutions, we choose the solution closest to our reference model m r e f as our final model (i.e., the solution makes the regularization term minimum). m r e f   represents the reference model described in 2.2. Unknown temperature anomalies are searched from −800 to 800 K with a search step length of 10 K. Unknown water contents are searched from 0 to 2 wt% H₂O with a search step length of 0.01 wt% H₂O in subduction zones, but for regions far away from subduction zones, we only search for solutions in the range of 0–1 wt% because the maximum water storage capacity of olivine polymorphs is 0.5–1 wt% along the normal geotherm (Ohtani et al., 2004). Harzburgite abundances are searched from 69.5 mol% to 92.5 mol% with a search step length of 1 mol%. We only keep the solutions with an error less than 0.2%. For the regions with a misfit greater than 0.2%, we cannot explain well, and other factors (e.g., the effect of dense hydrous magnesium silicates) should be involved in the future.

In this study, we use TX2019slab (Lu et al., 2019), a P and S wave joint tomography model based on a linearized inversion of the major P- and S- body wave phases, to constrain the water distribution of the MTZ. TX2019slab shares some common features with other P and S wave joint tomography models (Houser et al., 2008; Simmons et al., 2010; Tesoniero et al., 2015) and has a higher resolution in the MTZ. Except for the North America (NAM) subduction zone, high velocity anomalies (1.5%–4%) related to subducted slabs in 7 subduction zones (i.e., the Calabria (CAL), western Java (JAV), northern Honshu (HON), northern Mariana (MAR), Tonga (TON), Central America (CAM) and northern Chile (CHI) subduction zones.) are well resolved. Low velocity anomalies (−1% to −2%) are dominantly located in Southern Eurasia, South Pacific, and East Africa (Figure 4). Compared with other models, TX2019slab reveals a more refined slab geometry and higher amplitude of velocity anomalies within slabs since they incorporated a 3-D subducting slab model into their starting reference model in inversion. Lu and Grand (2016) pointed out that the incorrect imaging of slabs can bias tomography results (Lu et al., 2019), so a tomography model with more accurate slab structures can provide a more reliable water distribution in the MTZ. The above analysis shows TX2019slab is an appropriate model for revealing the water anomalies near subduction zones.

van Keken et al. (2011) used geodynamic calculation to model the thermal and petrological structures of the subduction zones and to predict water contents of downgoing slabs. Their calculation indicates that cold slabs, such as TON, MAR, HON, and JAV, dehydrate slowly and can transfer water into depths lower than 230 km. In contrast, warm slabs, such as NAM, become nearly anhydrous at shallow depths (Figures 6A,B), which is consistent with our results. Their study suggests that the CHI slab is warm and dehydrates significantly above 230 km depth, while our results show a clear hydrous signal in the MTZ close to the CHI slab. The possible reason for the discrepancy is that the subduction of the middle ocean ridge (Hu et al., 2016) at the CHI slab brings more water into the mantle than other slabs, so the slab remains hydrous in the MTZ, although it significantly dehydrates in the upper mantle.

Besides seismic velocity, electrical conductivity is also highly sensitive to water, and a small amount of water can significantly increase electrical conductivity (Karato, 1990). Sun et al. (2015) reported a high-resolution global electrical conductivity model, which shares some common features with our results (Figures 6C,D). Their model shows high electrical conductivity anomalies are mainly located in the HON, MAR, JAV, and CHI subduction zones (red areas in Figures 6C,D). Moderate electrical conductivity anomalies appear near the CAM subduction zone (white areas in Figures 6C,D). Low electrical conductivity anomalies appear close to the NAM subduction zone (blue areas in Figures 6C,D). Those characteristics are all consistent with our results (Figure 5). In addition, high conductivity anomalies are also located in central Eurasia and the northern Indian Ocean (red areas in Figures 6C,D), which show high water contents up to 0.6–1 wt% in the upper MTZ in our results (Figure 5). However, our results exhibit a hydrous MTZ beneath the TON subduction, which does not appear in the conductivity model. Since the electrical conductivity of minerals is not only affected by water but also by temperature and chemical composition (Huang et al., 2005; Yoshino et al., 2008; Yoshino and Katsura, 2009), we need to consider the combined effect of these factors to understand the high electrical conductivity in the MTZ.

According to our water distribution results and the water contents and water distributions of the upper mantle (Peslier et al., 2017; Ohtani, 2021), together with the information on deep mantle melting (Karato et al., 2020) and plate tectonic history (Brown et al., 2020), we propose a water circulation model in Figure 7.

Hydrous regions in the upper mantle are near subduction zones, where slabs transported 0.03–0.3 wt% water into the MTZ (Suetsugu et al., 2010; Peslier et al., 2017; Ohtani, 2021). Because plate tectonics began in the early Paleoproterozoic (2.2 Ga) (Brown et al., 2020), the subduction of oceanic plates may have been going on for hundreds of millions of years. Since the diffusion rate of hydrogen is slow (∼1 km for 100 My) but the motion of water-rich melts is fast (∼100–1,000 km for 100 My), the motion of water-rich melts reflects the redistribution of water in the Earth (Karato et al., 2020). The thermal structure of our inversion results demonstrates that it is difficult to produce intensive partial melting in the MTZ. However, when subducted slabs penetrate the 660-km discontinuity, dehydration and partial melting will occur due to the significant contrast of water storage capacity between the MTZ and the lower mantle (Litasov and Ohtani, 2002; Hirschmann et al., 2009; Inoue et al., 2010b; Litasov et al., 2013), and most released water enters the melt, resulting in water-rich melts below the 660-km discontinuity. Since the water-rich melt is lighter than the surrounding mantle below 660 km, it returns to the MTZ with upwelling flows (e.g., mantle plumes) (Sanloup et al., 2013; Sakamaki, 2017; Karato et al., 2020). Similarly, dehydration melting occurs when upwelling flows pierce the 410-km discontinuity. Because the hydrous melt above 410 km is heavier than the surrounding mantle, the water-rich melt will return to the MTZ with downwelling flows (e.g., subducted slabs) (Sanloup et al., 2013; Sakamaki, 2017; Karato et al., 2020). Similar to the plate tectonic, such redistribution of water may have been going on for hundreds of millions of years, so the water released from slabs can be accumulated in the MTZ, the upper MTZ, and the lower MTZ have accumulated 0.3–0.5 wt% and 0.15–0.2 wt% water, respectively (Figure 7), during the whole tectonic history.