Thermal State, Thickness, and Composition of the Lithospheric Mantle beneath the Upper Muna Kimberlite Field (Siberian Craton) Constrained by Clinopyroxene Xenocrysts and Comparison with Daldyn and Mirny Fields

To gain better insight into the thermal state and composition of the lithospheric mantle beneath the Upper Muna kimberlite field (Siberian craton), a suite of 323 clinopyroxene xenocrysts and 10 mantle xenoliths from the Komsomolskaya-Magnitnaya (KM) pipe have been studied. We selected 188 clinopyroxene grains suitable for precise pressure (P)-temperature (T) estimation using single-clinopyroxene thermobarometry. The majority of P-T points lie along a narrow, elongated field in P-T space with a cluster of high-T and high-P points above 1300 ◦C, which deviates from the main P-T trend. The latter points may record a thermal event associated with kimberlite magmatism (a “stepped” or “kinked” geotherm). In order to eliminate these factors, the steady-state mantle paleogeotherm for the KM pipe at the time of initiation of kimberlite magmatism (Late Devonian–Early Carboniferous) was constrained by numerical fitting of P-T points below T = 1200 ◦C. The obtained mantle paleogeotherm is similar to the one from the nearby Novinka pipe, corresponding to a ~34–35 mW/m2 surface heat flux, 225–230 km lithospheric thickness, and 110–120 thick “diamond window” for the Upper Muna field. Coarse peridotite xenoliths are consistent in their P-T estimates with the steady-state mantle paleogeotherm derived from clinopyroxene xenocrysts, whereas porphyroclastic ones plot within the cluster of high-T and high-P clinopyroxene xenocrysts. Discrimination using Cr2O3 demonstrates that peridotitic clinopyroxene xenocrysts are prevalent (89%) among all studied 323 xenocrysts, suggesting that the Upper Muna mantle is predominantly composed of peridotites. Clinopyroxene-poor or -free peridotitic rocks such as harzburgites and dunites may be evident at depths of 140–180 km in the Upper Muna mantle. Judging solely from the thermal considerations and the thickness of the lithosphere, the KM and Novinka pipes should have excellent diamond potential. However, all pipes in the Upper Muna field have low diamond grades (<0.9, in carats/ton), although the lithosphere thickness is almost similar to the values obtained for the high-grade Udachnaya and Mir pipes from the Daldyn and Mirny fields, respectively. Therefore, other factors have affected the diamond grade of the Upper Muna kimberlite field. Minerals 2020, 10, 549; doi:10.3390/min10060549 www.mdpi.com/journal/minerals Minerals 2020, 10, 549 2 of 20


Introduction
Mantle-derived xenoliths and xenocrysts from kimberlites provide direct information about the composition and thermal state of the lithospheric mantle underlying the ancient cratons at time of kimberlite eruption [1][2][3][4][5][6][7][8]. Constraints on mantle paleogeotherms (temperature change with depth) allow determination of the thermal gradient and thickness of the lithosphere [2,3,[9][10][11]. Besides the importance of understanding mantle dynamics and origin of cratons, knowledge of the thermal state and thickness of the subcratonic lithospheric mantle is also essential for predicting its diamond potential in the region of interest [12][13][14][15][16].
The common petrological approach to constraining the mantle paleogeotherm beneath cratons is based on the calculation of the P-T conditions of equilibration for mantle xenoliths using thermobarometry based on the exchange of components between multiple minerals [6,[17][18][19][20][21][22]. Another option is to use P-T data obtained from mantle xenocrysts, which represent single grains of disintegrated mantle xenoliths [2,15,16,23]. Single-mineral thermometers and barometers are well-developed for peridotitic clinopyroxenes and garnets [12,24]. The advantages of the second approach are as follows: (1) it can be used in cases where mantle xenoliths are absent or rare, e.g., due to secondary alteration; (2) xenocrysts are routinely obtained during diamond exploration; (3) it requires fewer analytical resources because each mineral grain is assumed to represent one mantle xenolith with an inferred mineralogy, for instance, single-clinopyroxene thermobarometry requires 3 or 4 times fewer microprobe analyses than conventional thermobarometry using mineral pairs; (4) its relative cheapness and efficiency (only one mineral is used in the P-T calculations) thus makes it possible to obtain a statistically meaningful number of P-T points, which is very important for robust estimation of the paleogeotherm.
Our knowledge of the thermal state of the mantle beneath the Siberian craton in Late Devonian-Early Carboniferous time comes mainly from numerous studies of mantle xenoliths from the Udachnaya kimberlite pipe, Daldyn field [1,6,10,17,[25][26][27][28]. Paleogeotherm constraints beneath the Upper Muna field were presented only in a few papers with slightly different results for the lithospheric thickness [2,23,29,30] because different approaches have been used.
In order to further evaluate the thermal state and structure of the lithospheric mantle beneath the Upper Muna kimberlite field at the time of kimberlite emplacement, we have analyzed a large suite of clinopyroxene xenocrysts from the Komsomolskaya-Magnitnaya (KM) pipe using the FITPLOT numerical fitting program (written by Dan McKenzie in 1988 and introduced by Mather et al. in 2011 [7], and based on a geotherm model proposed in [6,31]) and compare our results with the geotherm for the nearby Novinka pipe obtained in the same way [23]. To test the validity of single-mineral clinopyroxene thermobarometry, we compare our results with P-T data obtained from mantle xenoliths from the KM pipe. Discrimination of clinopyroxenes is based on Cr 2 O 3 and demonstrates that the mantle column beneath the KM pipe appears to be predominantly composed of peridotite. We also estimated paleogeotherms beneath the Daldyn and Mirny kimberlite field using P-T data for mantle xenoliths from literature and the FITPLOT program to compare the thermal state of the lithospheric mantle between different tectonic terranes of the Siberian craton. We used FITPLOT to calculate paleogeotherms because it has been widely used by petrologists in the last decade and allows comparison of the thermal state of the lithospheric mantle for different regions in a quantitative manner [7,13,15,23,32].

Geological Setting
The basement of the Siberian craton occupies~4 million km 2 and comprises Precambrian crust, mostly covered (~70%) by Riphean and Phanerozoic sediments and outcropping within the Anabar and
There are more than 20 known kimberlite fields in the Siberian craton located mostly within its northwestern part and grouped in the Yakutian kimberlite province [41][42][43]. Four distinct episodes of kimberlite magmatism are defined in the Siberian craton: Silurian-Early Devonian (420-400 Ma), Late Devonian-Early Carboniferous (350-380 Ma), Triassic (215-235 Ma), and Early Cretaceous-Jurassic (140-170 Ma) [44][45][46][47][48][49]. Diamond-rich kimberlites are overwhelmingly Late Devonian-Early Carboniferous in age. The Upper Muna, Daldyn, and Mirny fields are located within different terranes of the Anabar tectonic province: Daldyn, Markha, and Magan, respectively ( Figure 1). Within the Upper Muna field, 19 kimberlite bodies are known, and all are diamondiferous (3 dykes and 16 pipes, five of which are of commercial value and comprise the Verkhne-Munskoe deposit: Zapolarnaya, Deimos, Novinka, KM, and Poiskovaya) [43,50]. The Upper Muna kimberlites are located as a compact group in the Ulakh-Muna River valley, right tributary of the Muna River, and intruded into terrigenous-carbonate Middle to Upper Cambrian sediments. The KM pipe is situated close to the Novinka pipe (~100 m), studied in [23], and it is thought they have a common feeder dyke [43]. The KM pipe comprises three types of kimberlite rocks: kimberlite breccia, monticellite-bearing porphyritic kimberlite, and monticellite-free porphyritic kimberlite. According to terminology proposed in [43], kimberlite breccia is defined as kimberlite rocks containing >10% of country rocks, whereas porphyritic kimberlite hosts <10% of country rocks. There are no U-Pb perovskite or megacrystic-zircon age data for the KM pipe. Recent U-Pb perovskite and macrocrystic-zircon age determinations for six other kimberlite pipes from the Upper Muna field gave a narrow interval of kimberlite magmatism from 367 to 345 Ma [2,47,[51][52][53], corresponding to Late Devonian-Early Carboniferous. The emplacement age of KM is estimated as 334-382 Ma based on bulk groundmass K-Ar of kimberlite breccia [42,54]; Rb-Sr ages of phlogopite megacrysts from the pipe yielded 400-402 Ma [55]. These data are generally consistent with the Late Devonian-Early Carboniferous age of the Upper Muna field. The Novinka pipe has a U-Pb perovskite age of 355 ± 11 Ma [52].

Sample Descriptions
For this study, we used 323 clinopyroxene xenocrysts and 10 mantle xenoliths from the KM pipe. The heavy-mineral concentrates were obtained from material collected from the quarry left after bulk sampling inside the contour of the pipe at the present surface and represent a mixture of all exposed kimberlite types recognized at KM. The clinopyroxene grains were picked from the 0.5-4 mm fraction of heavy-mineral concentrates, without any preference to the size, color, or shape of the grains. The mantle xenolith collection is represented by one coarse spinel-garnet lherzolite, seven coarse garnet peridotites, and two porhyroclastic (sheared or high-T) peridotites.

Analytical Methods
Chemical analyses of clinopyroxene xenocrysts and minerals of mantle xenoliths from the KM pipe were obtained using a JEOL JXA-8100 microprobe at the Analytical Center for Multi-Elemental and Isotope Research of Siberian Branch of the Russian Academy of Science (V.S. Sobolev Institute of Geology and Mineralogy, Novosibirsk, Russia). Mineral compositions are given in (Tables S2 and S3).

Mineral Thermobarometry
For single-clinopyroxene thermobarometry, we used the combination of the enstatite-inclinopyroxene thermometer (T NT00 ) and Cr-in-clinopyroxene barometer (P NT00 ) proposed by Nimis and Taylor [24]. T NT00 is based on the exchange of the enstatite component between clinopyroxene and orthopyroxene. P NT00 is based on Cr exchange between clinopyroxene and garnet, where Cr in clinopyroxene is expressed as CaCr-Tschermak's component (CaCrAlSiO 6 ) [24]. Therefore, only clinopyroxene that was in equilibrium with orthopyroxene and garnet in the mantle can be suitable for T NT00 and P NT00 . Xenocrysts of such clinopyroxene are derived from garnet peridotites and pyroxenites and display specific compositional features (listed in Section 4.1) [23,24]. Moreover, P NT00 has intrinsic limitations related to the high sensitivity of pressure estimates to the activity of CaCr-Tschermak's components (a Cr = Cr−0.81(Na + K)·Cr/(Cr + Al) atoms per 6-oxygen formula unit), as low values of a Cr may provide large deviations, and to Cr# = Cr/(Cr + Al), because the barometer was calibrated on clinopyroxenes with Cr# ≤ 0.44. Ziberna et al. [23] revised the compositional filters to identify clinopyroxene compositions suitable for T NT00 and P NT00 , and we have used their protocol for optimum thermobarometry.
The P-T values for the studied mantle xenoliths from KM (Table S4, Section 4.5) and Mir pipes (Table S1, Section 4.5) were obtained using several combinations of thermometers and barometers, depending on the mineral pairs available, and the PTQuick software provided by Dmitry Dolivo-Dobrovolsky (Institute of Precambrian Geology and Geochronology, Russian Academy of Sciences, St. Petersburg).

Computation of the Paleogeotherm
Paleogeotherms were numerically fitted from the discrete P-T data using the FITPLOT program [6,7]. This program quantifies the lithosphere thickness, calculated as the intersection of the conductive geotherm with the adiabatic geotherm of the convecting mantle, and the surface heat flow. A detailed description of the way the geotherm is calculated can be found in Mather et al. [7].
The program requires several input parameters for the crust and mantle. We assumed a heat production rate in the upper mantle of 0.0 µW/m 3 and a potential temperature (T p ) for the asthenospheric isentrope of 1315 • C for all calculations. These input parameters are the same as Ziberna et al. [23] used for the nearby Novinka pipe. The FITPLOT program assumes a two-layer crust model that requires input of thickness and heat generation rates for the upper and lower crust; the values used in our calculations are shown in Table 1.

Clinopyroxene Thermobarometry
The sample selection for single clinopyroxene thermobarometry used the following protocol [23]: • We discarded grains for which at least one of three electron microprobe analyses do not satisfy the quality test, i.e., total cations calculated per 6 oxygen atoms fall outside 3.98-4.02 range (Tables S2  and S3).

•
Twenty-nine grains exhibited obvious compositional inhomogeneity, suggesting possible disequilibrium and were excluded from the study.

•
Thus, 323 homogeneous clinopyroxene xenocrysts were selected (Tables S2 and S3). For further evaluation, we used the average of three individual analyses per grain. This reduces uncertainties proportional to the square root of three.
• Equilibration with garnet was tested based on the Cr 2 O 3 vs. Al 2 O 3 discrimination diagram from [63], and 26 grains were discarded (Table S3).

•
All grains that were not excluded at previous stages satisfy the Cr# criterion, i.e., are within the 0.10-0.65 interval (Table S2).

•
A further 43 grains with aCr/Cr# ≤ 0.11 were discarded (Table S3), as aCr is the main building block in the P NT00 formulation.

•
Verification of clinopyroxene equilibration with orthopyroxene is difficult by means of simple compositional filters. However, it is recommended to use clinopyroxene with molar Ca/(Ca + Mg) ratio (Ca#) below 0.5 [23]. Six of the studied clinopyroxene xenocrysts have Ca# > 0.5, but they had already been excluded at previous steps.

•
It is also recommended to exclude samples with temperatures below 700 • C because T NT00 is not calibrated at <700 • C. Another reason to discard these points is that low estimates of temperature may indicate that clinopyroxene was not in equilibrium with orthopyroxene. We excluded one grain at this step.
Thus, according to the proposed filters [23], we accepted 188 clinopyroxene xenocrysts for T NT00 and P NT00 thermobarometry (Table S2). Judging by eye, the majority of accepted points lie along a narrow, elongated field in P-T space ( Figure 2a). However, there is a cluster of high-T points above 1300 • C that deviates from the main P-T-trend ( Figure 2a). This may be due to two reasons: (1) an artifact caused by underestimation of P NT00 at pressures above 6 GPa [64] because experimental data used for P NT00 calibration were only up to 6 GPa [24]; (2) possible inflection of the geotherm caused by a thermal event preceding kimberlite emplacement [16,65,66]. Our preferred interpretation is that this phenomenon reflects a heating event, which may have taken place just prior to entrainment in the kimberlite magma (details in Section 5.2).

Steady-State Paleogeotherm for the Mantle Beneath the KM Pipe
Regardless of the nature of the high-T cluster of points, we exclude P-T points above 1200 • C for calculation of the steady-state paleogeotherm for the lithospheric mantle beneath the KM pipe, i.e., we have calculated the geotherm that existed prior to the thermal perturbation that preceded kimberlite emplacement (Section 5.2). Choosing T = 1200 • C as the cut-off value not only excludes the cluster of high-T points but also eliminates samples with possibly underestimated pressures because P NT00 has not been calibrated to P > 6 GPa.
The thermal properties of the crust have been identified as a significant source of uncertainty when calculating cratonic paleogeotherms [7,69]. We attempted to use the best available local estimations for the structure of the crust and its heat production for our preferred model (model #1), we fixed the heat production rate in the upper and lower crust beneath the Upper Muna field at 0.76 and 0.076 µW/m 3 (Table 1) according to [35]. The thicknesses of the upper and lower crust were taken from [70] and they are 23 and 33 km, respectively ( Table 1). Fitting of the accepted P-T estimates (<1200 • C) yields a surface heat flow of 35.2 mW/m 2 and a lithosphere thickness of 229 km (Table 1, Figure 2a). The uncertainty of the geotherm calculation is plotted as thin lines (representing the 1-standard deviation misfit (Table 1)) parallel to the geotherm (Figure 2). The majority of the P-T points fall within this misfit.
temperature may indicate that clinopyroxene was not in equilibrium with orthopyroxene. We excluded one grain at this step.
Thus, according to the proposed filters [23], we accepted 188 clinopyroxene xenocrysts for TNT00 and PNT00 thermobarometry (Table S2). Judging by eye, the majority of accepted points lie along a narrow, elongated field in P-T space (Figure 2a). However, there is a cluster of high-T points above 1300 °C that deviates from the main P-T-trend (Figure 2a). This may be due to two reasons: (1) an artifact caused by underestimation of PNT00 at pressures above 6 GPa [64] because experimental data used for PNT00 calibration were only up to 6 GPa [24]; (2) possible inflection of the geotherm caused by a thermal event preceding kimberlite emplacement [16,65,66]. Our preferred interpretation is that this phenomenon reflects a heating event, which may have taken place just prior to entrainment in the kimberlite magma (details in Section 5.2).  Table 1). Diamond (D)-graphite (G) transition from [67].

Steady-State Paleogeotherm for the Mantle Beneath the KM Pipe
Regardless of the nature of the high-T cluster of points, we exclude P-T points above 1200 °C for calculation of the steady-state paleogeotherm for the lithospheric mantle beneath the KM pipe, i.e., we have calculated the geotherm that existed prior to the thermal perturbation that preceded kimberlite emplacement (Section 5.2). Choosing T = 1200 °C as the cut-off value not only excludes the Figure 2. P NT00 and T NT00 estimates for clinopyroxene xenocrysts from the Komsomolskaya-Magnitnaya pipe and four modeled mantle paleogeotherms calculated using the FITPLOT program (see Table 1). Diamond (D)-graphite (G) transition from [67]. McKenzie et al. [6], in their model for the Siberian craton (Udachnaya pipe), used higher values for the heat production rate both in the upper and lower crust-1.12 and 0.4 mW/m 2 . To demonstrate how the heat production rate in the crust affects geotherm parameters, we also calculated the paleogeotherm using these values (model #2, Table 1, Figure 2b). The resulting paleogeotherm requires a significantly higher surface heat flow (50.6 mW/m 2 ) and lithosphere thickness (253 km) than model #1 (Table 1).
We calculated two other paleogeotherm models (Table 1). Model #3 (Figure 2c) uses all 188 clinopyroxenes accepted for P NT00 and T NT00 . For model #4 (Figure 2c), we also included P-T estimates of 128 grains (avoiding seven samples showing crustal pressures) that were discredited by using the protocol (Section 3.1). Both models #3 and #4 yield surface heat flows and lithosphere thicknesses very similar to our preferred model #1, but with a higher misfit (Table 1). It is important to note that the vast majority of clinopyroxenes, filtered using the protocol, plot within misfit of model #4, whereas P-T points of clinopyroxene not recommended for P NT00 and T NT00 are significantly scattered (Figure 2d). This highlights the importance of the protocol proposed in [23] for the selection of appropriate clinopyroxene grains for single-clinopyroxene thermobarometry using P NT00 and T NT00 . The similarity of the geotherm parameters obtained in models #1 and #3 demonstrates that inclusion of the high-T samples (>1200 • C) does not affect the calculated steady-state geotherm and only increases the resulting misfit. This is because most of the P-T points lie below 1200 • C and are scattered along a pressure range 5-7 times wider than the high-T cluster. The consistency in output parameters of models #1 and #4 shows that, even if the criteria for sample selection for P NT00 and T NT00 are ignored, a large array of P-T data may constrain a reasonable paleogeotherm, but with greater uncertainties (Table 1).

Thermobarometry of Mantle Xenoliths from the Komsomolskaya-Magnitnaya Pipe
We used combinations of either the orthopyroxene-garnet barometer of Nickel and Green [71] (P NG85 ) with the clinopyroxene-orthopyroxene thermometer of Taylor [72] (T T98 ) for clinopyroxene-bearing peridotites, or P NG85 with the orthopyroxene-garnet thermometer of Nimis and Grütter [73] (T NG10 ) for clinopyroxene-free peridotites. These combinations have been recommended as the most robust for peridotites [73]. P-T estimates for the coarse peridotites are in very good agreement with the preferred clinopyroxene geotherm (model #1), whereas P-T points for the porphyroclastic peridotites are outside the misfit area and plot within the high-T cluster of clinopyroxene xenocrysts (Figure 3). The latter phenomenon, which can be related to the pre-kimberlite thermal event, will be discussed in Section 5.2.

Mantle Composition for the Komsomolskaya-Magnitnaya Pipe: Constraints from Clinopyroxene Xenocrysts
Here, we discuss the data set that includes only homogenous clinopyroxene xenocrysts with satisfactory microprobe analyses, i.e., 323 grains. The 188 (58%) clinopyroxene grains accepted for thermobarometry are assumed to be derived from orthopyroxene-bearing garnet peridotites. The 114 clinopyroxenes that were rejected also belong to the peridotite paragenesis according to the criterion Cr2O3 >0.5 wt% [63] (Figure 4). One hundred and nine (31%) of them plot in the field of garnet peridotites in the discriminant diagram ( Figure 4); among these, 109 grains show high Ca# >0.5 and might not have been in equilibrium with orthopyroxene. The remaining 5 (4%) grains are from spinel peridotites (Figure 4).

Mantle Composition for the Komsomolskaya-Magnitnaya Pipe: Constraints from Clinopyroxene Xenocrysts
Here, we discuss the data set that includes only homogenous clinopyroxene xenocrysts with satisfactory microprobe analyses, i.e., 323 grains. The 188 (58%) clinopyroxene grains accepted for thermobarometry are assumed to be derived from orthopyroxene-bearing garnet peridotites. The 114 clinopyroxenes that were rejected also belong to the peridotite paragenesis according to the criterion Cr 2 O 3 >0.5 wt% [63] (Figure 4). One hundred and nine (31%) of them plot in the field of garnet peridotites in the discriminant diagram ( Figure 4); among these, 109 grains show high Ca# > 0.5 and might not have been in equilibrium with orthopyroxene. The remaining 5 (4%) grains are from spinel peridotites (Figure 4). thermobarometry are assumed to be derived from orthopyroxene-bearing garnet peridotites. The 114 clinopyroxenes that were rejected also belong to the peridotite paragenesis according to the criterion Cr2O3 >0.5 wt% [63] (Figure 4). One hundred and nine (31%) of them plot in the field of garnet peridotites in the discriminant diagram ( Figure 4); among these, 109 grains show high Ca# >0.5 and might not have been in equilibrium with orthopyroxene. The remaining 5 (4%) grains are from spinel peridotites (Figure 4).  The other 21 clinopyroxene xenocrysts may belong to eclogite xenoliths or the megacryst suite or may represent kimberlitic phenocrysts. This non-peridotitic group comprises only 7%, and keeping in mind that sampling of the lithosphere column by the kimberlite magmas need not be representative [74], we can suggest that peridotites are the predominant rocks in the lithospheric mantle beneath the KM pipe. At least 5 grains from non-peridotitic clinopyroxenes, which are characterized by high Na 2 O, Al 2 O 3 , and low MgO, may be derived from eclogites ( Figure 5). The other 21 clinopyroxene xenocrysts may belong to eclogite xenoliths or the megacryst suite or may represent kimberlitic phenocrysts. This non-peridotitic group comprises only 7%, and keeping in mind that sampling of the lithosphere column by the kimberlite magmas need not be representative [74], we can suggest that peridotites are the predominant rocks in the lithospheric mantle beneath the KM pipe. At least 5 grains from non-peridotitic clinopyroxenes, which are characterized by high Na2O, Al2O3, and low MgO, may be derived from eclogites ( Figure 5).

Clinopyroxene Xenocrysts: Depth Distribution Profile for KM Pipe
The pressures estimated from PNT00 demonstrate that the clinopyroxene grains accepted for thermobarometry show a bimodal depth distribution with peaks at 130 and 180 km (Figure 6a). However, accepted clinopyroxenes represent only 58% of all grains. Although 113 clinopyroxenes equilibrated with orthopyroxene were not accepted for PNT00, their P-T conditions may be calculated by using a combination of TNT00 and the obtained geotherm (model #1) as an equation for the barometer. The depth distribution profile for all 221 pyroxenes that were in equilibrium with orthopyroxene is shown in Figure 6b and is similar to that in Figure 6a. This comprises 89% of the

Clinopyroxene Xenocrysts: Depth Distribution Profile for KM Pipe
The pressures estimated from P NT00 demonstrate that the clinopyroxene grains accepted for thermobarometry show a bimodal depth distribution with peaks at 130 and 180 km (Figure 6a). However, accepted clinopyroxenes represent only 58% of all grains. Although 113 clinopyroxenes equilibrated with orthopyroxene were not accepted for P NT00 , their P-T conditions may be calculated by using a combination of T NT00 and the obtained geotherm (model #1) as an equation for the barometer. The depth distribution profile for all 221 pyroxenes that were in equilibrium with orthopyroxene is shown in Figure 6b and is similar to that in Figure 6a. This comprises 89% of the grains and should yield a more representative depth distribution profile (Figure 6b).

The Mantle Paleogeotherm Beneath the Mir and Udachnaya Pipes
Input and output parameters for the Udachnaya and Mir geotherms are given in Table 1. We used original P-T values proposed by the authors for the mantle xenoliths and did not apply any filtering to the data. In the model #1 for both pipes, we used heat production rates in the upper and lower crust as in model #1 for the KM pipe. The Mir and Udachnaya pipes have similar thermal lithosphere thickness and surface heat flux (Figure 7, Table 1). In these models, the deepest (229 km) and the coldest lithosphere is observed beneath the KM pipe. However, the difference is not large and all three geotherms lie within their respective uncertainties. All paleogeotherms lie between the reference geotherms proposed by Hasterok and Chapman [77] with surface heat fluxes of 35 and 40 mW/m 2 .
For model #2, we used higher heat production rates in the upper and lower crust, and this yielded a higher thickness to slightly different degrees for all the pipes ( Table 1). The higher values of the input heat production rates in the crust drastically affect the resulting surface heat flux for the calculated geotherms (Table 1). For all three pipes, surface heat flux and lithosphere thickness increased by 12-15 mW/m 2 and 9-24 km in comparison to model #1. In Section 4.2, we proposed that the FITPLOT program fits the geotherms through the data array and is sensitive to the input crustal parameters. The higher values of the heat production in the crust (model #2) lead to a higher surface heat flux and a slightly oblique angle of the geotherm so it intersects the isentrope at greater depths (Section 5.1).  [2] (a) Only accepted for single crystal thermobarometry. P-T calculated by P NT00 and T NT00 . KM pipe-188 grains (58%). Novinka pipe-97 grains (56%) [23]. (b) All peridotitic xenocrysts assumed in equilibrium with orthopyroxene. P-T calculated by intersection of T NT00 with the paleogeotherm model #1. KM pipe-201 grains (89%). Novinka pipe-123 grains (71%) (data from [23]).

The Mantle Paleogeotherm Beneath the Mir and Udachnaya Pipes
Input and output parameters for the Udachnaya and Mir geotherms are given in Table 1. We used original P-T values proposed by the authors for the mantle xenoliths and did not apply any filtering to the data. In the model #1 for both pipes, we used heat production rates in the upper and lower crust as in model #1 for the KM pipe. The Mir and Udachnaya pipes have similar thermal lithosphere thickness and surface heat flux (Figure 7, Table 1). In these models, the deepest (229 km) and the coldest lithosphere is observed beneath the KM pipe. However, the difference is not large and all three geotherms lie within their respective uncertainties. All paleogeotherms lie between the reference geotherms proposed by Hasterok and Chapman [77] with surface heat fluxes of 35 and 40 mW/m 2 .
For model #2, we used higher heat production rates in the upper and lower crust, and this yielded a higher thickness to slightly different degrees for all the pipes ( Table 1). The higher values of the input heat production rates in the crust drastically affect the resulting surface heat flux for the calculated geotherms (Table 1). For all three pipes, surface heat flux and lithosphere thickness increased by 12-15 mW/m 2 and 9-24 km in comparison to model #1. In Section 4.2, we proposed that the FITPLOT program fits the geotherms through the data array and is sensitive to the input crustal parameters. The higher values of the heat production in the crust (model #2) lead to a higher surface heat flux and a slightly oblique angle of the geotherm so it intersects the isentrope at greater depths (Section 5.1).

A Mantle Paleogeotherm Beneath the Upper Muna Field
The common approach used to characterize the thermal state of the lithospheric mantle, used by petrologists for years, is a graphical comparison of the P-T data for mantle xenoliths and xenocrysts to the reference geotherms characterized in terms of surface heat flow [68] (hereafter PC77 geotherms). In this way, the best-fit geotherm is estimated qualitatively by eye. However, one of the main disadvantages of this method is that the resulting P-T data array often equally satisfies several reference geotherms. Other problems of this approach are discussed in detail in [7]. Currently, researchers both from scientific and diamond exploration communities use the FITPLOT program [10,13,15,32]. The FITPLOT program is based on an improved understanding of the thermal properties of the lithospheric mantle; it numerically fits obtained P-T data and allows the input of parameters for the crust estimated for the region of interest [7].
In order to compare our results with previous work, we nevertheless plotted our P-T estimates for clinopyroxene xenocrysts from the KM pipe against the reference PC77 geotherms (Figure 2d). Most of the KM samples lie between the 35 and 40 mW/m 2 PC77 geotherms and the high-T cluster plots at 45 mW/m 2 PC77 geotherm ( Figure 2) similar to clinopyroxene xenocrysts from the Novinka pipe [23]. P-T estimates made for the large dataset of various mantle xenocrysts and xenoliths from the KM, Novinka, Deimos, and Zapolarnay pipes in [29,30] are scattered between 35 and 40 mW/m 2 . Overall, a qualitative comparison concludes that our data are similar to previous results on the Upper Muna field. However, it is difficult to estimate the exact depth of the lithosphere if paleogeotherms are estimated in terms of the range of surface heat flux of the reference PC77 geotherms. It is possible to estimate only an interval for the thermal lithosphere-asthenosphere boundary in the case of the 40 and 45 mW/m 2 PC77 geotherms. Moreover, this is impossible in the case of PC77 geotherms with surface heat flux <40 mW/m 2 because their gradients increase dramatically at depth and they do not cross the isentrope at any point [7]. However, in these cases, applying data for the deepest xenoliths can be used as an alternative method to assess lithosphere thickness [7,22]. Thus, the deepest mantle xenolith from [29,30] yields a lithosphere thickness of ~250 km, whereas our deepest sample gives 210 km (± 20 km) for trimmed data and 235 km (± 25 km) for all data (uncertainty for PNT00 based on [64]). However, the deepest P-T estimates may represent an artifact rather than realistic values due to analytical error or/and uncertainties in the barometer and/or thermometer. Thus, the lithosphere

A Mantle Paleogeotherm Beneath the Upper Muna Field
The common approach used to characterize the thermal state of the lithospheric mantle, used by petrologists for years, is a graphical comparison of the P-T data for mantle xenoliths and xenocrysts to the reference geotherms characterized in terms of surface heat flow [68] (hereafter PC77 geotherms). In this way, the best-fit geotherm is estimated qualitatively by eye. However, one of the main disadvantages of this method is that the resulting P-T data array often equally satisfies several reference geotherms. Other problems of this approach are discussed in detail in [7]. Currently, researchers both from scientific and diamond exploration communities use the FITPLOT program [10,13,15,32]. The FITPLOT program is based on an improved understanding of the thermal properties of the lithospheric mantle; it numerically fits obtained P-T data and allows the input of parameters for the crust estimated for the region of interest [7].
In order to compare our results with previous work, we nevertheless plotted our P-T estimates for clinopyroxene xenocrysts from the KM pipe against the reference PC77 geotherms (Figure 2d). Most of the KM samples lie between the 35 and 40 mW/m 2 PC77 geotherms and the high-T cluster plots at 45 mW/m 2 PC77 geotherm ( Figure 2) similar to clinopyroxene xenocrysts from the Novinka pipe [23]. P-T estimates made for the large dataset of various mantle xenocrysts and xenoliths from the KM, Novinka, Deimos, and Zapolarnay pipes in [29,30] are scattered between 35 and 40 mW/m 2 . Overall, a qualitative comparison concludes that our data are similar to previous results on the Upper Muna field. However, it is difficult to estimate the exact depth of the lithosphere if paleogeotherms are estimated in terms of the range of surface heat flux of the reference PC77 geotherms. It is possible to estimate only an interval for the thermal lithosphere-asthenosphere boundary in the case of the 40 and 45 mW/m 2 PC77 geotherms. Moreover, this is impossible in the case of PC77 geotherms with surface heat flux <40 mW/m 2 because their gradients increase dramatically at depth and they do not cross the isentrope at any point [7]. However, in these cases, applying data for the deepest xenoliths can be used as an alternative method to assess lithosphere thickness [7,22]. Thus, the deepest mantle xenolith from [29,30] yields a lithosphere thickness of~250 km, whereas our deepest sample gives 210 km (±20 km) for trimmed data and 235 km (±25 km) for all data (uncertainty for P NT00 based on [64]). However, the deepest P-T estimates may represent an artifact rather than realistic values due to analytical error or/and uncertainties in the barometer and/or thermometer. Thus, the lithosphere thickness may be overestimated. On the contrary, if P-T data are limited or if the kimberlite magma did not sample the entire lithosphere column, it may be underestimated.
Griffin et al. [2] also used the reference PC77 geotherms, but they applied the algorithms proposed by Ryan et al. [12] to garnet xenocrysts from three pipes (Novinka, Zapolarnaya, and Zimnaya) to estimate the paleogeotherm beneath the Upper Muna field using the Ni-in-garnet (T R96 ) thermometer and Cr-in-garnet barometer (P R96 ). The P R96 barometer has a limitation, because reliable pressure estimates can only be done for Cr-saturated garnets that coexisted with spinel, otherwise only minimum pressures can be estimated [12]. Thus, the "garnet geotherm" was defined by the envelope of maximum pressure at each temperature (see Figure 4 in [2]), assuming those garnets coexisted with chrome spinel. The resultant paleogeotherm corresponds to the 38 mW/m 2 PC77 geotherm. However, above 1100 • C (~6 GPa), the locus of maximum P R96 no longer follows the 38 mW/m 2 geotherm and there are no garnets with P R96 > 6 GPa. This is probably because most high-temperature garnets are undersaturated in Cr, i.e., were not in equilibrium with chromite. The pressure of each garnet grain was derived by projecting its T R96 to the geotherm [2]. Griffin et al. [2] defined the lithosphere as depleted material with garnets containing <10 ppm Y, distinct from the convecting, presumably more "fertile" asthenosphere. Thus, the lithosphere-asthenosphere boundary was determined as the temperature limit of the distribution of Y-depleted garnets. This corresponds to a pressure of 6.5 GPa and a lithosphere thickness of~205 km. At temperatures above 1000 • C (P > 6.5 GPa), there are only Y-rich garnets. Griffin et al. [2] noted that high-temperature garnets have undepleted trace-element chemistry similar to the garnets of high-temperature porphyroclastic xenoliths, and, thus, the lithosphere-asthenosphere boundary defined by Griffin et al. [2] actually represents a "kink" in the lithospheric paleogeotherm (Section 5.2). Therefore, we suggest that the Griffin et al. [2] approach underestimated the lithosphere thickness beneath the Upper Muna field by ca 10%.
Ziberna et al. [23] constrained the mantle paleogeotherm beneath the Novinka pipe based on 173 clinopyroxene xenocrysts. The authors evaluated the paleogeotherm using the FITPLOT program indicating~34 mW/m 2 surface heat flow, a thermal lithosphere thickness of~225 km, and an over 100 km thick "diamond window" beneath the pipe at the time of kimberlite eruption.
Comparison of our models #1 and #2 geotherms for the KM pipe constrained by the FITPLOT program demonstrates that assumed heat production rates in the crust have drastic effects on the calculated paleogeotherm and lithosphere thickness. Heat production rates used in [6] were calculated based on the data from uplifted granulite terranes in Canada, whereas estimations in [35] are based on the upper and lower crustal xenoliths from Siberian kimberlites. We think that it is important to use the best available local estimates for the crustal parameters.
Thermal modeling has been used to propose that the present-day thermal lithosphere thickness varies from 300-350 km in the central part of the Siberian craton, where extremely low heat flow values (18-25 mW/m 2 ) are measured [69,78]. Thermo-petrological interpretations of seismic Vp-velocity model along the Peaceful Nuclear Explosions profiles confirm a 300 km [79] or 220-250 km [80] thick lithospheric mantle in this part of the Craton. Our model #2 gives a lithospheric thickness of more than 250 km and an extremely high surface heat flow (50.6 mW/m 2 ). Such inconsistency makes model #2 meaningless and further confirms that model #1 yields the preferred steady-state mantle paleogeotherm for the KM pipe. Model #1 excludes high-T points that may reflect a "kink" in the geotherm that was probably formed shortly before kimberlite emplacement (Section 5.2). If this geotherm inflection took place, then our model #1 paleogeotherm characterizes the thermal state of the lithospheric mantle beneath the Upper Muna field prior to kimberlite magmatic activity. Ziberna et al. [23] also did not include high-T points for their paleogeotherm calculation for the Novinka pipe and obtained similar results. Combining our model #1 for the KM pipe and results for the Novinka pipe [23] we can accept a~34-35 mW/m 2 surface heat flux and~225-230 km lithospheric thickness for the Upper Muna field. The consistency of P-T data for the coarse peridotites with the clinopyroxene geotherm model #1 (Figure 3) testifies that only a large P-T data set for clinopyroxene xenocrysts can be successfully used to robustly constrain the steady-state paleogeotherm. This confirms the conclusions of Mather et al. [7] that the carefully filtered xenocryst P-T data yield a paleogeotherm almost identical to that produced from the well-equilibrated xenoliths.
The differences in the present-day lithospheric thickness and surface heat flow [69,[78][79][80][81] and those obtained from the paleogeotherm in this study indicate thickening of the continental lithosphere with age in the central part of the Siberian craton. The progressive thickening of the continental lithosphere is typical for Archean cratonic mantle [69].

The "Kink" Problem of the Paleogeotherm in Lithosphere Mantle
In general, two textural types of peridotites xenoliths from kimberlites are defined: coarse and porphyroclastic [4,61,66]. Coarse peridotites are derived from the upper parts of the lithospheric mantle and their P-T estimates usually lie along a reference PC77 geotherm. Porphyroclastic peridotites are derived from the lower part of the lithosphere and have a significantly higher temperature range within a narrow range of pressures, demonstrating an inflection from the P-T array defined by granular peridotites. This "kink" (or in some cases a "step") in P-T trends is found in xenolith suites from many kimberlite localities worldwide. In the case of the KM pipe, also there is the cluster of high-T clinopyroxene xenocrysts located within a narrow pressure range ( Figure 3) and distinct from the calculated steady-state geotherm (model #1). Thus, these clinopyroxenes may represent fragments of porhyroclastic peridotites. It is thought that the deviation formed by porhyroclastic peridotites from the steady-state geotherm defined by coarse peridotites identifies a thermal event that affected the base of lithosphere shortly before or at the time of kimberlite emplacement, e.g., a mantle plume that may cause kimberlite magmatism [65,66,[82][83][84][85][86].
However, the kinked (or stepped) geotherm has been the subject of much debate, centered on whether it is an artifact of thermometers and barometers [65]. Most of the used thermometers and barometers are calibrated at pressures below 6 GPa, i.e., they do not cover the pressure range of porphyroclastic peridotites. Here, we studied two porphyroclastic peridotites from the KM pipe. Estimated P-T parameters for these xenoliths using combinations of different thermometers and barometers, recommended in [73], show scattering within a pressure range up to 1 GPa (Figure 3). Most of the estimates plot within the high-T cluster of clinopyroxene xenocrysts. It is remarkable that P and T estimated by P NT00 and T NT00 using clinopyroxenes from these xenoliths plots within the cluster formed by clinopyroxene xenocrysts (Figure 3). Recently, Brey et al. [87] recalibrated a barometer based on Al partitioning between orthopyroxene and garnet (Al-in-orthopyroxene barometer, P BBG08 ) at P > 6 and showed that P BBG08 adequately reproduces experimental pressures up to 8 GPa. P BBG08 yields pressures for these porphyroclastic peridotites somewhat lower than other barometers and enhances the "kink" effect ( Figure 3). This effect was also reported for other kimberlite localities [65]. Thus, pressures obtained for porphyroclastic peridotites using P BBG08 clearly show deviation from the steady-state geotherm and indicate the possibility of the existence of kinked geotherms.
Regardless of whether kinked geotherms are real or not, our model #1 constrains the steady-state geotherm before possible "kink" formation preceding kimberlite magmatism. Even if we include high-T points, which possibly reflect a "kink", their effect on the calculated parameters for the steady-state paleogeotherm is negligible.

Comparison of Lithospheric Thickness Beneath the Upper Muna, Daldyn, and Mirny Fields
The Udachnaya and Mir pipes are the most thoroughly characterized kimberlite localities within the Daldyn and Mirny fields. There are several studies with calculated paleogeotherms beneath the Daldyn field [6,10,30,88,89]. McKenzie et al. [6] calculated a paleogeotherm for the Udachnaya pipe using the FITPLOT program. The input parameters for crustal heat production rates as in our model #2. However, the lithosphere thickness (242 km) and surface heat flow (58.6 mW/m 2 ) obtained by McKenzie et al. [6] are significantly higher than in our model #2. This is because McKenzie et al. [6] used a much higher value for the thickness of the upper crust than predicted for the Daldyn field [70], and as a result, the overall heat production for the crust was larger than in our model #2. We prefer our model #1 for the paleogeotherm beneath the Udachnaya pipe as well as for the Mir and KM pipes because the input parameters are estimated for given kimberlite localities in the Siberian craton [35,70]. The paleogeotherm for the mantle beneath the Mirny field was calculated using the FITPLOT program for the first time.
Paleogeotherms constrained from garnets xenocrysts yielded different lithosphere thicknesses beneath the studied kimberlite fields: Daldyn (230 km) > Upper Muna (205 km) > Mirny (190 km) [2]. In contrast, both our models #1 and #2 demonstrate the similar lithosphere thickness for the Upper Muna, Daldyn, and Mirny fields. The "Garnet geotherms" may probably underestimate the lithospheric thickness and the Mirny and Upper Muna fields were studied based on a limited number of grains [2]. Our method yields more accurate results because it is based on the precise fitting of the large data suite and confirms that the lithospheric thickness at the time of kimberlite magmatism did not vary strongly across the central part of the Siberian craton.

Composition and Stratigraphy of Lithospheric Mantle Beneath the Upper Muna Field
Our results demonstrate that peridotitic clinopyroxenes are predominant (89%) among the studied xenocrysts from the KM pipe, suggesting that the sampled mantle is mainly composed of peridotites. Peridotitic clinopyroxene xenocrysts are also prevalent in the Novinka pipe [23]. Rare data on mantle xenoliths confirm that peridotites are the main rock types in the lithospheric mantle beneath the Upper Muna field, whereas eclogites are rare [90].
The depth distribution of xenocrysts of peridotitic clinopyroxene (Figure 6b) may be interpreted in two ways: (i) as the distribution of clinopyroxene-bearing peridotites in the mantle column, or (ii) as preferential sampling of rocks from different depths by the kimberlite magma. These possibilities are not mutually exclusive. The depth distribution profiles for the KM and Novinka pipes look slightly different (Figure 6b). If these two nearby pipes sampled the same lithospheric mantle, then the difference may be caused by the second scenario. Both profiles show a bimodal distribution for both pipes with the region of lower clinopyroxene concentration at 140-180 km between two peaks (Figure 6b). This may be evidence for a layer of clinopyroxene-poor or -free peridotitic rocks such as harzburgites and dunites beneath the Upper Muna field. This observation is consistent with garnet xenocryst data, which show that most of G10 (harzburgitic/dunitic) garnets from the Upper Muna kimberlites are also concentrated over this depth interval [2], as shown in Figure 6.

Implication for Diamond Potential of the Siberian Kimberlites
Both the KM and Novinka pipes are characterized by lower diamond grades (0.5 and 0.6, in carats/ton [91]) than the Udachnaya (2.1 [11]) and Mir (3.5 [91]) pipes. The diamond grade of kimberlites depends on both lithospheric mantle characteristics and processes occurring during kimberlite emplacement. The diamond potential of the lithospheric mantle depends on (i) its thickness and thermal state, which control the thickness of its volume located in the diamond P-T stability field, the so-called "diamond window", and (ii) balance between diamond-friendly and diamond-unfriendly events during its pre-kimberlite history. The "diamond window" is the area between the point at which the graphite-diamond transition intersects the paleogeotherm and the lithosphere-asthenosphere boundary (e.g., Figure 7). Many factors that take place during emplacement may affect the diamond potential of kimberlites, e.g., (i) rate of magma ascent and possible halts of rising magma at shallow depths, allowing resorption of diamond by the kimberlite magma, (ii) redox condition of the kimberlite melt, (iii) selective sampling of mantle rocks at different depths, (iv) dilution by country rock. Our model #1 demonstrates a thicker lithosphere and diamond window for the KM pipe than for the Udachnaya and Mir pipes (Figure 7). However, if the uncertainties are considered, this difference is not significant, and we can conclude that the lithospheric thickness for these three fields is similar. Therefore, the thermal state of the mantle is not the cause of the differences in diamond grade of kimberlites from the Upper Muna, Daldyn, and Mir pipes. Perhaps other factors described above cause differences in diamond grades. This requires further studies of both mantle samples and kimberlites of the Upper Muna field, and the KM pipe in particular.

•
The preferred steady-state mantle paleogeotherm was constrained by P-T estimates of selected clinopyroxene xenocrysts from the Komsomolskaya-Magnitnaya pipe, i.e., that (i) passed through the protocol proposed by Ziberna et al. [23], and (ii) showed T < 1200 • C. Clinopyroxenes with T > 1200 • C were excluded because they may represent a "kink" caused by a thermal event associated with kimberlite magmatism. • P-T estimates for the coarse peridotites from the KM pipe are consistent with the clinopyroxene geotherm, testifying that large P-T data sets for clinopyroxene xenocrysts alone can be successfully used to constrain the robust steady-state paleogeotherm.

•
Fitting of P-T data for clinopyroxenes from the KM pipe using the FITPLOT program yielded a mantle paleogeotherm similar to the nearby Novinka pipe [23]. Thus, we accept a~34-35 mW/m 2 surface heat flux, 225-230 km lithospheric thickness, and 110-120 thick "diamond window" for the Upper Muna field.

•
Peridotitic clinopyroxenes are predominant (89%) among xenocrysts studied from the KM pipe, testifying that mantle beneath the Upper Muna field is mainly composed of peridotites. This is consistent with previous studies.

•
The estimated lithosphere thickness beneath Upper Muna is slightly larger than for the Daldyn (218 km) and Mirny (210 km) fields. However, considering their mutual uncertainties, the lithosphere thicknesses for these three fields are similar.