Heat Transfer and Production in Cratonic Continental Crust: U‐Pb Thermochronology of Xenoliths From the Siberian Craton

Coupled U‐Pb and trace‐element analyses of accessory phases in crustal xenoliths from the Late Devonian Udachnaya kimberlite (Siberian craton, Russia) are used to constrain Moho temperature and crustal heat production at the time of kimberlite eruption. Rutile and apatite in lower‐crustal garnet granulites record U‐Pb dates that extend from 1.8 Ga to 360 Ma (timing of kimberlite eruption). This contrasts with upper‐crustal tonalites and amphibolites that contain solely Paleoproterozoic apatite. Depth profiling of rutile from the lower‐crustal xenoliths show that U‐Pb dates increase gradually from rim to core over μm‐scale distances, with slower‐diffusing elements (e.g., Al) increasing in concentration across similar length‐scales. The U‐Pb and trace element gradients in rutile are incompatible with partial Pb loss during slow cooling, but are consistent with neocrystallization and re‐heating of the lower crust for <1 Myr prior to eruption. Because Paleoproterozoic rutile and apatite dates are preserved, we infer that long‐term ambient lower‐crustal temperatures before this thermal perturbation were cooler than the Pb closure temperature of rutile and probably apatite (<400°C). The lower‐crustal temperature bounds from these data are consistent with pressure‐temperature arrays of Udachnaya peridotite xenoliths that suggest relatively cool geothermal gradients, signifying that the mantle xenoliths accurately capture the thermal state of the lithosphere prior to eruption. Combined, the xenolith data imply low crustal heat production for the Siberian craton (∼0.3 μW/m3). Nevertheless, such values produce surface heat flow values of 20–40 mW/m2, higher than measured around Udachnaya (average 19 mW/m2), suggesting that the surface heat flow measurements are inaccurate.

A tenet in developing thermal models for the lower crust based on U-Pb thermochronology is the notion that Pb undergoes thermally mediated volume diffusion at ambient lower-crustal temperatures. Indeed, xenolith transport in kimberlites and alkali basalts is too rapid (on the order of hours to days; Spera, 1984) to significantly disturb the U-Pb system in most accessory phases (Schmitz & Bowring, 2003). Xenoliths, however, may be heated during episodes of advective magmatic heating or fluid flow prior to their entrainment, and such heating could induce partial Pb loss in thermochronometers (e.g., Apen et al., 2020;Hacker et al., 2000;Jollands et al., 2018;Smit et al., 2016). Any inferences on lithospheric thermal structure based on xenoliths and their mineral constituents must therefore delineate primary signatures of crustal relaxation (on the order of 100-1,000 Myrs) from shorter-term heating by magmatic intrusion or metasomatism (<1 Myr). Spatially resolved analysis of multiple trace elements with different diffusivities can discriminate signatures of long-term cooling from shorter-term heating events in the lower crust: diffusive length scales should vary according to diffusivities, with faster diffusing elements showing more pronounced length scales of diffusion than slower diffusing ones (Apen et al., 2020;Garber et al., 2020;Holder et al., 2019;Stearns et al., 2016).
To test the fidelity of U-Pb thermochronology as a tool for constraining the long-term thermal history of the deep crust, we have undertaken laser ablation split stream (LASS) U-Pb and trace-element analyses of accessory phases in crustal xenoliths from the Udachnaya kimberlite located within the central Siberian craton (Figure 1). The xenolith suite covers a broad range of lithologies, including granitoids, amphibolites, and garnet-bearing granulites. This lithological diversity is interpreted as sampling different levels of the crustal column at Udachnaya. The new laser ablation data, combined with new thermobarometry and pseudosection modeling, are used to probe the thermal state of the crust prior to eruption and develop models of heat transfer and production within the Siberian craton.

Geologic Background and Previous Xenolith Work
Most of the Siberian craton is overlain by thick Phanerozoic sedimentary and volcanic cover , but information about crust formation and structure is gleaned from limited basement outcrops and kimberlite-hosted xenoliths (Figure 1). Based on geophysical data and outcrops in the Anabar shield, the craton is divided into five major terranes: Tungus, Magan, Daldyn-Markha, Olenek, and Aldan (see summaries in Rosen, 2002). Basement gneisses and granulites exposed in the Anabar shield reveal that the central craton is composed of Archean crust (ca. 3.6-2.7 Ga) reworked throughout the Archean and Proterozoic as well as juvenile Proterozoic crust (2.1-1.8 Ga) (Figure 1; Moyen et al., 2017;Paquette et al., 2017;Rosen, 1989;Smelov et al., 2012). Collision of the Tungus, Magan, Daldyn-Markha, and Olenek terranes to form the Siberian craton occurred 2.1-1.8 Ga based on the geochronology and geochemistry of Paleoproterozoic granitoids and metamorphic rocks of the Anabar shield (Rosen, 2002;Rosen et al., 2006). Re-Os dating of peridotite xenoliths yields similar age intervals for the formation and reworking of the underlying cratonic lithospheric mantle (Ionov et al., 2020;Moyen et al., 2017).
The crust in the vicinity of Udachnaya is 44-46 km thick Mackey et al., 1998). Surface heat flow data across the Siberian craton are sparse, but the six measurements within ∼100 km of Udachnaya yield an average value of 19 ± 3 mW/m 2 (1SD) (Fuchs et al., 2021, and references therein). As it stands, these data are among the lowest surface heat flow values observed for any craton (cf., Jaupart et al., 2016). If accurate, they require either exceedingly low heat production in the crust, a low mantle heat flux, or both (Jaupart et al., 2016;R. L. Rudnick et al., 1998). Alternatively, the surface heat flow measurements may be inaccurate, possibly compromised by the thick sedimentary/volcanic cover overlying the crystalline basement around Udachnaya , regional brine aquifers (Kitayama et al., 2017, and references therein), and/or Ice-Age glaciation (Birch, 1948;Jaupart & Mareschal, 2007). Low crustal heat production for the central Siberian craton is permissible based on regional magnetic Curie depths (580°C isotherm) estimated to be >35-km-deep (Gard & Hasterok, 2021).
All of the granulites exhibit granoblastic textures, and their mineral assemblages are well-preserved ( Figure 2). The xenoliths are generally fresh, by comparison with many of those from previous work, though some of the granulites show textural evidence for minor alteration, like pockets of very-fine grained material interpreted to be former melt (e.g., mafic garnet granulite 79-14) and biotite clusters interpreted to have formed by re-hydration (e.g., mafic garnet granulite 01-95) ( Figure 2). Rutile in garnet granulite xenoliths is typically 50-200 μm long and euhedral and occurs in the matrix-usually adjacent to garnet. Rutile is also found as inclusions within garnet and clinopyroxene (e.g., Figure 2c), indicating crystallization during peak granulite-facies metamorphism. Rutile appears free of ilmenite rims and exsolutions and does not contain obvious titanite or ilmenite overgrowths. Apatite grains are euhedral stubby prisms 50-200 μm long and 20-150 μm wide and are found both in thin exsolutions in their cores (Figure 2c).
We present results from thermobarometry and coupled U-Pb and trace-element analyses of accessory phases for the seven xenoliths. Zircon and apatite were analyzed for the tonalites and amphibolite (these rocks contain no other phases datable by U-Pb). Apatite and rutile were analyzed in all four of the granulites, whereas zircon was only found in two of them (01-95 and 36-14).

Electron Probe Microanalyses
Quantitative elemental analyses of major mineral phases were done using a Cameca SX-100 electron probe micro analyzer housed at the University of California, Santa Barbara. Measurements were made using the following beam conditions: an accelerating voltage of 20 kV, a beam current of 200 nA and a defocused beam of 2-5 μm diameter. A series of natural and synthetic standards were analyzed for calibration purposes.
We compared results from net-transfer and exchange reactions to pseudosection models developed using Perple_X (Connolly, 2009;Connolly & Kerrick, 1987). Models used the internally consistent thermodynamic data set of Holland and Powell (2011) and the metabasite set of solution models recommended by Green et al. (2016) and Palin et al. (2016). Whole-rock data used for pseudosection modeling are from Moyen et al. (2017); the whole-rock compositions of the newly analyzed xenoliths were determined using X-ray fluorescence at the Geoanalytical Lab at Washington State University. Details are reported in Supporting Information S1.

Laser Ablation Split-Stream Analyses
The xenoliths were processed using standard mineral separation techniques to recover whole grains (hand crushing, magnetic separation, and heavy liquid separation). Coupled U-Pb and trace element analyses on zircon, rutile, and apatite were done by LASS inductively coupled plasma mass spectrometry at UC Santa Barbara (Kylander-Clark, 2013, 2017. Analyses were conducted on polished grain interiors for conventional LASS analyses and on the unpolished exterior of whole grains for depth profiling. Analytical protocols are detailed in Supporting Information S1. The quoted uncertainties for all laser ablation data are 2σ and incorporate analytical uncertainty as well as additional uncertainties associated with reproducibility of secondary matrix-matched reference materials (Horstwood et al., 2016). Reported uncertainties for trace-element data are 2σ and only include analytical uncertainties. All U-Pb data are shown in Tera-Wasserburg concordia plots (Vermeesch, 2018) and are uncorrected for common-Pb. The quoted U-Pb dates for apatite and rutile are common-Pb corrected using a common-Pb intercept defined by a linear regression through the U-Pb data or assuming a Stacey and Kramers (1975) initial 207 Pb/ 206 Pb composition. The quoted dates for zircon are 207 Pb/ 206 Pb dates that are within 5% concordance of their respective 206 Pb/ 238 U date.
The sampling resolution of laser ablation spot analyses described above (50-μm-diameter spots for rutile) is too coarse to resolve age gradients. For this reason, depth profiling was done using continuous laser pulsing at a low frequency (e.g., Apen et al., 2020;Cottle et al., 2009;Garber et al., 2020;Holder et al., 2019;Smye & Stockli, 2014). Details regarding depth profiling are also reported in Supporting Information S1. Final pit depths were measured on the SEM and are on average ∼13 μm. The depth profile data are reported at 1 s intervals and the reported U-Pb uncertainties take into account the reproducibility of multiple reference rutiles (∼6% on 238 U/ 206 Pb and ∼4% on 207 Pb/ 206 Pb). In addition to U-Pb ratios, we also report Pb concentration profiles, 7 of 21 which we divided into radiogenic Pb (Pb*) and non-radiogenic (common) Pb (Pb c ) (see also Garber et al., 2020;Holder et al., 2019). The total concentrations of 207,206 Pb at each depth interval were calculated using the measured 238 U/ 206 Pb ratio, 207 Pb/ 206 Pb ratio, and U concentration and assuming an invariant 238 U/ 235 U ratio of 137.818 (Hiess et al., 2012). Concentrations of 207,206 Pb* were determined using common-Pb corrected dates and measured U concentrations. From this, we calculated the proportion of Pb c by subtracting 207,206 Pb* from the total 207,206 Pb concentrations. The final 207,206 Pb* profiles were further corrected for radiogenic Pb ingrowth following kimberlite eruption at 360 Ma.

Testing Thermal Histories With Depth Profiling
In order to quantify temperatures and timescales of diffusion using U, Pb, and other elements, we consider the inverse error functions of normalized elemental gradients (Crank, 1979), expressed as: where C x is the concentration at a given position (x) within the grain, C r is the concentration at the rim, C c is the concentration at the core (or in this study, the bottom of the ablation pit), D(T) is temperature-dependent diffusivity, and t is the time that the mineral spent at temperature T.
The length-scales of diffusion for elements in rutile are constrained by experiments. Diffusivities of Pb, Al, Si, Cr, Fe, Zr, Nb, Hf, and Ta in rutile vary by orders of magnitude, with relative diffusivities at a reference temperature of 700°C in the following order: (Cherniak, 2000;Cherniak & Watson, 2019;Cherniak et al., 2007;Marschall et al., 2013;Sasaki et al., 1985).
In addition to diffusion modeling, comparing U-Pb dates of different accessory phases aids in reconstructing thermal histories. Experimentally-and empirically derived Pb diffusion coefficients in rutile and apatite suggest that rutile has a higher Pb closure temperature than apatite (e.g., Smye et al., 2018). For simple linear cooling and similar grain sizes, one would then predict U-Pb dates from rutile to be exclusively older than apatite (Dodson, 1973). We discuss the thermal history of the Siberian craton lower crust within this framework.

Thermobarometry and Pseudosection Modeling
The stability field of the observed garnet-clinopyroxene-plagioclase-rutile assemblage in the mafic garnet granulites (samples 01-95, 79-14, and 01-34) encompasses a large P-T range (1.0-1.6 GPa and 600°C-900°C). These estimates are consistent with the intersections of GADS, AJQ, and Fe-Mg garnet-pyroxene reactions determined using THERMOCALC (0.9-1.5 GPa and 750°C-900°C), despite the absence of quartz in the xenolith's main mineral assemblage. The presence of relict orthopyroxene (occurring exclusively as cores within clinopyroxene) in mafic garnet granulite 01-95 documents crystallization at slightly lower pressures; assuming relict orthopyroxene was in equilibrium with garnet, Ca-in-orthopyroxene thermometry suggests earlier metamorphism at 0.9 GPa and 825°C. Temperatures based on the composition of garnet and clinopyroxene rims indicate cooler temperatures of 630°C-710°C. Intersections of garnet isopleths in Perple_X did not produce results compatible with the P-T data from endmember equilibria (see Supporting Information S1 for further details); this could reflect contamination of the xenolith's whole rock compositions by the kimberlite and associated fluids (e.g., M. Y. Koreshkova et al., 2011) or limitations in the application of activity models (e.g., pyroxene; Forshaw et al., 2019).
No whole-rock data exists for felsic granulite 36-14 so we rely on the positions and intersections of net-transfer and exchange reactions determined from THERMOCALC. The average P-T is 1.3 ± 0.4 GPa and 950°C ± 120°C; Fe-Mg temperatures from garnet and orthopyroxene rims are 700°C-750°C. Amphibole-plagioclase thermobarometry for both tonalites and amphibolites yields 0.3-0.5 GPa and 700°C-750°C. Only apatite from mafic garnet granulites 01-95 and 79-14 and felsic garnet granulite 36-14 has sufficient U concentrations to obtain geologically useful U-Pb dates (>0.5 ppm U). Apatite in felsic garnet granulite 36-14 yields U-Pb dates that scatter between 1.6 Ga and 600 Ma (all assuming a Stacey and Kramers common-Pb composition; Figures 3c and 3d). Most of the apatite in mafic garnet granulite 01-95 has relatively low U (<2 ppm U) and does not produce useful U-Pb dates, but 13 spots with >2 ppm U yield common-Pb corrected dates between 700 Ma and 300 Ma (Figures 3c and 3d). Mafic garnet granulite 79-14 contains apatite with common-Pb corrected U-Pb dates that range from 2.0 Ga and 400 Ma (Figures 3c and 3d). Many apatite grains from xenolith 01-95 contain monazite exsolution lamellae (confirmed by EDS and BSE images; see Figure 3c), and mixing between apatite and monazite is apparent in trace element versus U-Pb date patterns. Spot analyses older than 1.1 Ga define a linear array, with the 2.0 Ga endmember of this array characterized by greater (Gd/Yb) N ratios and higher Th/U ratios (Figures 3c and 3d). Rutile from the garnet granulite xenoliths shows significant age variability. In felsic granulite 36-14 rutile the oldest common-Pb corrected U-Pb date is ca. 1.6 Ga and the youngest date is ca. 760 Ma. Within each grain Zr-in-rutile temperatures are internally coherent, but temperatures vary among different grains (740°C-810°C, assuming a pressure of 1.3 GPa; Figure 3e). Mafic granulites 01-95 and 79-14 both have common-Pb corrected U-Pb dates that range from ca. 1.1 Ga to 360 Ma. Like the felsic garnet granulite, Zr-in-rutile temperatures in both samples are internally consistent but differ between grains (700°C-940°C at 1.2 GPa for sample 01-95 and 820°C-840°C at 0.8 GPa for sample 79-14). Rutile in mafic garnet granulite 01-34 has low U (<0.5 ppm) and produces common-Pb corrected dates with large uncertainties between 500 and 1,600 Ma ( Figure 3e). Zr-in-rutile temperatures for this sample are fairly uniform at 875°C-890°C (assuming a pressure of 1.5 GPa). Where present, coexisting apatite usually displays a similar range of U-Pb dates as observed in rutile (Figure 4).

U-Pb and Trace Element Data
Of the analyzed samples, granulites 36-14 and 01-95 contain rutile with enough U (<1 ppm) to produce robust U-Pb age depth profiles. In mafic garnet granulite 01-95, common-Pb corrected U-Pb date and Pb* concentration profiles increase linearly from 360 Ma at the rim to 1.2 Ga-600 Ma dates at the bottom of the pit (though two profiles are wholly 360 Ma throughout the profile; Figure 6). The outermost 1-2 μm portions of some rutile grains in mafic garnet granulite 01-95 show greater concentrations of U, Pb c , Zr, and Hf than in the interiors. Beyond the discrete rims, U, Al, and Zr gradually increase toward the grains' interior over a distance of 5-8 μm whereas all other elements show no significant variations with depth ( Figure 6). Rutile depth profiles for felsic garnet granulite 36-14 have common-Pb corrected U-Pb dates and Pb* concentrations that gradually increase from the rim toward the grain interior over a distance of 4-6 μm ( Figure 6); U-Pb dates at the outermost rim approach the 360 Ma eruption age. In felsic garnet granulite 36-14, rutile depth profiles show that Al concentrations also gradually increase toward the grain interior over 2-8 μm depths (Figure 6). Silicon and Fe also vary, showing relatively higher concentrations on the rims that decrease toward the core ( Figure 6). All other trace-elements are largely invariant with depth. 6. Discussion

Metamorphic Histories
In order to link the thermal histories of the xenoliths to heat fluxes in the lithosphere, it is necessary to evaluate what portions of the xenoliths' P-T history correspond to thermal conditions shortly prior to eruption (ca. 360 Ma) versus older episodes of metamorphism (e.g., during Proterozoic craton amalgamation). The tonalite and amphibolite xenoliths contain zircon with concordant Archean U-Pb dates and apatite with solely Paleoproterozoic U-Pb dates ( Figure 3). The data suggest that these rocks formed in the Archean and were later heated to temperatures high enough to fully reset the apatite (e.g., >700°C for >1 Myr) or recrystallized apatite during ca. 1.8 Ga amalgamation of the Siberian craton (Moyen et al., 2017;Paquette et al., 2017;Rosen et al., 2006). Of prime relevance, the presence of 1.8-1.7 Ga apatite in these rocks implies rapid cooling from 700°C to 750°C to below 400°C-450°C (nominal Pb closure temperature assuming 50-200 μm diameter grains and relevant cooling rates of 1°C-10°C/Myr) at 0.4-0.5 GPa, or 15-20 km depth. These rocks remained below 400°C-450°C and did not experience significant re-heating since 1.8-1.7 Ga (Figure 5b).
Establishing P-T paths for high-pressure granulites like those investigated in this study is difficult given the propensity of cation-exchange thermometers to be more affected by retrogression than net-transfer reactions (e.g., Frost & Chacko, 1989;Pattison & Bégin, 1994). The 0.8-1.5 GPa equilibration pressures (or 30-60 km depth) determined from pseudosection modeling and GADS/GAHS thermobarometry likely reflect (near) peak metamorphic conditions at 1.9-1.  Fe-Mg exchange thermometry in garnet and pyroxene rims indicate cooling that may correspond to some amount of decompression/exhumation. The cooler rim temperatures, however, are also unlikely to represent ambient temperatures at 360 Ma, but instead reflect closure of diffusional Fe-Mg exchange between garnet and pyroxene (Harley, 1989). In this regard, the U-Pb system in rutile and apatite, which undergo diffusional closure at significantly cooler temperature (Smye et al., 2018), provide further insights into the thermal evolution of the Udachnaya lower crust.
In garnet granulite xenoliths, rutile and apatite record dates that are significantly younger than the well-documented 2.0-1.8 Ga period of craton amalgamation recorded in zircon (Figure 4). In particular, the oldest U-Pb rutile dates from garnet granulite xenoliths 79-14, 01-95, and 36-14 are 1.6-1.1 Ga (Figure 4). These dates could be interpreted as the time when the rocks cooled below rutile Pb closure following slow cooling from 1.8 Ga or a discrete re-heating event at 1.6-1.1 Ga. If the rocks underwent slow cooling through the Pb closure temperature of rutile and apatite since 1.8 Ga, one would predict that apatite dates would be as young or younger than rutile dates (calculated Pb closure temperatures for rutile are 500°C-600°C and for apatite are 350°C-450°C using relevant 50-500-μm-radius grains and cooling rates of 0.01°C-0.1°C/Myr). In mafic garnet granulite 01-95, apatite U-Pb dates are consistently younger than the oldest rutile dates, and in felsic garnet granulite 36-14, the oldest apatite U-Pb dates are within uncertainty of the oldest 1.6 Ga rutile date (Figure 4). The exception is mafic garnet granulite 79-14, where about half of the apatite dates are significantly older than rutile dates (up to 2.0-1.8 Ga; Figure 4). The presence of monazite exsolution lamellae in apatite from this sample complicates comparisons between apatite and rutile U-Pb dates; because Pb diffusion in monazite occurs at significantly higher temperatures than apatite (Cherniak et al., 2004), any spot analyses incorporating monazite will produce a spread of U-Pb dates unrelated to thermally mediated volume diffusion. Based on trends in U-Pb date versus Gd/Yb and Th/U ratios, we conclude that ca. 2.0 Ga U-Pb dates represent the timing of monazite exsolution formation; apatite spot dates >1.1 Ga are due to mixed monazite/apatite analyses and are therefore not geologically meaningful. Accounting for this, the oldest apatite U-Pb date in mafic garnet granulite 79-14 is younger than the oldest rutile U-Pb date.
Based on apatite and rutile U-Pb dates, we suggest that the ca. 1.6-1.1 Ga upper age intercepts represent Mesoproterozoic re-heating of the lower crust by a discrete thermal event, and are not due to slow monotonic cooling through rutile/apatite Pb closure since 1.8 Ga. In support of this interpretation, we note that the oldest apatite and rutile U-Pb dates in felsic garnet granulite 36-14 are both 1.6 Ga, requiring rapid cooling through the nominal Pb closure temperatures of both phases (Figure 5b). No 1.6-1.1 Ga crystallization dates have been documented previously in outcrops or in xenoliths from the central Siberian craton (M. Koreshkova et al., 2009;Paquette et al., 2017;Priyatkina et al., 2016), with the exception of rift-related dike swarms in the Anabar shield and Olenek uplift that record similar ages (see Cherepanova et al., 2013;Gladkochub et al., 2012; and references therein). We speculate that a similar dike emplacement event occurred around Udachnaya at 1.6-1.1 Ga, creating a thermal pulse that heated the lower crust above the Pb closure temperatures of rutile and apatite (akin to resetting of rutile in Slave craton xenoliths by dike swarm intrusion; Davis, 1997). The spread of rutile U-Pb dates down to 360 Ma provides further insight into the thermal evolution of the deep crust following the Mesoproterozoic.

Re-Heating the Lower Crust?
The objective of multi-element depth profiling is to ascertain whether age and geochemical gradients within minerals are a result of thermally mediated volume diffusion or re-/neo-crystallization. If Pb was diffusing out of rutile during slow cooling in the lower crust over billion year timescales, this would imply residence at 500°C-600°C (temperatures of partial Pb retention in rutile) prior to eruption (path I in Figure 5b) (e.g., Schmitz & Bowring, 2003;T. Blackburn et al., 2011). On the other hand, Pb-loss could have been induced by a thermal event associated with kimberlite magmatism, which would imply long-term ambient lower crustal temperatures cooler than the nominal Pb closure temperature of rutile (path ii in Figure 5b). A third alternative is that gradients reflect re-/neo-crystallization in the deep crust during or shortly before kimberlite magmatism.
In order to evaluate the different scenarios discussed above, we calculated the inverse error functions of the normalized Pb and trace elements gradients in the 22 rutile depth profiles in mafic garnet granulite 01-95 and 18 rutile depth profiles in felsic garnet granulite 36-14. Gradients produced by thermally mediated volume diffusion are expected to conform to an error function-reflected in the linearity of the inverse error function-and the slopes of the gradients are predicted to scale by 1 √ 4 × ( ) × (Crank, 1979; see also Section 4.4). Figure 7 illustrates representative depth profiles of rutile in mafic garnet granulite 01-95 (profile 13) and felsic garnet granulite 36-14 (profile 9). In profile 13 of mafic garnet granulite 01-95, a 1-μm-thick rim is followed by monotonous increases in Pb*, Al, and Zr toward the grain interior ( Figure 7a); the latter portion of the profiles conform to an error function, but the corresponding slopes yield values that all overlap within uncertainty (Figure 7b). Similarly, rutile depth profile 9 from felsic garnet granulite 36-14 shows Pb* and Al gradients that also follow an error function with indistinguishable slopes (Figure 7d). The available experimentally determined diffusivities for Al, Pb, and Zr predict decoupling of these elements across relevant crustal temperatures (D Al ≈ 1E−9*D Pb* and D Al ≈ 1E−9*D Zr at 600°C, e.g., Cherniak and Watson (2019) and references therein). Therefore, if the gradients captured in the rutile interiors all formed by the same process, they cannot be explained by diffusion over billions of years as this would produce systematic variability in the slopes in Figures 7b and 7d. While it is possible that the Al profiles developed during an earlier episode of intense heating and the Pb* profile during later slow cooling, coincidentally producing gradients show nearly identical inverse error functions, we argue that the depth profiles capture rutile crystallization and/or shorter-term heating near the time of entrainment in the kimberlite.
We propose that partial U-Pb resetting of rutile from the Udachnaya lower crustal xenolith suite reflects transient heat advection from and interactions with melts/fluids prior to kimberlite eruption. The observation that different rutile depth profiles exhibit different U-Pb dates and elemental concentrations at the rim suggests flux-limited element transport at the grain boundary (Kohn et al., 2016;Smye et al., 2018). The development of elemental gradients could have been controlled by adjacent phases (garnet, pyroxene, feldspar) or a fluid phase that imposed different boundary conditions across at the grain edge. The involvement of fluids is supported texturally by the presence of pockets of fine-grained material and secondary biotite in the xenoliths, which may be remnants of a melt/fluid (Figure 2). There is also ample evidence for metasomatism in the Udachnaya peridotite xenolith suite (Agashev et al., 2013;Boyd et al., 1997;Golovin et al., 2018;Ionov et al., 2010;M. Y. Koreshkova et al., 2011;V. Shatsky et al., 2008). Agashev et al. (2013) identified multiple episodes of metasomatism in the peridotite xenoliths: carbonatite metasomatism (addition of Ca, Al, and LREE); silicate melt metasomatism (enrichment of REE, Y, and Zr resulting in formation of garnet and clinopyroxene in otherwise highly depleted peridotites; see also Boyd et al. (1997)), including enriched-basaltic melt metasomatism (local enrichment of Fe and Ti). In addition, Golovin et al. (2018) described Na-K-Ca-rich melt inclusions in olivine from sheared peridotite xenoliths that resulted from interactions with primitive kimberlite magmas in the mantle prior to eruption. Regarding Udachnaya granulite xenoliths, Koreshkova et al. (2011) noted that kimberlite melt interactions are required to account for bulk enrichments of LREE, Th, and U (with contributions from kimberlite-related fluids; Kamenetsky et al., 2007).  Figure 6). Beyond this rim, there is monotonous increases in Al, Zr, and Pb* toward the grain interior that could be related to thermally mediated volume diffusion. This latter portion of the profiles conform to an error function, but the corresponding slopes yield values that all overlap within uncertainty (panel b). (c and d) Similarly, rutile depth profile 9 from felsic garnet granulite 36-14 (panels (c and d)) shows Pb* and Al gradients also follow an error function form and produce indistinguishable slope values. Data included in the dashed area were used to calculate slopes. Thermally mediated volume diffusion over billions of years would be expected to produce differences in the slopes of the inverse error functions that vary by . The profiles therefore correspond to transient heating or (re)crystallization immediately prior to eruption. The involvement of one or more of the putative metasomatic media discussed above could be responsible for the development of the complex age and elemental patterns observed in the rutile depth profiles. Trace-element rutile-melt partition coefficients-especially the high field strength elements like Zr, Hf, Nb, Ta-depend greatly on melt composition (Foley et al., 2000;Klemme et al., 2005; and references therein). For example, Klemme et al. (2005) note the compatibility of high field strength elements in rutile over silicate melt with an andesitic or rhyolitic composition. Neither of these melt compositions would be wholly relevant to the Udachnaya xenoliths as they are too silicic and viscous, but we speculate that the enriched U, Hf, and Zr at the rims relative to the rutile interior could have resulted from a combination of the melt compositions documented in other Udachnaya xenoliths. An additional consideration is that portions of the rutile depth profiles record diffusion (see Figures 7a and 7c) in response to heating at relatively high temperatures and over short time scales.
Thin rims (1-2 μm) observed on rutile in mafic garnet granulite 01-95 are enriched in U, Pb c , Hf, and Zr. We interpret this feature as reflecting neo-or re-crystallization of rutile associated with precursory kimberlitic magmatism. This process is distinct from the eruption process itself because: (a) timescales of transport to the surface are too short to cause significant Pb diffusion (Doucet et al., 2014;Schmitz & Bowring, 2003;Spera, 1984); and (b) shallower tonalite and amphibolite xenoliths from the same kimberlite pipe have apatite with undisturbed Paleoproterozoic U-Pb isochrons (Figures 3c and 3d). The elevated Zr concentrations in the 01-95 rutile rims indicate heating to high temperatures prior to eruption (e.g., 9,000 ppm Zr in profile 13 equates to >1000°C, and higher still in other profiles; Figure 6). Thus, in addition to metasomatism, it is probable that such extreme temperatures prior to eruption induced Pb loss. Assuming that a portion of the rutile Pb profiles does reflect thermally mediated volume diffusion (Figure 7), bounds can be placed on the duration of heating. The observed Pb* gradients in rutile from mafic garnet granulite 01-95 are consistent with average heating time-scales of <18, <3, and <0.6 Myr at 900°C, 1000°C, and 1100°C, respectively; time-scales of heating based on Pb* profiles in rutile from felsic garnet granulite 36-14 are <4, <0.6, and <0.1 Myr at 900°C, 1000°C, and 1100°C, respectively. Because slower-diffusing Al shows similar gradients to Pb* in the rutile profiles ( Figure 6), the time-scales of heating calculated from Pb gradients are maxima. Nonetheless, the transient nature of this event is consistent with timescales of advective heating prior to eruption proposed for other xenolith suites, such as those from the North Atlantic craton (Smit et al., 2016), Kaapvaal craton (Jollands et al., 2018) and Tibetan Plateau (Hacker et al., 2000).
If deep crustal heating was depth-dependent, then samples with greater U-Pb resetting resided at greater depths relative to those with minimal overprinting (e.g., the mafic granulites were deeper than the felsic granulites; Figure 4). However, there is no correlation between the extent of rutile U-Pb resetting and peak pressure, and advective heating may have been heterogeneous. This inference is supported by data on Udachnaya mantle xenoliths that suggest similar heterogeneous heating in the mantle lithosphere (e.g., Liu et al., 2022).
Finally, if the U-Pb heterogeneity was caused by processes that operated immediately preceding kimberlite eruption, the preservation of Proterozoic rutile and apatite dates indicates that the deep crust of the Siberian craton generally resided at temperatures cool enough that Pb was not actively diffusing out of the rutile and apatite crystal lattice prior to this late-stage heating. For the typical 50-200 μm-wide apatite grains this translates to temperatures below 400°C for ∼1 Gyr (Figure 8).

Crustal Heat Production Models for the Siberian Craton
The data discussed above provide new temperature constraints that can be used to calculate crustal heat production. Previous approaches to determining crustal heat production include using measured abundances of U, Th, and K in lower-crustal granulite terranes exposed at the surface or reconstructed from xenoliths (e.g., Gruber et al., 2021;Nicolaysen et al., 1981;. Each of these approaches has advantages and disadvantages. The former is accessible and enables extensive sampling, but granulites terranes may have only resided In order to preserve 1.1 Ga rutile and apatite dates, lower-crustal temperatures were below 400°C (down arrow at 10 3 Myr). Rutile Pb diffusivity data from Cherniak (2000) and apatite Pb diffusivity data from Cherniak et al. (1991). transiently in the mid-deep crust and may not be representative of the lower crust (e.g., Bohlen & Mezger, 1989). The latter provides in situ samples of deep crust, but suffers from unquantifiable uncertainties because it is unknown how representative a xenolith suite is of a crustal column in general (i.e., whether the whole column was sampled and if it was sampled in direct proportion that it exists within the crust). The approach we take exploits the temperature bounds imposed by U-Pb thermochronology to constrain heat fluxes near the Moho. That is, our data require long-term Moho temperatures below 400°C (temperatures cool enough to preserve Mesoproterozoic rutile and apatite U-Pb dates).
As a secondary check, we compare permissible geothermal gradients to P-T data from garnet peridotite xenoliths also entrained by the Udachnaya kimberlite(e.g., Liu et al., 2022;R. L. Rudnick et al., 1998). If the garnet peridotite xenoliths were in equilibrium at ambient mantle lithosphere conditions, any geotherm encompassing the mantle xenolith data would also yield Moho temperatures cooler than ∼400°C (i.e., temperatures for the lower crust and Moho should be consistent with closed-system behavior of Pb in rutile and apatite).
We assess different heat production models for the lithospheric column at Udachnaya using steady-state one-dimensional conductive geotherms (following Hasterok and Chapman (2011) and Furlong and Chapman (2013)). Our goal is not to determine the exact crustal heat production value for Udachnaya, but rather to evaluate general trends in geotherms that match the available xenolith data. We evaluate three different distribution models of heat production within the crust: (a) uniform crustal heat production; (b) stepwise decrease in crustal heat production with depth; and (c) exponential decrease in heat production with depth ( Figure 9).
In model A, we chose a uniform average Archean crustal heat production of 0.65 μW/m 3 (Jaupart & Mareschal, 2014). In model B, we calculate heat production values based on whole-rock measurements of U, Th, and K in crustal xenoliths (Moyen et al., 2017) and adopt a stepwise change in heat production assuming two layers: an upper crust (0-15 km) consisting of tonalite only and lower crust (15-45 km). The average upper crustal tonalite heat-production is 0.36 μW/m 3 , whereas the lower crustal average heat production of mafic granulites is 0.26 μW/m 3 for a total crustal heat production of 0.29 μW/m 3 . These heat production values are based on bulk whole-rock measurements, which we note are likely biased toward higher values given the likelihood of contamination of the whole rocks by kimberlite magma or related fluids (Golovin et al., 2018Kamenetsky Figure 9. Conductive geotherm models. Model A assumes uniform crustal heat production (average Archean heat production value of Jaupart and Mareschal (2014)). Model B assumes a stepwise change in heat production at 15 km. Model C assumes an exponential decrease in crustal heat production. All models employ heat production value of 0.006 μW/m 3 for the mantle lithosphere (after McIntyre et al., 2021). Moho is at 45 km (after Cherepanova et al., 2013). Garnet peridotite P-T data are shown in gray, with darker gray spots representing sheared peridotites (data from Liu et al. (2022)). In models (B and C), the surface heat flow values of 20-35 mW/m 2 are consistent with preservation of Proterozoic apatite and rutile U-Pb dates (<400°C, represented by the gray bar at 0.8-1.4 GPa) and the peridotite P-T data. These surface heat flow values are generally higher than the average surface heat flow measurements around Udachnaya (19 ± 3 mW/m 2 ). et al., 2007; M. Koreshkova et al., 2009). Model C assumes an exponential decrease in heat production, starting with the average tonalite heat-production at the surface and grading to an average mafic granulite value at 45 km depth for an integrated crustal heat production of 0.30 μW/m 3 . All of the models utilize a mantle lithosphere heat production value of 0.006 μW/m 3 (after McIntyre et al., 2021), a Moho depth of ∼45 km , and are calculated with fixed surface heat flow values between 20 and 50 mW/m 2 (encompassing the range observed in cratons globally; Jaupart & Mareschal, 2014).
Some interesting observations derive from these modeled geotherms. First, the geotherms encapsulating the peridotite P-T data suggest Moho temperatures that are cooler than the temperature bounds suggested by apatite and rutile U-Pb dates (<400°C; Figure 8). Notably, the cool geotherms are consistent with relatively deep Curie depth estimates for the Siberian craton (>35 km; Gard & Hasterok, 2021), highlighting potential magnetization below the Moho (Ferré et al., 2013). Second, in all three models the Udachnaya peridotite P-T data require geotherms that are uniformly higher than the measured surface heat flow value of 19 ± 3 mW/m 2 (represented approximately by the 20 mW/m 2 geotherms in Figure 9). The peridotite xenolith P-T data do indeed reflect some lithospheric heating prior to eruption (see Liu et al., 2022), such that the ambient geotherm at 360 Ma should lie closer to the cooler geotherms ( Figure 9). Even so, reconciling the Udachnaya peridotite P-T data with the observed low surface heat flow around Udachnaya would require substantial cooling of the lower crust since kimberlite emplacement 360 Ma; the secular decay of radioactive elements alone could not account for this (∼100°C/Gyr; Jaupart et al., 2016). The Siberian xenolith data are, however, consistent with the global average surface heat flow range of 20-50 mW/m 2 for cratons (Jaupart et al., 2016). Given these considerations, we assert that surface heat flow measurements made in the central Siberian craton are inaccurate, having been potentially compromised by the thick and permeable sedimentary and volcanic cover , largely infiltrated by brines (Alexeev et al., 2007(Alexeev et al., , 2022Kitayama et al., 2021), or the chilling effect of Pleistocene glaciations (Birch, 1948).
Of the three crustal models, the constraints from both lower crustal temperatures and mantle xenolith P-T arrays are best matched by models B and C (Figure 9), which have extremely low total crustal heat production in the Siberian craton, on the order of ∼0.3 μW/m 3 . This could reflect a bulk mafic crust composition (e.g., R. L. Rudnick & Gao, 2014) or extensive depletion of heat-producing elements during high-temperature metamorphism and melt extraction associated with cratonization (e.g., Carlson et al., 2005). The presence of lower-crustal felsic granulites at Udachnaya and nearby sites require some amount of silicic material at depth (e.g., Moyen et al., 2017;M. Y. Koreshkova et al., 2011;Shatsky et al., 2019). Previous studies have shown that accessory phases in high-grade felsic rocks retain Th even after undergoing extensive melt depletion (e.g., Alessio et al., 2018;Bea, 1996; R. L. Rudnick & Fountain, 1995), such that high heat-producing felsic lithologies may be subordinate to mafic lithologies in the deep crust of the Siberian craton.

Conclusions
New laser ablation split-stream U-Pb and trace-element data from accessory phases in xenoliths from the ca. 360 Ma Udachnaya kimberlite, Siberian craton, show that the U-Pb system is sensitive to transient heating of the lower crust and not just slow cooling. Length-scales of U-Pb and trace element gradients revealed by depth profiling of rutile suggest that depth-dependent heating of the deep crust occurred for <1 Myr before eruption. Preserved Mesoproterozoic dates in both apatite and rutile long-term lower-crustal temperatures were cool enough that Pb was not actively diffusing out of the both phases (<400°C). The lower-crustal temperature bounds imparted by these data are consistent with P-T arrays of Udachnaya peridotite xenoliths that suggest relatively cool geothermal gradients, signifying that the mantle xenoliths accurately capture the thermal state of the lithosphere prior to eruption. Combined, the xenolith data imply low crustal heat production for the Siberian craton (∼0.3 μW/m 3 ). Although exceedingly low, such values still produce surface heat flow values higher than those measured around Udachnaya (average 19 mW/m 2 ), suggesting that the surface heat flow measurements are too low (possible compromised by thick sedimentary cover and deep regional aquifers). Collectively, the xenolith data support low crustal heat production in the Siberian craton, implying a relatively mafic bulk crust composition or a HPE-depleted felsic bulk crust.

Data Availability Statement
The data are archived in the DRYAD data repository and the DOI will be made publicly accessible upon publication of this manuscript (https://doi.org/10.25349/D94K79).