1 Introduction
Most if not all continents are cored by Archaean cratonic nuclei surrounded by Proterozoic belts (Artemieva and Mooney, 2001; Bleeker, 2003, Fig. 1). These nuclei comprise greenstone belts delineating crustal-scale domes cored by magmatic/migmatitic rocks of felsic composition (Anhaeusser et al., 1969; Bouhallier et al., 1995; Collins et al., 1998; Condie, 1981, 1994; De Ronde and de Wit, 1994; De Wit, 1998; Hickman, 1984; Kröner, 1985:Van Kranendonk et al., 2007a). The mafic/ultramafic complexes forming the greenstone belts, including komatiites, imply partial melting in a mantle with a ambient temperature 200–300 °C higher than the present-day one, which is consistent with the 1600 °C mantle adiabat shown in Fig. 2 (Herzberg et al., 2010; Pollack, 1986). The widespread occurrence of high-temperature metamorphic mineral assemblages and of migmatites attests to a high crustal geothermal gradient during the Archaean (Brown, 2008; Brown and Johnson, 2018, Fig. 2). The magmatic record of the oldest cratonic nuclei, exemplified by the Pilbara carton and the Barberton terrane of the Kaapvaal craton, respectively in Australia and South Africa, typically spans over one billion years (Nd model ages ranging from 4.0 to 3.0 Ga calculated according to DePaolo (1981) using data from the Georoc database and references therein) and U–Pb ages on zircon ranging from 3.7 to 2.7 Ga (Georoc database and references therein and Zeh et al., 2007, 2009, 2011; Fig. 3). The Pilabara and Barberton are probably the most emblematic and best studied Archaean terranes and might have been part of a single Archaean craton (Bleeker, 2003; Cheney, 1996) representing a surface of more than 500 000 km2, which is similar to the Altiplano. This peculiar prolonged record may reflect a succession of continental growth/reworking events associated with lithospheric scale dynamics (Bédard, 2006; Bleeker and Hall, 2007; Condie, 2005; Hickman and Van Kranendonk, 2012; O'Neil et al., 2013; Percival et al., 2001; Zeh et al., 2009). In contrast, geophysical data combined with petrologic analysis of xenoliths from the subcontinental mantle suggest the formation of a depleted and buoyant lithospheric mantle root, allowing for the stabilization of these cratons at the early stage of the Archaean crust accretion (Griffin et al., 2003; Jordan, 1978; Pollack, 1986; Poupinet et al., 2003).

Schematic geologic map showing the distribution of Precambrian cratons composed of Archaean nuclei in red and Proterozoic belts in pink (modified from Artemieva and Mooney, 2001; Bleeker, 2003).

Geological and experimental constraints on the Archaean geothermal gradient. The red dots correspond to P–T conditions in the crust recorded by Archaean metamorphic rocks extracted from the dataset compiled in Brown (2008) updated in Brown and Johnson (2018) and references therein. The green line corresponds to the 1600 °C adiabat proposed for the Archaean mantle geotherm based on the conditions required for the generation of komatiites (Herzberg et al., 2010; Pollack, 1986). The pink shaded area corresponds to the positions of likely Archaean continental geotherms in consistence with these constraints.

Geochronological record of the Pilbara craton in Australia (bottom) and of the Barberton terrane part of the Kaavaal craton in South Africa (top) extracted from the GEOROC database with additional data from Zeh et al., 2007, 2009, 2011. The Pilbara and Barberton probably corresponded to a single Archaean cratonic nuclei (Cheney, 1996) representing a surface of more than 500 000 km2. Nd model ages (TDM) are calculated according to DePaolo (1981). Nd model ages suggest a continuous transfer of magmas from the mantle to the crust from the Hadean to the MesoArchaean. U–Pb ages point to a continuous record of magmatic differentiation and/or crustal reworking for close to 1 Ga from the EoArchaean to the transition with the Proterozoic.
These features feed an exciting ongoing debate regarding the early Earth dynamics and the processes leading to continental crust formation and evolution (Stern, 2008, Fig. 4). A wide range of models has been proposed regarding the nature of the ill-defined “primitive crust” with proposed compositions ranging from mafic to felsic (Albarède, 1998; Condie, 1998; Hawkesworth et al., 2017; Kemp et al., 2010; Stein and Hofmann, 1994; Wilde et al., 2001). Most authors agree that about 50% of the present-day continental crust had been extracted by 3.5 Ga from the mantle (Belousova et al., 2010; Dhuime et al., 2012). Some authors invoke a specific dynamical regime for Archaean time dominated by upwelling of mantle plumes responsible for the emplacement of large igneous complexes, later reworked to form the early continental crust (Bédard, 2006, 2018; Campbell and Hill, 1988; Condie, 2005; Stein and Hofmann, 1994; Van Kranendonk et al., 2007b; Zegers and van Keken, 2001). Others propose that, even in a hotter Earth, Archaean crust formation and evolution were related to plate tectonics in a way similar to the present-day situation (Arculus, 1999; Bleeker, 2003; Martin, 1986; Moyen and van Hunen, 2012; Rapp et al., 1991; Smithies et al., 2003). The widespread long-lived magmatic and tectonic activity typically recorded by Archaean cratonic nuclei, with concordant U–Pb zircon ages spanning for tens to hundreds of My, is attributed either (i) to a succession of magmatic-tectonic events associated with recurrent mantle plumes (Cawood et al., 2006; Condie, 2005; Van Kranendonk et al., 2007b), best exemplified with the Superior Province (Bédard, 2006) or (ii) to repeated subduction of slabs with a shallow dip to account for the widespread magmatic activity over large regions (Cawood et al., 2006; Percival et al., 2001) as illustrated by the Slave Province (Bleeker and Hall, 2007). More recently, Harris and Bédard (2014) and Bédard (2018) revived the old idea of continental drift in a model of periodically destabilized stagnant lid crust in which magmatism and deformation of cratonic nuclei are related to nascent continental collision controlled by the motion of continental blocks dragged by mantle flow. The investigation of this scenario by a 2D coupled petrological–thermomechanical tectono-magmatic numerical model indicates that it is prone to account for the petrogenesis of the continental crust during the Eo- and MesoArchaean (Sizova et al., 2015).

Schematic representations of proposed models for Archaean geodynamics (modified from Stern, 2008). Left: Stagnant lithospheric lid with accretion of oceanic crust above mantle plumes, some of which lead to the formation of oceanic plateaux that differentiate into continental nuclei. Gravitational instabilities insure recycling of the base of the lithospheric mantle and of the mafic lower crust. Right: Nascent plate tectonics with (i) accretion of oceanic crust along mid-oceanic ridges, (ii) accretion of continental crust through formation of magmatic arcs along subduction zones, and (iii) recycling of continental crust along convergent plate boundaries.
In any case, the early Earth was potentially characterized by a much higher geothermal gradient owing to the initial heat produced by the Earth accretion-differentiation and to a more intense decay of heat producing isotopes (HPI). Both of these heat sources have decreased significantly since the formation of the Earth (Bédard and Birch, 1965; Jaupart and Mareschal, 1999; Mareschal and Jaupart, 2013; Schubert et al., 1980; Stern, 2008; Turcotte, 1980; Wasserburg et al., 1964). In this contribution, we use the present-day thermal and lithologic structure of Archaean cratons provided by a synthesis of geological and geophysical data (Artemieva and Mooney, 2001; Mooney et al., 1998) to define a generic Archaean craton and calculate its thermal evolution in one dimension considering solely the effect of secular exponential decrease of HPI (238U, 235U, 232Th, and 40K) with time. A sensitivity analysis accounting for the variable range of relevant parameters (present crustal thickness, lithosphere thickness and surface heat flow) was achieved to provide plausible bounds of the crustal temperature profile during the Archaean. We then place these results in the context of geological observations and discuss their implications in terms of Archaean dynamics and crustal differentiation.
2 Structure and composition of Archaean cratonic nuclei
By definition, cratons are segments of continental crust that did not experience any major tectonic event since the end of the Precambrian. In the following, we explore available data on the present-day composition, structure and thermal characteristic of the oldest parts of Archaean cratons (Table 1) in order to define generic Archaean cratonic nuclei. The present-day thickness of the Archaean cratonic crust typically ranges from 30 to 45 km with an average of 40 km (Mooney et al., 1998; Schubert et al., 1980). Thinner cratonic crustal sections tend to be close to continent/ocean transitions along passive margins such as the Dharwar craton in South India. In these regions, the cratons margins have probably been affected by Phanerozoic lithospheric scale deformation and the present-day structure is not representative of the initial thickness of Archaean cratonic nuclei. Geological mapping of these nuclei indicate that they are typically structured in crustal-scale domes cored by felsic rocks (migmatites also designated by the term “grey gneisses” and plutonic rocks forming the so-called TTG suite) surrounded by greenstone belts formed of metasedimentary rocks alternating with mafic to ultramafic metavolcanics (Bédard, 2006; Hickman and Van Kranendonk, 2012; Moyen and Martin, 2013; Van Kranendonk et al., 2007b). Greenstones and felsic migmatites have been affected by high-temperature/low-pressure metamorphism ranging from greenschist to granulite facies (Brown, 2008). The solidus of these type of rocks (amphibolites, TTG gneisses) ranges from 700 to 900 °C (see compilation in Moyen and Martin, 2013). The presence of these rocks at the present-day surface implies that the initial Archaean crust was 5–15 km thicker than currently observed. We consider a thickness of 40 km for the present-day generic Archaean crust, which is likely a lower bound, leading to a conservative (lower bound) estimate of the total heat production within the crust. Rare exposed sections of the mid to lower Archaean cratonic crust and xenoliths point to a predominance of felsic to intermediate rocks, with an average granodioritic composition associated with variable proportions of mafic rocks (Percival and West, 1994; Rudnick et al., 1998; Rudnick and Gao, 2003; Wedepohl, 1995). This lithological assemblage is corroborated by geophysical data (Mooney et al., 1998). For the sake of simplicity and given uncertainties regarding the composition of the Archaean lower crust, we consider a uniform distribution of heat producing isotopes with depth throughout the crust and we neglect heat production in the lithospheric mantle. The present-day concentration of HPI in the crust is inferred from surface heat flow, crustal thickness, and depth of the thermal lithosphere.
Compilation of data available for Archaean cratons.
Names of the cratons | Surface heat flow | Crustal thickness | Lithospheric thickness | Ages | |
mW·m−2 | km | km | Technique | Ga | |
Barberton, Kaapvaal | 40–50 | 35–40 | 200–250 | XD | 3.64–3.2 |
Labrador, North American Craton | 40–60 | 35 | 200–250 | XD | 3.9–3.4 |
Slave, North American Craton | 40–50 | 35 | 180–220 | XD | 3.96–3.1 |
Nuvvuagittuk, Canadian shield north | 40 | 35 | 200–250 | XD | >3.35 |
Pilbara, Western Australia | 40 | 35 | 200–300 | RST | 3.5–3.0 |
Aldan, Siberian craton | 40–50 | 35–40 | 200–250 | XD | 3.4–3.2 |
Central Brazil shield | 40–50 | 40 | 200 | RST | 3.2–3.0 |
Southern Granulite terrane, Indian craton | 40 | 25 | 100–200 | GST–XD | 3.4–3.0 |
Zimbabwe | 40–50 | 40 | 200–250 | XD | 3.64–3.2 |
Ukrainian shield, East-European craton | 40 | 35–40 | >200 | E | 3.6–3.0 |
The present-day thickness of the thermal lithosphere is constrained by (i) pressure–temperature conditions retrieved from subcratonic mantle xenolith carried back to the surface in recent lavas, and (ii) the interpretation of tomographic data (Poupinet and Shapiro, 2009; Poupinet et al., 2003; Roy et al., 1968). Xenolith data are consistent with a present-day thickness of the thermal lithosphere between 200 and 250 km for most cratons (Artemieva and Mooney, 2001). Tomographic data corroborate these estimates at first order but tend to yield a thicker lithosphere (Table 1), with the notable exception of the China and Siberia cratons that are underlain by a thinner lithospheric root (Poupinet and Shapiro, 2009). We consider xenolith data to be more representative of the long-term thermal lithosphere and an average of 225 km has been chosen for the present-day thickness of the thermal lithosphere beneath the Archaean generic craton.
The present-day surface heat flow of Archaean cratons ranges from 40 to 60 mW m−2 (Lachenbruch, 1970; Stern, 2008), and we have chosen the typical value of 45 mW m−2 for the present-day generic Archaean craton. Indeed, zones of high surface heat flow mostly correspond to the edges of cratons that have been recently (during the Phanerozoic) affected by rifting, and it seems reasonable to take a value that is closer to the minimum one. The thermal conductivity of rocks is inversely correlated with temperature, resulting in a variation from about 2 to 4 W m−1·K−1 for the continental crust and from about 3 to 6 W m−1·K−1 for the mantle between 600 and 1300 °C (Hofmeister et al., 2007; Whittington et al., 2009). As we do not want to explore the effect of conductivity on the thermal evolution of the Archaean crust, constant values of 2.5 W m−1·K−1 and of 3.5 W m−1·K−1are attributed to the continental crust and to the lithospheric mantle, respectively.
3 A model for the thermal evolution of Archaean cratonic nuclei
In this section, we present a simple model for the thermal evolution of a generic Archaean cratonic nuclei (Fig. 5). The present-day geothermal gradient of cratonic regions in the lithosphere, in the absence of significant heat advection, is defined by a conductive equilibrium (Carslaw and Jaeger, 1959; England and Thompson, 1984; Richardson and Oxburgh, 1978) expressed by:
(1) |
(2) |
(3) |

Geometric characteristics and boundary conditions used in the model for the computation of the thermal evolution of a generic Archaean craton. a) Present-day characteristics of a generic Archaean craton and the boundary conditions that serve as a basis to calculate the present-day crustal radiogenic heat production, ), and the mantle heat flow, . The temperature is fixed at 0 °C at the surface and at 1300 °C at the base of the lithosphere. The surface heat flow, , is fixed at 45 mW m−2. The thickness of the crust, zc, is fixed at 40 km, and the thickness of the lithosphere, zl, is fixed at 225 km. The thermal conductivities of the crust, kc, and of the mantle, km, are constant and fixed at 2.5 W m−1·K−1 and at 3.5 W m−1·K−1, respectively. The heat production is nil in the mantle. b) Geometric characteristics and boundary conditions that serve as a basis for calculating the thermal evolution with time of a generic Archaean craton. The mantle heat flow, , is considered to be the same as the present-day one and the heat production at a given time, , is calculated taking into account the natural radioactive decay of heat producing isotopes, 238U, 235U, 232Th, 40K. The thickness of the thermal lithosphere, , and the surface heat flow, , underlined in grey, are left free.
This yields the expressions of T(z), the temperature as a function of depth (z) in the crust and in the mantle:
(4) |
(5) |
(6) |
A value for the present-day crustal radiogenic heat production is obtained by fixing the base of the conductive lithosphere as the 1300 °C isotherm, the surface temperature as the 0 °C isotherm (Fig. 3).
The present-day crustal concentration in radioactive isotopes (U238, U235, Th232, K40) is then inferred from the present-day crustal radiogenic heat production (). The contribution of the various heat producing isotope is given by (Turcotte, 1980; Turcotte and Schubert, 1982):
(7) |
The past heat production, , is calculated as a function of time accounting for the secular decrease of the concentration of heat producing isotopes, which is given by:
(8) |
(9) |
(10) |
Considering the assumptions and boundary conditions presented above, the present-day geothermal gradient calculated for the Archaean generic craton is characterized by a Moho temperature of 484 °C and a radiogenic crustal heat production of 0.74 mW m−3 (Fig. 6) yielding a crustal radiogenic heat flow of 29.6 mW m−2. These values are at first order consistent with the ones determined today by geophysical data and experimental means, namely between 400 and 600 °C for the Moho temperature and a range of values between 0.29 and 1.00 mW m−3 for the radiogenic crustal heat production in Archaean regions (Artemieva and Mooney, 2001; Rudnick et al., 1998; Stern, 2008; Turcotte and Schubert, 1982). This consistency gives some credit to our underlying assumptions and thus to the thermal evolution that we model for the generic Archaean craton.

Generic model for the thermal evolution of an Archaean craton. In this model, the thermal evolution of the generic Archaean craton is solely a function of the secular decay of its initial concentration in heat producing isotopes since the Archaean. The present-day and Archaean values for the concentration in HPI have been computed based on the 1D conductive heat equation using values chosen after a critical analysis of the present-day structural and thermal characteristics of Archaean cratons (Table 1) and boundary conditions discussed in the text. The successive geotherms are calculated considering the secular decay of HPI (Table 2). The 700–900 °C temperature range corresponds to typical solidus temperatures of continental acidic to mafic rocks. The zone affected by partial melting from 3.5 to 2.5 Ga is coloured in yellow and indicated by the arrow. The 1300 °C isotherm is underlined as a proxy for the base of the thermal lithosphere.
Taking into account the decay of HPI, the calculated geothermal gradient for the generic Archaean craton at 3.5 Ga is characterised by a Moho temperature of 900 °C, a thermal lithosphere thickness of 132 km (Fig. 6) and a crustal heat production of 2.04 mW m−3 yielding a crustal radiogenic heat flow of 81.4 mW m−2.
The Archaean Moho temperature of the generic craton is significantly higher than the present-day one and this is only caused by the higher radioactive heat production in Archaean time (2.04 mW m−3 instead of 0.74 mW m−3). This higher heat production also implies a higher surface heat flux and a thinner thermal lithosphere in Archaean time. Secular decrease of the HPI content with time is responsible for exponential cooling of the Archaean crust at a pace that is solely controlled by the decay constant of radioactive isotopes, corresponding typically to a half-life in the order of the billion years. For example, according to this model, the temperature equivalent to the continental solidus (taken as ranging between 700 and 900 °C) reaches the Moho at ca. 2.6 Ga, about one billion years after the accretion of the Archaean crust, taken as 3.5 Ga. In other words, the main outcome of this model is that the Archaean cratonic crust is potentially partially molten for a billion years after its accretion. In this case, the actual temperature of the molten region will be buffered by the latent heat of melting and by melt migration and the value calculated here will not be reached (Depine et al., 2008; Hodges et al., 1988; Stüwe, 1995).
The limitations of this model are related to the major unknowns regarding the thermal evolution of the Archaean lithosphere, namely (i) the evolution of the mantle heat flux through time, (ii) the temperature-dependence of conductivity, and (iii) the distribution of HPI within the crust. In this generic model, by considering a homogenous distribution of HPI and fixed thermal parameters, we neglect the impact of these parameters. Higher heat production of the mantle, more vigorous convection and release of accretion heat may be responsible for a higher mantle heat flow during the Archaean than the one assumed here. Accordingly, the mantle heat flow value considered in our model is rather conservative, and taking a higher mantle heat flux will increase the temperature in the crust and strengthen the conclusion that the crust had been likely partially molten for about one billion year. Temperature dependence of conductivity may modify the exact shape of the geotherm, but is likely of second order considering the other sources of uncertainties of our calculation. The choice of a homogeneous distribution of HPI in the crust is the main limitation of our model. Some authors have identified an empirical relationship between heat production of surface rocks and surface heat flow at the scale of “heat flow provinces” (Lachenbruch, 1970; Roy et al., 1968), although it has been challenged in some cases (Alessio et al., 2018; Bea, 2012; McLennan and Taylor, 1996). This relationship implies a higher concentration of HPI in the upper crust, which would cause an increase in the thermal gradient close to the surface, but a reduction of the Moho temperature. A more thorough investigation of this parameter is beyond the scope of this paper, but will be addressed in future research.
The thermal evolution of the early Earth is subject to debate, but all authors agree that it is dominated by secular cooling since at least 3.0 Ga (Herzberg et al., 2010; Korenaga, 2006; Labrosse and Jaupart, 2007). Nevertheless, none of these thermal parameters affect the secular decay of HPI with time and thus the global cooling pattern of Archaean cratons. Accordingly, despite its exponential form, secular decay of HPI with time results in slow cooling of the Archaean crust owing to the half-life of radioactive isotopes (238U, 235U, 232Th, and 40K), that is, in the order of a billion years. For a present-day geotherm of the Archaean generic craton characterized by a Moho temperature of about 480 °C, it is calculated that the cooling rate at Moho depth was 0.12 °C/Myr, with an Archaean Moho temperature of about 900 °C at 3.5 Ga. This implies that the temperature in the lower crust should have stayed at a temperature above the solidus temperature for a billion year (Fig. 6).
4 Sensitivity analysis of the 1D thermal model
Given the range of possible values for the parameters influencing the thermal history of the generic craton, it is essential to quantify the uncertainties associated with the predicted temperature profiles back in time, and more precisely the probability density function of temperature profiles. For this purpose, the range of crustal thicknesses, lithosphere thicknesses and surface heat flows at present time given in Tables 1 and 2 may be represented by probability density functions. For each of these parameters, we assume a Gaussian probability function with a mean and standard deviation chosen to cover within 2σ (95.45% of probability) the range of inferred values for most of the craton list in Table 1 (see Table 3). A Monte Carlo method is then used to propagate the probability density functions of each parameter and estimate the probability density functions associated with the temperature profiles and temperatures of the Moho as a function of time (Figs. 7 and 8).
List of the parameters and of their values used in the models.
Parameters | Value(s) | Unit | |
a) Present day | |||
Equations 4 and 5 | |||
Surface heat flux | 45 | mW·m−2 | |
Crustal thickness | 40 | km | |
Lithospheric thickness | 225 | km | |
Temperature at the base of the lithosphere | 1573 | K | |
Temperature at the surface | 273 | K | |
Thermal conductivity of the crust | 2.5 | W·m−1·K−1 | |
Thermal conductivity of the lithospheric mantle | 3.5 | W·m−1·K−1 | |
b) Thermal evolution with time | |||
Equations 6 and 7 | |||
Mean heat production | |||
238U | 9.37E-05 | W·kg−1 | |
235U | 5.69E-04 | W·kg−1 | |
232Th | 2.69E-05 | W·kg−1 | |
40K | 2.769E-05 | W·kg−1 | |
Decay constant | |||
238U | 1.55E-10 | yr−1 | |
235U | 9.85E-10 | yr−1 | |
232Th | 4.81E-11 | yr−1 | |
40K | 5.54E-10 | yr−1 | |
Density of the crust | 2.7 | ||
Eqs.(9) and (10) | |||
Crustal thickness | 40 | km | |
Thermal conductivity of the crust | 2.5 | W·m−1.K−1 | |
Thermal conductivity of the lithospheric mantle | 3.5 | W·m−1.K−1 | |
Temperature at the surface | 273 | K | |
Mantle heat flow | 15.4 | mW·m−2 | |
Heat production (after Eqs. (3), (5), and (6)) | |||
Present day | 0.74 | μW·m−3 | |
0.5 Ga | 0.81 | μW·m−3 | |
1 Ga | 0.90 | μW·m−3 | |
1.5 Ga | 1.02 | μW·m−3 | |
2 Ga | 1.17 | μW·m−3 | |
2.5 Ga | 1.38 | μW·m−3 | |
3 Ga | 1.65 | μW·m−3 | |
3.5 Ga | 2.03 | μW·m−3 |
Parameters of the (normal) probability density function used in the Monte Carlo simulation. The range of values within 2σ corresponds to a probability of 95.45%.
Parameter | Mean | Standard deviation (σ) | Range of values within 2σ |
Crustal thickness | 40 km | 2.5 km | 35–45 km |
Present lithosphere thickness | 225 km | 37.5 | 150–300 km |
Present surface heat flow | 45 mW m−2 | 2.5 mW m−2 | 40–50 mW·m−2 |

Temperature profile within the crust and the mantle at 3.5 Ga. The red curve is the average value of the probability density function of the temperature profile. The pink areas delineate the 1σ and 2σ values, corresponding respectively to a probability of 68.27% and 95.45%.

Evolution of the Moho temperature (taken at 40 km) from 3.5 Ga up to the present day. The red curve is the average value of the probability density function. The pink areas correspond to 1σ and 2σ values. The probability that the Moho temperature was greater than 700 °C as a function of age can be inferred from a Monte Carlo simulation (See Table 4).
This analysis indicates that it was highly probable (94.74%) that the Moho temperature of the Archaean generic craton was higher than 700 °C (a proxy for the continental solidus) at 3.5 Ga. It is also very likely that the Moho temperature remained above the continental solidus for 500 Myr (76.55%) and the probability for a temperature above the solidus for 1 Ga is about 39% Table 4).
Probability that the Moho temperature was above the solidus temperature of the continental crust (700 °C) as a function of time.
Probability that T > 700 °C at 1.5 Ga | 1.32% |
Probability that T > 700 °C at 2.0 Ga | 10.34% |
Probability that T > 700 °C at 2.5 Ga | 39.09% |
Probability that T > 700 °C at 3.0 Ga | 76.55% |
Probability that T > 700 °C at 3.5 Ga | 94.74% |
4.1 Implication for Archaean crustal evolution and tectonics
The thermal model presented above, proposing that the Archaean crust remained partially molten for about a billion years after its accretion, is consistent with the presence of migmatites coring crustal-scale domes in most if not all Archaean nuclei, attesting to widespread partial melting of the Archaean crust. This proposition provides a new key to interpret the geological and geophysical record of Archaean rocks and opens new perspectives on the understanding of Archaean crustal dynamics.
First of all, it is important to mention that this model considers a felsic-to-intermediate dominant composition for the Archaean crust, which was not significantly modified since its initial accretion. It is clearly an end-member neglecting the input of juvenile crust after initial accretion. This end-member is advocated by some authors based on isotopic signatures of Hadean zircon grains with ages as old as ca. 4.4 Ga (Kemp et al., 2010; Wilde et al., 2001). Besides zircons, the oldest rocks preserved in Archaean cratonic nuclei are TTG gneisses yielding U–Pb zircon ages as old as ca. 4.0 Ga in the Slave Province (Bowring et al., 1989; Iizuka et al., 2007) and as ca. 3.9 Ga in the North Atlantic craton of Labrador and Greenland (Nutman and Collerson, 1991; Nutman et al., 1997). The ubiquitous magmatic complexes forming Archaean cratonic nuclei display a similar evolution from early emplacement of large volumes of TTG suites attributed to partial melting of meta-igneous mafic rocks during the Meso and NeoArchaean, followed by emplacement of more differentiated granites (high-K metaluminous and peraluminous granites) considered as reflecting the contribution of metasediments (Laurent et al., 2014). This evolution is consistent with a single differentiation trend for cratonic nuclei during the entire Archaean period, as opposed to a succession of events of juvenile inputs followed by periods of crustal differentiation. Moreover, the lower crustal refractory mafic residue left after extraction of felsic melts might be recycled back into the mantle by gravitational instabilities shifting the average composition of Archaean crust toward the felsic pole (Johnson et al., 2014).
Secondly, it should be mentioned that a similar thermal model has been used by Rey and Coltice (2008) to explore the secular evolution of the crustal thickness with the idea that gravity-driven flow precludes the persistence of a thick low-viscosity partially molten layer at the base of the Archaean crust (Flament et al., 2008). This idea is seducing, but considers that the stability of the crust is controlled by its integrated strength in the vertical dimension (Kusznir and Park, 1987; Sandiford and Powell, 1990), which is not applicable, for example, to the Himalaya–Tibet region or the Andean plateau, the thickest continental crust on the present-day Earth displaying a substantial partially molten orogenic root (Nelson, 1996; Schilling et al., 1997; Schilling and Partzsch, 2001; Vanderhaeghe and Teyssier, 2001). This example demonstrates that discussing the stability of a weak partially molten orogenic root requires a 3D approach and that considering the 1D integrated strength is too simplistic. Indeed, the presence of cool and strong boundaries around thermally relaxed orogenic cores might prevent the hot and low-viscosity orogenic root to flow laterally, which is marked by the formation of high continental plateaux (Vanderhaeghe et al., 2003). Moreover, the structural and metamorphic record of Archaean rocks is marked by pervasive deformation and crustal-scale domes and basins dissected by thick shear zones with vertical foliations active under amphibolite to granulite facies indicative of distributed crustal thickening coeval with gravitational instabilities within a crustal-scale partially molten layer (Chardon et al., 2009; Collins et al., 1998; Martelat et al., 2000; Van Kranendonk et al., 2007b; Wiemer et al., 2018). The development of these domes is associated with the emplacement of plutonic and volcanic rocks that are for the most part segregated from this partially molten layer and display ages spanning the entire Archaean period (Ketchum et al., 2004; Laurent et al., 2014; Zeh et al., 2011).
Models interpreting these structural, metamorphic, and magmatic features as reflecting either (i) repeated upwelling of mantle plumes associated with polyphased metamorphism and associated crustal differentiation (Bédard, 2006), or (ii) successive tectonic accretion of magmatic arcs to the margin of older cratonic nuclei along prolonged convergent plate boundary marked by shallow subduction (Bleeker and Hall, 2007; Helmstaedt, 2009; Smithies et al., 2003), imply that the tectonic history of Archaean cratons is related to lithospheric-scale dynamics. This contrasts with the models proposed for the stabilization of the Archaean cratons that invoke the development of a stable, buoyant and resistant lithospheric root just after the accretion of the Archaean crust (Griffin et al., 2003; Pollack, 1986; Poupinet et al., 2003). As an alternative, following Bédard (2018), we propose to interpret the geological record of the oldest part of Archaean cratons in the frame of an initially largely partially molten Archaean crust affected by widespread deformation of the crust, above a roughly stabilized lithospheric keel. During the EoArchaean, the thick partially molten layer will promote convection as recently proposed for such types of Archaean cratons (Wiemer et al., 2018), but also for Phanerozoic thermally relaxed orogenic roots (Riel et al., 2016; Vanderhaeghe et al., 2018, Fig. 9a). In this scenario, the sequence of U–Pb ages recorded by zircon grains at the scale of grains, rocks, and cratons is not attributed to a succession of tectonic–magmatic events, but is rather interpreted as reflecting repeated dissolution–crystallization cycles of zircon grains as they are entrained in convection cells affecting the partially molten crust. The vigour of convective overturn decreases with progressive cooling during the Archaean, which is associated with fractional crystallisation of magmas generated by partial melting of an ultramafic/mafic primitive crust leading to the differentiation of the tonalite–trondhjemite–granodiorite suite. The last stage of cooling and crystallization of the Archaean crust is marked by the development of crustal-scale diapirs and by the segregation and emplacement of highly differentiated high-K granites (Fig. 9).

Model for the dynamics of a partially molten Archaean crust. a) The EoArchaean mafic crust is affected by widespread partial melting and the development of gravitational instabilities associated with the segregation of magmas leading to the emplacement of plutonic rocks of the TTG suite. b) The NeoArchaean crust records the end of crystallization of the partially molten layer, leading to the formation of TTG migmatites (grey gneisses) and to the emplacement of highly differentiated high-K granites.
5 Conclusion
Modelling of the thermal evolution of a generic Archaean craton indicates that its Moho may have experienced a temperature close to 900 °C after crustal accretion associated with a thickness of the thermal lithosphere of about 130 km. This initial stage is followed by exponential cooling owing to secular decay of heat producing elements. The main outcome of this model is that the Archaean crust might have remained partially molten for a billion years before the solidus temperature reached the Moho at the Archaean–Proterozoic transition. Considering that the Archaean crust cooled slowly from a partially molten state opens new perspectives on the way to interpret the geological record of Archaean cratons. In particular, it provides new perspectives to take into account the large range of geochronological ages and widespread homogeneous deformation marked by domes and gravitational instabilities that characterize the core of the continent as well as the associated single progressive differentiation trends of cratons during the Archaean Eon.
Acknowledgments
This paper is the result of a ten years long scientific journey that involved a number of Master students, including Séverine Blouin, Julie Garry, Paul Grauwin, Sylvain Michel, Sophie Audonneau, Luca Guillaumot, and Madeline Olivier. Bruno Scaillet is warmly thanked for smooth editorial handling. The current version of this paper benefited from constructive reviews by Laurent Jolivet and an anonymous reviewer. The scientific demonstration, although similar at first order, was significantly improved to respond to the critical comments of Patrice Rey and of an anonymous reviewer on a previous version of this paper. This contribution is part of the URCO project supported by the NEEDS CNRS program.