Radiogenic heat production provides a thermal threshold for Archean cratonization process

The processes by which cratonization of Archean crust occurred remain enigmatic (e.g., Cawood et al., 2013; Vanderhaeghe et al., 2019; Hawkesworth and Jaupart, 2021). A potentially undervalued role is that of radiogenic crustal heat production, particularly in priming the mid- and deep crust for large-volume partial melting and the ultimate transfer to the shallower crust. In the Yilgarn craton (Australia), cratonization was a prolonged process, with nearly continuous magmatism between 2830 and 2600 Ma. Final cratonization was characterized by a widespread late “bloom” of high-K granitic magmatism that peaked between 2660 and 2640 Ma with very high radiogenic heat production (interquartile range of 6.6–11.8 µW/m³). Such prolonged magmatism, and particularly the late granite bloom, has been inferred to have been mantle driven (e.g., Bodorkos and Sandiford, 2006; Mole et al., 2019), although it largely occurred over a time frame with negligible mantle input. However, previous studies likely underestimated the magnitude of radiogenic heat production in the Yilgarn craton and therefore ruled it out as being a major contributor to cratonization, instead favoring external drivers (e.g., Bodorkos and Sandiford, 2006). A comprehensive new geochemistry data set from the Yilgarn craton (Smithies et al., 2023, 2024) allows for a thorough assessment of the role of radiogenic heat production in facilitating magmatism and long-term stability, without the need for external (mantle) drivers. In this study, we aimed to address whether the radiogenic heat production caused a thermal regime in the mid- and deep crust near, at, or above the solidus, setting up a “threshold” scenario primed for protracted large-volume melting, followed by long-term cooling. To achieve this, we utilized large-scale data sets collected by the Geological Survey of Western Australia, including geochemistry and robust metamorphic constraints, integrated with a detailed stratigraphic and event framework (e.g., Mole et al., 2019; Ivanic et al., 2022) and crustal architecture, to perform one-dimensional (1-D) thermal modeling of the Archean evolution of the Yilgarn craton. More generally, protracted magmatism and late granite “blooms” recognized in other Archean cratons (e.g., Superior Province, Slave craton, Pilbara craton) may reflect the signature of the role of radiogenic heat production in cratonization.

The Yilgarn craton is dominated by Archean granite-greenstone terranes formed between ca. 3100 and 2600 Ma, including the Youanmi, Narryer, and South West terranes of the West Yilgarn craton, and the Eastern Goldfields superterrane made up of the Kalgoorlie, Kurnalpi, Burtville, and Yamarna terranes (Fig. 1). Although cratonwide magmatism, including the 2660–2640 Ma late granite bloom, affected all the terranes from ca. 2720 to 2600 Ma, the Eastern Goldfields superterrane is younger and more juvenile than the West Yilgarn craton (e.g., Champion and Cassidy, 2007; Wyche et al., 2012; Mole et al., 2019).

An extensive, high-quality geochemical data set of 5545 unaltered felsic igneous rocks (Smithies et al., 2023, 2024; Table S1¹) includes 522 dated (U-Pb zircon) samples, >90% of which have crystallization ages between ca. 2832 and 2600 Ma (Fig. S2). These can be divided (see Supplemental Material) into groups broadly reflecting:

1. Melting of mafic crust at deep-crustal (high-Ca, high-Sr/Y granites) and mid-crustal (high-Ca, low-Sr/Y granites) depths, equivalent to the sodic Archean tonalite-trondhjemite-granodiorite (TTG) series;

2. Melting, at a range of crustal depths, of variably reworked crustal sources (low-Ca granites), equivalent to typical Archean high-K granites but including a charnockitic variant (high-Ca [high-Ti] granites); and

3. Melting of a dominantly lithospheric mantle source (diorites, high field strength element granites, sanukitoids, syenites).

Additional geochemical data sets include lithologies that make up greenstone belts in the Youanmi terrane (Table S2; including banded iron formation, n = 52; felsic volcanic rocks, n = 419; mafic rocks, n = 1835; siliciclastic metasedimentary rocks, n = 219; ultramafic rocks, n = 329), and mafic igneous rocks (gabbro, n = 30; metagabbro, n = 36) from the Narryer terrane (Table S3), as a proxy for the mafic lower crust.

The West Yilgarn metamorphic data set comprises 30 samples (Fig. 1; Table S4). The methodology is outlined in the Supplemental Material. There are 14 samples with metamorphic ages ranging from ca. 2730 to 2622 Ma (Fig. 1) and 16 undated samples inferred to record Neoarchean metamorphism. The quantitative pressure-temperature (P-T) data reflect peak Neoarchean metamorphism at amphibolite- to granulite-facies conditions, consistent with previous work (e.g., Gee et al., 1981; Goscombe et al., 2019). The median peak metamorphic T ranges from 558 °C to 890 °C (Fig. S1A), and the median peak P ranges from 2.85 to 8.05 kbar (two outliers excluded; Fig. S1B). The median peak thermobarometric ratio ranges from 68 to 225 °C/kbar (Fig. 1). This data set provides P-T constraints through which model thermal gradients must pass.

Calculated heat production for the key lithological groups used in the 1-D models is provided in Figure 2, at present day (A₀) and at 2640 Ma (A₂₆₄₀). The TTG group has median A₀ = 1.80 µW/m³, and high-K granite groups have median A₀ = 5.15 µW/m³, which is significantly higher than both the global average granite (= 2.5 µW/m³; Haenel et al., 1988) and the mean felsic value previously reported for the Yilgarn craton (= 1.88 µW/m³; Hawkesworth and Jaupart, 2021; Fig. 2A). The heat production map (Fig. 3A) shows that vast tracts of the exposed surface have A₀, which is significantly higher than global average granite. High heat production at the surface coupled with low to moderate surface heat flow (Fig. 3A) imply that modern-day heat production at depth must be lower, consistent with extraction of heat-producing elements (HPEs; K, Th, U) from deeper sources to produce these upper-crustal radiogenic granites. Median A₂₆₄₀ for the felsic groups ranges from 3.5 to 8.8 µW/m³ (Fig. 2B), and the heat production map for 2640 Ma (Fig. 3B) shows that the lithologies presently exposed would have had considerable heat production. Approximately 89% of the dated granites have calculated heat production at the age of crystallization (A_Cryst) > 2.5 µW/m³ (Fig. S2), with a gradual increase until the highest A_Cryst values in the high-K granite (= low-Ca granite groups) at 2680–2600 Ma.

To investigate the effects of radiogenic heat production on the Archean thermal evolution of the Yilgarn crust, 1-D numerical models were conducted. Model setups are shown in Figure 4A and explained in the Supplemental Material. For most lithologies, heat production was derived from age-corrected values calculated from the geochemical data set. However, for melt products derived from known sources, a proportion of the vertically integrated heat production was removed from each of the sources. Steady-state models were run at present day and at three different time intervals, broadly capturing key stages when the crustal column had significant changes to its configuration: 2900–2760 Ma (model 1), 2760–2650 Ma (model 2), and 2650–2600 Ma (model 3). The predicted thermal gradients are intended to represent the end of each model, although the absolute age at which the system reaches steady state could be different, especially for model 3. There are two results for model 3: one based only on heat production as a thermal source (solid red line, Fig. 4B), and one that also includes the advective thermal effects of high-K granite emplacement (dashed red line, Fig. 4B). Other thermal perturbations, such as mafic underplating or intrusion, are not considered in any of the models, nor is the thermal dependence on conductivity considered; thus, the predicted thermal gradients can be treated as minimum values.

The models attempt to reconstruct a simple but realistic crustal evolution based on the existing stratigraphic and magmatic framework, with the sampled rocks providing constraints on crustal compositions. There are three important premises in the model setups: First, radiogenic heat production of many of the crustal components is well above global averages (Fig. 2). Second, the models are internally consistent, apart from the transitional TTG layer in the present-day setup that has an increase in A₀ relative to model 3 (Fig. 4A). Possible explanations are addressed in the Supplemental Material. Third, the high-K granitic rocks that characterize the late granite bloom have the highest radiogenic heat production and must have played a fundamental role in cratonization. Sandiford et al. (2001) and Sandiford and McLaren (2006) recognized that the thickness of high-heat-producing layers is a critical factor in the magnitude of their contribution to crustal heat flow and the consequent thermal regime, with thicker layers having greater influence. The 3-km present-day thickness of these radiogenic granites represents a minimum due to erosion, and therefore a minimum contribution to the total crustal heat production. Significantly, when extrapolating back in time, and to source regions for these granites, this thick layer implies that the Yilgarn crust contained elevated rates of heat production at depth prior to the late granite bloom. This agrees with recent work by Hartnady et al. (2024), which showed that fractionation of Th and U in the Yilgarn craton accompanied crust formation and differentiation as early as 3.3 Ga, suggesting that nearly the entire budget of these elements was part of the crustal profile by the Neoarchean, and that high-temperature metamorphism and crustal anatexis further redistributed HPEs.

Numerous significant outcomes arise from the minimum thermal gradients in the 1-D models (Fig. 4). First, the ca. 2760 (model 1) and ca. 2650 Ma (model 2) thermal gradients are essentially the same (Fig. 4B), which depict elevated thermal conditions that intersect the tonalite wet solidus at ~20 km depth, and both dehydration curves at 25 km (Fig. 4C). This implies that for at least 100 m.y., the thermal regime of the mid- to deep crust was near to, at, or above conditions for partial melting. This protracted magmatism is evident in the geological record (e.g., Fig. S2), and it suggests that prolonged melting did not sufficiently redistribute HPEs to reduce the thermal gradients until the next major crustal reconfiguration (model 3). This invites the possibility that the final granite bloom of high-K granites in model 3 (at ca. 2650–2600 Ma) occurred because the crust was thermally primed and subsequently crossed the tipping point. Second, after 2650 Ma (model 3), the crust was characterized by a cooling thermal gradient (Fig. 4B) due to the redistribution of significant radiogenic heat production at depths that were sufficiently shallow to lessen their thermal impact. The Nd isotope record for high-K granites in particular shows that older sources are more refractory than younger sources (Smithies et al., 2023), strongly suggesting that the deeper crust became progressively more refractory with time, as the thermal gradient progressively cooled. These factors would make it more difficult for the crust to melt, ultimately resulting in the cessation of melting and final cratonization by ca. 2600 Ma. Third, compared to the thermal gradients from Bodorkos and Sandiford (2006) for the Eastern Goldfields superterrane, there is a significant difference in thermal state (Fig. 4B). For example, in our models, the thermal gradient intersects the tonalite wet solidus at 20 km (at 2760–2650 Ma), whereas in Bodorkos and Sandiford (2006), it is 300 °C cooler at ca. 2715 Ma and 150 °C cooler at 2650 Ma, far below the solidus. Therefore, we consider that Bodorkos and Sandiford (2006) underestimated the contribution of radiogenic heat production to promoting large-volume melting and chemical differentiation of Yilgarn craton crust. Our revised finding highlights the value of the more comprehensive data sets available to us that have revealed the greater magnitude of an internal driver for cratonization, without the necessity for external thermal inputs.

The predicted thermal gradients at 2760–2650 Ma capture many of the metamorphic data points over that interval (Fig. 4B). However, the metamorphic record at 2641–2622 Ma appears, at first glance, to be at odds with the 1-D models. There is no realistic scenario for the magnitude and vertical distribution of radiogenic heat production that can explain this part of the metamorphic data set. However, as the metamorphic age data overlap with the ages of late “bloom” granites, and indeed many of the sampled metamorphic rocks are from among such granites, the metamorphic data may instead be explained by a thermal regime involving (significant) advective heat (dashed red line, Fig. 4B). That is, these granites caused an advective thermal regime superimposed on that caused by radiogenic heating alone and resulted in regional contact metamorphism, with Neoarchean thermobarometric ratios that lie at the extremes of the high-T/P regime of Brown and Johnson (2018).

Our results show that the thick layer of high-K, highly radiogenic granitic rocks currently exposed at the surface of the Yilgarn craton has significant implications for the thermal character and evolution of the craton. When integrated back in time to source rock types, radiogenic heat production was elevated in large parts of the Archean crustal column for protracted periods of time. This set up an internally driven thermal regime whereby the mid- and deep crust was at temperatures near or at the solidus from at least 2760 Ma, primed for (large-volume) melting, without the need to invoke external drivers for magmatism. This protracted melting and upward transfer of magma and heat-producing elements caused the crust to become more refractory over time, requiring progressively higher temperatures to continue melting and driving melting to shallower sources. This would have been a self-fulfilling process toward cratonization after 2650 Ma because upward migration of heat production combined with increasingly refractory deeper crust would have resulted in cooler and stronger crust (e.g., Sandiford and McLaren, 2006). Our models show that radiogenic heat internally primed the crust for cratonization—neither mantle heat nor other specific tectonic events were needed—and may account for the Archean trend from sodic to potassic granites preserved in the global record (e.g., Laurent et al., 2014).

We are grateful to two anonymous reviewers and C. Kinney for constructive reviews on an earlier version of the manuscript and to T. Rushmer for editorial handling. The authors publish with the permission of the executive director, Geological Survey of Western Australia, and acknowledge funding from the Government of Western Australia Exploration Incentive Scheme.