## Abstract

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 ^{40}Ar/^{36}Ar. 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.

## INTRODUCTION

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., (*4*–*8*)]. 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 ^{40}Ar/^{36}Ar 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 ^{40}Ar added from the sample, the estimates provided by Pujol *et al.* (*2*) are better to be seen as an upper bound on the ^{40}Ar/^{36}Ar 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 ^{40}Ar/^{36}Ar 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 ^{40}Ar/^{36}Ar 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 ^{40}Ar/^{36}Ar 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.

## MODEL FORMULATION AND RESULTS

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 (*t*_{s}), the initial and present-day crustal recycling rates (*R*_{s} and *R*_{p}), and the decay constants of crustal recycling and growth rates (κ_{r} and κ_{g}). The extent of crustal reworking is controlled by the reworking factor (*f*_{rw}), 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 ^{40}Ar/^{36}Ar in the Archean atmosphere. This stage is the most time-consuming, and from 8 × 10^{6} iterations of Monte Carlo sampling, we have collected a total of ~3 × 10^{3} successful solutions that exhibit good agreements with all of the observational constraints (Fig. 1).

The a posteriori distributions of the parameters *R*_{s}, *R*_{p}, *t*_{s}, κ_{r}, κ_{g}, and *f*_{rw} characterize important features of crustal evolution models selected at different stages (Fig. 2). On the basis of our previous work (*8*), *R*_{s} and *R*_{p} have a priori ranges of 0 to 10 × 10^{22} kg Ga^{−1} and 0 to 2 × 10^{22} 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 *t*_{s} 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 *f*_{rw} 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 *R*_{s} larger than 2.5 × 10^{22} kg Ga^{−1} (Fig. 2A), *R*_{p} smaller than 9.5 × 10^{21} 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 *f*_{rw} (Fig. 2F). Both recycling and reworking are responsible for vigorous crustal degassing in the early Earth and, consequently, the rise of atmospheric ^{40}Ar/^{36}Ar 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 ^{40}Ar/^{36}Ar to ~100 at the beginning of the Archean, which is important for matching the Archean constraints (Fig. 1A). The distribution of *t*_{s} is rather uniform with some preference to later onset (Fig. 2E), which gives time for ^{40}K to decay and stock terrestrial inventories of ^{40}Ar. 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.

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 ^{40}Ar/^{36}Ar 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 ^{40}Ar/^{36}Ar (Fig. 1A).

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 ^{40}Ar from the solid Earth ranges from 7.1 × 10^{7} to 3.6 × 10^{8} mol year^{−1}, which is broadly consistent with the estimate of Bender *et al.* (*19*) (1.1 ±0.1 × 10^{8} 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 (*K*_{mc}), oceanic crust (*K*_{mo}), and hot spot islands (*K*_{mp}). Among them, *K*_{mc} 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 *K*_{mo} and *K*_{mp} are inferred from the thermal history of Earth, with the possibility of incomplete degassing taken into consideration (Table 2). With *K*_{mc} being the largest contributor to the atmospheric ^{40}Ar 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 ^{40}Ar in the early Archean through recycling (*K*_{rc}) and reworking (*K*_{rw}) (Fig. 4A), thanks to abundant parent isotope ^{40}K 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.

## DISCUSSION

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 ^{40}Ar/^{36}Ar 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 ^{40}Ar (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 ^{40}Ar. 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*, *22*–*24*)], 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 ^{40}Ar/^{36}Ar 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.

## METHODS

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 rates

(1)

$${K}_{\text{rc}}(t)={R}_{\mathrm{s}}+\frac{{R}_{\mathrm{p}}-{R}_{\mathrm{s}}}{1-{e}^{-{\mathrm{\kappa}}_{\mathrm{r}}({t}_{\mathrm{p}}-{t}_{\mathrm{s}})}}(1-{e}^{-{\mathrm{\kappa}}_{\mathrm{r}}(t-{t}_{\mathrm{s}})})$$(2)

$$\frac{{\mathit{dM}}_{\text{cc}}(t)}{\mathit{dt}}={K}_{\text{cc}}(t)-{K}_{\text{rc}}(t)$$(3)where *M*_{cc}(*t*) is the mass of continental crust at time *t*. As shown in Eq. 3, the time derivative of *M*_{cc}(*t*) is equal to the difference between the crustal generation rate, *K*_{cc}(*t*), and the crustal recycling rate, *K*_{rc}(*t*). The present-day crustal mass *M*_{cc}(*t*_{p}) is set to be 2.09 × 10^{22} kg. The term *t*_{s} denotes the onset of crustal generation and recycling; *R*_{s} and *R*_{p} are the rates of crustal recycling at *t*_{s} and *t*_{p}, respectively; and κ_{g} and κ_{r} are the decay constants for *K*_{cc} and *K*_{rc}, 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 *K*_{cc}(*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,

(*t*), we use the complementary nature between the crust and the mantle

(4)where *M*_{M}(*t*) is the mass of the mantle at time *t* and

is the number of ^{40}K 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

${K}_{\text{mc}}^{\text{prod}}(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, *K*_{mc}, is calculated as follows

(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 *f*_{rw}, which varies between 0.1 and 0.8 in our model. The crustal reworking rate is therefore calculated as

(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 *f*_{rw} 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 *t*

(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 as

(8)where δ(*t*) is the Dirac delta function. The present-day cumulative formation age distribution, CFD(τ), may then be calculated as

(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 as

(10)

The present-day cumulative surface age distribution, CSD(τ), may then be calculated as the present-day surface age distribution integrated over τ

$$\text{CSD}(\mathrm{\tau})=\frac{1}{{M}_{\text{CC}}({t}_{\mathrm{p}})}{\displaystyle {\int}_{0}^{\mathrm{\tau}}}s({t}_{\mathrm{p}},{\mathrm{\tau}}^{\prime})d{\mathrm{\tau}}^{\prime}$$(11)

### Thermal evolution

The thermal evolution of Earth constrains the mantle processing rates to generate oceanic crust, *K*_{mo}, and hot spot islands, *K*_{mp}. First, we calculate *K*_{mo} 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 *K*_{mo} can be constrained on the basis of plate velocity, *V*, and the initial depth of mantle melting, *Z*, as

(12)where *K*_{mo}(*t*_{p}) is the present-day mantle processing rate at mid-ocean ridges, which is estimated to be 6.7 × 10^{14} kg year^{−1} (*33*).

The initial depth of melting *Z*(*t*) is controlled by mantle potential temperature, *T*_{p}, which can be calculated as (*32*)

(13)where *g* is the gravitational acceleration (9.8 m s^{−2}), ρ_{m} is the mantle density (3300 kg m^{−3}), and

is 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 as

$$V(t)=V({t}_{\mathrm{p}}){\left(\frac{Q(t)}{Q({t}_{\mathrm{p}})}\frac{{T}_{\mathrm{p}}({t}_{\mathrm{p}})}{{T}_{\mathrm{p}}(t)}\right)}^{2}$$(14)where the present-day plate velocity *V*(*t*_{p}) is set to 5 cm year^{−1} (*34*), and the present-day mantle heat flux *Q*(*t*_{p}) 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, *H*_{cc}(*t*_{p}), as

(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*)

(16)where ε_{2} is another random variable, whose range is between −1 and 1.

The evolution of mantle potential temperature, *T*_{p}(*t*), can be tracked backward in time using the heat production, *H*(*t*), the mantle heat flux, *Q*(*t*), and the core heat flux, *Q*_{c}(*t*), according to the following global energy balance (*31*)

(17)where *C*_{m} is the heat capacity of the whole mantle (4.97 × 10^{27} J K^{−1}).

In our model, we take the present-day core heat flux *Q*_{c}(*t*_{p}) as a free parameter, which can vary from 5 to 15 TW. To track the evolution of core heat flux *Q*_{c}, we consider *Q*_{c} changes linearly with time

(18)where Δ*Q*_{c} 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 (^{238}U, ^{235}U, ^{232}Th, and ^{40}K) as follows

(19)where *c _{i}* and

*p*are the present-day relative concentration and the heat generation rate of the isotope in interest, λ

_{i}*is the decay constant, and*

_{i}*H*(

*t*

_{p}) 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

*H*

_{cc}(

*t*

_{p})

(20)where ε_{3} is a random variable, whose range is between −1 and 1. Because the uncertainties of the terrestrial heat flux *Q*(*t*_{p}), the present-day continental crust heat production *H*_{cc}(*t*_{p}), and the present-day mantle heat production *H*(*t*_{p}) 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*(*t*_{p}) over the entire Earth history following Korenaga (*1*). This assumption results in mantle potential temperature *T*_{p} 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*) ≈ α*T*_{p}(*t*)^{β}, results in *T*_{p} 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 *Q*_{c}(*t*), *Q*(*t*), and *H*(*t*), we can integrate Eq. 17 backward in time to obtain the mantle potential temperature *T*_{p}(*t*) and then calculate the mantle processing rate to generate oceanic crust *K*_{mo} with Eq. 12.

Second, by assuming a linear relation between the mantle processing rate to generate hot spot islands *K*_{mp} and the core heat flux *Q*_{c}, we calculate *K*_{mp} as follows

(21)where *K*_{mp}(*t*_{p}) is the present-day rate of plume mass flux, which can be estimated from the present-day plume buoyancy flux, *f _{pb}*(

*t*

_{p}), as

(22)where *f*_{pb}(*t*_{p}) is considered to be 55 × 10^{3} 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 ^{40}Ar/^{36}Ar values (>10,000) reported for mantle-derived ultramafic rocks suggest a substantial degassing of primordial ^{36}Ar 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 *t*_{d}. At *t*_{d}, the mantle loses a fraction of primordial argon, *F*_{d}, through sudden degassing to atmosphere, and after *t*_{d}, 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 *t*_{d} and the degassing fraction *F*_{d} are both treated as free parameters in our model, with *t*_{d} varying from 0.05 to 0.1 Ga and *F*_{d} varying from 10 to 80%, respectively.

From the beginning of the Solar System, *t*_{0}, to the time of sudden degassing *t*_{d}, the parent isotope ^{40}K in the BSE gradually decays through time, whereas the amount of daughter isotope ^{40}Ar in the BSE increases accordingly. Meanwhile, ^{40}Ar and ^{36}Ar in the BSE experience sudden degassing at *t*_{d}

In the BSE domain

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{40}\mathrm{K}}^{\text{BSE}}(t)=-\mathrm{\lambda}{N}_{{}_{}{}^{40}\mathrm{K}}^{\text{BSE}}(t)$$(23)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{BSE}}(t)={\mathrm{\lambda}}_{\mathit{e}}{N}_{{}_{}{}^{40}\mathrm{K}}^{\text{BSE}}(t)-{F}_{\mathrm{d}}\mathrm{\delta}(t-{t}_{\mathrm{d}}){N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{BSE}}(t)$$(24)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\text{BSE}}(t)=-{F}_{\mathrm{d}}\mathrm{\delta}(t-{t}_{\mathrm{d}}){N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\text{BSE}}(t)$$(25)

In the atmosphere domain (Atm)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{Atm}}(t)={F}_{\mathrm{d}}\mathrm{\delta}(t-{t}_{\mathrm{d}}){N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{BSE}}(t)$$(26)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\text{Atm}}(t)={F}_{\mathrm{d}}\mathrm{\delta}(t-{t}_{\mathrm{d}}){N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\text{BSE}}(t)$$(27)where λ and λ_{e} are the total decay constant and the branch decay constant of ^{40}K, 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 ^{40}K in the continental crust (CC) is denoted by

, which can also be expressed as

$${N}_{{}_{}{}^{40}\mathrm{K}}^{\text{CC}}(t)=\frac{{C}_{{}_{}{}^{40}\mathrm{K}}^{\text{CC}}(t){M}_{\text{cc}}(t)}{{m}_{{}_{}{}^{40}\mathrm{K}}}$$(28)where *C* is the concentration, *M* is the reservoir mass, and *m* is the atomic mass of the isotope in interest.

From *t*_{d} to the present-day *t*_{p}, argon degassing is modeled as follows. The amount of ^{40}K 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 ^{40}K 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 *K*_{mc}, oceanic crust *K*_{mo}, and hot spot islands *K*_{mp} and the latter parameterized by the crustal recycling rate *K*_{rc} and the reworking rate *K*_{rw}.

In the BSE domain

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{40}\mathrm{K}}^{\text{BSE}}(t)=-\mathrm{\lambda}{N}_{{}_{}{}^{40}\mathrm{K}}^{\text{BSE}}(t)$$(29)

In the continental crust domain (CC)

$${N}_{{}_{}{}^{40}\mathrm{K}}^{\text{cc}}(t)={N}_{{}_{}{}^{40}\mathrm{K}}^{\text{cc}}({t}_{\mathrm{p}}){e}^{\mathrm{\lambda}t}\frac{{M}_{\text{cc}}(t)}{{M}_{\text{cc}}({t}_{\mathrm{p}})}$$(30)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{CC}}(t)={\mathrm{\lambda}}_{\mathrm{e}}{N}_{{}_{}{}^{40}\mathrm{K}}^{\text{CC}}(t)-[{K}_{\text{rw}}(t)+{K}_{\text{rc}}(t)]{N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{CC}}(t)$$(31)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\text{CC}}(t)=-[{K}_{\text{rw}}(t)+{K}_{\text{rc}}(t)]{N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\text{CC}}(t)$$(32)

In the mantle domain (M)

$${N}_{{}_{}{}^{40}\mathrm{K}}^{\mathrm{M}}(t)={N}_{{}_{}{}^{40}\mathrm{K}}^{\text{BSE}}(t)-{N}_{{}_{}{}^{40}\mathrm{K}}^{\text{CC}}(t)$$(33)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\mathrm{M}}(t)={\mathrm{\lambda}}_{\mathrm{e}}{N}_{{}_{}{}^{40}\mathrm{K}}^{\mathrm{M}}(t)-{K}_{\mathrm{m}}(t){N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\mathrm{M}}(t)$$(34)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\mathrm{M}}(t)=-{K}_{\mathrm{m}}(t){N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\mathrm{M}}(t)$$(35)

In the atmosphere domain

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{Atm}}(t)={K}_{\mathrm{m}}(t){N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\mathrm{M}}(t)+[{K}_{\text{rw}}(t)+{K}_{\text{rc}}(t)]{N}_{{}_{}{}^{40}\mathrm{A}\mathrm{r}}^{\text{CC}}(t)$$(36)

$$\frac{d}{\mathit{dt}}{N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\mathrm{M}}(t)={K}_{\mathrm{m}}(t){N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\mathrm{M}}(t)+[{K}_{\text{rw}}(t)+{K}_{\text{rc}}(t)]{N}_{{}_{}{}^{36}\mathrm{A}\mathrm{r}}^{\text{CC}}(t)$$(37)where *K*_{m} is the combination of *K*_{mc} and the effective parts of *K*_{mo} and *K*_{mp}. 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 *K*_{mo} and *K*_{mp} are modeled by free parameters

and

${f}_{\text{eff}}^{{K}_{\text{mp}}}$, 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 *K*_{m} is calculated as follows

(38)

In the above mass transfer equations, the five mass transfer rates *K*_{mo}, *K*_{mc}, *K*_{mp}, *K*_{rc}, and *K*_{rw} 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.

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.

## REFERENCES AND NOTES

- ↵
- ↵
- ↵
S. R. Taylor, S. M. McLennan,

*The Continental Crust: Its Composition and Evolution*(Blackwell, 1985). - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
D. R. Hilton, D. Porcelli, Noble gases as mantle tracers, in

*Treatise on Geochemistry*, H. D. Holland, K. K. Turekian, Eds. (Elsevier, 2003), vol. 2, pp. 277–318. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
R. L. Rudnick, S. Gao, Composition of the continental crust, in

*Treatise on Geochemistry*, H. D. Holland, K. K. Turekian, Eds. (Elsevier, 2003), vol. 3, pp. 1–64. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
C. Jaupart, S. Labrosse, J. C. Mareschal, Temperatures, heat and energy in the mantle of the Earth, in

*Treatise on Geophysics*, G. Schubert, Ed. (Elsevier, 2007), vol. 7, pp. 253–303. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵

**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.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).