Heterogeneous distribution of water in the mantle transition zone inferred from wavefield imaging

Heterogeneous distribution of water in the mantle transition zone inferred from wavefield imaging

Earth and Planetary Science Letters 505 (2019) 42–50 Contents lists available at ScienceDirect Earth and Planetary Science Letters www.elsevier.com/...

4MB Sizes 1 Downloads 38 Views

Earth and Planetary Science Letters 505 (2019) 42–50

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

Heterogeneous distribution of water in the mantle transition zone inferred from wavefield imaging Yinzhi Wang a,∗ , Gary L. Pavlis a , Mingming Li b a b

Department of Earth & Atmospheric Sciences, Indiana University, 1001 E 10th Street, Bloomington, IN 47405, USA School of Earth and Space Exploration, Arizona State University, PO Box 876004, Tempe, AZ 85287, USA

a r t i c l e

i n f o

Article history: Received 18 January 2018 Received in revised form 12 September 2018 Accepted 8 October 2018 Available online 24 October 2018 Editor: J.P. Avouac Keywords: mantle transition zone receiver functions imaging mantle convection

a b s t r a c t The amount of water in the Earth’s deep mantle is critical for the evolution of the Earth. Mineral physics studies revealed that the mantle transition zone can store several times the volume of water in an ocean. However, the actual water distribution in the transition zone remains enigmatic. We used the highest resolution images produced of scatterers in the North-American transition zone derived from teleseismic data recorded by the Earthscope Transportable Array. We find that the transition zone is filled with previously unrecognized small-scale heterogeneities that produce pervasive, negative polarity signals. Simulations demonstrated the images can be explained by low-velocity bodies shaped as distributed blobs or near vertical structures. We suggest these low-velocity bodies may be heterogeneously distributed, water-enriched subducted harzburgites produced through long-term accumulation. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Water in the Earth’s mantle significantly reduce mantle viscosity and the melting temperature of mantle rocks. Water could have a very important role in the Earth’s evolution by affecting mantle dynamics, partial melting processes, and the degassing history of the Earth. Mineral physics experiments have shown that the mantle transition zone is capable of storing a large amount of water in the crystal structure of wadsleyite and ringwoodite, which are the high-pressure polymorphs of olivine in the transition zone (Bell and Rossman, 1992; Bercovici and Karato, 2003; Hirschmann, 2006; Pearson et al., 2014; Ringwood and Major, 1967). Although some geophysical observations (Emry et al., 2015; Fei et al., 2017; Huang et al., 2005; Smyth and Jacobsen, 2006; van der Lee et al., 2008) and direct evidence from inclusion of diamond samples (Pearson et al., 2014) suggest that the transition zone could be hydrous, the distribution of water in the transition zone is insufficiently mapped (Houser, 2016; Jacobsen and van der Lee, 2006; Karato, 2011; Suetsugu et al., 2006). Seismic observations have provided one way to constrain the distribution of water in the mantle transition zone (Suetsugu et al., 2006; van der Meijde et al., 2003). Shen et al. (1998) and Gao et al. (2002) appear to be the first to document the existence of neg-

*

Corresponding author. E-mail address: [email protected] (Y. Wang).

https://doi.org/10.1016/j.epsl.2018.10.010 0012-821X/© 2018 Elsevier B.V. All rights reserved.

ative polarity P to S scattering phases within the transition zone. After several regional receiver function studies confirmed the presence of similar features at different depths (Bagley et al., 2013; Bonatto et al., 2015; Gao et al., 2010; Jasbinsek et al., 2010; Schmandt, 2012; Shen et al., 2008; Shen and Blum, 2003), Shen et al. (2014) suggest that a negative phase within the transition zone is a global feature. Tauzin et al. (2013) used the common conversion point (CCP) stacking method to image the fine structure of the transition zone. They identified negative phases in the transition zone over the entire western US. Tauzin et al. (2017) observed similar features that suggested a direct link to active subduction in the western Pacific. However, CCP stacking has the tendency to produce artificial horizontal structures (Pavlis and Wang, 2015; Wang and Pavlis, 2016a). It is not surprising that these studies interpreted their results as several minor discontinuities at varying depth ranges and of variable size. Our results suggest that interpretation is likely due to an artifact of horizontal appearance in CCP stacking results. Using the data from complete coverage of the USArray, this paper provides a new understanding of the negative phases observed in these previous studies within the transition zone. Our new results are possible because of the increase in resolution made possible by the plane wave migration (PWMIG) method (Pavlis, 2011a; Poppeliers and Pavlis, 2003a, 2003b), which was shown in previous papers to provide higher resolution of small-scale heterogeneities in the transition zone than the CCP stacking method used in previous work (Liu and Pavlis, 2013; Pavlis and Wang, 2015; Wang and Pavlis, 2016a). It is not sub-

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

ject to the horizontal layering assumption in CCP stacking, and the result are limited mainly by diffraction. PWMIG produces an image volume with superior resolution of scattering from irregular smallscale anomalies. We also validate the results using point source simulation as a tool for error appraisal, which is lacking in previous work on imaging the scattering features in the transition zone. The error appraisal demonstrates the limitations of the data and strengthens the interpretation.

43

Liu and Pavlis (2013). That simulator uses a point-source Born approximation to forward model the scattering of small objects (see Supplementary Text 1.1). By superimposing a large number of single scattering, point response functions, this approach can also be used to model any continuous features of a range of geometries in the subsurface when the spacing of the points is smaller than the dimensions of the Fresnel zone. The supplement discusses a theoretical ambiguity in that model and how we have validated the results against other simulation methods.

2. Data and method 3. Results The foundation of this paper is a 3D image we constructed from 141,080 receiver function estimates from the EarthScope Automated Receiver Survey (EARS) (Crotwell and Owens, 2005) recorded between May 1990 and May 2015. These receiver functions were computed by the IRIS facility using the iterative deconvolution method (Ligorría and Ammon, 1999). Following the same workflow of Wang and Pavlis (2016a), we reprocessed this dataset with the generalized iterative deconvolution method (Pavlis and Wang, 2015; Wang and Pavlis, 2016b) to reshape the results to a 0.1 Hz Ricker wavelet instead of the Gaussian wavelet used by the EARS. Pavlis and Wang (2015) show the Generalized Radon Transform imaging method used by PWMIG has poor resolution with a Gaussian wavelet due to the implicit spatial smoothing of the operator. Dramatic improvements were achieved by reshaping the data to a Ricker wavelet before imaging. PWMIG is a fully 3D vector imaging method. This imaging technique is a variant of the technology of “migration” that is the foundation of all active source seismic imaging in the oil and gas industry today. Our approach is currently the only fully threedimensional migration technique successfully applied to Earthscope TA data. While the details for the imaging workflow are discussed in Wang and Pavlis (2016a), it is worth noting that previous simulations showed that PWMIG has superior resolution (∼140 km in the transition zone) for non-horizontal structures than the simpler and more common CCP stacking method. By considering signals in P to S receiver function as the result of forward scattering process, PWMIG can migrate all scattered S wave energy back to the scatterer to image structures more complex than horizontal layers assumed by CCP stacking. See Supplementary Text 1.1 for a discussion on the theoretical foundation of the scattering process. The images were migrated to true spherical geometry using a selection of 3D tomography models we obtained from IRIS EMC (IRIS DMC, 2011) including DNA13, USSL14, MIT15, US16 and GyPSuM models. We found USSL14 (Schmandt and Lin, 2014) was superior for imaging the transition zone because it did not underestimate depth of the 410 and 660 discontinuities (d410 and d660) east to the Great Plains. Previous work (Wang and Pavlis, 2016a) argued that these offsets were related to inadequate estimates of high velocities in the upper mantle under the craton. On the other hand, the results of this paper focusing on the amplitude features within the transition zone do not depend on the choice of models, so we used the updated model determined by Burdick et al. (2017) for consistency with previous work. Simulation validation is a significant component to address strength and limitation of any scattered wave imaging methods. In this study, two simulation methods are used to provide a better perspective to interpret the imaging result. The first is a radially symmetric layered model simulator (Pavlis, 2011a) that simulates only the primary conversions from flat layer boundaries. Although it simplifies the scattering to be of constant amplitude and angle independent, we have found it a useful tool to evaluate how random noise influences the image in the transition zone when used in combination with random noise added to simulated data. In addition, we use the point-source simulation method developed by

Our result reveals the transition zone to be a scattering region with strong preference towards negative scattering objects. Fig. 1A shows example cross-sections that illustrate how the amplitude we recover for the transition zone is negative everywhere we can see. For the full perspective see Movies S1 and S2. The only exceptions are a few points at the edges of the imaging domain where the results are not reliable. We imported these data into a modern 3D seismic reflection interpretation package. We experimented with both automated methods and careful manual picking. With this careful look at the data we recognized that we could find no structure with a consistent dimension larger than ∼140 km horizontally. Earlier work with the same data demonstrated structure at that scale is near the diffraction limit (Wang and Pavlis, 2016a). By ‘consistent’ we mean that despite their appearance, they do not define coherent, continuous horizons spanning significant distances in all directions. Fig. 1 illustrates an example with manual picking. For this example, we attempted to create surfaces starting from the largest negative amplitudes in line 38. We manually tracked each horizon north and south till we reached a near zero point. Each surface truncates a short distance north and south of the initial points and cannot be tracked reliably in 3D (Fig. 1C, Movie S3). We find the same thing happens for all the strong negative horizons in the transition zone. Fig. 2A shows a simple metric describing this bulk property for the whole transition zone. We plot the amplitude averaged vertically between d410 and d660. A question of concern is whether the seismic anomalies in Fig. 2A are artifacts of noise modulated by migration artifacts linked to the d410 and d660. To test this hypothesis, we constructed a model with flat d410 and d660 surfaces similar to previous work (Pavlis, 2011b). The simulation matches exactly the coverage of the original dataset. To each synthetic seismogram, we added Gaussian white noise with a rms amplitude −6 dB below the largest amplitude. We then processed the simulation data with the same workflow as the real data. Fig. 2B shows a map of average amplitude computed from this simulation. Unlike the real data, substantial portions of the region show a positive average amplitude. The histogram shows that amplitude is distributed across the positive and negative ranges with a negative shift. We attribute this shift to the negative side lobe of Ricker wavelet signal from d410 and d660 surfaces smeared across the transition zone by the migration operator as predicted by our hypothesis. However, this simulation shows that random noise and smearing by the migration operator cannot explain the observed strong preference of negative returns from the transition zone. The preference for negative amplitudes in the data is much stronger than that of the simulation. Another unique feature of PWMIG is that it is comparable to “prestack migration” methods used in reflections seismic imaging. That means we separately image stacks of data from different source regions. We use the coherence metric defined in an earlier paper (Pavlis, 2011a), which measures the similarity in volume sections of dimension 140 × 140 × 20 km down-sampled from the original grid at 35 × 35 × 2 km. A key result from the coherence analysis is that many of the strong negative bodies within

44

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

Fig. 1. Interpreted negative scattering objects. The original image grid is defined in the GCLgrid (Fan et al., 2006) that maps the grid with Mercator projection with origin at 41◦ latitude and −109.05◦ longitude. The grid showing here is distorted to cuboid from the original grid when imported into the seismic interpretation package. They are screen shots from the seismic interpretation package we used which treated these data like 3D seismic reflection images. Line 40 shown here is the great circle path going through the origin of the original grid. (A) Cross-sections through part of the image volume. The negative scattering objects marked with lines of different colors in the cross-sections is manually interpreted following all of the strong and significant negative amplitude from Line 38. The y-axis is depth in km. These sections are conformal maps from the true spherical section to a rectangle. (B) Map showing the locations of these cross-sections as red lines. The spacing between each cross-section is 70 km and all are 4515 km long. (C) A snapshot from Movie S3 showing the interpreted segments as surfaces in 3D. The east-to-west cross-section is Line 38. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

the transition zone have coherence comparable to that of d410 and d660 but are more localized (Movie S4). Localized high coherence is expected even with pure noise, but the simulation results in Fig. 2C suggest the coherence of the data is much higher than expected from random noise. Fig. 2C shows a median stack of coherence profiles from all grid points at constant depth and characterize the uncertainty by interquartile. The much higher value of median coherence of the data compared to the simulations argues that this is not a pure artifact from noise. The upper quartiles show an even stronger difference. That fact is a quantitative validation of our observation that many areas show coherence comparable to d410 and d660. Fig. 2C also shows that the signal from the transition zone has a higher average coherence than signals from immediately above d410 and below d660. Finally, we note that the coherence stack in the lower part of the transition zone (∼520 km to 660 km) has higher median coherence than the upper transition zone (410 km to ∼520 km). Taking 520 km as a dividing line we produced Fig. 3 by dividing the transition zone into two different vertically averaged amplitude maps. Negative scattering still dominates both the upper and lower portions. The lower portion, however, has greater negative amplitudes in most of the region. Notice that the negative scattering in the western US is weaker in the upper transition

zone and stronger in the lower transition zone. In the central US, both the upper and lower portions are relatively weak with the upper transition zone showing only slightly larger negative amplitudes. Negative scattering is consistently largest at all depths under the eastern US. We also note a linear grain in both amplitude maps orienting from northwest to southeast. The direction of this grain is intriguing as it is approximately perpendicular to Farallon–North America plate motion over the past 70 Ma. That geometry is illustrated in Fig. 3 by drawing one of the flow lines computed by Pavlis et al. (2012) to about 30 Ma. Simulation with a model of regularly distributed point scatterers showed that amplitude from irregular coverage of the data cannot produce patterns of such scale (Fig. S1). This observation hints that the processes that are producing the negative scattering are being aligned by the mantle flow field. Earlier work (Wang and Pavlis, 2016a) showed a similar pattern for the fabric of roughness observed on the d410 and d660. This is particularly compelling in light of recent work by Tauzin et al. (2017) who directly link negative scattering bodies in the transition zone to active subduction in the western Pacific as well as the work by van der Lee et al. (2008) showing evidence of hydrous upwellings originated from past episodes of subduction of the Farallon Plate.

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

45

Fig. 2. Average amplitude analysis. (A) Map view of the average radial component amplitude between d410 and d660 surfaces derived from the image volume. The histogram is the distribution of the average amplitude values in (A). (B) The same average amplitude map from the radially symmetric layer simulation with added white noise. (C) A comparison of the stacking coherence from the 3D seismic image (red) and the radial symmetric layer simulation (blue). The thick lines are the median stack of coherence, and the shaded area comprises one quartile in each direction from the median.

Fig. 3. Map view of the average amplitude of the upper transition zone (left) and lower transition zone (right). The upper and lower portions are divided by a constant depth of 520 km. The thick dashed line is from the work by Pavlis et al. (2012), which is a flow line from the interpreted Farallon slab interface based on prior plate motions and tomography models with an origin at the Pacific–North America–Farallon (Gorda) triple junction near Cape Mendocino, California. It is used as a reference for relative plate motion direction.

4. Model simulations Previous studies all considered the negative phases seen in the transition zone as minor discontinuities that appear at certain depths. Our results suggest this is likely an artifact of the CCP stacking method’s tendency to produce horizontal layers from form of the migration impulse response (Pavlis and Wang, 2015) for that

method. Resolution tests in previous studies show that our method is capable of resolving point heterogeneities and dipping structures to the scale of the resolution limit imposed by diffraction (Liu and Pavlis, 2013; Wang and Pavlis, 2016a). As noted above multiple methods we applied to follow these negative horizons indicated they are always irregular shapes in 3D. Since there are no continuous horizons that exceed the resolution limit, it suggests that

46

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

Fig. 4. Models of heterogeneities distributed in the mantle transition zone. (A) Model of randomly distributed blobs treating the low-velocity bodies as point scatterers. (B) Model of near vertical flows treating the observation as artifact due to limited dipping resolution.

previous attempts of interpreting them as multiple minor discontinuities may not be correct. Scattered wave imaging methods do a poor job of retrieving seismic velocities compared to tomography methods. However, given that the negative amplitudes we observe are as strong as the d410 and d660 signal, it implies these negative bodies could have up to ∼4% lower S wave velocity (V S ) than the surrounding mantle, as that is the scale of change in velocity at d410 and d660 in all standard earth models. One could hypothesize that such low-velocity bodies within the transition zone could result from either a compositional or thermal heterogeneity. Seismic tomography studies have argued for both, but tomography methods cannot distinguish the two different processes (Karato, 2011; Meier et al., 2009; Suetsugu et al., 2006). Our result based on the higher resolution scattered wave imaging gives a different perspective. The scattering process as well as the 0.1 Hz shaping wavelet from this dataset restricted these low-velocity bodies to be much smaller than ∼100 km in size (details in Supplementary Text 1.2). Note that this scale length is fundamentally different from resolution controls on seismic tomography models. Tomography models are always regularized by some form of smoothing to limit resolution to the order of several hundred kilometers. If structures were truly of the geometry defined in tomography models they would be invisible to scattered wave images where the scattering strength scales more closely with the gradient in velocities. We infer that these low-velocity bodies are relatively small, and the velocity gradient at their boundaries must be sharp. Thus, the lowvelocity bodies are unlikely to be a pure thermal anomaly due to the difficulty in maintaining sharp thermal gradient within a small volume. Fig. 4 illustrates two alternative models we considered to explain the geometry of the compositional heterogeneity imaged from the USArray. The simplest is a set of randomly distributed low-velocity blobs of slight difference in size (Fig. 4A). That conceptual model is founded on the “plum pudding” model of preserved heterogeneity proposed in previous work (Allègre and Turcotte, 1986; Carlson, 1988). The second model (Fig. 4B) is based on previous geodynamic modeling work (Motoki and Ballmer, 2015; Richard and Iwamori, 2010; van der Lee et al., 2008) that argue we can expect diapir structures in the transition zone created by density instability in water-rich heterogeneities.

Fig. 5. Map view of the average amplitude within the transition zone from the randomly distributed point simulation.

4.1. Scattering from distributed heterogeneities The concept illustrated in Fig. 4A is that the negative scattering we observe is created by widely-distributed, small, low-velocity bodies that scatter P waves to S waves but are too small to be visible by tomographic inversions. Here, we conducted experiments to test if this could explain the observations with a random scattering simulation. We built a model with a random distribution of point scatterers within the transition zone. These scatterers are randomly sampled with a probability following uniform distribution on the imaging grid that has a nominal horizontal spacing of 35 km and a vertical spacing of 2 km. Note that due to the spherical geometry of the imaging grid, the actual horizontal spacing is depth-dependent and less than the nominal spacing defined at the surface. Also, the scattering strengths of the scatterers are normally distributed and defined by the mean and standard deviation of amplitudes measured from the real data (Fig. 2A). The average amplitude map recovered from this simulation shows a pattern similar to the original image but considerably smoother (Fig. 5). This suggests the real data may have a pattern not captured by a random simulation. This result strengthens our observation that an alignment of higher amplitudes perpendicular to North America–Farallon plate motion could be a marker of the flow field.

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

47

Fig. 6. 3D view of the simulation models used to test the dip resolution of PWMIG in the transition zone. The densely-spaced spheres illustrated form lines colored with five different gray shades. The dips of the lines are 15◦ , 30◦ , 45◦ , 60◦ and 90◦ and scale progressively from dark to light gray with increasing dip. The thin black lines are state boundaries at the surface and the thin gray lines are the outline of the image grid. This is a three-dimensional view from above at a point south-southeast of the U.S.

According to the scattering classification in Aki and Richards (1980, Fig. 13.11), the scale length of scatterer imaged at the depth of the transition zone can be as small as ∼25 km (details in Supplementary Text 1.2). Due to the increasing possibility of multiple scattering at length scales smaller than the wavelength, however, the Born Approximation on which both the imaging and simulation methods are founded may no longer be valid. Random medium theory or multiple scattering is generally considered to properly describe such a medium. We suggest, however, that this does not imply that the negative amplitudes from our result are artifacts. Signals being scattered multiple times within a random medium should average out its preference on the polarity. The distinct preference for negative scatterers within the transition zone thus makes it less likely a case due to multiple scattering. In fact, the highly smoothed pattern from our random point simulation suggests that by increasing the spatial density of seismic stations as well as the frequency of the data, it is likely to approximate our observation with this “plum pudding” model that is limited by diffraction. Therefore, the persistent negative amplitude is consistent with a set of small low-velocity bodies randomly distributed in the transition zone. 4.2. Scattering from dipping structures Fig. 4B illustrates an alternative geometry that has been suggested for small-scale structure in the transition zone (Motoki and Ballmer, 2015). The ability to distinguish that geometry in our 3D images depends on our ability to resolve steeply dipping structures. Imaging a diapir, in fact, is comparable to the challenge of imaging a salt done in seismic reflection data. To understand the limitations of our technique we built a set of models with line structures of varying dips distributed within the transition zone. These models simulate a set of equally spaced dipping lines of low shear wave velocity (Fig. 6). The line structures are constructed from the superposition of point scatterers with a vertical spacing of 2 km extending from 410 km to 660 km. Each simulation has eighteen linear objects of the same dip angle. We ran five different

simulations with 15◦ , 30◦ , 45◦ , 60◦ and 90◦ dips. Each of the five resulting data sets were processed exactly like the real data. It is clear from these simulations that dips up to 30◦ can be well resolved (Fig. 7). The 45◦ dipping structures are resolvable but with serious deficiencies. This demonstrates that the current upper limit of our dip resolution is between 30◦ and 45◦ . We also find the resolution is more susceptible to data coverage variations at higher dip angles. The reason this happens is well known from seismic reflection imaging (e.g. Yilmaz, 2001). Increasing dip resolution puts increasing demands on how raw data are mapped to the image to reconstruct the scatterer geometry. Furthermore, a key to imaging steeply dipping structures at these depths is proper handling of turning rays. That capability is lacking in the current implementation of PWMIG yielding this dip resolution limit. A key point of these simulations is that in all the steeply dipping (≥45◦ ) simulations (Fig. 7), the top and bottom of the structure is accentuated and the region between is smeared across a wider area. The results become indistinguishable from a pair of point scatterers at the top and bottom of each line segment. This implies that with these data and this method we cannot distinguish steeply dipping structure from point scatterers. Slightly deflected or vertical flow as illustrated in Fig. 4B is a feasible model for the heterogeneities we are observing. 5. Heterogeneous distribution of water in the mantle transition zone Accumulation of subducted oceanic crust has been proposed previously to explain compositional heterogeneity in the transition zone (Shen and Blum, 2003). Tauzin et al. (2017) observed pervasive low-velocity scatters in the Northwest Pacific region where a stagnant slab is well constrained by tomography models and seismicity. The transition zone may contain some sporadically distributed subducted oceanic crust if small-scale convection occurs due to slab stagnation (Motoki and Ballmer, 2015). The oceanic crust has been found to be compositionally less dense than the background pyrolitic mantle at depths of ∼100 km into the lower

48

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

Fig. 7. Cross-sections of the dipping resolution test at 15◦ , 30◦ , 45◦ , 60◦ and 90◦ of dip. The cross-sections are cutting through the center of the image volume from west to east. The wiggle traces are filled on negative.

mantle below the transition zone (Ringwood, 1990), and as a result, a significant amount of crustal materials could be trapped in the uppermost lower mantle (Nakagawa et al., 2010). In addition to the subducted oceanic crust, slab subduction brings ∼10 times the volume of depleted lithosphere (e.g., harzburgites) to the deep mantle. It is thus hard to imagine a scenario that could trap significantly more volume of oceanic crust than the depleted lithospheric materials throughout the transition zone. Therefore, the low-velocity structures that universally presented in the transition zone beneath the US (Fig. 2) are difficult to explain by subducted oceanic crust alone. Here, we propose that the wide distribution of low-velocity bodies in the transition zone may be explained by water-enriched subducted harzburgites accumulated in the transition zone over Earth’s history and being influenced by the current flow field. Note that because the theoretical ambiguity discussed in Supplementary Text 1.1 that the plane wave boundary theory could also apply to our observation, the response from a low-velocity body could also appear as a negative phase followed by a positive. In that case, the diapir structures shown in Fig. 4B is preferred as that structure is likely to create a smoother velocity gradient at the bottom suppressing the positive phase in the image. To demonstrate the interpretation, we computed a geodynamic model to illustrate a mechanism for the formation of compositional heterogeneities through subduction of oceanic crust and depleted oceanic lithosphere. The details about the model can be found in the Supplementary Materials. The modeling results show that after cold slabs are subducted to the lowermost mantle, a significant amount of material

is entrained by upwelling mantle flows and the parent slab material is later stirred into the background mantle (Movie S5). Fig. 8 shows a snapshot of the temperature field and the composition field. A key point the model demonstrates for this paper is that at small scales the compositional field is much more heterogeneous than the temperature field. Compositional heterogeneities dominate the transition zone with the volume of harzburgites much larger than that of the basalts (Fig. 8D) due to the larger fraction of harzburgites (depleted mantle) in the slabs. The geometry of the accumulated harzburgites in the transition zone is time-dependent, and is a complex mixture of distributed blobs and folds (Fig. 8D, Movie S5), which is not unlike the small scattering objects illustrated in Fig. 4. These basic results are not expected to change by adding more complexities to the geodynamic model, such as the effects of phase transitions and compositionally-dependence of viscosity. Here, we show that the compositional heterogeneity in the transition zone can be formed by entraining the subducted oceanic crust and depleted lithosphere from the lowermost mantle to the transition zone and stirring them with the background mantle material. Similar compositional heterogeneities in the transition zone can also be produced by other mechanisms, such as through small-scale convection in the transition zone due to slab stagnation (Motoki and Ballmer, 2015). Although the transition zone could contain a significant amount of randomly distributed harzburgites being subducted, the shear wave velocity of dry harzburgites would be predicted to be higher than pyrolite, which is commonly used to define average back-

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

49

Fig. 8. Formation of compositional heterogeneity. (A) Snapshot of temperature field at t = 0.006. (B) Snapshot of the composition field at t = 0.006. (C) Zoomed-in of the temperature field in regions marked by green box in (A). (D) Zoomed-in of the compositional field in regions marked by green box in (B). The background mantle material, basalts, and harzburgites are shown by white, red and blue color, respectively. The horizontal black dashed lines in (A) and (B) mark the depth of 410 km and 660 km.

ground mantle. However, if the harzburgites are enriched in water, their shear wave velocity could be significantly reduced given its high olivine content (Baker and Beckett, 1999; Xu et al., 2008). It has been well established that the high-pressure polymorphs of olivine in the transition zone have high water capacity (Bell and Rossman, 1992; Bercovici and Karato, 2003; Hirschmann, 2006; Pearson et al., 2014; Ringwood and Major, 1967). It is possible for water to be carried into the transition zone with mechanisms such as the superhydrous phase B (Inoue et al., 2006; Shieh et al., 1998). A locally hydrous transition zone is also supported by the direct evidence from ringwoodite in diamond samples with ∼1 wt.% H2 O (Pearson et al., 2014) as well as a number of compelling indirect evidences (Jacobsen and van der Lee, 2006; van der Meijde et al., 2003). Because harzburgite has larger portion of olivine content than pyrolite (Stixrude and Lithgow-Bertelloni, 2012), it could be more water enriched than the pyrolite in the transition zone. Single-crystal elastic moduli measurement of hydrous wadsleyite and ringwoodite under transition zone temperature and pressure conditions (Mao et al., 2012, 2011) indicate that V S of wadsleyite is 3.5% lower with 1.93 wt.% H2 O and the V S of ringwoodite is 4.5% lower with 1 wt.% H2 O. These measurements are consistent with our estimation of ∼4% velocity decrease indicating the extensive presence of pockets of wadsleyite or ringwoodite with at least 1 wt.% H2 O in the transition zone. The greater negative amplitudes in the lower portion of the transition zone agree with the inference of wave velocities at 520 to 660-km depth being more sensitive to hydration than the upper portion (Mao et al., 2012). Thus, we propose that the observed low-velocity bodies in the entire transition zone are heterogeneously distributed, water-enriched subducted harzburgites. These are likely preserved compositional heterogeneity from all of Earth history with emplacement and hydration both having a long history. Such heterogeneities are most likely created by accumulation in a long-term process spanning Earth’s history and being influenced by the current flow field driven by the Farallon slab. Acknowledgements We sincerely thank Suzan Van der Lee and an anonymous reviewer whose comments helped improve this paper significantly. We acknowledge Schlumberger for providing access to the Petrel seismic interpretation software package. Data from the TA network were made freely available as part of the EarthScope USArray

facility, operated by Incorporated Research Institutions for Seismology (IRIS) and supported by the National Science Foundation, under Cooperative Agreements EAR-1261681. The cyberinfrastructure of this research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative at IU is also supported in part by Lilly Endowment, Inc. This material is based upon work supported by the National Science Foundation under Grant No. CNS-0521433. Support for this work came from the U.S. National Science Foundation Earthscope Program under award EAR-1358149. Appendix A. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.epsl.2018.10.010. References Aki, K., Richards, P.G., 1980. Quantitative Seismology: Theory and Methods, vol. II. W.H. Freeman & Co, San Francisco. Allègre, C.J., Turcotte, D.L., 1986. Implications of a two-component marble-cake mantle. Nature 323, 123–127. https://doi.org/10.1038/323123a0. Bagley, B., Courtier, A.M., Revenaugh, J., 2013. Seismic shear wave reflectivity structure of the mantle beneath northeast China and the northwest Pacific. J. Geophys. Res., Solid Earth 118, 5417–5427. https://doi.org/10.1002/jgrb.50385. Baker, M.B., Beckett, J.R., 1999. The origin of abyssal peridotites: a reinterpretation of constraints based on primary bulk compositions. Earth Planet. Sci. Lett. 171, 49–61. https://doi.org/10.1016/S0012-821X(99)00130-2. Bell, D.R., Rossman, G.R., 1992. Water in Earth’s mantle: the role of nominally anhydrous minerals. Science 255, 1391–1397. https://doi.org/10.1126/science.255. 5050.1391. Bercovici, D., Karato, S., 2003. Whole-mantle convection and the transition-zone water filter. Nature 425, 39–44. https://doi.org/10.1038/nature01918. Bonatto, L., Schimmel, M., Gallart, J., Morales, J., 2015. The upper-mantle transition zone beneath the Ibero-Maghrebian region as seen by teleseismic Pds phases. In: Special Issue on Iberia Geodynamics: An Integrative Approach from the Topo-Iberia Framework. Tectonophysics 663, 212–224. https://doi.org/10.1016/ j.tecto.2015.02.002. Burdick, S., Vernon, F.L., Martynov, V., Eakins, J., Cox, T., Tytell, J., Mulder, T., White, M.C., Astiz, L., Pavlis, G.L., van der Hilst, R.D., 2017. Model update May 2016: upper-mantle heterogeneity beneath North America from travel-time tomography with global and USArray data. Seismol. Res. Lett. 88, 319–325. https:// doi.org/10.1785/0220160186. Carlson, R.W., 1988. Layer cake or plum pudding? Nature 334, 380–381. https:// doi.org/10.1038/334380a0. Crotwell, H.P., Owens, T.J., 2005. Automated receiver function processing. Seismol. Res. Lett. 76, 702–709. https://doi.org/10.1785/gssrl.76.6.702.

50

Y. Wang et al. / Earth and Planetary Science Letters 505 (2019) 42–50

Emry, E.L., Nyblade, A.A., Julià, J., Anandakrishnan, S., Aster, R.C., Wiens, D.A., Huerta, A.D., Wilson, T.J., 2015. The mantle transition zone beneath West Antarctica: seismic evidence for hydration and thermal upwellings. Geochem. Geophys. Geosyst. 16, 40–58. https://doi.org/10.1002/2014GC005588. Fan, C., Pavlis, G.L., Tuncay, K., 2006. GCLgrid: a three-dimensional geographical curvilinear grid library for computational seismology. Comput. Geosci. 32, 371–381. https://doi.org/10.1016/j.cageo.2005.07.001. Fei, H., Yamazaki, D., Sakurai, M., Miyajima, N., Ohfuji, H., Katsura, T., Yamamoto, T., 2017. A nearly water-saturated mantle transition zone inferred from mineral viscosity. Sci. Adv. 3, e1603024. https://doi.org/10.1126/sciadv.1603024. Gao, S.S., Silver, P.G., Liu, K.H., Kaapvaal Seismic Group, 2002. Mantle discontinuities beneath Southern Africa. Geophys. Res. Lett. 29, 129-1. https://doi.org/10.1029/ 2001GL013834. Gao, Y., Suetsugu, D., Fukao, Y., Obayashi, M., Shi, Y., Liu, R., 2010. Seismic discontinuities in the mantle transition zone and at the top of the lower mantle beneath eastern China and Korea: influence of the stagnant Pacific slab. In: Special Issue on Deep Slab and Mantle Dynamics. Phys. Earth Planet. Inter. 183, 288–295. https://doi.org/10.1016/j.pepi.2010.03.009. Hirschmann, M.M., 2006. Water, melting, and the deep Earth H2 O cycle. Annu. Rev. Earth Planet. Sci. 34, 629–653. https://doi.org/10.1146/annurev.earth.34.031405. 125211. Houser, C., 2016. Global seismic data reveal little water in the mantle transition zone. Earth Planet. Sci. Lett. 448, 94–101. https://doi.org/10.1016/j.epsl.2016.04. 018. Huang, X., Xu, Y., Karato, S., 2005. Water content in the transition zone from electrical conductivity of wadsleyite and ringwoodite. Nature 434, 746–749. https:// doi.org/10.1038/nature03426. Inoue, T., Ueda, T., Higo, Y., Yamada, A., Irifune, T., Funakoshi, K., 2006. High-pressure and high-temperature stability and equation of state of superhydrous phase B. In: Jacobsen, S.D., van der Lee, S. (Eds.), Geophysical Monograph Series. American Geophysical Union, Washington, DC, pp. 147–157. https://doi.org/10.1029/ 168gm11. IRIS DMC, 2011. Data services products: EMC. https://doi.org/10.17611/DP/EMC.1. Jacobsen, S.D., van der Lee, S. (Eds.), 2006. Earth’s Deep Water Cycle, Geophysical Monograph Series. American Geophysical Union, Washington, DC. https:// doi.org/10.1029/gm168. Jasbinsek, J.J., Dueker, K.G., Hansen, S.M., 2010. Characterizing the 410 km discontinuity low-velocity layer beneath the LA RISTRA array in the North American Southwest. Geochem. Geophys. Geosyst. 11, Q03008. https://doi.org/10.1029/ 2009GC002836. Karato, S., 2011. Water distribution across the mantle transition zone and its implications for global material circulation. Earth Planet. Sci. Lett. 301, 413–423. https://doi.org/10.1016/j.epsl.2010.11.038. Ligorría, J.P., Ammon, C.J., 1999. Iterative deconvolution and receiver-function estimation. Bull. Seismol. Soc. Am. 89, 1395–1400. Liu, X., Pavlis, G.L., 2013. Appraising the reliability of converted wavefield imaging: application to USArray imaging of the 410-km discontinuity. Geophys. J. Int. 192, 1240–1254. https://doi.org/10.1093/gji/ggs088. Mao, Z., Jacobsen, S.D., Frost, D.J., McCammon, C.A., Hauri, E.H., Duffy, T.S., 2011. Effect of hydration on the single-crystal elasticity of Fe-bearing wadsleyite to 12 GPa. Am. Mineral. 96, 1606–1612. https://doi.org/10.2138/am.2011.3807. Mao, Z., Lin, J.-F., Jacobsen, S.D., Duffy, T.S., Chang, Y.-Y., Smyth, J.R., Frost, D.J., Hauri, E.H., Prakapenka, V.B., 2012. Sound velocities of hydrous ringwoodite to 16 GPa and 673 K. Earth Planet. Sci. Lett. 331–332, 112–119. https://doi.org/10.1016/j. epsl.2012.03.001. Meier, U., Trampert, J., Curtis, A., 2009. Global variations of temperature and water content in the mantle transition zone from higher mode surface waves. Earth Planet. Sci. Lett. 282, 91–101. https://doi.org/10.1016/j.epsl.2009.03.004. Motoki, M.H., Ballmer, M.D., 2015. Intraplate volcanism due to convective instability of stagnant slabs in the mantle transition zone. Geochem. Geophys. Geosyst. 16, 538–551. https://doi.org/10.1002/2014GC005608. Nakagawa, T., Tackley, P.J., Deschamps, F., Connolly, J.A.D., 2010. The influence of MORB and harzburgite composition on thermo-chemical mantle convection in a 3-D spherical shell with self-consistently calculated mineral physics. Earth Planet. Sci. Lett. 296, 403–412. https://doi.org/10.1016/j.epsl.2010.05.026. Pavlis, G.L., 2011a. Three-dimensional, wavefield imaging of broadband seismic array data. Comput. Geosci. 37, 1054–1066. Pavlis, G.L., 2011b. Three-dimensional wavefield imaging of data from the USArray: new constraints on the geometry of the Farallon slab. Geosphere 7, 785–801. Pavlis, G.L., Sigloch, K., Burdick, S., Fouch, M.J., Vernon, F.L., 2012. Unraveling the geometry of the Farallon plate: synthesis of three-dimensional imaging results from USArray. Tectonophysics 532–535, 82–102. https://doi.org/10.1016/j.tecto. 2012.02.008. Pavlis, G.L., Wang, Y., 2015. Shaping wavelet effects in scattered wave imaging of P to S conversion data. Geophys. J. Int. 203, 373–383. https://doi.org/10.1093/gji/ ggv163. Pearson, D.G., Brenker, F.E., Nestola, F., McNeill, J., Nasdala, L., Hutchison, M.T., Matveev, S., Mather, K., Silversmit, G., Schmitz, S., Vekemans, B., Vincze, L., 2014.

Hydrous mantle transition zone indicated by ringwoodite included within diamond. Nature 507, 221–224. https://doi.org/10.1038/nature13080. Poppeliers, C., Pavlis, G.L., 2003a. Three-dimensional, prestack, plane wave migration of teleseismic P -to-S converted phases: 1. Theory. J. Geophys. Res. 108. https:// doi.org/10.1029/2001JB000216. Poppeliers, C., Pavlis, G.L., 2003b. Three-dimensional, prestack, plane wave migration of teleseismic P -to-S converted phases: 2. Stacking multiple events. J. Geophys. Res. 108. https://doi.org/10.1029/2001JB001583. Richard, G.C., Iwamori, H., 2010. Stagnant slab, wet plumes and Cenozoic volcanism in East Asia. In: Special Issue on Deep Slab and Mantle Dynamics. Phys. Earth Planet. Inter. 183, 280–287. https://doi.org/10.1016/j.pepi.2010.02.009. Ringwood, A.E., 1990. Slab-mantle interactions: 3. Petrogenesis of intraplate magmas and structure of the upper mantle. Chem. Geol. 82, 187–207. https://doi.org/10. 1016/0009-2541(90)90081-H. Ringwood, A.E., Major, A., 1967. High-pressure reconnaissance investigations in the system Mg2 SiO4 –MgO–H2 O. Earth Planet. Sci. Lett. 2, 130–133. https://doi.org/ 10.1016/0012-821X(67)90114-8. Schmandt, B., 2012. Mantle transition zone shear velocity gradients beneath USArray. Earth Planet. Sci. Lett. 355–356, 119–130. https://doi.org/10.1016/j.epsl. 2012.08.031. Schmandt, B., Lin, F.-C., 2014. P and S wave tomography of the mantle beneath the United States. Geophys. Res. Lett. 41, 6342–6349. https://doi.org/10.1002/ 2014GL061231. Shen, X., Zhou, H., Kawakatsu, H., 2008. Mapping the upper mantle discontinuities beneath China with teleseismic receiver functions. Earth Planets Space 60, 713–719. https://doi.org/10.1186/BF03352819. Shen, X., Yuan, X., Li, X., 2014. A ubiquitous low-velocity layer at the base of the mantle transition zone. Geophys. Res. Lett. 41, 2013GL058918. https://doi.org/ 10.1002/2013GL058918. Shen, Y., Blum, J., 2003. Seismic evidence for accumulated oceanic crust above the 660-km discontinuity beneath southern Africa. Geophys. Res. Lett. 30, 1925. https://doi.org/10.1029/2003GL017991. Shen, Y., Sheehan, A.F., Dueker, K.G., de Groot-Hedlin, C., Gilbert, H., 1998. Mantle disco‘ntinuity structure beneath the southern East Pacific Rise from P-to-S converted phases. Science 280, 1232–1235. https://doi.org/10.1126/science.280. 5367.1232. Shieh, S.R., Mao, H., Hemley, R.J., Ming, L.C., 1998. Decomposition of phase D in the lower mantle and the fate of dense hydrous silicates in subducting slabs. Earth Planet. Sci. Lett. 159, 13–23. https://doi.org/10.1016/S0012-821X(98)00062-4. Smyth, J.R., Jacobsen, S.D., 2006. Nominally anhydrous minerals and Earth’s deep water cycle. In: Jacobsen, S.D., van der Lee, S. (Eds.), Geophysical Monograph Series. American Geophysical Union, Washington, DC, pp. 1–11. https://doi.org/ 10.1029/168gm02. Stixrude, L., Lithgow-Bertelloni, C., 2012. Geophysics of chemical heterogeneity in the mantle. Annu. Rev. Earth Planet. Sci. 40, 569–595. https://doi.org/10.1146/ annurev.earth.36.031207.124244. Suetsugu, D., Inoue, T., Yamada, A., Zhao, D., Obayashi, M., 2006. Towards mapping the three-dimensional distribution of water in the transition zone from P-velocity tomography and 660-km discontinuity depths. In: Jacobsen, S.D., Lee, S.V.D. (Eds.), Earth’s Deep Water Cycle. American Geophysical Union, pp. 237–249. https://doi.org/10.1029/168gm18. Tauzin, B., van der Hilst, R.D., Wittlinger, G., Ricard, Y., 2013. Multiple transition zone seismic discontinuities and low velocity layers below western United States. J. Geophys. Res., Solid Earth 118, 2307–2322. https://doi.org/10.1002/jgrb.50182. Tauzin, B., Kim, S., Kennett, B.L.N., 2017. Pervasive seismic low-velocity zones within stagnant plates in the mantle transition zone: thermal or compositional origin? Earth Planet. Sci. Lett. 477, 1–13. https://doi.org/10.1016/j.epsl.2017.08.006. van der Lee, S., Regenauer-Lieb, K., Yuen, D.A., 2008. The role of water in connecting past and future episodes of subduction. Earth Planet. Sci. Lett. 273, 15–27. https://doi.org/10.1016/j.epsl.2008.04.041. van der Meijde, M., Marone, F., Giardini, D., van der Lee, S., 2003. Seismic evidence for water deep in Earth’s upper mantle. Science 300, 1556–1558. https://doi.org/ 10.1126/science.1083636. Wang, Y., Pavlis, G.L., 2016a. Roughness of the mantle transition zone discontinuities revealed by high-resolution wavefield imaging. J. Geophys. Res., Solid Earth 121, 6757–6778. https://doi.org/10.1002/2016JB013205. Wang, Y., Pavlis, G.L., 2016b. Generalized iterative deconvolution for receiver function estimation. Geophys. J. Int. 204, 1086–1099. https://doi.org/10.1093/gji/ ggv503. Xu, W., Lithgow-Bertelloni, C., Stixrude, L., Ritsema, J., 2008. The effect of bulk composition and temperature on mantle seismic structure. Earth Planet. Sci. Lett. 275, 70–79. https://doi.org/10.1016/j.epsl.2008.08.012. Yilmaz, Ö., 2001. Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data, 2nd ed. Investigations in Geophysics. Society of Exploration Geophysicists, Tulsa, OK. https://doi.org/10.1190/1.9781560801580.