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

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.

. 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  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. 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;) 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, 2017Zhan, , 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 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., 2014Palot et al., , 2016Pearson et al., 2014;Shirey et al., 2019;Smith et al., 2016Smith et al., , 2018Thomson et al., 2014;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 a l., 2016Smith et a l., , 2017Smith et a l., , 2018Smith et a l., , 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).

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., 2016Smith et al., , 2018Stachel, et al., 2000a;Stachel, et al., 2000b;Stachel et al., 2005;Thomson et al., 2014;Walter et al., 2008Walter et al., , 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 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;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 SHIREY ET AL.
10.1029/2020AV000304 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) Table 1 for diamond type definitions) diamond (left) with tomographic image of ringwoodite and breyite (interpreted to be retrograded Caperovskite; 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 CH 4 or H 2 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 CH 4 and H 2 that show equilibration of metallic liquid with a hydrous environment. Figures and interpretation from Smith et al. (2016).  Tschauner et al. (2018) a Type 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). b Number of inclusion-bearing diamonds as published is given. Number of sublithospheric diamonds examined to find inclusions for study is in parentheses. c Melt 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). d CLIPPIR 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. (2016Smith et al. ( , 2017. Depths estimated from at-pressure analog silicate mineral assemblages extrapolated to entrapment depths. e IIb diamonds are inferred to crystallize from aqueous fluids based on their boron content (associated with seawater and serpentinization) and the hydrogenrich jackets (CH 4 ± H 2 ) that surround the silicate minerals (Smith et al., 2018). f Estimated 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). g Pressure estimated from thermdynamic modeling of phase diagram for magnesioferrite and ferropericlase (Palot et al., 2016). h Pressure estimate derived from stability field for ringwoodite. Ringwoodite contained 1.5 wt% water . i Pressure estimate derived from stability field for AlSiO 3 (OH) (also known as "Phase Egg") under mantle conditions. Fluid estimate based on high porosity at room pressure (Wirth et al., 2007). j Pressure estimates made using residual pressures, equation of state of Ice VII and H 2 O-based fluid intersection with mantle adiabat (Tschauner et al., 2018).
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;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 (

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., 2016Smith et al., , 2018Thomson 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  indicates reactions between the mantle and diamond-forming fluids, suggesting that mantle rocks are hosting mobile fluids in veins.

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 3 phases (breyite plus perovskite) and are thought to crystallize when carbonatitic melts derived from the basaltic portion of the slab react with the mantle 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 (δ 18 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 (δ 13 C as low as −28.7‰) far below that of the average ambient mantle composition in δ 13 C of around −5‰ (Table 1; Bulanova et al., 2010;Palot et al., 2014Palot et al., , 2017Shirey et al., 2019;Walter et al., 2008Walter et al., , 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., , 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 (δ 13 C as low as −26.9‰ and −20.8‰, respectively) that ranges to far below the average ambient mantle (Table 1; Smith et al., 2016Smith et al., , 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 (δ 56 Fe = + 0.79 to +0.90) and significantly less radiogenic ( 187 Os/ 188 Os = 0.1109 to 0.1115) than the "normal" oceanic mantle ranges (i.e. δ 56 Fe = −0.2 to +0.2; 187 Os/ 188 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).
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  Table S1 for naming and precise locations.
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 2 O) and much less hydrous (<2 wt% H 2 O) or nominally anhydrous assemblages.

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).
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 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 <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-SiO 2 -H 2 O 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 H 2 O 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% CO 2 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. 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 2 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.
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;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  is intersected at ∼20 GPa to produce carbonatitic melt and majorite, Ti-rich CaSiO 3 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. 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% H 2 O 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 CO 2 " and "more CO 2 " 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 , 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).

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  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(Billen, , 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.
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. SHIREY ET AL.  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.

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;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 δ 13 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 δ 13 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;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., 2008Walter et al., , 2011 whereas aqueous sublithospheric diamonds are dominated by mantle mineral assemblages (e.g., Smith et al., 2016Smith et al., , 2018Smith et al., , 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.
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., 2016Smith et al., , 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 SHIREY ET AL.  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, 2017Zhan, , 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.