Research ArticleGEOCHEMISTRY

Argon constraints on the early growth of felsic continental crust

See allHide authors and affiliations

Science Advances  20 May 2020:
Vol. 6, no. 21, eaaz6234
DOI: 10.1126/sciadv.aaz6234


The continental crust is a major geochemical reservoir, the evolution of which has shaped the surface environment of Earth. In this study, we present a new model of coupled crust-mantle-atmosphere evolution to constrain the growth of continental crust with atmospheric 40Ar/36Ar. Our model is the first to combine argon degassing with the thermal evolution of Earth in a self-consistent manner and to incorporate the effect of crustal recycling and reworking using the distributions of crustal formation and surface ages. Our results suggest that the history of argon degassing favors rapid crustal growth during the early Earth. The mass of continental crust, highly enriched in potassium, is estimated to have already reached >80% of the present-day level during the early Archean. The presence of such potassium-rich, likely felsic, crust has important implications for tectonics, surface environment, and the regime of mantle convection in the early Earth.


The continental crust covers ~40% of Earth’s surface and accounts for ~80% of the total volume of Earth’s crust. Its relatively low density and great thickness allow its surface to be exposed above sea level, providing a unique environment for life. The evolution of continental crust is directly linked with a number of important issues in Earth sciences, such as the origin of plate tectonics (1), the history of mantle degassing (2), and the distribution of heat-producing elements within the solid Earth (3). During crustal formation, incompatible elements in the parental mantle, which include rare earth elements and heat-producing elements such as K, U, and Th, are preferentially partitioned to the crust, and noble gases are degassed to the atmosphere. In addition, crustal recycling continuously returns incompatible elements to the mantle, and along with crustal reworking, they redistribute these elements among silicate reservoirs. Consequently, even though the continental crust comprises only ~0.5% of the mass of the bulk silicate Earth (BSE), it is a major geochemical reservoir, whose composition places important constraints on the chemical budget of other terrestrial reservoirs. Thus, understanding how the mass of continental crust and its composition evolve helps unravel the workings of the Earth system as a whole.

A number of Earth scientists have tried to constrain crust-mantle evolution [e.g., (48)]. Whereas the existing models of continental growth exhibit great diversity, it has recently been shown that such diversity is mostly superficial (9). Following the classification scheme of Korenaga (9), we can group them into three categories: crust-based, mantle-based, and others. Crust-based models are estimates on the present-day distribution of the formation age of continental crust [e.g., (5, 10, 11)]. They do not contain information of the crust that has been recycled into the mantle, thereby serving as the lower bound on net crustal growth. Mantle-based models aim at the net growth of continental crust by using the complementary nature of the depleted mantle and the crust [e.g., (7, 8)]. Because of crustal recycling, the appearance of these two types of models can be very different, but this does not necessarily indicate inconsistency. The latest crust-based model (11) and the latest mantle-based model (8) are in agreement with each other, once the effect of recycling is taken into account.

In this study, we aim to constrain continental evolution using the history of argon degassing. This approach belongs to the third category of crustal growth models [e.g., (2, 12)], because it aims at net growth of continental crust using less direct, albeit quantitative, constraints than the mantle-based approach. Nevertheless, using the degassing history of Earth allows us to understand crustal growth in the broader framework of the Earth system. Inert noble gases residing in the atmosphere are an integrated result of degassing from Earth’s interior, and as such, they can potentially provide important insights into the evolution of terrestrial reservoirs [e.g., (13)]. The relation between crustal growth and argon degassing has been explored most recently by Pujol et al. (2); they provided estimates on the 40Ar/36Ar of the Archean atmosphere derived from the measurements of Archean hydrothermal quartz and used the box model of Hamano and Ozima (14) to infer crustal growth. Given the possibility of radiogenic 40Ar added from the sample, the estimates provided by Pujol et al. (2) are better to be seen as an upper bound on the 40Ar/36Ar of the Archean atmosphere. Even with this limitation, their estimates still present a new exciting opportunity for this field, because degassing models have long suffered from the lack of observational constraints on the ancient atmosphere.

However, as in a number of previous studies, Pujol et al. (2) used directly the net growth of continental crust in their model. Such simplistic treatment of continental growth does not acknowledge that net growth is a product of a dynamic balance between new crust generation and crustal recycling, and as a result, the effect of crustal growth on argon degassing is severely underestimated in their study. Also, the other aspects of their model, such as the timing of sudden degassing and the mode of thermal evolution, are difficult to justify. For example, a sudden degassing event, which plays an important role in the early phase of argon degassing, represents the collective effect of giant impacts, and its timing corresponds to the Moon-forming giant impact. While the details of giant impacts remain debated, all of the solutions of Pujol et al. (2) begin with an atmospheric 40Ar/36Ar of 50 at 170 million years (Ma) after the beginning of the Solar System. With an earlier timing of the last giant impact, which is more consistent with recent estimates on the age of the Moon [e.g., (15)], the sudden degassing event cannot effectively elevate the atmospheric 40Ar/36Ar as high as 50. Similarly, their assumption on the history of mantle degassing, which is inherited from the model of Hamano and Ozima (14), is grossly outdated, given the recent development on the thermal evolution of Earth [e.g., (1, 9)].

In light of the above, we have built a more comprehensive model of argon degassing to better understand its constraints on crustal growth. Our model is the first to investigate the full effects of crustal evolution, including both crustal recycling and reworking, on the degassing history of Earth. Here, the term “crustal recycling” is used to denote the loss of continental crust to the mantle through subduction and delamination, whereas the term “crustal reworking” refers to the processes that change isotopic compositions and reset the apparent age of the preexisting crust (e.g., erosion, sedimentation, and partial melting). Crustal recycling affects both the continental crust and the mantle, while crustal reworking takes place within the continental crust. Also, the net growth of continental crust is defined as the evolution of continental mass present on Earth’s surface, i.e., with the effect of crustal recycling already taken into account. Our model also incorporates geophysical and geochemical constraints on the thermal evolution of Earth and uses several robust geological observations to depict a coherent story of continental evolution. Most previous crustal growth models are based on one kind of observation from either geochemical, geophysical, or petrological aspect. However, these observations provide different and complementary constraints on crust formation. In our new model, we use the 40Ar/36Ar of the present-day atmosphere (16) and the Archean atmosphere (2) to constrain argon degassing, the distributions of formation age (11) and surface age (17) of continental crust to constrain the extent of crustal recycling and reworking, and the Archean and Proterozoic mantle potential temperatures (18) to constrain the thermal history of Earth. In what follows, we first present a brief description of our modeling strategy, then summarize model results, and discuss how the history of argon degassing constrains crustal evolution. The details of our box modeling are given in Methods.


We have built a box model to simulate the degassing history of argon using the following three terrestrial reservoirs: the mantle, the continental crust, and the atmosphere. As in traditional degassing models, the oceanic crust is not considered as a major reservoir in our model because of its relatively short residence time. However, degassing caused by oceanic crust formation is tracked throughout Earth history and is part of the total rate of mantle degassing. Among the three reservoirs, the mantle degasses argon into the atmosphere through magmatism, and the crust does so through crustal recycling and reworking. Potassium-40, which decays to argon-40, is transferred to the continental crust from the mantle during partial melting, and part of it is recycled back into the mantle through subduction. To understand how the model is constrained by different observations, we conduct our modeling in three stages, first with crustal evolution only, then with thermal evolution, and lastly with the history of argon degassing. We choose to use simple Monte Carlo sampling for all three stages of our modeling. By comparing the a posteriori distributions of model parameters with their a priori distribution, this simple random sampling allows us to better quantify how each stage constrains crustal evolution. A total of 16 independent variables (Table 1) are used to control the details of these three stages, and 6 dependent variables (Table 2) are used to ensure internal consistency between the thermal and chemical evolution of Earth. Following Rosas and Korenaga (8), crustal evolution is parameterized using the onset of crustal growth (ts), the initial and present-day crustal recycling rates (Rs and Rp), and the decay constants of crustal recycling and growth rates (κr and κg). The extent of crustal reworking is controlled by the reworking factor (frw), which is the fraction of the initial reworking rate with respect to the initial recycling rate. A wide range of crustal evolution can be modeled by varying these six parameters, and it is relatively easy to find various combinations of those parameters that can satisfy the observed crustal formation and surface age distributions. Even with simple Monte Carlo sampling, the sampling efficiency is about 10%. Next, we couple the accepted models of crustal evolution with different models for the thermal evolution of Earth by taking into account the uncertainties of heat productions and heat fluxes of terrestrial reservoirs. Because these uncertainties are relatively small, it is also easy to find adequate thermal budgets to satisfy the history of mantle cooling during the Archean and Proterozoic; the sampling efficiency at the second stage is as high as about 50%. Last, we couple the successful combinations of crustal evolution and thermal evolution with different scenarios of argon degassing by testing the impact of early sudden degassing as well as incomplete degassing during mantle magmatism. Successful solutions at the third stage are chosen on the basis of the argon isotopic abundances of the present-day atmosphere and 40Ar/36Ar in the Archean atmosphere. This stage is the most time-consuming, and from 8 × 106 iterations of Monte Carlo sampling, we have collected a total of ~3 × 103 successful solutions that exhibit good agreements with all of the observational constraints (Fig. 1).

Table 1 List of independent variables in the argon degassing model.

View this table:
Table 2 List of dependent variables in the argon degassing model.

View this table:
Fig. 1 Observational constraints and the distribution of successful model solutions.

Observational constraints used to select successful solutions are shown in orange. The 50th and 90th percentiles of successful solutions are shown in dark green and light green, respectively. (A) 40Ar/36Ar in the atmosphere. Orange bars are based on 40Ar/36Ar measurements of the Archean hydrothermal quartz by Pujol et al. (2), and the orange dot is the present-day 40Ar/36Ar from Lee et al. (16). (B) Mantle potential temperature. Orange dots are the Archean and Proterozoic mantle potential temperatures from Herzberg et al. (18). (C) Present-day cumulative distribution of crustal formation age. The orange line is from Korenaga (11). Ga, billion years. (D) Present-day cumulative distribution of crustal surface age. The orange line is based on the distribution of zircon U-Th crystallization ages (17).

The a posteriori distributions of the parameters Rs, Rp, ts, κr, κg, and frw characterize important features of crustal evolution models selected at different stages (Fig. 2). On the basis of our previous work (8), Rs and Rp have a priori ranges of 0 to 10 × 1022 kg Ga−1 and 0 to 2 × 1022 kg Ga−1, respectively, and κr is varied from −3 to 3 Ga−1 to cover from convex to concave evolution patterns for crustal recycling, whereas κg is varied from −1 to 30 Ga−1 to cover from late-stage to instantaneous crustal growth. The onset time ts is varied from 0.057 to 0.567 Ga, after the beginning of the Solar System, i.e., from the likely timing of the Moon-forming giant impact (15) to the end of Hadean. The temporal evolution of crustal reworking is assumed to follow that of crustal recycling, and different reworking rates are tested by varying frw within the range of 0.1 to 0.8. Compared to stage 3, stages 1 and 2 do not provide tight constraints on crustal evolution; most parameters exhibit nearly uniform distributions (Fig. 2, D to F), with a slight preference to low recycling rates (Fig. 2, A to C). However, all successful solutions at stage 3 have Rs larger than 2.5 × 1022 kg Ga−1 (Fig. 2A), Rp smaller than 9.5 × 1021 kg Ga−1 (Fig. 2B), and κr larger than 0.2 Ga−1 (Fig. 2C), favoring crustal evolution with intense recycling at the beginning followed by rapid decrease. The rate of crustal reworking is similarly high to that of crustal recycling, as indicated by the distribution of frw (Fig. 2F). Both recycling and reworking are responsible for vigorous crustal degassing in the early Earth and, consequently, the rise of atmospheric 40Ar/36Ar in the Archean (Fig. 1A). Net crustal growth models are characterized by κg, and ~80% of the successful solutions have κg larger than 3 Ga−1, i.e., rapid crustal growth in the early Earth (Fig. 2D). Such rapid crustal generation contributes to elevating the atmospheric 40Ar/36Ar to ~100 at the beginning of the Archean, which is important for matching the Archean constraints (Fig. 1A). The distribution of ts is rather uniform with some preference to later onset (Fig. 2E), which gives time for 40K to decay and stock terrestrial inventories of 40Ar. The evident contrast between the distributions from stage 3 and those from previous stages indicates that degassing history is sensitive to different models of crustal formation and thus can place useful constraints on continental growth. For the a posteriori distributions of other model parameters, see fig. S1.

Fig. 2 A posteriori distributions of crustal evolution parameters, based on ~2 × 104, ~2.5 × 104, and ~3 × 103 successful Monte Carlo solutions from stages 1, 2, and 3, respectively.

Distributions from stages 1, 2, and 3 are shown in green, blue, and red, respectively. (A) Initial recycling rate, (B) present-day recycling rate, (C) decay constant for crustal recycling, (D) decay constant for crustal generation, (E) onset time for crustal formation, and (F) reworking factor.

The preferred models of crustal evolution, as determined by the parameters mentioned above, are visualized in Fig. 3. About 80% of the net crustal growth display rapid formation during the early Archean, with the crustal mass being comparable to the present-day value (Fig. 3A). Such early formation challenges the popular notion of less than 30% of continental crust in the Archean (2, 7, 12), but it is essential to match the atmospheric 40Ar/36Ar in the Archean (Fig. 1A). All successful solutions exhibit high crustal generation during the early Earth, which requires vigorous mantle melting and results in substantial mantle degassing (Fig. 3B). The rapid crustal growth generates a substantial amount of Hadean and Archean continental crust, but as constrained by the distributions of crustal formation and surface ages, intense early recycling and reworking are also necessary to erase and reset the records of the old crust (Fig. 3, C and D), with only ~1 and ~10%, respectively, preserved at present day from the Hadean and the early Archean (Fig. 1, C and D). As a consequence, vigorous crustal degassing during the early Archean is achieved through recycling and reworking, and along with mantle degassing, they have elevated the atmospheric 40Ar/36Ar (Fig. 1A).

Fig. 3 Net crustal growth and crustal generation and destruction rates based on ~3 × 103 successful Monte Carlo solutions.

The middle 50 and 90% of our successful solutions are shown in dark green and light green, respectively. The medians of the successful solutions are shown in black dotted lines. (A) Net crust growth, (B) crustal generation rate, (C) crustal recycling rate, and (D) crustal reworking rate.

To understand how different degassing processes contribute to the atmospheric argon abundance, their instantaneous and cumulative contributions are compared in Fig. 4. In our model, the modern degassing rate of 40Ar from the solid Earth ranges from 7.1 × 107 to 3.6 × 108 mol year−1, which is broadly consistent with the estimate of Bender et al. (19) (1.1 ±0.1 × 108 mol year−1). It is noted that our model can reproduce an acceptable present-day degassing rate without using such a constraint. The mantle degasses argon at different rates while generating continental crust (Kmc), oceanic crust (Kmo), and hot spot islands (Kmp). Among them, Kmc is calculated from the rate of crustal generation, with the consideration of continental crust being the secondary product of oceanic crust (Eqs. 4 and 5), while Kmo and Kmp are inferred from the thermal history of Earth, with the possibility of incomplete degassing taken into consideration (Table 2). With Kmc being the largest contributor to the atmospheric 40Ar through Earth history (Fig. 4A), approximately 50% of the present-day abundance originates in the generation of continental crust (Fig. 4B), which demonstrates that continental growth greatly affects the degassing history of Earth. The continental crust can release a substantial amount of 40Ar in the early Archean through recycling (Krc) and reworking (Krw) (Fig. 4A), thanks to abundant parent isotope 40K in the crust. Whereas the major element composition of continental crust is not directly considered in our modeling, the crust that is as enriched in potassium as the present-day crust should be closer to felsic than mafic. A large quantity of such felsic-like crust in the early Earth [e.g., (20)] conflicts with the prevailing notion of little felsic crust in the early Archean [e.g., (3, 12, 21)], requiring a careful rethinking of tectonics, surface environment, and the style of mantle convection during the early Earth.

Fig. 4 Instantaneous and cumulative contributions of individual mass transfer rates to atmospheric 40Ar based on the median of ~3 × 103 successful Monte Carlo solutions.

Mantle processing rates to generate continental crust (Kmc), oceanic crust (Kmo), and hot spot islands (Kmp) are shown in red, dark blue, and green, respectively. Crustal recycling (Krc) and reworking (Krw) rates are shown in orange and purple, respectively. (A) Instantaneous contributions and (B) cumulative contributions of individual mass transfer rates to atmospheric 40Ar.


Our results show that the degassing history of argon is sufficiently sensitive to different modes of net crustal growth, and all of our successful solutions are consistently characterized by rapid crustal generation with intense recycling and reworking during the early Earth (Fig. 3). The preference for such rapid crustal evolution is largely guided by the high 40Ar/36Ar of the Archean atmosphere. In our model, the potassium content in the continental crust is first tracked backward in time using its decay constant and present-day concentration and then scaled to be proportional to the growth of continental crust. As a consequence, the early crust is assumed to be as enriched in potassium as the present-day crust. The considerable contribution of crustal degassing to the Archean atmospheric 40Ar (Fig. 4) indicates that such a large amount of potassium-enriched, thus possibly felsic crust during the early Earth was essential. The substantial amount of potassium-enriched crust does not necessarily indicate that the continental crust like today already existed in the early Archean, but the equivalent amount of potassium in the crust is required to explain the atmospheric 40Ar. In other words, our model constrains the evolution of crustal mass on the basis of the assumption that the crust is as enriched in potassium as the present-day crust through Earth history; if the early continental crust is not as felsic as assumed [e.g., (3)], our estimate on net crustal growth should serve as a lower bound. Because early crustal growth is already very rapid in our model, however, we suspect that the true crustal growth would not substantially deviate from this lower bound. The secular evolution of crustal composition has been studied with a variety of approaches using the geochemistry of sedimentary and igneous rocks, the weighted average of stratigraphic successions, crustal xenoliths, and seismic crustal structure [e.g., (3, 2224)], but the composition of early continental crust remains controversial [e.g., (12, 20, 21)]. Our study provides a new constraint from a different perspective.

One important feature of our model is the simultaneous application of multiple observational constraints to ensure the internal consistency among the thermal evolution, crustal evolution, and degassing history of Earth, which also allows us to quantitatively investigate the effects of recycling and reworking on crustal degassing. This feature is one of the important differences between this study and Pujol et al. (2). Instead of mechanistically relating various crust-mantle differentiation processes with mantle degassing, Pujol et al. (2) modeled mantle degassing in an abstract manner, with only one parameter. As one can see from Fig. 4, however, the mantle degasses argon through three types of magmatism, each of which follows a different evolutionary trend. To make matters worse, their modeling of mantle degassing assumes an exponential decrease in the vigor of mantle convection. Recent progress on the thermal evolution of Earth, however, suggests more sluggish mantle convection in the past [e.g., (1, 9)]. Moreover, Pujol et al. (2) adopted the approach of Hamano and Ozima (14) to model the rate of crustal degassing, which was inferred from the difference between K-Ar mineral ages and Rb-Sr whole-rock ages; i.e., they only considered the contribution of reworking to crustal degassing. Without taking into account the effect of crustal recycling, the total crustal contribution to degassing history is underestimated (Fig. 4). With such low crustal degassing, Pujol et al. (2) were able to match the Archean atmospheric constraint only by introducing a late onset of sudden degassing, which effectively elevated their atmospheric 40Ar/36Ar to 50 (see their fig. 2a). As mentioned in Introduction, sudden degassing at 170 Ma after the beginning of the Solar System is probably too late to be consistent with what geochronological studies suggest [e.g., (15)]. Pujol et al. (2) provided an invaluable observational constraint on the ancient atmosphere, but for the aforementioned reasons, we reached different conclusions using the same Archean data.

The substantial differences between Pujol et al. (2) and this study underline the importance of treating crustal evolution properly in degassing models, although this issue is not widely appreciated. For example, a recent study of mantle xenon isotopes (25) suggests that the deep volatile cycles shifted from a net degassing to a net regassing regime around 2.5 Ga ago. They used, however, the net growth of continental crust to calculate the mantle degassing rate corresponding to continental growth, which is the same treatment used by Pujol et al. (2). A recent study on Archean komatiites (26) suggests that the subduction of water was initiated before 3.3 Ga ago, which is more consistent with our modeling results. For these reasons, the degassing models that do not fully acknowledge the effect of continental growth [e.g., (25, 27)] deserve to be revisited carefully.

As mentioned in Introduction, our model of crustal growth belongs to the third category, because it is inferred indirectly from the degassing history of Earth. Nevertheless, this model is in good agreement with the recent mantle-based model by Rosas and Korenaga (8), thereby reinforcing the notion of rapid crustal growth during the Hadean and the early Archean (Fig. 3A and their fig. 1c). On the basis of the samarium-neodymium isotope systems, Rosas and Korenaga (8) were able to place a tighter bond on crustal evolution because the mantle-based approach is more direct. Given that potassium is much more incompatible than samarium and neodymium, however, our model can better constrain the compositional evolution of continental crust. Also, as a heavy noble gas, the atmospheric argon keeps an integrated degassing history of the bulk Earth, thus suffering less from preservation issues than mantle-based models, which rely critically on preserved rock samples. Crust-based models are essentially the present-day distribution of crustal formation ages, and with consideration of crustal recycling, both our study and Rosas and Korenaga (8) are in good agreement with the recent crust-based model of Korenaga (11), so the different approaches appear to be converging. The large offset between the rapid crustal growth (Fig. 3A) and the nearly linear distribution of formation ages (Fig. 1C) can be explained by a time-integrated effect of crustal recycling (8). This model of crustal evolution, characterized with rapid crustal growth and efficient crustal recycling in the early Earth, resembles closely what Armstrong suggested almost 40 years ago (28).

The net growth of continental crust is controlled by a dynamic balance between crustal generation and recycling. The approximately zero net crustal growth after the Hadean, but with nonzero crustal recycling, indicates that new crust must have been continuously formed at the same rate of recycling (Fig. 3). Such persistent crustal formation and destruction through Earth history are consistent with the onset of plate tectonics in the very early Earth [e.g., (29)]. Here, the term “plate tectonics” is defined in a broad sense, a mode of mantle convection with the continuous, wholesale recycling of surface layer, as opposed to stagnant lid convection. Our result of rapid crustal growth with efficient early recycling suggests vigorous mantle convection in the early Hadean, followed by a gradual decrease to the present-day level. Such decline in crustal generation may be explained by the secular cooling of the mantle (18), because a cooler mantle yields less voluminous melting. The corresponding decrease in crustal recycling may be attributed to an increasing preservation potential of continental crust (9). The positive net water flux from the surface into the mantle, as inferred from deep water cycle and continental freeboard (30), suggests the gradual hydration of the convecting mantle through time, and as a consequence, the continental mantle lithosphere, which remains dry owing to its generation mechanism, becomes stronger with respect to the convecting mantle, making crustal recycling less efficient (1). The mode of mantle convection in the early Earth is still controversial [e.g., (17, 21, 29)], and a more careful geodynamical work is warranted (31) to discuss whether mantle convection had switched from a different mode to modern plate tectonics under early Earth conditions.


As briefly described in the “Model formulation and results” section, the terrestrial reservoirs considered in our box model are the mantle, the continental crust, and the atmosphere. The mantle degasses argon to the atmosphere through mantle magmatism, and the crust does so through crustal recycling and reworking. The potassium is transferred between the mantle and crust during partial melting and subduction. We modeled the crustal evolution first to constrain the mass transfer rates between mantle and crust, then simulated the thermal evolution of Earth to infer the rates of mantle magmatism, and calculated the degassing history of argon using the degassing rates acquired from the above two stages. The model formulation for each stage is described in detail below.

Continental crust growth

As mentioned in the “Model formulation and results” section, we follow Rosas and Korenaga (8) for the parameterization of crustal growth and recycling ratesMcc(t)=Mcc(tp)1eκg(tpts)(1eκg(tts))(1)Krc(t)=Rs+RpRs1eκr(tpts)(1eκr(tts))(2)dMcc(t)dt=Kcc(t)Krc(t)(3)where Mcc(t) is the mass of continental crust at time t. As shown in Eq. 3, the time derivative of Mcc(t) is equal to the difference between the crustal generation rate, Kcc(t), and the crustal recycling rate, Krc(t). The present-day crustal mass Mcc(tp) is set to be 2.09 × 1022 kg. The term ts denotes the onset of crustal generation and recycling; Rs and Rp are the rates of crustal recycling at ts and tp, respectively; and κg and κr are the decay constants for Kcc and Krc, respectively. A wide variety of crustal growth patterns can be tested in our model by varying these parameters, covering from late growth to nearly instantaneous growth.

The crustal generation rate Kcc(t) describes how much mass has been added to the crust with time. To calculate the mantle processing rate corresponding to the generation of the continental crust, Kmcprod(t), we use the complementary nature between the crust and the mantleKmcprod(t)=Kcc(t)NK40CC(t)Mcc(t)MM(t)NK40M(t)(4)where MM(t) is the mass of the mantle at time t and NK40M(t) is the number of 40K atoms in the mantle at time t. Both of them can be inferred from mass balance among the mantle, the continental crust, and the BSE.

However, the continental crust does not result from the single-stage melting of the mantle. Considering that at least part of the continental crust is likely to be produced through the secondary melting of oceanic crust, treating all of Kmcprod(t) as the mantle processed rate to generate the continental crust can overestimate the total mantle processing rate. Following Padhi et al. (32), therefore, the mantle processing rate solely responsible for the generation of continental crust, Kmc, is calculated as followsKmc(t)=max(Kmcprod(t)Kmo(t),0)(5)

Crustal reworking is likely to be in sync with recycling, as the former is necessary to break down the crust into subductable sediments, so we set the secular evolution of reworking to be similar to that of recycling, with a factor of frw, which varies between 0.1 and 0.8 in our model. The crustal reworking rate is therefore calculated asKrw(t)=Krc(t)frw(6)

Taking into account that reworking is responsible for the difference between the distributions of crustal formation and surface ages (11), we choose the acceptable reworking factor frw by calculating these distributions and comparing to observational constraints. To do so, we first follow Rosas and Korenaga (8) to model the formation age distribution of continental crust, m(t,τ), where t is the time and τ represents the formation age. The summation of the m(t,τ) over time τ gives the total crustal mass at time tMCC(t)=0tm(t,τ)dτ(7)

In our model, the continental crust is considered to be a homogeneous reservoir at all times, so recycling uniformly affects the crustal parts that are formed at different times; i.e., crustal recycling is independent of formation age. The evolution of m(t,τ) with such age-independent recycling may be expressed asm(t,τ)t=Kmc(t)δ(tτ)Krc(t)MCC(t)m(t,τ)(8)where δ(t) is the Dirac delta function. The present-day cumulative formation age distribution, CFD(τ), may then be calculated asCFD(τ)=1MCC(tp)0τm(tp,τ)dτ(9)

Second, we model the cumulative surface age distribution (CSD) considering the combined effects of crustal generation, recycling, and reworking to the surface age distribution s(t,τ) of continental crust. Crustal reworking resets the apparent age of older crust, and because we consider the continental crust as a single reservoir, reworking uniformly affects older crust. That is, the surface age distribution, s(t,τ), may be calculated ass(t,τ)t=Kmc(t)δ(tτ)Krc(t)MCC(t)s(t,τ)+Krw(t)δ(tτ)Krw(t)MCC(t)s(t,τ)(1δ(tτ))(10)

The present-day cumulative surface age distribution, CSD(τ), may then be calculated as the present-day surface age distribution integrated over τCSD(τ)=1MCC(tp)0τs(tp,τ)dτ(11)

Thermal evolution

The thermal evolution of Earth constrains the mantle processing rates to generate oceanic crust, Kmo, and hot spot islands, Kmp. First, we calculate Kmo considering that the oceanic crust is generated through decompressional melting beneath mid-ocean ridges, so its production rate can be directly linked with the thermal history of Earth. To be more concrete, the mantle processing rate for oceanic crust Kmo can be constrained on the basis of plate velocity, V, and the initial depth of mantle melting, Z, asKmo(t)=Kmo(tp)Z(t) V(t)Z(tp)V(tp)(12)where Kmo(tp) is the present-day mantle processing rate at mid-ocean ridges, which is estimated to be 6.7 × 1014 kg year−1 (33).

The initial depth of melting Z(t) is controlled by mantle potential temperature, Tp, which can be calculated as (32)Z(t)=Tp(t)1150 m(1.2×107(dTdP)S)(13)where g is the gravitational acceleration (9.8 m s−2), ρm is the mantle density (3300 kg m−3), and (dTdP)Sis the adiabatic gradient in the mantle (1.54 × 10−8 K Pa−1).

Then, using the relationship between surface heat flux and plate velocity, the temporal evolution of plate velocity can be calculated asV(t)=V(tp)(Q(t)Q(tp)Tp(tp)Tp(t))2(14)where the present-day plate velocity V(tp) is set to 5 cm year−1 (34), and the present-day mantle heat flux Q(tp) is calculated as the difference between the present-day total terrestrial heat flux (46 ± 3 TW) (35) and the present-day continental crust heat production, Hcc(tp), asQ(tp)=(46+3ε1)HCC(tp)(15)where ε1 is a random variable, which can vary between −1 and 1.

The present-day continental crust heat production is considered to be (24)HCC(tp)=7.5+2.5ε2(16)where ε2 is another random variable, whose range is between −1 and 1.

The evolution of mantle potential temperature, Tp(t), can be tracked backward in time using the heat production, H(t), the mantle heat flux, Q(t), and the core heat flux, Qc(t), according to the following global energy balance (31)CmdTp(t)dt=H(t)Q(t)+Qc(t)(17)where Cm is the heat capacity of the whole mantle (4.97 × 1027 J K−1).

In our model, we take the present-day core heat flux Qc(tp) as a free parameter, which can vary from 5 to 15 TW. To track the evolution of core heat flux Qc, we consider Qc changes linearly with timeQc(t)=ΔQc(tpt)/tp+Qc(tp)(18)where ΔQc is the difference between initial and present-day core heat flux, which can vary between 2 and 5 TW (36).

To integrate Eq. 17, we need to express H and Q as functions of time. Following Korenaga (33), H can be expressed as a function of time using the decay constants and heat production rates of major heat-producing elements within Earth (238U, 235U, 232Th, and 40K) as followsH(t)=H(tp)i=14cipieλiti=14cipi(19)where ci and pi are the present-day relative concentration and the heat generation rate of the isotope in interest, λi is the decay constant, and H(tp) is the present-day mantle heat production, which is calculated as the total BSE heat production of 16 ± 3 TW (37) minus the present-day continental crust heat production Hcc(tp)H(tp)=(16+3ε3)HCC(tp)(20)where ε3 is a random variable, whose range is between −1 and 1. Because the uncertainties of the terrestrial heat flux Q(tp), the present-day continental crust heat production Hcc(tp), and the present-day mantle heat production H(tp) are not related to each other, the three random variables, ε1, ε2, and ε3, vary independently to each other.

As for the mantle heat flux Q(t), we assume it to be constant (i.e., same as Q(tp) over the entire Earth history following Korenaga (1). This assumption results in mantle potential temperature Tp rising gradually from the present-day value of 1350°C to approximately 1700°C during Earth history. The classical scaling law, which can be expressed as Q(t) ≈ αTp(t)β, results in Tp quickly rising and reaching infinity, and such phenomenon is known as thermal catastrophe. Compared with the classical scaling law, using a constant Q is more comparable with the evolution of the mantle potential temperature [e.g., (18)] as well as the geochemical model of Earth’s composition. For more details of the two heat flux scaling laws, see Korenaga (1, 33). This validity of relatively constant heat flux becomes uncertain before ~3.5 Ga ago, as mantle potential temperature in such a deep time is not constrained observationally, and it is probably unrealistic to represent the thermal state of the very early Earth. However, the impact of this early phase of thermal evolution on argon degassing is limited thanks to our parameterization of mantle degassing rates (Eq. 5).

Knowing Qc(t), Q(t), and H(t), we can integrate Eq. 17 backward in time to obtain the mantle potential temperature Tp(t) and then calculate the mantle processing rate to generate oceanic crust Kmo with Eq. 12.

Second, by assuming a linear relation between the mantle processing rate to generate hot spot islands Kmp and the core heat flux Qc, we calculate Kmp as followsKmp(t)=Kmp(tp)Qc(t)Qc(tp)(21)where Kmp(tp) is the present-day rate of plume mass flux, which can be estimated from the present-day plume buoyancy flux, fpb(tp), asKmp(tp)=fpb(tp)αΔT(22)where fpb(tp) is considered to be 55 × 103 kg s−1 (38), α is the thermal expansivity (4 × 10−5 K−1), and ΔT is the temperature anomaly associated with mantle plumes (200 K).

The history of argon degassing

The different modes of degassing are considered before and after a sudden degassing event, which is assumed to have occurred in the early Earth [e.g., (14)]. The high 40Ar/36Ar values (>10,000) reported for mantle-derived ultramafic rocks suggest a substantial degassing of primordial 36Ar in the early Earth. Such catastrophic degassing is likely to have resulted from the highly energetic phase of planetary accretion. We denote the end of such an intense degassing phase by td. At td, the mantle loses a fraction of primordial argon, Fd, through sudden degassing to atmosphere, and after td, the mantle and the crust lose argon to the atmosphere in a gradual manner, corresponding to continuous crustal formation and destruction through Earth history. Because the detailed information of giant impacts is unknown, the timing of sudden degassing td and the degassing fraction Fd are both treated as free parameters in our model, with td varying from 0.05 to 0.1 Ga and Fd varying from 10 to 80%, respectively.

From the beginning of the Solar System, t0, to the time of sudden degassing td, the parent isotope 40K in the BSE gradually decays through time, whereas the amount of daughter isotope 40Ar in the BSE increases accordingly. Meanwhile, 40Ar and 36Ar in the BSE experience sudden degassing at td

In the BSE domainddtNK40BSE(t)=λNK40BSE(t)(23)ddtNA40rBSE(t)=λeNK40BSE(t)Fdδ(ttd)NA40rBSE(t)(24)ddtNA36rBSE(t)=Fdδ(ttd)NA36rBSE(t)(25)

In the atmosphere domain (Atm)ddtNA40rAtm(t)=Fdδ(ttd)NA40rBSE(t)(26)ddtNA36rAtm(t)=Fdδ(ttd)NA36rBSE(t)(27)where λ and λe are the total decay constant and the branch decay constant of 40K, respectively, and N is the number of atoms of the isotope denoted in the subscript contained in the reservoir denoted in the superscript. For example, the number of 40K in the continental crust (CC) is denoted by  NK40CC(t), which can also be expressed asNK40CC(t)=CK40CC(t)Mcc(t)mK40(28)where C is the concentration, M is the reservoir mass, and m is the atomic mass of the isotope in interest.

From td to the present-day tp, argon degassing is modeled as follows. The amount of 40K in the continental crust is tracked backward in time according to its present-day abundance and then scaled to be proportional to the growth of the continental crust. The abundance of 40K in the mantle is calculated from mass balance among the BSE, the mantle, and the continental crust. Argon degassing is considered to take place during mantle magmatism as well as crustal recycling and reworking, with the former parameterized by the mantle processing rates to generate continental crust Kmc, oceanic crust Kmo, and hot spot islands Kmp and the latter parameterized by the crustal recycling rate Krc and the reworking rate Krw.

In the BSE domainddtNK40BSE(t)=λ NK40BSE(t)(29)

In the continental crust domain (CC)NK40cc(t)=NK40cc(tp)eλtMcc(t)Mcc(tp)(30)ddtNA40rCC(t)=λeNK40CC(t)[Krw(t)+Krc(t)]NA40rCC(t)(31)ddtNA36rCC(t)=[Krw(t)+Krc(t)]NA36rCC(t)(32)

In the mantle domain (M)NK40M(t)=NK40BSE(t)NK40CC(t)(33)ddtNA40rM(t)=λeNK40M(t)Km(t)NA40rM(t)(34)ddtNA36rM(t)=Km(t)NA36rM(t)(35)

In the atmosphere domainddtNA40rAtm(t)=Km(t)NA40rM(t)+[Krw(t)+Krc(t)]NA40rCC(t)(36)ddtNA36rM(t)=Km(t)NA36rM(t)+[Krw(t)+Krc(t)]NA36rCC(t)(37)where Km is the combination of Kmc and the effective parts of Kmo and Kmp. Given that argon may not degas completely from basalts erupted in deep water, we consider the possibility of incomplete degassing at mid-ocean ridges (39). We also allow incomplete degassing for hot spot magmatism for three reasons: (i) The plume buoyancy flux estimate of Sleep (38) is subject to large uncertainties (40); (ii) flux estimates based on swell topography are likely to be the upper bound of the buoyancy flux (41); and (iii) thick lithosphere may hinder some of plume magma to reach the surface. The effective parts of Kmo and Kmp are modeled by free parameters  feffKmo and  feffKmp, which can vary from 50 to 100% and 10 to 100 % ,respectively. We assume 100% efficiency during all of other degassing processes. The total mantle degassing rate Km is calculated as followsKm(t)=Kmc(t)+Kmp(t) feffKmp+Kmo(t) feffKmo(38)

In the above mass transfer equations, the five mass transfer rates Kmo, Kmc, Kmp, Krc, and Krw are the critical unknowns in our model. Constraining them with crustal growth and thermal evolution assures us to conduct the model in a self-consistent manner.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank three anonymous reviewers for constructive comments and suggestions. Funding: This article was based on work supported by the NSF under grant EAR-1753916. Author contributions: M.G. performed the modeling and wrote the manuscript. J.K. designed the project, discussed the results, and commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. MATLAB scripts used for our modeling are also included in the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article