Slab Transport of Fluids to Deep Focus Earthquake Depths—Thermal Modeling Constraints and Evidence From Diamonds

Abstract

The nature and cause of deep earthquakes remain enduring unknowns in the field of seismology. We present new models of thermal structures of subducted slabs traced to mantle transition zone depths that permit a detailed comparison between slab pressure/temperature (P/T) paths and hydrated/carbonated mineral phase relations. We find a remarkable correlation between slabs capable of transporting water to transition zone depths in dense hydrous magnesium silicates with slabs that produce seismicity below ∼300‐km depth, primarily between 500 and 700 km. This depth range also coincides with the P/T conditions at which oceanic crustal lithologies in cold slabs are predicted to intersect the carbonate‐bearing basalt solidus to produce carbonatitic melts. Both forms of fluid evolution are well represented by sublithospheric diamonds whose inclusions record the existence of melts, fluids, or supercritical liquids derived from hydrated or carbonate‐bearing slabs at depths (∼300–700 km) generally coincident with deep‐focus earthquakes. We propose that the hydrous and carbonated fluids released from subducted slabs at these depths lead to fluid‐triggered seismicity, fluid migration, diamond precipitation, and inclusion crystallization. Deep focus earthquake hypocenters could track the general region of deep fluid release, migration, and diamond formation in the mantle. The thermal modeling of slabs in the mantle and the correlation between sublithospheric diamonds, deep focus earthquakes, and slabs at depth demonstrate a deep subduction pathway to the mantle transition zone for carbon and volatiles that bypasses shallower decarbonation and dehydration processes.

Introduction

Deep (300–700 km) earthquakes occur at pressure (P) and temperature (T) conditions that are well beyond those optimal for brittle failure, making their existence a long-standing problem in solid-Earth geophysics. The vast majority of these events occur within the downgoing oceanic lithosphere in subduction zones. While hypocentral depth distributions within individual slabs vary (Figures 1a and S1–S5), there is a global bimodal depth distribution characterized by a decreasing number of events between 70 and 300 km depths and a distinct increase of seismicity in the lower transition zone between 500 and 700 km depths (Figure 1b). In this study, we focus on this increase in seismicity especially in the lower half of this bimodal distribution (i.e., earthquakes below ∼450-km depth), to which we will refer simply as “deep earthquakes,” in contrast to the shallower (but still deep) “intermediate-depth seismicity” that occurs between 70 and 300 km depths.

Figure 1. Seismicity as a function of depth. (a) Thermal model of the Flores Sea subducting plate with earthquake hypocenters showing the bimodal earthquake distribution on a single subducting plate transect (also shown as transect 17 of Figure S4). “x-sec distance” is the distance from the trench. (b) Histogram of seismicity of 28 years of global seismic data from the NEIC catalog (1990–2018) that includes the Flores Sea plate data, the other 22 transects shown in Figures S1–S5 and all other subducting slabs. The comparison of the seismicity of a single plate and the histogram of global data is presented to highlight that the global depth distribution is an aggregate of the thermal properties and composition of all subduction plates.

The causes of intermediate-depth earthquakes and the deep earthquakes considered here are not necessarily the same. However, given the experimental challenges of studying brittle fracture processes at the pressures and temperatures typical for deep earthquakes, we consider mechanisms for both categories of earthquakes. Several reviews have been written on the characteristics and causes of deep earthquakes (e.g., Frohlich, 1989; Green & Houston, 1995; Houston, 2015; Zhan, 2020). Early theories about the physics behind deep earthquakes often considered volume-changing phase transformations for their genesis (Benioff, 1963; Dziewonski & Gilbert, 1974; Evison, 1967; Leith & Sharpe, 1936). This explanation is appealing due to the known presence of phase transformations at transition-zone depths where “normal” brittle faulting is excluded due to the very high confining pressures. However, deep earthquakes repeatedly have been characterized dominantly by double-couple focal mechanisms (i.e., simple slip on a planar fault surface) similar to those observed for shallow earthquakes (e.g., Kawakatsu, 1991; Kikuchi & Kanamori, 1994; Leith & Sharpe, 1936; Stimpson & Pearce, 1987), suggesting that large volumetric changes (i.e., collapse mechanisms) are not responsible for these events.

Modern theories of how brittle failure is possible at conditions present in slabs below 70-km depth fall roughly into three non-mutually exclusive categories. The first category broadly spans those mechanisms that rely on phase transformations that result in seismicity through volume decreases and/or the formation of nanocrystalline reaction products. For intermediate-depth earthquakes, this typically involves hydrous phases like lawsonite eclogite (Jung et al., 2004) or antigorite (Gasc et al., 2017; Incel et al., 2017; Plümper et al., 2017), though these mechanisms should not be confused with other fluid-related mechanisms such as “dehydration embrittlement” discussed below. For deep earthquakes, most of the focus has been on the development of a metastable olivine wedge within the transition zone (e.g., Green, 2007; Green & Burnley, 1989; Green & Houston, 1995; Kirby et al., 1996; Schubnel et al., 2013; Stein & Rubie, 1999; Y. Wang et al., 2017). These theories focus either on the volume decrease that is caused by the transformation of olivine which leads to collapse structures known as “anti-cracks,” or on the formation of interconnected segments of nanocrystalline spinel reaction products (e.g., Schubnel et al., 2013; Y. Wang et al., 2017) whose weakness can accommodate slip despite the tremendous confining pressures in the lower transition zone.

The second category involves shear heating and thermal run-away mechanisms that enable significant fault rupture due to high temperatures generated along a propagating fault tip, some of which also incorporate reduced grain size along discrete planes of weakness. These mechanisms are thought to operate both for intermediate-depth earthquakes (e.g., John et al., 2009; Kelemen & Hirth, 2007; Ohuchi et al., 2017; Prieto et al., 2013; Thielmann, 2018) and for deep earthquakes (e.g., Ogawa, 1987; Wiens, 2001). We note that other mechanisms may act as triggers that then allow for the propagation of faulting through this shear heating and thermal instability model (e.g., Wiens, 2001; Zhan, 2017, 2020).

The third category we will generally call “fluid-triggered deep seismicity.” The fate of water and other volatiles including carbon into subduction zones is a long-standing, heavily studied, but still insufficiently understood the problem (see Faccenda, 2014 and references therein). The potential role of these volatiles in generating intermediate-depth or deep seismicity is similarly ill-constrained. For the purpose of this study, we include in this category the commonly cited mechanism of “dehydration embrittlement,” in which the release of water from dehydration reactions results in an increase in pore pressure that is sufficient to reduce the normal stresses to the point where brittle failure is possible (e.g., Meade & Jeanloz, 1991; Raleigh & Paterson, 1965). However, even at intermediate depths, dehydration reactions do not necessarily produce net increases in volume (e.g., Dobson et al., 2002; Gasc et al., 2017; Jung et al., 2004). We therefore include other possible mechanisms related to the presence of fluids, whether generated by local reactions or by fluids channelized to updip segments of the slab (e.g., Shiina et al., 2017; Wilson et al., 2014). We note that there is overlap here with phase transformation mechanisms, in that some fluid-triggered intermediate-depth seismicity may form due to dehydration reactions that, in addition to releasing fluids, produce nanocrystalline reaction products that can lead to extreme weakening and the transfer of stress to the surrounding rock matrix, resulting in brittle failure (e.g., Ferrand, 2019; Ferrand et al., 2017; Jung et al., 2004). Similar mechanisms that involve both dehydration reactions and weakening due to nanocrystalline reaction products have not, to our knowledge, been explored specifically for deep earthquakes, though our results may warrant such an investigation.

The ability of any fluid-triggered seismicity mechanisms to explain deep earthquakes hinges first on the ability of subducting slabs to transport hydrous phases to sufficient depth. This depends on the pressure/temperature path of the subducting plate which, in turn, depends on the temperature distribution in the slab before it began subducting (and therefore the age of the slab at the trench) and its tectonic history before reaching a given depth (which is a function of, among others, convergence rate and dip angle). Thermal modeling has suggested that, at least in the cases of the warmest slabs (young slabs with shallow dips and/or low convergence velocities), dehydration is largely complete above about 100-km depth (van Keken et al., 2011). In most subduction zones, however, only the crust is expected to dehydrate by those depths (Komabayashi et al., 2004; Maruyama & Okamoto, 2007; Ohtani, 2015; Ohtani et al., 2004; Omori et al., 2004; Thompson, 1992). The amount of water present in the mantle lithosphere of subducting slabs is a topic of ongoing research, but a number of recent studies highlight potentially tens of kilometers of hydration below the oceanic Moho made possible by extensional outer-rise faulting as the plate bends to enter the trench (Bloch et al., 2018; Boneh et al., 2019; Cai et al., 2018; Faccenda, 2014; Fromm et al., 2006; Grevemeyer et al., 2018; Ranero et al., 2003; Wagner et al., 2020). The fate of this deeper reservoir of water in descending plates has been the topic of significant research (Barcheck et al., 2012; Hacker, 2008; Syracuse et al., 2010; van Keken et al., 2011) though often primarily at depths <300 km.

Discoveries in diamond petrology have provided insights into these long-standing questions. Studies of sublithospheric diamonds (that crystallize at ∼300–700 km depths) demonstrate the ability of slabs to transport and release aqueous fluids, metallic liquids, and/or carbonatitic melt at the depths of deep earthquakes (Burnham et al., 2015; Palot et al., 2014, 2016; Pearson et al., 2014; Shirey et al., 2019; Smith et al., 2016, 2018; Thomson et al., 2014; Thomson, Walter, et al., 2016; Walter et al., 2008). Mineral inclusions in some sublithospheric diamonds have been shown to be compositionally compatible with the high-pressure phases expected to be present in subducted mid-ocean ridge basalts (Stachel, 2001a; Thomson et al., 2014; Walter et al., 2011), indicating that these diamonds crystallized within or near the oceanic crust of subducted slabs. Mineral inclusions in other diamonds suggest an origin from the slab mantle (Regier et al., 2020; Smith et al., 2016, 2017, 2018, 2021). In many cases, reconstructed high-pressure phase assemblages of sublithospheric inclusions indicate that they were encapsulated by the diamond within the mantle transition zone and near the upper/mantle lower mantle boundary (e.g., Harte, 2010; Harte et al., 1999; Hutchison et al., 2001; Stachel et al., 2005). This has provided a way forward to use mineralogy to trace Earth's deep water and carbon cycles (Harte, 2010; Plank & Manning, 2019; Walter et al., 2008). These studies have provided clearer estimates of the types of fluids involved, the depths of diamond crystallization and fluid release, and the host rocks from which the fluids were derived (Figures 2a–2c; Table 1).

Figure 2. Examples of sublithospheric diamonds from the transition zone, mineral inclusions, and evidence for aqueous and other fluids. (a) Type IIa (see footnote a in Table 1 for diamond type definitions) diamond (left) with tomographic image of ringwoodite and breyite (interpreted to be retrograded Ca-perovskite; center) and infrared spectra demonstrating ∼1.5 wt% water in the ringwoodite from Pearson et al. (2014). (b) Type IIb diamond (left) with a scanning electron microscope image (SEM; center) of a composite inclusion inclusion of jeffbenite and NaAl-pyroxene (Smith et al., 2018) that we interpret to be a deep mantle mineral such as bridgmanite (see Armstrong & Walter, 2012; Hutchison et al., 2001). Other inclusions (right) are routinely jacketed in fluid envelopes of CH4 or H2 visible with Raman spectroscopy suggesting hydrogen diffusion out of inclusions that were equilibrated with a hydrous environment. Figures from Smith et al. (2018). (c) Type IIa diamond (left) containing black metal inclusions which SEM images (center) reveal to be exsolved assemblages of cohenite, pyrrhotite, and Fe-Ni metal. As in (b), Raman spectroscopy of indicates fluid jackets of CH4 and H2 that show equilibration of metallic liquid with a hydrous environment. Figures and interpretation from Smith et al. (2016).

Table 1 Summary of Sublithospheric Diamond Types, Mineral Inclusions, Depths, and Isotopic Compositions

Fluid or melt type	Diamond typea	No. of diamondsb	Inclusions at Earth’s surface	Inclusions at mantle depths	Estimated depth (km)	Diamond δ13C	References
Carbonatiticc	IIa IaB	90 (500)	garnet+cpx breyite+perovskite spinel+nepheline spinel+kalsilite opx+jeffbenite opx/spinel+olivine coesite+kyanite	majorite Ca-perovskite CF phase NAL phase Al-bridgmanite bridgmanite stishovite	300–700	−3.5‰ to −28.7‰	Brenker et al. (2005) Bulanova et al. (2010) Burnham et al. (2015) Harte et al. (1999) Harte (2010) Hutchison et al. (2001) Thomson et al. (2014) Thomson, Walter, et al. (2016) Walter et al. (2008, 2011) Zedgenizov et al. (2014)
Metallicd	CLIPPIR IIa rare IaB	83 (∼8,000)<	Fe3C+Fe-Ni+FeS breyite garnet+cpx	metallic liquid Ca-perovskite majorite	360–750	−3.8‰ to −26.9‰	Smith et al. (2016)
Aqueous	IIb (blue)e	46 (∼2,000)	breyite+larnite opx+olivine coesite+kyanite ferropericlase opx+jeffbenite FeS+Fe-carbide+FeO	Ca-perovskite bridgmanite stishovite ferropericlase Al-bridgmanite metallic alloy	360–750+f	−1.8‰ to −20.8‰	Smith et al. (2018)
	IIa	1 (6)	ferropericlase magnesioferrite brucite	ferropericlase	300–700g		Palot et al. (2016)
	IIa	1 (400)	ringwoodite breyite	ringwoodite Ca-perovskite	300–660h		Pearson et al. (2014)
	IIa	1	AlSiO3(OH) stishovite fluid	AlSiO3(OH)	320–660i		Wirth et al. (2007)
	IaAB	13	Ice VII (Na,K)Cl ilmenite Fe, Ni, C	aqueous fluid ilmenite Fe, Ni, C	330–650j		Tschauner et al. (2018)

^aType II diamonds are very low nitrogen or nitrogen-free, designated as “IIa” if they are boron-free and “IIb” if they are boron bearing. CLIPPIR diamonds are a subset of type IIa diamonds as defined by Smith et al. (2016). IaA diamonds have nitrogen substituting in pairs and IaB diamonds have nitrogen substituting in clusters of 4 nitrogen atoms with a vacancy. IaAB diamonds have both substitution mechanisms (e.g. Breeding and Shigley, 2009).

^bNumber of inclusion-bearing diamonds as published is given. Number of sublithospheric diamonds examined to find inclusions for study is in parentheses.

^cMelt compositions for the carbonatitic melt types are inferred petrologically from mineralogy of the inclusion assemblage and their trace element compositions where available (Bulanova et al., 2010; Walter et al., 2008). Depths estimated from at-pressure analog silicate mineral assemblages, exsolved majorite garnet compositions, and high-pressure majorite breakdown depths. Significant water content has been measured on CaSiTi-perovskite (2,300 ppm) and majorite (900 ppm) by SIMS (Hauri, Thomson, Shirey, Walter unpublished data).

^dCLIPPIR diamonds carry frozen samples of metallic liquids allowing the metallic liquid compositions to be analyzed directly. Total number of specimens from two studies of Smith et al. (2016, 2017). Depths estimated from at-pressure analog silicate mineral assemblages extrapolated to entrapment depths.

^eIIb diamonds are inferred to crystallize from aqueous fluids based on their boron content (associated with seawater and serpentinization) and the hydrogen-rich jackets (CH₄ ± H₂) that surround the silicate minerals (Smith et al., 2018).

^fEstimated depths from inclusion mineralogies extrapolated to entrapment depths and the greatest minimum calculated remnant pressures of 14.1 GPa in ferropericlase based on X-ray diffraction measurement of ferropericlase unit cell and a mantle temperature at entrapment of up to 1,727°C (Angel et al., 2017; Smith et al., 2018).

^gPressure estimated from thermdynamic modeling of phase diagram for magnesioferrite and ferropericlase (Palot et al., 2016).

^hPressure estimate derived from stability field for ringwoodite. Ringwoodite contained 1.5 wt% water (Pearson et al., 2014).

ⁱPressure estimate derived from stability field for AlSiO₃(OH) (also known as “Phase Egg”) under mantle conditions. Fluid estimate based on high porosity at room pressure (Wirth et al., 2007).

^jPressure estimates made using residual pressures, equation of state of Ice VII and H₂O-based fluid intersection with mantle adiabat (Tschauner et al., 2018).

Fluids in the Transition Zone—Evidence From Sublithospheric Diamonds

Sublithospheric diamonds are brought to Earth's surface by kimberlitic magmas sourced from near mantle transition zone depths. Diamond carbon composition, nitrogen abundance variations, and mineral inclusion chemistry record important features from the mantle transition zone region such as the depth of diamond crystallization, growth from migrating fluids, and the composition of these fluids in relation to host rocks in the slab.

Depths of Origin

The depths of formation of sublithospheric diamonds (Table 1) are estimated either directly from majoritic garnet mineral barometry or based on pressure-temperature phase relations applied to specific inclusion chemical stoichiometries and/or mineral assemblages occurring in single diamonds (e.g., Bulanova, 1995; Harte & Cayzer, 2007; Harte et al., 1999; Hayman et al., 2005; Hutchison et al., 2001; Kaminsky et al., 2001; Smith et al., 2016, 2018; Stachel, et al., 2000a; Stachel, et al., 2000b; Stachel et al., 2005; Thomson et al., 2014; Walter et al., 2008, 2011). Majoritic garnet inclusions have provided the most numerous estimates of the depths of sublithospheric diamond and inclusion formation based on empirical barometers that predict pressure as a function of mineral composition (Beyer & Frost, 2017; Collerson et al., 2010; Harte & Cayzer, 2007; Stachel, et al., 2000a; Thomson et al., 2021).

Harte and Cayzer (2007) documented clinopyroxene exsolution from majoritic garnet inclusions in superdeep diamonds and used reconstructed bulk compositions to estimate inclusions pressures, recognizing that using the uncorrected garnet composition yields minimum pressures (i.e., more pyroxene component dissolves into garnet at high pressure). Thomson et al. (2021) noted that most reported majoritic garnet inclusion compositions do not account for clinopyroxene exsolution (see inset to Figure 3), thereby causing a systematic underestimation of pressure that can be quantified. In their case study, Thomson et al. (2021) used reconstructions of eight majorite-cpx inclusions from the Juina region of Brazil (four from Harte and Cayzer [2007]) and identified a systematic underestimation in pressure ranging from ∼2 to 8 GPa with a mean of about 4 GPa. Based on a regression to the observed pressure deviations, a pressure-dependent model correction was obtained.

Figure 3. Histograms (bottom) showing pressures calculated for majoritic garnets in a global data set of 221 inclusions using the machine-learning barometer of Thomson et al. (2021). The blue histogram shows results for reported compositions. The green histogram shows a model pressure-dependent correction for exsolution based on the the results presented in Thomson et al. (2021). The inverted histogram on the top shows subduction zone earthquakes from Figure 1b in 50 km bins with orange-brown bars. Inset images show two examples of pyroxene exsolution in majoritic garnet from the Juina region, Brazil (Thomson et al., 2014; see also Harte & Cayzer, 2007).

Figure 3 shows pressures calculated for a global database of majoritic garnet inclusion compositions (n = 221) based on the recent machine-learning barometer calibrated against more than 500 experimental majoritic garnet compositions (Thomson et al., 2021). The reported (uncorrected for exsolution) garnet compositions record pressures of about 6–22 GPa, with two distinct modes in the distribution at around 8 and 13 GPa, with the higher pressure mode predominant. When the Thomson et al. (2021) model exsolution correction is applied to the database, pressures shift upwards to about 13–23 GPa (∼400–600 km), with modes around 14 and 18 GPa. If these corrections are correct, then majoritic garnet inclusions apparently formed almost exclusively within the transition zone. We also note that nearly all inclusions occupying the higher pressure mode in the global database are low-Cr, low-Mg# inclusions with an eclogitic paragenesis, which based on major and trace element characteristics have been attributed to equilibration with a melt or fluid derived from subducted oceanic crust (Bulanova et al., 2010; Thomson, Kohn, et al., 2016; Thomson, Walter, et al., 2016; Walter et al., 2008).

If the original high-pressure mineral structure of an inclusion is preserved to the surface, it can also provide direct information about the depth of diamond formation based on mineral phase relations. For example, a specimen of pristine Mg-rich ringwoodite places its host diamond unquestionably in the lower half of the mantle transition zone when the inclusion was incorporated by its diamond host (e.g., Pearson et al., 2014). Conversely, most inclusions have experienced retrograde transitions (Table 1) either to lower pressure polymorphs of the same composition and stoichiometry such as MgSiO₃, CaSiO₃, (Mg_xFe_1−x)O (Burnham et al., 2016; Davies et al., 2004; Harte et al., 1999; Hutchison et al., 2001; Stachel, et al., 2000b) or to unmixed-assemblages (e.g., garnet-pyroxene, perovskite-breyite, enstatite-jeffbenite) depending on composition and phase relations (Brenker et al., 2002; Bulanova et al., 2010; Harte & Cayzer, 2007; Hayman et al., 2005; Kaminsky et al., 2001; Nestola et al., 2016; Walter et al., 2008, 2011; Zedgenizov et al., 2014). Phase associations in single diamonds such as the co-occurrence of MgSiO₃ as low-Ni enstatite (former bridgmanite) and CaSiO₃ as breyite (former Ca-perovskite), jeffbenite, and ferropericlase have been interpreted as indicative of a derivation at deep transition zone and upper lower mantle depths (Burnham et al., 2016; Davies et al., 2004; Harte et al., 1999; Hutchison et al., 2001; Stachel, et al., 2000b). While the sole occurrence of some of these phases, such as breyite, could be interpreted to indicate a deep asthenosphere origin (e.g., Woodland et al., 2020), the co-existence of such phases with ringwoodite (Pearson et al., 2014), MgSiO₃, and ferropericlase is suggestive of a higher pressure origin for many breyite inclusions.

Other sublithospheric diamonds carry composite assemblages recrystallized from original high-pressure mineralogy to lower pressure conditions. Composite inclusions with the bulk stoichiometry of Ti-rich Ca-perovskite have major and trace elements having chemical compositions associating them with subducted oceanic crustal components, and carbonated melts derived from oceanic crust in particular, although phase relations do not provide tight constraints on formation depths (Brenker et al., 2005; Bulanova et al., 2010; Thomson et al., 2014; Thomson, Walter, et al., 2016; Walter et al., 2008). Composite inclusions with the bulk stoichiometry of aluminous bridgmanite, CF-phase, and NAL-phase have been interpreted to have formed in relation to subducted ocean crust at conditions near the upper-mantle/lower-mantle boundary (Thomson et al., 2014; Walter et al., 2011). Overall, the recrystallized mineral inclusions carried by many sublithospheric diamonds, when ascribed to their original mineralogy, indicate depths of diamond crystallization between ∼300 and ∼700 km (Table 1).

Diamond Growth: Evidence for Fluids and Fluid Movement

Observations and experiments indicate that both lithospheric and sublithospheric diamonds form by crystallization from fluids rather than by solid-state transformation from graphite. Internal growth textures revealed by cathodoluminescence imaging show that many lithospheric diamonds have thin, intricate, concentric growth zoning, as well as ample evidence of resorption and recrystallization (e.g., Bulanova, 1995; Shirey et al., 2013). Both textures indicate crystallization from a fluid (e.g., Sunagawa, 1984). Fractured diamonds can be re-cemented with more diamond crystallization, which is also indicative of multiple pulses of fluid-mediated growth (Czas et al., 2018). The changing carbon and nitrogen isotopic composition and nitrogen abundance of diamond from core to rim is chemical evidence for progressive growth from a fluid (e.g., Boyd et al., 1987; Smart et al., 2011).

Sublithospheric diamonds have been observed to possess intricate growth layers similar to lithospheric diamonds (e.g., Bulanova et al., 2010; Hayman et al., 2005; Palot et al., 2017; Thomson et al., 2014) that are evidence of fluid-mediated growth. Just as in lithospheric diamonds, coupled small-scale C isotope and N abundance measurements in sublithospheric diamonds can be explained by minor isotopic fractionation during diamond growth from fluids or melts (Palot et al., 2017). In sublithospheric diamonds, zoning may be disturbed or overprinted by subsequent deformation textures (Shirey et al., 2013; Smith et al., 2016, 2018; Thomson et al., 2014).

Textural and chemical evidence shows that the continued growth of diamonds is the product of mobile fluids flowing through localized channels. Tomographic imaging of some large eclogitic mantle xenoliths from the continental lithosphere shows diamonds positioned along grain boundaries, that are sometimes clearly aligned with intergranular planes that likely served as fluid pathways (Anand et al., 2004; Czas et al., 2018; Ishikawa et al., 2008; Liu et al., 2009; Spetsius & Taylor, 2003). While no diamondiferous mantle xenoliths from the below the lithosphere are known, evidence from majorite inclusions with compositions intermediate between basalt and peridotitic endmembers (Thomson, Walter, et al., 2016) indicates reactions between the mantle and diamond-forming fluids, suggesting that mantle rocks are hosting mobile fluids in veins.

From the study of sublithospheric diamonds and their mineral inclusions, three different types of diamond-forming fluids have been postulated (Figure 2; Table 1): (1) carbonatitic melts (e.g., Walter et al., 2008); (2) C-O-H supercritical fluids containing water and methane (e.g., Smit et al., 2016; Stachel & Luth, 2015); and (3) Fe-Ni-S-C metallic liquids (e.g., Smith et al., 2016). Diamond precipitation can only occur if the species of carbon in the fluids (such as CO₂, CH₄, CO₃²⁻, and HCO) encounters changes in temperature (Luth & Stachel, 2014; Stachel & Luth, 2015), pH (Sverjensky & Huang, 2015), redox conditions (e.g., Frost & McCammon, 2008; Rohrbach & Schmidt, 2011; Stagno et al., 2013; Thomson, Walter, et al., 2016), or some combination thereof. Therefore, for diamonds to nucleate and grow, fluids must be able to permeate into the surrounding solid rock matrix. In summary, from textures, mineral inclusion compositions, and experimental-theoretical studies of fluid compositions, the sublithospheric diamond itself is perhaps the best tangible piece of physical evidence that freely migrating fluids are present between ∼300 and 700 km depths.

A Subducted Oceanic Slab Source for Sublithospheric Diamond-Forming Fluids

Diamond and inclusion compositions provide not only constraints on the depth at which the diamonds form, but the origin of the fluids from which they are grown. These fluids, in turn, vary depending on the type of diamond in question. The fluids from which carbonatitic sublithospheric diamonds are derived are most closely associated with carbonate-bearing oceanic crust. These diamonds contain mineral inclusions such as Ca-rich, Cr-poor majorite, and Ti-rich CaSiO₃ phases (breyite plus perovskite) and are thought to crystallize when carbonatitic melts derived from the basaltic portion of the slab react with the mantle (Thomson, Walter, et al., 2016; Walter et al., 2008). This is consistent with trace element patterns (Bulanova et al., 2010; Stachel, 2001; Stachel et al., 2005; Thomson et al., 2014; Walter et al., 2008) and heavy oxygen isotopic compositions exhibited by inclusions (Burnham et al., 2015; Ickert et al., 2013; Regier et al., 2020) that support derivation from altered oceanic crust. Indeed, such very heavy oxygen isotopic compositions (δ¹⁸O from +6‰ to +10‰) are at the extreme end of the seawater altered basalt compositional range and could be derived directly from carbonate melts released from altered oceanic crust (Regier et al., 2020). These sublithospheric diamonds are characterized by isotopically variable, typically light carbon with a composition (δ¹³C as low as −28.7‰) far below that of the average ambient mantle composition in δ¹³C of around −5‰ (Table 1; Bulanova et al., 2010; Palot et al., 2014, 2017; Shirey et al., 2019; Thomson, Walter, et al., 2016; Walter et al., 2008, 2011). This isotopically light carbon is best explained by the presence of carbonate formed by the oxidation of organic carbon in the deep, altered oceanic crust (e.g., Li et al., 2019).

The fluids from which metallic and aqueous sublithospheric diamonds are derived are compositionally and isotopically linked to oceanic lithospheric mantle that was serpentinized before subduction. These diamonds contain mantle minerals such as former bridgmanite and Ca-perovskite, ferropericlase, ringwoodite (Harte, 2010; Nestola et al., 2018; Pearson et al., 2014; Smith et al., 2018), and hydrated equivalents such as brucite (Palot et al., 2016). The high content (ppm levels) of the crustal element, boron, in blue-colored, aqueous sublithospheric diamonds is a consequence of the high partition coefficient for boron in serpentine during seawater alteration on the ocean floor. Similarly, the carapace of methane commonly found around mineral inclusions in aqueous diamonds (Smith et al., 2016, 2018) results from the abundant methane produced in the highly reducing environment of serpentinized peridotite. Sublithospheric diamonds derived from metallic and aqueous growth media are also characterized by isotopically light carbon (δ¹³C as low as −26.9‰ and −20.8‰, respectively) that ranges to far below the average ambient mantle (Table 1; Smith et al., 2016, 2018). This isotopically light carbon is likely derived from organic compounds in the altered peridotite (e.g., Ménez et al., 2018). Recent iron and osmium isotopic analyses show that metallic inclusions in CLIPPIR diamonds are much heavier (δ⁵⁶Fe = + 0.79 to +0.90) and significantly less radiogenic (¹⁸⁷Os/¹⁸⁸Os = 0.1109 to 0.1115) than the "normal" oceanic mantle ranges (i.e. δ⁵⁶Fe = −0.2 to +0.2; ¹⁸⁷Os/¹⁸⁸Os = 0.124 to 0.129), which is explained by derivation from magnetite and awaruite formed by low-temperature, seafloor serpentinization of oceanic lithospheric peridotite (Smith et al., 2021). Thus, the overwhelming evidence from mineral inclusions and stable isotopic compositions of sublithospheric diamonds show that the fluids from which they formed originate from the crust and mantle portions of subducted oceanic slabs.

Methods

Significant study has been done to model the temperatures of subducting slabs over a range of different spatial scales to investigate fluid migration and the causes of intermediate-depth and deep seismicity. Previous study to predict subduction zone thermal properties worldwide focused on a series of transects with laterally averaged slab geometries designed to study properties of slabs in the uppermost 250 km of the Earth, averaged over 500-km wide swaths along strike (e.g., Barcheck et al., 2012; Syracuse & Abers, 2006; Syracuse et al., 2010). We follow a different approach here by focusing on transects that best allow us to compare the thermal properties of slabs to the depth distribution of deep slab seismicity. Here, we select transects with the following attributes: (1) a clear connection in seismicity and/or tomographic imaging between the slab entering the trench and the slab passing into/through the mantle transition zone (i.e., no torn or orphaned slabs); and (2) an absence of significant variability in slab geometry within 100 km along strike of the trench, which is important for our two-dimensional (2D) thermal modeling and for the projection of seismicity along strike onto a 2D profile. Uncertainties in thermal models still exist due to numerous assumptions, including that the age of the slab at the trench and the convergence rates have not changed over the time (∼5–23 Myr) required for the slabs to reach their maximum modeled depths.

Selection of Subduction Zone Profiles

The locations of the slab profiles we use in this study are shown in Figures 4 and S1–S5 and Table S1. Our selection of profile locations is inspired by those in Syracuse and Abers (2006) and Syracuse et al. (2010) but modified to better suit the needs for this study. In some cases, transects were omitted or moved along trench strike to ensure that all of our modeled transects complied with the two aforementioned constraints. For example, transects along the Philippine Trench appear to show deep earthquakes, but seismic tomography studies indicate that these may be due to the completely subducted Molucca Sea slab, not the currently subducting Philippine Sea plate (e.g., Fan & Zhao, 2018). Further to the north, where the Philippine plate subducts along the Ryukyu trench, seismic imaging and plate reconstructions suggest the down-going slab collides with a stagnated portion of the Pacific plate within the transition zone, making the attribution of deep seismicity to one slab or the other difficult (Wu et al., 2016). To the west, we do include a transect through the Bonin subduction zone for completeness, though recent work highlights both the complexity of the geodynamics produced by the en-echelon arrangement of subducting slabs and the presence of a nearby slab tear (Faccenna et al., 2018; Zhang et al., 2019). The northern end of the Tonga-Kermadec trench comprises a folded slab that possibly collides with detached portions of the Australian slab subducting beneath the Vanuatu trench (e.g., Richards et al., 2011). In this case, we have included one transect just south of where the putative interaction between these slabs is located, but the complexity of the region may not be captured in our 2D models. Further to the west, the subducting slab along the Solomon trench beneath the Bismark Sea is heavily contorted and likely torn in the lower part of the upper mantle (Holm & Richards, 2013) which adds three-dimensional (3D) complexity that is well beyond the scope of this study. In some cases, such as the westernmost end of the Sumatra arc, the slab subducts normally in the shallow upper mantle, but is heavily deformed in the lower upper mantle and transition zone (Pesicek et al., 2008).

Figure 4. Locations of the 23 subduction zone profiles used in this study. For detail of profiles 1–22, see Figures S1–S5. See main text for selection criteria and Table S1 for naming and precise locations.

A challenging subduction zone to model with 2D profiles and with assumptions on the invariance of plate age, convergence rate, and geometry lies along the western margin of South America. Here, the down-going Nazca plate changes from flat slab geometry to normally dipping slab geometry multiple times along strike. In addition, past episodes of flat slab subduction have been inferred along most of the South American Margin since the early Miocene (see review by Ramos and Folguera [2009]). However, South America is also the one region where a comparatively young slab appears to generate deep seismicity along a significant portion of its margin. We have selected three transects for which the current slab geometry does not vary within 100 km along the strike of the trench. Only one of these (Bolivia) samples the band of deep seismicity that extends from the southern margin of the Peruvian flat slab to the northern margin of the Pampean flat slab in central Chile/Argentina. We note that the central Andes at the latitude of our transect likely experienced an episode of flat subduction from 18 to 12 Ma (James & Sacks, 1999; Ramos & Folguera, 2009). The slab that is currently at the bottom of the transition zone is likely subducted >15 Ma given the current convergence rate and slab geometry. The part of the slab above that would have subducted more recently, meaning that the Nazca plate currently located in and above the mantle transition zone was probably originally part of a flat slab that has since foundered and is now subducting normally. This violates the assumptions of our 2D thermal modeling, but we include this transect for completeness given the extensive research that has been done on deep earthquakes beneath the South American margin.

Subduction Zone Thermal Modeling

Modeling of the thermal profiles of 23 subduction transects (Figures 4 and S1–S5; Table S1) allows us to investigate how transition zone slab seismicity varies with slab temperature and the stability of hydrous or carbonate minerals. Transects are oriented perpendicular to the trench and the convergence speed is obtained by projecting the convergence velocity onto the transect, which provides the most accurate 2D representation of 3D structure (Bengtson & van Keken, 2012). Slab geometry along each transect is determined using Slab2 (Hayes et al., 2018). For the plate motions of each of our transects, we use MORVEL2010 (DeMets et al., 2010) from which convergence direction and velocity are obtained. The one exception is for our transect through the Mariana's subduction zone, for which we account for back-arc spreading in order to more accurately reflect the rate of subduction of the downgoing Pacific plate (Syracuse et al., 2010). Specifically, we use the Pacific-Eurasia pole of rotation from Sella et al. (2002) together with the Philippine Sea—Eurasia and Philippine Sea—Northern Mariana's poles of rotation from Kato et al. (2003) to calculate the convergence rate at the trench along our transect. Subducting plate ages (Müller et al., 2008) are taken 50 km oceanward of the trench. For those transects where the slab geometry in Slab2 exceeds vertical orientation (folded-over slabs), we assume a near-vertical slab geometry for simplicity.

High-resolution finite element modeling (following Syracuse et al. [2010], van Keken et al. [2011], and Wei et al. [2017]) is used to predict the thermal structures of the subducting slab for each transect. The model uses a kinematic slab that follows the geometries described above, and a dynamic mantle wedge with temperature- and stress-dependent rheology based on dislocation creep in dry olivine. The models are based on the Boussinesq assumption but consider the effect of the exothermic 410-km discontinuity as well as choice of adiabat (either 0.3 or 0.5°C/km) by a posteriori addition. A concise description of the modeling with equations and solution method is in the supplement to Wei et al. (2017).

The effects of varying the adiabat and mantle potential temperature on our calculated slab trajectories are shown in Figures S8a and S8b (also see Figure S8 and Table S1 for subduction zone index). Supporting Information includes a description of the directories and the text files for temperatures at the top of the sediments, at the top of the oceanic crust, at the Moho (6.5 km below the top of the crust), and at the coldest place for each subduction transect for the four thermal models considered (mantle potential temperatures of 1,350°C or 1,422°C, and adiabats of 0.3 or 0.5°C/km). The actual data set (http://doi.org/10.5281/zenodo.4386270) is found in the ZENODO data repository.

Mapping Slab Seismicity to Slab Geotherms and Pressures

The earthquake data set is based on the U.S. Geological Survey-National Earthquake Information Center (NEIC) catalog (https://earthquake.usgs.gov/earthquakes/search) from 1990 to 2018 with magnitudes >4, depths >75 km, and depth errors <25 km. For each transect (Figures 4 and S1–S5), we project all events within 100 km of the transect onto each cross-section (Figure S6a). These projected locations, along with the original event latitude/longitudes, can be found in the “Events” directory of the supplementary data set.

To plot these events on a phase diagram, we must determine a location within the downgoing plate that these events occur. Given the large errors in the locations (and particularly in the depths) of these events, we cannot simply map the earthquake locations to the thermal models which all feature steep thermal gradients in the slab as a function of distance from the slab surface. Instead, we assume that each earthquake could have occurred at one of four locations within the slab. These four locations are meant to provide temperature bounds for the different lithologies in question. We map each earthquake (via closest distance) to (1) the top surface of the slab (from Slab2; Hayes et al., 2018) which we take as the top of the sedimentary layer, (2) the top of the basaltic oceanic crust 2 km deeper within the slab from the top surface, (3) the slab Moho (6.5 km below the top of the basaltic crust), and (4) a surface defined by the coldest point within the slab at a given distance from the trench (see Figure S6b for an example of earthquakes projected to the slab Moho). For each earthquake and for each slab layer, we use the projected location to assign the temperature and depth of that event. Depth is then converted to pressure using a linear interpolation of PREM (Dziewonski & Anderson, 1981). These projections can be found in the “PTeq” directory in the supplementary data set.

For comparing phase relations to these temperatures, we compare sedimentary lithologies to the temperatures at the top of the slab and the top of the oceanic crust, basaltic lithologies to the top of the oceanic crust to the top of the Moho, and mantle lithologies to the slab Moho and coldest slab temperatures. In doing so, the temperatures in each case are the highest and lowest temperatures possible for that composition.

Phase Relations in the Slab: Carbonate-Bearing Crust and Hydrated Mantle

Representative phase diagrams have been constructed from the literature to accompany the slab thermal modeling, specifically focusing on the crustal and uppermost mantle portions of the slab. These phase diagrams are generalized but allow for a comparison of slab P/T paths with mineralogy showing the potential for generating fluid/melt. Conditions of fluid/melt production can then be related to the presence or absence of deep-focus earthquakes at the appropriate P/T conditions. In each phase diagram, the critical feature is either the presence of the solidus in the hydrated or carbonate-bearing crust or the phase boundaries between hydrous mineral-bearing assemblages with significant water carrying capacity (e.g., ∼4 to >10 wt% H₂O) and much less hydrous (<2 wt% H₂O) or nominally anhydrous assemblages.

Results Water-Saturated Phase Diagrams, Slab Geotherms, and Deep Earthquakes

The oceanic lithosphere has three major settings near Earth's surface where it can acquire aqueous fluids: at the mid-ocean ridge, in the lithosphere as it cools and ages with its travel away from the ridge, and in the fore-arc before subduction (e.g., Faccenda, 2014). Sediments accumulated during slab residence on the ocean floor hold water and sit directly on the hydrothermally altered and hydrated basaltic crustal layer at the surface of the slab (Alt et al., 2013; Boschi et al., 2013; Früh-Green et al., 2004; Karson et al., 2015). Plate-bending at the outer rise produces lithospheric fractures that allow seawater to circulate deeply (see the summary in Faccenda, 2014) and serpentinize mantle lithosphere and low seismic velocities in the mantle portion of the slab (e.g., Ranero & Sallarès, 2004). Geophysical modeling (e.g., Faccenda, 2014), as well as a growing number of seismic studies (Bloch et al., 2018; Boneh et al., 2019; Cai et al., 2018; Faccenda, 2014; Fromm et al., 2006; Grevemeyer et al., 2018; Hatakeyama et al., 2017; Ranero et al., 2003; Wagner et al., 2020) indicate that fracture systems can result in hydration some 20–30 km into the slab. This is significant, not only for the total volume of water subducted with the slab, but also because plate interiors remain cold to greater depths, delaying the release of water from these parts of the slab.

The fate of water present in the oceanic lithospheric mantle as it subducts into the transition zone is assessed using a representative phase diagram for hydrous peridotite (Figure 5a). These are the important features of the phase diagram that dictate the amount of water that can be retained by the slab lithosphere into the deeper mantle: (1) the hydrous-phase dehydration curve (thick solid line in Figure 5a), (2) the antigorite breakdown curve (thick dashed line in Figure 5a), and (3) the field between total dehydration and antigorite breakdown where the hydrous, 10 Å phase (±Mg-sursassite) is stable. The maximum water storage capacity in antigorite-bearing slab peridotite is about 4–5 wt% (Iwamori, 2004; Komabayashi & Omori, 2006), whereas 10 Å phase-bearing (±Mg-sursassite) peridotitic assemblages have a lower water storage capacity of ∼1–2 wt% (Fumagalli & Poli, 2005; Iwamori, 2004; Pawley et al., 2011). In contrast, at higher pressures dense hydrous magnesium silicate phases (DHMS; Frost, 1999) in peridotitic bulk compositions can store at least 5 wt% water and potentially much more (>10 wt%), such that the amount of water transported to the deeper mantle is likely fixed by the water storage capacity of either antigorite or 10 Å phase-bearing peridotitic assemblages in cooler slabs (Figure 5a). In order for slab lithosphere to transport the maximum amount of water (e.g., 4–5 wt%) into the deeper mantle, hydrated lithosphere must remain below the maximum temperature at which antigorite-bearing peridotite can transform to phase A-bearing peridotite without loss of bulk water, sometimes referred to as a “choke-point” (e.g., Iwamori, 2004), which we show at about 6.5 GPa and 580°C on Figure 5a (white star). We recognize that phase relations will depend in detail on compositional parameters like water activity, oxygen fugacity, silica, and alumina contents, and that reactions are generally continuous and not univariant, such that the position of the dehydration boundaries will shift. The depiction in Figure 5a is drawn to be compatible with previous renditions of hydrous peridotite phase relations based on experiments and thermodynamic modeling (e.g., Fumagalli & Poli, 2005; Iwamori, 2004; Poli & Schmidt, 2002; Ulmer & Trommsdorff, 1995; van Keken et al., 2011).

Figure 5. Phase relations of hydrated mantle peridotite and basaltic crust relative to P/T paths with projected earthquakes. (a) Pressure versus temperature plot for hydrous peridotite. Earthquakes are projected onto temperature-pressure paths calculated for 23 modern subduction zone transects at the slab Moho for a mantle potential temperature of 1350ºC and an adiabatic gradient of 0.5°C/km. Red symbols show earthquakes in slabs that do not have seismicity below ∼350 km depth. Yellow symbols show Bolivia (transect 8) and Bonin (transect 19) which are included for completeness but whose modeled temperatures are unlikely to accurately reflect the true slab temperature (see discussion in text). Blue symbols show all other subduction zones with seismicity below 300 km depth that can accurately be modeled in two dimensions. Also shown are generalized phase relations for model hydrous peridotite. Phase boundaries at [less than]10 GPa are constructed after Poli and Schmidt (2002), Fumagalli and Poli (2005), Pawley et al. (2011) and Iwamori (2004) and references therein. At P>10 GPa the fields and dehydration curves reflect phase relations in the MgO-SiO2-H2O system based on the experimental study of Komabayashi and Omori (2006) and are similar in topology to those shown in the studies of Kawamoto (2004) and Iwamori (2004). The white star shows an estimate of the maximum temperature at which antigorite-bearing peridotite can transform to phase A-bearing peridotite without loss of bulk water. Phase labels designate generalized regions of stable hydrous magnesian silicates: Ant, antigorite; Chl, chlorite; 10Å, 10 angstrom phase; PhA, phase A; PhE, phase E; Bru, brucite; PhD, phase D; SHB, superhydrous phase B; Wad, wadsleyite; Ring, ringwoodite; Brg, bridgmanite. Estimated storage capacities for hydrous peridotite are shown as H2O wt% (after Iwamori, 2004; Komabayashi & Omori, 2006). See the supplement and Figures S1–S8 for further details about slab thermal models, earthquake depths, subduction zone transects and uncertainties, variations in phase relations and results of other thermal models. (b) Pressure versus temperature plot for hydrous basalt and deep earthquakes. Earthquake data for all transects are projected onto the top of the crust (purple and pink) versus the mantle at the Moho (dark green and light green). Purple and dark green indicate slabs with deep seismicity; pink and light green indicate slabs with no deep seismicity. These paths are superposed onto fluid-saturated phase relations for hydrous basalt (after Poli & Schmidt, 2002) which depicts the crustal mineralogy. Heavy black line shows the boundary between higher water-carrying modal mineralogy (medium blue), lower water carrying mineralogy (light green), anhydrous mineralogy (yellow), and melt (light pink). K-holl = potassium-rich hollandite. Red lines are from Kerrick and Connolly (2001) and are contours in wt% CO2 as carbonate remaining for carbonate-bearing basalt that is dehydrating during subduction. Note that many thermal paths cross the melt/dehydration front and that the cooler paths could retain significant carbonate.

Recognizing these caveats, we compare the P/T path for slab Moho temperatures traced by our projected earthquakes for a mantle potential temperature of 1,350°C and an adiabatic gradient of 0.5°C/km with the phase diagram in Figure 5a. Slabs that have deep earthquakes versus slabs that do not have deep earthquakes are separated by color (blue vs. red, respectively). The slab Moho P/T paths show that hydrated mantle lithosphere in most hotter slabs is predicted to dehydrate from chlorite or 10 Å phase-bearing assemblages into nominally anhydrous assemblages of olivine, orthopyroxene, and clinopyroxene at less than ∼200-km depth. However, some slabs without deep earthquakes pass below the total dehydration curve but above the choke point before the transition to a phase A-bearing assemblage. We note that in these few cases, earthquakes cease at <350 km. Conversely, hydrated peridotite in many colder slabs that have deep-focus earthquakes (blue) remain in the antigorite-bearing peridotite field before the transition into phase A-bearing peridotite and traverse below the choke point. In these cases, there is no expectation for water release upon antigorite breakdown or at higher pressures due to DHMS phase transitions (e.g., from phase A to phase E to phase D to superhydrous phase B) because the water content should always remain below the DHMS-bearing peridotite storage capacity. Deep water retention enabled by DHMS phases has been postulated in a number of previous studies that have proposed that for the coldest slab geotherms, DHMS can be hosts for deep water transport (Frost, 1999; Harte, 2010; Iwamori, 2004; Maurice et al., 2018; Ohtani, 2015; Ohtani et al., 2004; Omori et al., 2004; Poli & Schmidt, 2002; Thompson, 1992; van Keken et al., 2011).

The thermal profiles in Figure 5a are especially dependent on choices for the mantle potential temperature and the adiabatic gradient and we show a range of models in Figure S8. For example, choosing a mantle potential temperature of 1,422°C and an adiabat of 0.5°C/km (the warmest slab temperatures we considered) moves most of the blue curves above the choke point but below the dehydration reaction where 10 Å phase (±Mg-sursassite) is stable. In this case, much less water would be expected to be transported to deeper levels in peridotite (∼1 wt% water capacity). In contrast, a model with a 1,350°C potential temperature and 0.3°C/km has all of the slabs with deep earthquakes passing below the choke point. We also note that the temperatures in these models are calculated at the Moho, the hottest portion of slab lithospheric mantle, whereas hydration may extend tens of kilometers below the Moho, where temperatures can remain below the choke point even in the hottest slab models (see Figure S8). When considering all these caveats and uncertainties with respect to the phase relations and the thermal modeling, we suggest that the following observations remain robust: (1) there is a clear distinction between colder P/T paths that do have deep earthquakes and warmer ones that do not; (2) that the earthquakes in warmer slabs generally terminate in a region close to the dehydration trough; and (3) that along P/T paths that can plausibly allow for DHMS phases to form without prior antigorite dehydration (i.e., remain below the choke point), earthquakes occur at depths extending to the bottom of the mantle transition zone.

The two transects that we include for completeness where our 2D modeling may not fully reflect the complexities experienced along them are highlighted in yellow (Figure 5a). The first is the transect through Bolivia, where the likelihood of radical changes in slab geometry along the transect over the course of the past 20 Ma (roughly the time since the material now in the transition zone first entered the trench) suggests that our predicted temperature may not accurately reflect the slab's true thermal state or evolution (see Section 5.3 for discussion). The second is the Bonin transect where the recent work indicates a slab tear may be present near our profile (Faccenna et al., 2018; Zhang et al., 2019). In this case, we cannot rule out that the temperatures could be higher than modeled if flow through the slab tear serves to increase the temperatures above the downgoing slab. However, if a slab tear is not present, or if it does not affect the thermal structure along our profile, then we note that this slab, too, is roughly consistent with our other observations of cold slabs with deep seismicity. Other than these two profiles, included here for completeness, the observed pattern is robust for all subduction zones that can accurately be studied with 2D thermal models.

The fate of water and hydrous melt released by the oceanic crust as it subducts is shown by comparing temperatures at the top of the oceanic crust and at the Moho traced by our projected earthquakes compared to the generalized phase diagram for hydrous basalt (Figure 5b). We consider the crust for completeness here and to highlight how our earthquake projections show fluid and melt release versus mineralogy of the crust at depths less than 300 km (Hacker, 2008; Maruyama & Okamoto, 2007; Poli & Schmidt, 2002; van Keken et al., 2011) where they are thought to be directly related to intermediate-depth seismicity. As in Figure 5a, slabs that have deep earthquakes versus slabs that do not have deep earthquakes are separated by color (pink and light green: no deep earthquakes; purple and dark green: with deep earthquakes) to illustrate potential differences in dehydration behavior. At the top of the oceanic crust (pink/purple), slabs with and without deep earthquakes heat up at low pressure (2–3 GPa) and begin to dehydrate around 500°C as chlorite, talc, lawsonite, and amphibole schists that could contain up to ∼5 wt% water convert to zoisite eclogite that likely contains <1 wt% water (Schmidt & Poli, 1998). The remaining water present in zoisite eclogite means that at ∼700°C and 3 GPa some paths cross the water-saturated solidus for hydrous eclogite and produce hydrous melt. Temperatures at the base of the crust (Moho, light green/dark green) show that lawsonite eclogite can remain stable and transport crustal water to greater depths. However, the phase transformation to K-hollandite at ∼300-km depth marks the point at which all slab crustal sections are predicted to be dominantly dehydrated (e.g., Maruyama & Okamoto, 2007), and for this reason we do not consider the crust to be a likely source of hydrous fluids deeper in the mantle.

Melting Relations of Carbonate-Bearing Slab Lithologies, Slab Geotherms, and Deep Earthquakes

Carbon is found throughout the oceanic lithosphere (see review by Plank and Manning [2019]) including sediments, crust, and mantle lithosphere. The sedimentary carbon budget, comprised of seawater-precipitated carbonate and organic carbon, is readily released into the mantle wedge and effectively removed through the volcanic arc (Mason et al., 2017; Plank & Manning, 2019). Recent analyses of deep crustal gabbros (Li et al., 2019) reveal that they have a carbonaceous component whose isotopic composition suggests organic compounds are oxidized to carbonate in situ and occur deep within the oceanic crust. Carbonate can also be found in the peridotite portion of the oceanic lithosphere, especially in samples emplaced at fracture zones (Li et al., 2019). Abyssal peridotite at the ridge carries organic carbon (Ménez et al., 2018) that is likely augmented by the ambient flux of CO₂ to the base of the lithosphere as it ages in the ocean basins (e.g., Clerc et al., 2018). Thus both the oceanic crust and mantle portion of the lithosphere can contain significant carbon.

The fate of carbon as carbonate in the downgoing oceanic crust is assessed by plotting the modeled slab surface pressures and temperatures relative to the melting curves of carbonate-bearing sediments (Figure 6a) and carbonate-containing basaltic oceanic crust (Figure 6b). In the warmer slabs, carbonate in metasedimentary rocks of pelitic composition will be lost by generating carbonatitic melt below ∼8 GPa (Grassi & Schmidt, 2011) in the presence of residual water. This is because the carbonated melt solidus has a low-temperature melting cusp similar to the dehydration cusp in hydrated peridotites (Figure 5a). Note that carbonate-bearing metasedimentary crust in most warm slabs with residual water will intersect this cusp which also coincides with the cessation of seismicity in warm slabs. Conversely, sediments from colder slabs are not predicted to form melts from carbonate at any depth above the lower mantle—indeed, the solidus curves and thermal profiles diverge with increasing depth.

Figure 6. Phase relations of carbonate-bearing crust relative to P/T earthquake paths. (a) Pressure versus temperature plot for carbonate-bearing sediments at the top of the crust. The field of melt is in pink and the key subsolidus mineralogies and their changes with pressure are shown in shades of light blue, yellow, and green: phengite, white mica; cpx, clinopyroxene; Na carb, Na-rich carbonate. Earthquake projection paths for top of the sediments are for a potential temperature of 1,350°C and an adiabatic gradient of 0.5°C/km and are colored as in Figure 5a. These paths are compared to the melting curves for dry carbonate-bearing metasedimentary rock (pelite; brown solid) and metasedimentary rock with 1 wt% H2O from Grassi and Schmidt (2011; blue dashed). All metasedimentary rocks are at the top of the slab and will dehydrate similar to the paths shown in Figure 5b. The truncation of the red earthquake paths at the melting cusp at 9 GPa and 950°C shows that warm slabs with carbonate will dehydrate and produce melt whereas the cold slabs can retain their carbonate into the transition zone whether they have dehydrated or not (e.g., Gorman et al., 2006; Kerrick & Connolly, 2001) as indicated by the “less CO2” and “more CO2” labels along the hot slab versus cold slab trajectory, respectively. Note that for any path, dry sediment will retain all its carbonate. (b) Pressure versus temperature plot for carbonate-bearing basalt/gabbro at the top of the oceanic crust with paths (1,350°C/0.5°C/km) colored as in Figure 5a. The melting curves are from experimental studies: green dashed (Kiseeva et al., 2013), purple solid (Thomson, Walter, et al., 2016), and blue dashed (Zhang et al., 2020). Dol, dolomite; mag, magnesite; arag, aragonite; Na carb, Na-rich carbonate. P/T paths at this depth in cold slabs show that the basaltic crust will retain its carbonate in the form of magnesite ± aragonite until it intersects carbonate-bearing basalt solidus around 20 GPa for these model parameters. Note that the addition of any unreleased water from below in the lithosphere will drop the temperature of the solidus at a given pressure (downward black arrows).

For subduction of altered oceanic crust, the survivability of carbonate beyond the arc volcanic front is a critical function of crustal dehydration. As illustrated schematically by the red lines on Figure 5b, modeling shows that while some carbonate in the warmest slabs will likely be dissolved by aqueous fluids and removed during dehydration, cooler slabs and carbonate deeper in the crust will survive past the volcanic front (Alt et al., 2013; Gorman et al., 2006; Kerrick & Connolly, 2001). Beyond about 250 km, carbonate takes the form of magnesite ± aragonite until the depressed carbonate-bearing basalt solidus is intersected at about 15–20 GPa or greater, 1,150°C–1,250°C (Kiseeva et al., 2013; Thomson, Walter, et al., 2016; Zhang et al., 2020; Figure 6b). For the P/T paths with deep earthquakes shown, with a potential temperature of 1,350°C and an adiabatic gradient of 0.5°C/km, the solidus of Thomson, Walter, et al. (2016) is intersected at ∼20 GPa to produce carbonatitic melt and majorite, Ti-rich CaSiO₃ perovskite, and stishovite. We note that the position of the melting curve depends on the original composition of the eclogite (see Figure 6.2 and Table 6.1 of Yaxley et al. [2019]), and that variations in the carbonate-bearing basalt solidus among recent experimental studies shown in Figure 6b primarily reflect variations in the bulk Ca/Mg ratio. Any water left in the crust or introduced from below by the dehydration of the slab mantle at transition zone depths could also have the effect of lowering the carbonate-bearing basalt solidus significantly (Dasgupta et al., 2007). Removal of carbonate induced by such infiltration from below could be particularly effective in warm slabs (e.g., Gorman et al., 2006).

The net effect in subducting carbonate-bearing crustal lithologies is that carbonate not removed by the more labile fluids and melting beneath the volcanic front will be retained in the subducting crust until around 500 km or deeper where the carbonate-bearing eclogite solidus is intersected and a melt is produced (Figure 6b). In the case of cold subduction zones, the pressures and temperatures at which deep earthquakes occur are similar to the conditions required for melting of carbonate-bearing oceanic crust. For hotter subduction zones, no earthquakes are seen near carbonate-bearing basaltic solidus conditions.

Discussion Correlation Between Deep Slab Hydration, Deep Earthquakes, and Deep Diamonds

Our slab thermal models show a strikingly strong correlation between slabs with deep earthquakes and slabs that could retain water to transition zone depths in their hydrated mantle lithosphere. The choke point near 6.5 GPa and 580°C produced by the higher temperature stability of antigorite and the incoming stability of phase A (white star, Figure 5a) provides a plausible explanation for the separation between slabs that have seismicity below ∼350 km depth and those that do not. The modeling shows that within cold subducting slabs, hydrous fluid retention into the mantle transition zone is likely if not expected. However, our modeling alone provides no mechanism for fluid release, a problem recognized by Green et al. (2010). Missing from their analysis and indeed every summary paper on deep focus seismicity before or after (e.g., Frohlich, 1989; Green & Houston, 1995; Houston, 2015; Zhan, 2020) is an understanding that fluids can be released from the slab at transition zone depths (e.g., Richard et al., 2007) and that these fluids can move through surrounding mantle (e.g., Sun et al., 2020). The necessary evidence for fluids and their movement is provided by sublithospheric diamonds as the discovery of a ringwoodite inclusion in diamond with ∼1.5 wt% water (Pearson et al., 2014) demonstrates. Consistent with this evidence are conclusions from many studies of sublithospheric diamonds (Figures 2a–2c; Table 1) where inclusion mineralogy constrains the depths to which fluids can be transported, and inclusion geochemistry ties the fluids to deeply subducted slabs.

The DHMS phases that carry the water to the transition zone are stable throughout the transition zone at the temperatures expected for cold slabs in our models, which raises the question of how a dehydration mechanism could cause earthquakes if no dehydration is predicted. Omori et al. (2004) considered that phase transitions among the DHMS phases can release water throughout the deep upper mantle and transition zone, possibly accounting for deep seismicity. We find this mechanism unlikely. The storage capacity in DHMS-bearing peridotite is higher than the expected water contents in subducted peridotite lithosphere at depths beyond antigorite breakdown such that no fluid release would be expected during phase transitions (Iwamori, 2004; Komabayashi & Omori, 2006). We suggest the answer lies in the dynamics of slab subduction and the potential for slabs to stall and heat up in the transition zone (Bina, 1997). There is significant evidence from seismic tomography (Li et al., 2008) and dynamical modeling (Billen, 2008, 2020) that slabs tend to bend, thicken, and stagnate in the transition zone as they approach the higher viscosity and slower convection of the lower mantle (Mao & Zhong, 2018). Variations in strain rate caused by slab deformation at these depths may play an important contributing role to the detailed spatial distribution of deep seismicity, independent of triggering mechanism (Billen, 2020). Slower slab movement means that warming is inevitable as the ∼1,000°C mantle of the cool slab interior at 15 GPa (Figure 7) sits in an ambient mantle of ∼1,500°C.

Figure 7. Dehydration model for earthquake generation at 500–700 km depths from the crust to the mantle in the lithosphere. The pressure-temperature diagram depicting the phase relations for hydrous peridotite is generalized from Figure 5a and the carbonate-bearing basalt solidi (red solid and dashed) are from Figure 6b. H2O contents for the slab reflect maximum estimates expected to make it past antigorite breakdown in subducted peridotitic lithosphere. Other water contents are generalized values for maximum storage capacity in peridotite lithologies at higher temperatures (Iwamori, 2004; Komabayashi & Omori, 2006). The histogram of subduction zone earthquakes from Figure 1b is shown in 50 km bins with light brown bars. P/T paths with projected earthquakes are for the top of the oceanic crust (purple and pink) or for the peridotitic slab Moho (dark green and light green). Purple and dark green = slabs with deep seismicity, pink and light green = slabs with no deep seismicity. The temperatures at the top of the oceanic crust indicate the warmest temperatures in this model for oceanic crust that could cross the carbonate-bearing basalt solidi. The Moho temperatures indicate the warmest temperatures in this model for peridotite compositions that may include DHMS phases. In warm slabs both the slab crust and mantle extensively dehydrate by depths of 250 km. Note the alignment of the end of the light pink and light green seismicity and the right edge of the sharply decreasing intermediate-depth earthquake frequency on the left of the diagram illustrating this point. In cool slabs the dense hydrous magnesium silicate (DHMS) phases will be stable and transport water deep into the transition zone. Once the cool slabs stall and warm up in the transition zone (gray arrows), breakdown of the DHMS phases produces fluids that initiate the peak in earthquakes at 550–650 km and provide the fluids from which aqueous sublithospheric diamonds and their inclusions precipitate. Note that the cessation of earthquakes at ∼700 km depth corresponds closely to the stability region where super hydrous phase B coexists with bridgmanite and can store more water (∼2%) than in wadsleyite and ringwoodite dominated lithologies (∼1%). The top of the crust is warmer, fully dehydrated (e.g., Figure 5b) and intersects the carbonate-bearing basalt solidi at about the same depth as the proposed dehydration of the peridotite. Because little or no warming of the crust is needed, production of carbonatitic melt from the altered oceanic crust would occur before the cooler slab mantle has time to warm up.

The amount of heating needed to reach the DHMS dehydration boundary at ∼1,200°C will depend on the specific slab, mantle potential temperature, and adiabatic gradient, and is likely around (∼100°C–300°C; Figures 7 and S8). The time it takes for the Moho temperature to reach 1,200°C by diffusive heating of the slab from ∼1,000°C after stalling at 600-km depth is expected to be between 10 and 30 Myr. In a general sense, the time scale for slab Moho warming through 100°C interval is roughly the same as the time taken for the slab to reach 600-km depth (e.g., for the Tohoku location, it is 11 Myr). As the stalled slab warms, temperatures will cross the DHMS dehydration boundaries. Phase E-, super-hydrous phase B-, and phase D-bearing peridotite assemblages with high water storage capacities of 5 wt% or greater will transform to much less hydrous assemblages containing wadsleyite and ringwoodite in modal proportions that would produce bulk peridotite with ∼1 wt% water storage capacity (Iwamori, 2004; Komabayashi & Omori, 2006). This residual amount of water is considerably less than the ∼4–5 wt% that can be transferred to DHMS phases upon antigorite breakdown in hydrous peridotite, and we hypothesize that it produces a diamond forming region in the phase diagram at temperatures above ∼1,100°C–1,250°C and between ∼13 and 23 GPa (blue boundary, Figure 7).

Intriguingly, the cessation of seismicity in hot slabs and the cessation of seismicity in cold slabs bear a close relationship to the mineralogy of the slab mantle. As noted above, the cessation of seismicity in hot slabs and the drop-off in the intermediate-depth earthquake histogram (Figure 7) coincides roughly with the hydrous peridotite dehydration curve at around 5–7 GPa (Figures 5a and 7). Similarly, the abrupt cessation of seismicity in cold slabs at ∼700 km shown in Figures 1, 5a, and 7 corresponds generally with the appearance of a superhydrous phase B + bridgmanite assemblage with higher water carrying capacity than a ringwoodite-bearing assemblage. It is unclear whether dehydration is simply suppressed at this point by reincorporation into this assemblage or whether the slab has simply become fully dehydrated. We also note that the seismic velocity anomalies identified in Schmandt et al. (2014) and interpreted to represent hydrous melt occur at ∼700 km and may be typical of the phase-change-induced fluid release of diamond-forming fluids that we ascribe to subducting slabs.

Correlation Between Carbonate Melting, Deep Earthquakes, and Deep Diamonds

Carbonatitic sub-lithospheric diamonds (Table 1) provide direct evidence of the presence of slab-derived carbonatitic melts at the depths of deep earthquakes (Bulanova et al., 2010; Thomson et al., 2014; Thomson, Walter, et al., 2016; Walter et al., 2008). It is important to note that we are not referring to carbonatitic melts that could have been produced by the abundant carbonate in the sedimentary layer of the crust. These carbonates would likely have been removed shallower than 250 km during slab dehydration in the forearc (Plank & Manning, 2019) and by deeper melting such as that at the cusp in the water-bearing carbonate-bearing sediment solidus (Figure 6a; see Grassi & Schmidt, 2011). Additionally, the carbon isotopic composition of carbonatitic sublithospheric diamonds, which have an average δ¹³C of −25 and range from −28‰ to −3‰ (Bulanova et al., 2010; Thomson et al., 2014), is much too light to have been derived from carbonate in sediments with a heavy average δ¹³C of 0 and a range from −3‰ to +5‰ (Cartigny, 2005). Rather, in this examination of the role of carbonate, we are only considering the role of cooler slabs and carbonate in the deeper, altered igneous oceanic crust (Li et al., 2019), especially since it is only the cooler slabs that produce deep seismicity regardless of mechanism (Figures 5a, 6a, and 6b). This association with oceanic crust is cemented further by oxygen isotope measurements on majorite inclusions in sublithospheric diamonds (Burnham et al., 2015; Ickert et al., 2013; Regier et al., 2020). Just as with water in altered peridotite, carbonate in the deeper igneous oceanic crust of hot slabs will be removed shallower than 300 km depths (Figure 5b; e.g., Kerrick & Connolly, 2001), especially if water is infiltrating from cooler, deeper crustal levels (e.g., Gorman et al., 2006; Martin & Hermann, 2018).

The melting curve of carbonate-bearing oceanic crust (Kiseeva et al., 2013; Thomson, Walter, et al., 2016; Zhang et al., 2020) is crossed by the P/T path of cold slab crust in the transition zone (Figures 6b and 7). The melting conditions for carbonate in oceanic crust are similar to the dehydration phase boundary for hydrated peridotite (Figure 7), but since the oceanic crust is always hotter than the oceanic Moho in a subducting plate, no additional heating scenario (i.e., stalling of the slab in the transition zone) is required to produce carbonatitic melt in contrast to the case for dehydration of peridotite. The presence of any residual water or water derived from dehydrating DHMS is likely to lower carbonate-bearing crust melting temperatures and facilitate the process. The overlapping depth ranges of dehydration in the slab mantle lithosphere and carbonate melting of the slab crust suggest that it is not possible, using available phase equilibria only, to differentiate their relative contributions toward the genesis deep earthquakes. These uncertainties, along with the 3D complexity of slab thermal evolution, make a simple or precise comparison between slab geotherms, carbonate melting curves, and slab seismicity difficult to interpret without additional constraints.

However, examination of sublithospheric diamond petrology in the context of our slab P/T paths and slab-related seismicity provides some clue as to how carbonatitic melts and aqueous fluids might be related at depth. On the scale of individual diamonds, carbonatitic diamonds are dominated by oceanic crustal mineral assemblages (e.g., Bulanova et al., 2010; Thomson et al., 2014; Walter et al., 2008, 2011) whereas aqueous sublithospheric diamonds are dominated by mantle mineral assemblages (e.g., Smith et al., 2016, 2018, 2021). Neither paragenesis fully excludes the other (Table 1) which gives the impression that the fluids generated have not been fully confined to their mantle or crustal hosts and may have migrated (i.e., Regier et al., 2020). On the diamond locality or slab domain scale, collections of sublithospheric diamonds are usually not mixtures of similar proportions of carbonatitic and aqueous types found together. For example, diamondiferous kimberlites from Juina, Brazil are known for erupting sublithospheric diamonds whose parental melts are derived chiefly from the carbonate-bearing basaltic crust whereas the sublithospheric diamond suite of kimberlites from Letseng, Premier, and Karowe, in southern Africa, are sourced from fluids derived from hydrated mantle peridotite. These petrological features of the sublithospheric diamond suite are readily explained if, in a cooler slab that is stagnating and warming, the carbonate-bearing oceanic crust starts melting first because its solidus is closer to melting than the slab mantle is to dehydrating. Dehydration of the hydrated slab peridotite might lag behind by as much as 10 Myr because of the time taken to warm up by conduction. This time lag would be sufficient to place carbonatitic and aqueous sublithospheric diamonds in different groups in relation to their age and physical location such that they could be sampled separately by the discrete and geologically pinpointed nature of kimberlitic volcanism (e.g., Giuliani & Pearson, 2019).

Fluid-Triggered Mechanisms for Deep Seismicity

A full description of a fluid-triggered mechanism for deep earthquakes remains for future study. However, the evidence we present for the transport of both water and carbonate to transition zone depths and the release of hydrous fluids and/or carbonate melting at these depths is in many ways easy to integrate with the existing dominant theories for both intermediate-depth and deep earthquakes, and therefore does not necessarily require the formulation of an entirely new mechanism. For deep earthquakes, the most common mechanism cited relies on anti-cracking caused by volumetric reductions associated with the phase transformation of olivine, metastable beyond its normal pressure range, to more densely packed but weak nanocrystalline structures common in the transition zone. These are the key elements necessary: (1) the production of anti-cracks (small voids) caused by the reduction in volume of small pockets of recently transformed material that result in a transfer of accumulated stress onto the surrounding untransformed lattice; and (2) the formation of nanoparticulate reaction products that are weak and allow for the evolution of a fault surface as anti-crack regions begin to link up (e.g., Green & Houston, 1995). In many ways, the mineral involved in the phase transformation is secondary. The dehydration of DHMS phases may also produce a net volume reduction and weak nanocrystalline reaction products. Melting along the depression in the solidus of carbonate-bearing oceanic crust may similarly produce a volumetric decrease, as is required by the negative Clapyron slope given a positive entropy of melting (Figure 6b). In this case, the weak melt would substitute for weak nanocrystalline reaction products. Either of these mechanisms (the dehydration of DHMS phases or the formation of carbonate melt) may be capable of producing small weak zones that coalesce into a more substantial fault plane. Eventually, movement on this nascent fault can generate enough heat to both increase the phase transformation process and ultimately propagate into adjacent regions through a thermal runaway shear heating mechanism. Where Green and Houston (1995) and others describe the transformation of metastable olivine producing “anti-cracks” and nanocrystalline weak reaction products to form deep earthquakes, Ferrand (2019) describes the dehydration of antigorite resulting in small distributed volume decreases and weak reaction products to form intermediate-depth earthquakes. We are simply proposing that DHMS phases and/or carbonate can play the same role for deep earthquakes as antigorite does for intermediate-depth earthquakes under this model.

As described by Zhan (2017), the size of faults associated with deep earthquakes is much larger than any proposed metastable olivine wedge or, for that matter, any plausible region of slab hydration. While slabs are now believed to be hydrated to depths of tens of kilometers (Bloch et al., 2018; Boneh et al., 2019; Cai et al., 2018; Faccenda, 2014; Fromm et al., 2006; Grevemeyer et al., 2018; Ranero et al., 2003; Wagner et al., 2020), the very largest seismic events likely require a secondary mechanism such as shear heating to propagate into surrounding “stable” areas. However, smaller seismicity should be able to remain localized within the region of phase transformation, and so may provide information on the size of that region. For the metastable olivine wedge theory, warm slabs are predicted to be able to generate deep seismicity only within a very narrow (3–6 km) span of the coldest part of the subducting slab. We would argue that warm slabs do not have seismicity and that hydrous phases are preserved in the slab mantle from the slab Moho to the depth of the coldest portion of the slab. This could explain the presence of double seismic zones illuminated in the transition zone (i.e., Wiens et al., 1993), and the overall breadth of deep seismicity zones, though this may also be due to significant errors in location, particularly as a function of depth.

The Bolivian slab appears to be an important exception to the cold slab–fluid transport relationship in that it is currently a slab with apparently moderate thermal parameters while having deep seismicity. We speculate that since subduction was faster and delivered older oceanic lithosphere at the trench in the Miocene (Sdrolias & Müller, 2006; Somoza & Ghidella, 2012), the deep lithosphere might be associated with this older, colder subduction that was able to avoid dehydration and was therefore capable of supporting seismicity. This will be tested in the future study. The Bolivian slab example serves to remind us that slabs of irregular geometry and varying subduction parameters (age, convergence rate, and dip angle), which we purposely did not consider, have seismicity whose genesis cannot be assessed without much more complex four-dimensional thermal modeling that can account for 3D geometries that change over time along with varying convergence rates and subducted plate ages.

Implications for Deep Focus Earthquakes and Sublithospheric Diamond Growth

Our observations relating slab thermal profiles, dehydration and decarbonation melting phase equilibria, deep earthquake depth distribution, and sublithospheric diamond petrology suggest that deep focus earthquake hypocenters could mark, in a general sense, the release and migration of fluids that lead directly to sublithospheric diamond growth. We can only speculate as to the exact mechanism by which release of hydrated and/or carbonated fluids at transition zone depths induces seismicity (e.g., dehydration embrittlement), but we suggest there is a link between sublithospheric diamonds that are now forming in the mantle (Figure 8) and where deep-focus earthquakes occur. If deep subduction with slab stalling is a product of the post-Archean Earth (e.g., Klein et al., 2017), then the potentially younger ages of sublithospheric diamonds implied by the current sparse data (Bulanova et al., 2010; Harte & Richardson, 2012; Hutchison et al., 2012) compared with the more ancient lithospheric diamonds (e.g., Howell et al., 2020) is consistent with such a growth environment.

Figure 8. Blocky and irregular type IIa diamonds thought to be crystallized from the type of fluids that lead to earthquakes in the mantle transition zone. The typical textures that suggest a mode of growth that requires the coalescence of a significant amount of diamond-forming fluid in a mantle fracture network. The names and sizes of diamonds as follows: (a) the “Constellation” (813 ct); (b) an unnamed stone (124 ct); (c) the “Queen of Kalahari” (342 ct); (d) an unnamed stone (472 ct). All diamonds are from the Karowe Mine, Botswana. (Photos courtesy of John Armstrong, Lucara Mining Corp.; used with permission).

The sublithospheric diamond suite contains the largest diamonds made within Earth (Figure 8). They have unusual, non-octahedral shapes compared to lithospheric diamonds and/or display textures or dense dislocation networks associated with intense deformation (Bulanova et al., 2010; Shirey et al., 2013; Smith et al., 2016, 2018). These very large diamonds and their anhedral, often slab-like morphologies (Figures 8a–8d), have posed an enigma in the diamond community since their discovery. While it seems too speculative to assert that the dilatancy in mantle host rock caused by earthquake rupture could make space for anomalously large diamonds to grow, it is clear that fluids must be locally abundant enough to enable the growth of diamonds to sizes exceeding 500 carats and in rare cases >1,000 carats. Large size, low inclusion density, and slab-like morphologies diamonds demonstrate fluid coalescence in quasi-planar voids in the mantle large enough to accommodate their sizes.

Conclusions

Our results show clearly that both hydrous fluids and carbonate melts are not only possible but expected to form from cold subducting slabs in the mantle transition zone. Given the cumulative evidence from diamonds and the thermal modeling presented here, we suggest that mechanisms involving these fluids are likely involved in the genesis of deep earthquakes in part because fluids are so clearly related to intermediate-depth seismicity (Shiina et al., 2017; Wiens, 2001; Zhan, 2017, 2020). The mechanisms need not be the same for hydrous fluids and carbonatitic melts. In the case of the former, it is worth considering whether the dehydration of DHMS phases could produce similar nanocrystalline particles and volumetric changes to those that have been postulated for the metastable olivine wedge. For the latter, localized melting could produce weakening and/or volumetric changes that could trigger seismicity and/or promote localized shear instabilities. Significant work needs to be done to test these possibilities, both experimentally and seismically. Given the diversity of observed characteristics in deep seismicity (see Zhan, 2020), it is possible that both mechanisms play a role at different stages of a subducting slab's evolution.

Acknowledgments

Helen Janiszewski reminded us of the state of the deep earthquake problem at an Earth and Planets Laboratory geophysics reading group which led to this study. Evan Smith, Karen Smit, David James, and Alan Linde read and critiqued an early draft of this manuscript. The authors thank them and German Prieto and Andy Thomson for discussions. The authors thank Doug Wiens, Oliver Tschauner, Ben Harte, and an anonymous reviewer for their comments that helped us improve the presentation. The authors thank the Carnegie Institution for Science and the University of Alberta for continuing support. Diamond research is supported by the Deep Carbon Observatory and NSF grant EAR-2025779 to Steven B. Shirey and Michael J. Walter.

Conflict of Interest

The authors declare no competing interests or affiliations.

Data Availability Statement

No new analyses on minerals were carried out for this study and all data on natural diamond samples and experiments used in the study come from the literature as cited in this study. The earthquake data set is based on the U.S. Geological Survey-National Earthquake Information Center (NEIC) catalog (https://earthquake.usgs.gov/earthquakes/search). All data produced by modeling and synthesis in this study are available in the main manuscript, the supporting information, or uploaded data files. The uploaded data files for each slab subduction transect that support the findings of this study are available in the supplementary data set for this study (https://doi.org/10.5281/zenodo.4386270) in the ZENODO data repository.