# Differentiated Archean oceanic crust: Its thermal structure, mechanical stability and a test of the sagduction hypothesis

Alberto Roman, Nicholas Arndt

► To cite this version:

Alberto Roman, Nicholas Arndt. Differentiated Archean oceanic crust: Its thermal structure, mechanical stability and a test of the sagduction hypothesis. Geochimica et Cosmochimica Acta, 2020, 278, pp.65-77. 10.1016/j.gca.2019.07.009. insu-03584767

HAL Id: insu-03584767

https://insu.hal.science/insu-03584767v1

Submitted on 5 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

**Differentiated Archean oceanic crust: its thermal structure, mechanical stability and a test of the sagduction hypothesis**

Alberto Roman, Institut de Physique du Globe de Paris, 1 Rue Jussieu, 75005 Paris,
Now at Jet Propulsion Laboratory, 4800 Oak Grove Dr, Pasadena, CA 91109

Nicholas Arndt, ISTerre, University Grenoble Alpes, 38400 Grenoble, France

**Abstract**

Many recent studies conclude that plate tectonics started about 3 billion years ago in the mid Archean. The transition from a pre-subduction regime to modern plate tectonics is reported to be marked by changes in trace element ratios or isotopic compositions that monitor the rate of growth of the continental crust, the appearance of eclogitic inclusions in diamonds, or an apparent change in the composition of the upper crust. Behind most of these arguments is the hypothesis that, early in Earth history when the mantle was hotter, subduction was intermittent or impossible. If so, a mechanism other than subduction must have created the granitoids that dominate Archean continental crust. One alternative, commonly referred to as sagduction, proposes that the base of thick oceanic crust founders and partially melts to generate granitic magma. Here we evaluate the sagduction process, starting by discussing two crucial concepts: 1) thick oceanic crust is internally differentiated, with hydrated basalt being restricted to the uppermost layers, 2) the generation of granitic magma requires that water and basalt is present in the lower part of the crust or is taken deep into the mantle. We present the results of numerical modelling that demonstrates that when intrusion is taken into account, the lower portion of the crust is well above dehydration temperatures and therefore essentially dry. We show that any deformation within thick, differentiated crust is restricted to the lowermost layers of dry, infertile mafic-ultramafic cumulates that lack the ingredients essential for the generation of granitic magma. Given the implausibility of the sagduction process, we suggest that subduction was the main mechanism that generated granitoid magmas, in the Archean as today.

32

33

**1. Introduction**

Granitoids are abundant in all preserved segments of old continental crust (Bowring & Williams, 1999, Moyen & Martin, 2012). They constitute >60% of most outcrops and in many geophysical surveys are seen to dominate Archean lower crust (Abbott et al., 2013; Artemieva & Thybo, 2013). The presence of hydrous minerals like amphibole and mica provides clear evidence that these granitoids crystallized from hydrous magmas, and numerous experimental studies (e.g. Beard & Lofgren, 1991, Moyen & Stevens, 2006, Zhang et al., 2013) show that voluminous granitic magma is generated only if the source contains several percent of water. Most Archean granitoids belong to the tonalite-trondhjemite-granodiorite (TTG) suite. Depletion of heavy rare earth elements in these rocks is attributed to garnet extraction, either at the source or during crystallization (Jahn et al., 1981, Moyen & Martin, 2012). In mafic assemblages, garnet is stable only at depths >30-40 km (Ito & C. Kennedy, 1971, Johnson et al., 2013), so, to generate Archean continental crust requires that both basalt and water are taken to lower/sub crustal depths. When these conditions are met, partial melting of the basalt produces hydrous magma of intermediate composition that subsequently is transformed by

34

35

36

37

38

39

40

41

42

43

44

45

46

47
48 fractional crystallization, commonly accompanied by wall rock assimilation, into more
49 evolved granitoid magmas.

50 In many models for the formation of Archean granitoid crust (e.g. Bédard, 2006; Johnson et
51 al., 2013; Sizova et al., 2015; Jain et al. 2019, Piccolo et al, 2019), the base of thick oceanic
52 crust or an oceanic plateau either melts in place, or in the “sagduction model”, sinks into hot
53 underlying mantle, then partially melts to generate granitic magma. It is important here to
54 distinguish "crust-forming" sagduction, the process that generates new continental crust, from
55 what might be called "intra-crustal" sagduction, the vertical movement of diapirs that produces
56 the dome-and-basin structure of Archean cratons (Bouhallier et al., 1983, Choukroune et al.,
57 1995, Robin & Bailey, 2009, Thébaud & Rey, 2013, Van Kranendonk et al., 2004) We
58 recognise that vertical displacements of mafic and felsic rocks were common in the Archean
59 but believe that this process merely redistributes material within the crust, rather than being a
60 process that generates continental crust.

61 In this paper we challenge the hypothesis that sagduction, drip-melting, and other processes
62 that act within thick mafic crust, generated Archean continental crust. We base our approach
63 on two main arguments: 1) thick Archean oceanic crust like that generated from hot Archean
64 mantle (McKenzie & Bickle, 1988, Sleep & Windley, 1982) does not have a uniform basaltic
65 composition from top to base, as assumed in most models of Archean geodynamic processes
66 (Bédard, 2006, Johnson et al., 2013, Sizova et al., 2015, van Hunen & van den Berg, 2008,
67 Nebel et al., 2019). Instead this crust was strongly differentiated, both in terms of lithology
68 and water content; 2) the lower part of thick Archean oceanic crust consisted of ultramafic
69 cumulates that were essentially anhydrous and contained neither the basaltic component nor
70 the water needed to generate large volumes of granitic magma.

## 2. The Differentiated Structure of Thick Oceanic Crust

### 2.1 Lithological variations

74 Many observations, and solid petrological theory, demonstrate that thick piles of mafic
75 material such as oceanic crust or oceanic plateaus are not vertically homogeneous: instead,
76 they are differentiated, with basaltic uppermost layers, gabbroic intermediate layers, and lower
77 layers that consist of mafic-ultramafic cumulates.

78 Observational evidence of strong differentiation is found in intrusions and igneous complexes
79 ranging from simple units like 1-10m thick komatiite lava flows (Arndt et al., 2008) or the
80 300m thick Palisades sill (Shirley, 1987) through to the 40 km thick Ontong Java plateau
81 (Fitton et al., 2004), all of which are internally differentiated. Likewise, the lithology of the 9
82 km thick Bushveld Complex changes downwards from diorite at the top through pyroxene
83 and pyroxene-plagioclase cumulates to basal olivine cumulates (Cawthorn, 1996; Wilson
84 2012).

85 In making reference to intrusions in crustal settings, we do not suggest that they are direct
86 analogues of thick Archean oceanic crust. Rather we use them to illustrate a basic petrological
87 principle, namely that as mafic magma intrudes, because of gravitational settling or in-situ
88 growth of early-formed mafic minerals, it differentiates into lower mafic-ultramafic lower
89 portions and more evolved upper portions. Although the mechanism that causes this
90 differentiation may be subject to debate, there is no doubt that the process operates in all but
91 the smallest, most rapidly cooled intrusions.

92 These processes act at a larger scale to differentiate thick oceanic crust and oceanic plateaus.
93 To generate a large thickness of mafic material requires a high degree of mantle melting which
94 produces a large volume of picritic magma (McKenzie & Bickle, 1988). As schematically
illustrated in Figure 1a, two separate processes then act in concert to produce lower ultramafic layers and upper basaltic layers. First, the primary picritic magma is relatively dense and when it enters the crust, it tends to be trapped as intrusions near the crust-mantle boundary (Cox, 1980). Lithospheric mantle is absent during the growth of oceanic crust and temperatures at the Moho are close to the liquidus temperature. Most oceanic plateaus form at or near spreading centres (Neal et al. 2015) and lithosphere will also be very thin to absent. Incoming magmas therefore do not freeze but cool slowly. Fractional crystallization produce a less dense, evolved basalt, which ascends towards the surface, and a portion of mafic-ultramafic cumulates, which remain at or near the base of the crust. Second, as the crust thickens, density influences the level at which incoming magmas intrude. Dense, olivine-rich, picritic magmas are trapped at lower levels while less dense, more evolved basaltic magmas ascend to higher levels.

We illustrate these processes using Herzberg et al.'s (2010) compilation of the compositions of Archean mafic volcanic rocks which has been widely used in previous studies focusing on Archean geodynamic regimes (e.g. Johnson et al. 2013). The Precambrian tholeiitic basalts in the data-set are not primary melts, with MgO contents between 10 and 14 wt%. The compositions of their primary magmas is found by progressively adding olivine until equilibrium with a fertile mantle peridotite is reached (see Herzberg and Asimow, 2008 for details on the method). Using this approach, we calculated that up to 35 wt% olivine has to be added to the original samples in order to obtain compositions in equilibrium with the mantle (Fig. 1b). Moreover, there is a strong correlation between the MgO content in the primary melts and the amount of olivines that has to be added. Whatever the mechanism producing this trend, this observation shows that picritic magmas undergo large amounts of fractionation before eruption.

The method of Herzberg et al. (2010) allows us also to estimate the compositions of the fractionated olivines and that of the residual mantle (that is, the composition of the mantle left behind the extraction of the primary liquid). These compositions can then be used to calculate the densities of both materials, as illustrated in Figure 1c. The minerals in the cumulates crystallized from melts that contained more Fe than the primary melt, and their densities (3330-3350 kg m⁻³) are slightly higher than that of the residual mantle (3310-3350 kg m⁻³). The ultramafic cumulates in the lower crust contain very little Al and therefore very little garnet, the mineral that imparts high density to eclogite. For this reason, the cumulate densities are significantly less than that of eclogite (> 3370 kg m⁻³) which is often considered to drive sagduction (e.g. Sizova et al., 2015). The relatively small density contrast at the Moho is an important parameter for the stability of the whole crustal column.

Herzberg et al. (2010) proposed that mantle potential temperatures reached a maximum (Tp = 1600°C) towards the end of the Archean and that partial melting in ambient mantle produced picritic magmas containing 20% wt% MgO. The observations outlined above suggest that fractional crystallization of this magma would yield 1) evolved magmas containing 10-14 wt% MgO that erupt as lava at the surface, 2) 30% of olivine cumulate that resides in the lower portion, and 3) an intermediate portion of gabbroic composition. We used this conceptual model of Archean crust, shown in Figure 1a, when testing sagduction in our numerical experiments.
## 2.2 Variations in water content in thick differentiated crust

Not only the lithology but also the water content can be predicted to change downwards. Even in relatively thin modern ocean crust, strongly hydrated basaltic lava is restricted to the uppermost layer, and underlying gabbros and mafic-ultramafic cumulates have low water contents (Hacker et al., 2003a; Jarrard, 2003). Water is rare to absent for two reasons. First, the lower part of the crust is heated by conduction from underlying mantle and advection of heat in incoming magmas. This results in temperatures above about 500-700°C, the levels where the dominant water-bearing phases such as serpentine, lawsonite and chlorite break down (Faccenda, 2014). In thicker hotter Archean oceanic crust, a much larger portion will be above the dehydration temperatures, as discussed later. Second, seawater normally penetrates only into the uppermost part of the crust. At magma-poor slow-spreading ridges, active faulting allows water to circulate down to depths of several kilometers, but in magma-rich fast-spreading ridges, the penetration is minimal (Faccenda, 2014). Most of the water that is released during subduction enters the crust and lithosphere only when the plate bends and fractures as it enters the subduction zone. Thick, extant oceanic crust and oceanic plateaus of the type considered in sagduction models most closely resemble fast-spreading ridges and hydrothermal circulation will affect only their uppermost levels.

Many numerical and petrological models that simulate Archean geodynamic processes (e.g. (Bédard, 2006, Johnson et al., 2013, Moore & Webb, 2013, van Hunen & van den Berg, 2008) consider that the entire Archean oceanic crust is made up of hydrated basalt. To arrive at this structure, the authors propose that basaltic lava becomes hydrated through interaction with seawater at the sea floor, then is transported to the base of the crust by burial under subsequent flows. We test this hypothesis in our modelling.

## 3. Formation and Instability of the Archean Crust

The sagduction mechanism is based on the assumption that the lower part of >40 km thick mafic crust becomes unstable, founders into the underlying mantle, and then partially melts to generate granitic material. The instability is related to the high density of the lower crust, either because it converts to eclogite, or because of an elevated Fe content.

In our study we modelled the thermal evolution of a growing crustal pile in which both intrusion and eruption occur. This allows us to constrain both the water content at depth and to draw some inferences on the density and viscosity profile of the crust. The latter parameter exerts an important control on the stability of the lower crust and on the coupling between the upper and lower crustal layers. Then we evaluated the development of instabilities and the onset of foundering within differentiated mafic-ultramafic crust using finite element code. Reducing the study of the sagduction process to two simpler subproblems (a thermal problem and a mechanical one) is motivated by the fact that the thermal structure at the onset of the instability is set by its previous accretion history. Even though this strategy is undoubtedly a simplification of the complex and poorly understood processes of Archean crust formation, it allows us to clearly illustrate the limits of the sagduction model.

### 3.1 Thermal models: methods

We model the formation of Archean crust assuming a scenario similar oceanic plateaus, where the crust is constructed by repetitive eruptions and injections of sills. The total flux of magma from the mantle is partitioned between intrusive rocks and erupted lavas according to the ratio $IE = \text{volume intrusive rocks}/\text{volume erupted lavas}$.
Global compilations of volcanic output (Crisp, 1984, White et al., 2006) indicate that during the growth of both continental and oceanic crust, most magma is intruded, not erupted. Estimates of the proportion of intrusive to extrusive rocks range from about 5:1 in thin ocean crust, to 10:1 in thicker sequences such as the Kerguelen plateau. Seismic tomography has revealed that even in regions of continental break-up, such as the rifted North Atlantic margin, the intrusive to extrusive ratio is at least equal to 1 (White et al. 2006).

We assume that heat transport in the crust is controlled by diffusion and neglect the effects of hydrothermal circulation for the reasons discussed in section 2.2. We further assume that diffusion proceeds mainly in the vertical direction, which is a good approximation if the intrusions forming the crust have low aspect ratios (i.e. are thin sill-like intrusions).

The 1-D heat diffusion equation is:

$$ \rho C_p \frac{\partial T}{\partial t} + \rho L \frac{\partial M}{\partial t} = k \frac{\partial^2 T}{\partial z^2} \quad (1) $$

where ρ is density, C<sub>p</sub> is the heat capacity, T is the temperature, L is the latent heat of crystallization, M is the melt fraction, t is time, and k is the thermal conductivity. By using the chain rule, this equation can be rewritten as

$$ \rho C_p \frac{\partial T}{\partial t} + \rho L \frac{\partial M}{\partial t} \frac{\partial T}{\partial t} = k \frac{\partial^2 T}{\partial z^2} \quad (2) $$

which requires the definition of M(T). The details of the dependence of melt fraction on temperature vary with pressure and water content (Müntener et al. 2001), and are important for quantifying the compositional evolution of the system. However, they have only a minor effect on the bulk energy balance (Annen, 2017) so that we can safely assume a simple linear form such that:

$$ \begin{array}{ll} \text{212} & \begin{array}{l} \text{M=0,} \\ \text{M=(T-T_s) / (T_l-T_s)} \end{array} \\ \text{213} & \begin{array}{l} T<T_s \\ T_l<T<T_s \end{array} \\ \text{214} & \begin{array}{l} \text{M=1} \end{array} \\ \text{215} & \begin{array}{l} T>T_l \end{array} \end{array} \quad (3) $$

Equations 2 and 3 are solved with an implicit finite differences scheme which has been benchmarked with solutions of Furlong et al. (1991) for cooling of a single magma batch with uniform background temperature.

As magma rises to the surface, it is trapped at density or rheological interfaces where it cools, differentiates, and eventually resumes its ascent. This process is likely to repeat itself at many crustal levels, so that during accretion many storage zones are simultaneously active at different depths. A complete exploration of these thermal histories would require that we consider a prohibitively large number of polybaric crystallization paths and, for the sake of simplicity, we consider only rapid separation of liquid from cumulates in a single chamber, so that intrusion and eruption of lava are contemporaneous. At each accretion event, part of the magma coming from the mantle is emplaced as an intrusion with a prescribed thickness at liquidus temperature, while the rest erupts and is instantaneously cooled to zero degrees. The respective thickness of the intrusion and of the erupted material are determined by the IE ratio. For instance, if IE=2 and each accretion event add 300 m of mantle derived material to the crust, 100 m are emplaced as lavas at the surface, and 200 m as a sill at depth. Grid points are displaced below the intrusion by an amount corresponding to the total thickness of the newly
233 accreted material, so that the size of the numerical grid increases as the crust grows. From the 234 density considerations outlined above, magma is expected to stall at the Moho, so we assume 235 that intrusions are emplaced at the crust-mantle interface, whose depth increases as the crust 236 grows. We stop the calculations when the Moho reaches a depth of 40 km, which corresponds 237 to the supposed thickness of the Archean crust.

238 We apply constant temperature boundary conditions at the top and bottom of the domain, with 239 the value at the bottom corresponding to the temperature at the base of the lithosphere. We 240 varied this value between 1300-1600 °C which encompasses the range between present day 241 mantle potential temperature and the Herzberg et al.'s (2010) estimate for the Archean mantle.

242 A self-consistent evaluation of lithospheric thickness and temperature would require that we 243 also account for mantle convection. Here, we consider that the crust grows rapidly (typically 244 a few millions of years), as observed in oceanic plateaus, so the temperature at the bottom can 245 be taken as fixed during the accretion process. Mantle convection controls the temperatures 246 variations of the lithosphere on longer time-scales, but the water content of the crust is set 247 during accretion, during which hydrated lavas are buried by subsequent flows. The purpose 248 of this section is to constrain the water content at depth, so that we use a constant temperature 249 boundary condition and neglect convection. Using this approach, the effects of different 250 temperatures at the base of the lithosphere and variations in its initial thickness can be studied 251 by varying these parameters.

252 Similar considerations apply to radioactive decay. Due to the short time scales considered, we 253 neglect radioactive processes in our calculations. However, it is clear that their effect is to 254 increase the estimated temperatures.

255 The numerical values of all parameters are listed in Table 1.

256

<table><thead><tr><th>PARAMETER</th><th>SYMBOL</th><th>VALUE</th></tr></thead><tbody><tr><td>THERMAL CONDUCTIVITY</td><td>k</td><td>2.5, 2.0* (W K<sup>-1</sup> m<sup>-1</sup>)</td></tr><tr><td>DENSITY</td><td>ρ</td><td>2800, 3200* (kg m<sup>-3</sup>)</td></tr><tr><td>LATENT HEAT</td><td>L</td><td>4.2 10<sup>5</sup> (J kg<sup>-1</sup>)</td></tr><tr><td>HEAT CAPACITY</td><td>C<sub>p</sub></td><td>1.5 10<sup>-3</sup> (J kg<sup>-1</sup> K<sup>-1</sup>)</td></tr><tr><td>SOLIDUS TEMPERATURE</td><td>T<sub>s</sub></td><td>1100 (°C)</td></tr><tr><td>LIQUIDUS TEMPERATURE</td><td>T<sub>l</sub></td><td>1300 (°C)</td></tr><tr><td>SURFACE TEMPERATURE</td><td>T<sub>t</sub></td><td>0 (°C)</td></tr></tbody></table>

257

258 **Table 1:** Fixed parameters used in the calculations. * correspond to the value assigned to mantle rocks. Rock physical properties values are derived from Chapman, 1992, and Bohrson and Spera, 2001.

259

260

261

262

263

## 3.2 Results

### 3.2.1 Temperature variation within a growing pile of lava flows
In the first step we calculate the thermal evolution when all magma is erupted (Figure 2), which corresponds to a “heat-pipe Earth” (e.g. Moore et al., 2013). In this case, the heat flux associated to magma transport is equal to zero, because lavas are emplaced at the surface at 0°C. Lava eruption buries previous flows and cold material is displaced to the base of the crust. As shown in Figure 2 this displacement outpaces diffusion, so that cold lavas reach 40 km depth with minor diffusive heating. As a consequence, temperatures throughout the crust are extremely low, everywhere less than 400°C. Higher emplacement rates (1 Myr in the figure) result in even lower temperatures.

### 3.2.2 The effect of intrusion into the lava pile

In Figure 3 we illustrate the effect of intrusion. The approach is similar to that of Annen and Sparks (2002), except that we allow part of the magma to escape and erupt at the surface, contributing to crustal accretion. While intrusion acts to increase the temperature in the lower crust, extrusion acts in the opposite way via downward advection of cold material. We study the effects of these two competing process by varying the intrusive to extrusive ratio. Michaut and Jaupart (2006) have shown that the main controls on the cooling rate of an intrusion is the heat flux associated with magma injection, which scales as the ratio between intrusion thickness and injection time interval. Thus, we fix the intrusion thickness and vary the period over which the crust is formed, resulting in different magmatic fluxes.

Figure 3 shows the temperature evolution in a crustal column that formed over a time interval of 1 Myr, and in which the magmatic flux is partitioned equally between intrusive and extrusive rocks (IE ratio=1). In this example, the initial crust thickness is set to 1 km with an initial lithospheric thickness of 40 km, and a bottom temperature kept at 1400 °C. The intrusions are emplaced at the Moho, whose depth increases with time as the crust grows. At the beginning the crust is thin, so that the Moho is shallow and cold. In this phase, the heat brought in by the first pulses of magma is removed efficiently due to the high temperature contrast. However, as the crust thickens, the thermal boundary layer overlying the injection zone grows as $\delta=(\kappa t)^{1/2}$ (where $\kappa$ is the thermal diffusivity and t is time), so that heat diffuses away less and less rapidly. Thus, temperatures increase progressively both at the intrusion level and in the overlying crust. As it can be seen from the red curves, temperatures at depths greater than 20 km in the lower part of the crust are well above 500-700 °C the limit of stability of water-rich phases such as chlorite, lawsonite or serpentine. Water-poor minerals like amphibole or mica persist to higher temperatures, but experimental studies have shown that the amount of melt generated during the breakdown of these minerals in very small (e.g. Aranovich et al. 2014).

Figure 4a shows the final thermal profiles for calculations in which we varied the magmatic flux such that a 40 km crust formed in 1, 5 and 10 Myr. When the crust formed rapidly, a 15 km thick region of the lower crust has near-constant temperatures close to that of injection. On the other hand, if the crust formed in 10 Myr, temperatures at the Moho reach a maximum value slightly below the solidus and decline almost linearly with decreasing depth. However, even in this case, temperatures in the lower 15 km of the crust remain above the dehydration temperatures. This is readily apparent in Figure 4b where Moho temperature is plotted as a function of depth.

An interesting feature, which could be expected on physical grounds, is that the situation is the opposite in the upper crust. High fluxes correspond to rapid extrusion of lavas at the
312 surface, so that cold flows are rapidly advected downward, reducing temperatures in the upper
313 10 km. This results in a very steep gradient between 10 and 30 km. Conversely, when the crust
314 forms over ten million years, flows are buried at a slower pace, and diffusion of heat from
315 below acts to increase the temperature in the upper 15-20 km.

316 Figure 5 shows the effect of changing the intrusive to extrusive ratio for crust that formed in
317 5 Myr. Doubling the intrusive to extrusive ratio results in temperatures that are about 200
318 degrees higher over the whole crustal column. The small amount of lavas erupted at the
319 surface is not sufficient to decrease temperatures in the upper crust, while the increased
320 volume of intrusive rocks results in a higher heat flux due to injections.

321 In the calculations shown here the bottom of the grid is maintained at a temperature of 1400
322 °C. The effect of this parameter and the thickness of the initial lithosphere were studied in few
323 calculations (not shown) and, in the range of parameters explored, were found not to play an
324 important role. This is because the heat flux associated with magmatic intrusions is several
325 orders of magnitudes higher than the steady state conductive heat flux from the mantle for
326 crust forming in time interval < 10-20 Myr. A crust that forms over longer periods of time
327 will have geotherms controlled by the depth and temperature of the bottom of the lithosphere.

328 In such a case, temperatures of less than 600 °C at the Moho are obtained if the lithospheric
329 thickness is more than 100 km. However, such values of lithospheric thickness are
330 inconsistent with the high magmatic fluxes needed to generate a thick crust.

331 In summary, in all the models explored in this paper, temperatures in the lower crust are
332 significantly higher than the upper limit of breakdown of water-rich hydrous phases. We
333 conclude that if Archean crust grew by a combination of intrusion and extrusion, the lower
334 30-20 km were dry.

335
336
337 **3.3 Gravitational instabilities and the generation of felsic crust**

338 Two important conclusions stem from the discussion and calculations presented in the
339 previous sections. The first is that crust with a bulk composition similar to that of a picritic
340 melt will have a lower portion composed of ultramafic cumulates (i.e. olivine ± pyroxene) and
341 mafic lithologies containing plagioclase will be restricted to the upper portion of the crust.

342 The second is that high temperatures associated with intrusive activity will cause dehydration
343 and drive off water from the rocks in the lower crust, so that water will also be confined to the
344 upper crustal layers. Given that plagioclase and water are necessary ingredients to produce
345 voluminous granitoid magma, the growth of continental crust requires that rocks that normally
346 reside in the upper part of the oceanic crust are brought to the base of the crust or into the
347 underlying mantle by some physical process. As it can be seen in Figure 1c, the ultramafic
348 cumulates of primary magmas have densities that are greater than the residual mantle, with a
349 density contrast that is exacerbated for picritic liquids. These cumulates are therefore
350 potentially unstable and may founder into the underlying mantle. If this dense cumulate layer
351 generates gravitational instabilities, it is possible that the upper wet, basaltic layer could be
352 entrained to the mantle, thus meeting the requirements for continental crust formation.

353
354 In our simplified mechanical model of differentiated mafic-ultramafic crust, an upper gabbroic
355 layer overlies a lower ultramafic cumulate layer. The crust is underlain by lithosphere whose
356 density is lower than that of the cumulates, due to its lower iron content. This density contrast
357 drives the gravitational instability. The upper crustal layer is less dense than the cumulates
due to the presence of plagioclase (we ignore at this stage the presence of even less dense
hydrous phases).

Our analysis closely follows previous studies (Conrad & Molnar, 1997, Houseman & Molnar, 1997, Jull & Kelemen, 2001) focusing on the development of Rayleigh-Taylor instabilities in the crust and lithosphere. The main difference here is that our target is not to estimate the timescales over which an instability develops, but rather to understand under which conditions foundering of ultramafic cumulates may entrain an upper buoyant layer down in the mantle. Because coupling between different layers is expected to become efficient during the late non-linear phase of instability growth, we use a finite element code outlined in the Appendix.

A key rheological feature of differentiated Archean crust is the continuity of the effective viscosity at the Moho, because both the lower crust and upper mantle are composed mainly of olivine and the lower crust lacks plagioclase. The effective viscosity is strongly temperature-dependent which requires that we define the Archean geotherm. However, when using the conventional exponential formulation for the variation of viscosity with depth, the effects of the thermal structure of the crust can be studied by varying a single dimensionless parameter γ. If this parameter is large (>>1) all the layers have the same viscosity; if small, viscosities increase rapidly above the Moho. The results are then discussed a posteriori in relation with the findings of Section. 3.2. The other dimensionless parameter controlling the dynamics of the system is the buoyancy number B, which expresses the ratio of the density differences between the three layers. Models with high B correspond to the case in which the upper gabbroic layer is strongly buoyant compared to the density contrast driving the instability. Explicit definitions of both dimensionless parameters are given in the Appendix.

To illustrate the effects of buoyancy alone, Figure 6a shows the time-lapse evolution of an experiment in which all layers have the same viscosity (γ >> 1). A buoyancy number B of 3.0 corresponds to the case in which densities are 3150 kg m⁻³, 3350 kg m⁻³, 3300 kg m⁻³, for the upper crust, cumulates and mantle lithosphere respectively. The cumulate layer starts to thicken on the left side where an initial perturbation is imposed. The upper layer, which initially is coupled to the underlying unstable layer, thickens to the left and thins to the right. However, once the instability is fully formed, the upper and lower layers rapidly decouple and, after foundering of cumulates, the original thickness of the upper layer is restored to almost its initial value.

Figure 6b tracks the evolution of the different interfaces. As shown by the violet solid line, the perturbation of the cumulate layer initially grows very slowly and then abruptly accelerates when the displacement is only a small fraction of the cumulate layer thickness (when dimensionless time is about 384). The instability of the lower crustal layer starts to significantly affect the upper layer when the perturbation amplitude is 0.15-0.20. The upper crust (black solid line) thickens quickly to slightly more than 0.25 until the instability detaches, then declines rapidly to less than 0.1. Tracking the deformation in the thinned portion of the domain (on the right side of the computational domain, marked as dashed dotted line) shows that maximum thinning is about half the maximum thickening, suggesting that the thickened portion of the crust is very small and localized. We conducted a series of experiments in which we varied only the buoyancy number B (by increasing the density of the cumulate layer) between 3 and 6, and kept all the other parameters constant. For high buoyancy numbers B >3, the maximum thickening of the upper crust is <0.8 and decreases
rapidly between B= 3 and B = 4 (Figure 7). Figure 7 shows that if the upper crustal layer is strongly buoyant (B>5), its deflection will be of the order of 20% of the cumulate layer.

In Figure 8 we show the effect of the exponential decrease in viscosity in the vicinity of the Moho. We set the rheological parameter γ at 0.65, which corresponds to a 20 km thick cumulate layer and a characteristic length-scale of viscosity decay of 13 km. The dominant wavelength is smaller than when the viscosity is constant: previous studies (Conrad & Molnar, 1997), and our calculations, show that when γ < 1, the characteristic wavelength and the timescale of instability development scales with characteristic length-scale of viscosity decay and not the thickness of the unstable layer. Consequently, only the lower part of the cumulate layer becomes unstable and the upper part of the crust deforms only slightly (as seen in the light/dark grey stripes in Figure 8). The calculation was performed for the limiting case in which B=1; i.e. the upper crustal layer and the mantle have the same density. This shows that if γ < 0.65, deformation in the upper crust is negligible, independent of rock density.

## 4. Discussion

### 4.1 Application to the Archean crust

The results of thermal and dynamical calculations place some constraints on the sagduction process in the Archean. The thermal structure and mechanical behavior of the crust are linked by viscosity, which is controlled by temperature and which exerts a dominant control on the coupling between the lower and upper crust.

Figure 4 illustrates two end-member scenarios. In the first, temperatures are nearly constant in the whole lower crust, which corresponds to the case in which the crust formed over a short period of time (for instance 1 Myr, blue line in Figure 4a). In this situation we can use the results for a crust with constant viscosity (Figures 6 and 7) to evaluate whether foundering of ultramafic cumulates at the base of the crust is a viable mechanism to generate felsic continental crust. In this case, the main force that prevents entrainment of the upper hydrated and plagioclase rich layer is buoyancy, as expressed by the buoyancy number B.

The maximum deformation in the upper layer is small (< 25% of the lower crust thickness) when the B number is higher than 4. In Figure 1b, the density of the cumulate layer is between 3350-3330 kg m⁻³, while that of the residual mantle is only slightly less, at 3310-3340 kg m⁻³, resulting in a maximum density difference driving the gravitational instability of 40 kg m⁻³. Under these conditions, upper crustal rocks whose density is less than 3190 kg m⁻³ will experience only minor deformation during foundering of the lower crust. The density that we assumed for the upper crust is higher than the density of rocks containing plagioclase (<200 kg m⁻³) and much higher than typical values of fully hydrated basaltic rocks (2500 kg m⁻³).

Now consider the case in which the crust forms over longer periods of time (5-10 Myr in Figure 4). In this situation, temperature in the crust is nearly linear and the viscosity variation is exponential. Temperatures at the Moho are between 1000-1200°C, with corresponding temperature gradients between 25-30°C. In these conditions, the characteristic length-scale of viscosity decay L is less than 5 km, which is a fraction of the thickness of the unstable lower crust, so that γ << 1 (see the Appendix for the details of the values used for the estimation).

In these conditions the two layers completely decouple and the upper hydrated basalts are left undeformed.
We note that foundering of lower crustal cumulates could cause upwelling of the asthenosphere and that upward advection of heat could lead to minor melting in the crust. However, given the absence of water, the amount of melt would be low. In addition, this melt forms at shallow levels and low pressure where garnet is not stable. Melt produced under these conditions would not have the depletion of heavy rare earth elements that characterizes Archean granitoids.

## 4.2 Limitations of the approach

The thermal models shown here have been focusing on a scenario in which the Archean crust is formed by vertical stacking of sills and flows, in a similar fashion to oceanic plateaus. In this situation sagduction is a viable mechanism only if no intrusion occurs. However, the possibility remains to be explored that 40 km Archean crust formed by lateral accretion of dikes, such as in modern oceanic ridges. In this context, it seems likely that the depth at which hydrated phases are present is set by hydrothermal circulation. In fact, in this scenario there is no mechanism advecting hydrated material from the surface to deep levels (contrary to our models in which hydrated flows are buried by subsequent eruption). Thus, on one hand, it is possible that in the case of extremely fractured Archean crust, water could have been present at large depths. On the other hand, as noted in Section 2.2 such crust is likely characterized by high temperatures as in modern fast spreading ridges (Fontaine et al., 2017). High temperatures would promote viscous flow of its lower portion which, in turn, would close cracks and pores. It seems therefore more likely that water would be, even in this case, confined to the uppermost levels of the crust. Nevertheless, given the critical role played by water in the generation of granites, it is clear that such process must be studied in a quantitative way using numerical techniques that model hydrothermal circulation and viscous deformation of the lower crust in consistent way, as it is done for modern ocean ridges (e.g. Theissen-Krah et al., 2016, Fontaine et al., 2017).

## 4.3 Implications

We have shown that the proposed sagduction models fail to produce large volumes of felsic magma with the composition of Archean granitoids. What of the alternatives? Most models for the generation of Archean granitoids assume that the entire crust had a basaltic composition like that of the typical Archean tholeiites. This is inconsistent with arguments that have been proposed and accepted by petrologists for many decades. These tholeiites have evolved compositions that cannot be generated by melting of normal mantle peridotite and are very different from the picrites generated in hot Archean mantle. For the reasons outlined in section 2, we believe that rocks with basaltic compositions would be confined to the upper portions of the crust. Basalts also have low densities, particularly if hydrated. Due to their position at the top of the crust and their low densities, they are most unlikely to spontaneously founder or drip into the mantle.

It remains the possibility that down-going portions of long-lived mantle convection cells could drag the crust down into the mantle. This process might affect to lower part of the crust or the underlying lithosphere, but because of the decoupling between lower and upper crust that appears in all our models, it will not affect the upper layers. However, if down-dragging of the lower crust persisted, the upper crust would thicken and it could heat up to reduce its viscosity. In this case it might be dragged down by foundering cumulates. One would expect that during its downward displacement, the increase in temperature would destabilise the
497 hydrous minerals and the foundered crust would become relatively dry. This complex
498 interplay must be addressed by future studies.

499 Piccolo et al. (2019) have suggested that heat advected into the crust in magmas from the
500 mantle might lead to basalt melting and the generation of granitoid magma. This model also
501 runs into problems if we accept our concept of a differentiated crust in which hydrated basalt
502 is limited to the uppermost levels. If the magma intrudes into dry cumulates deep in the crust,
503 little to no melting will take place; and if the magma were able to ascend to shallow levels
504 where hydrous minerals were stable, low temperatures would restrict the amount of melting,
505 and any melt would form at low pressure. In the absence of garnet, the melt would not have
506 the HREE depletion that characterizes Archean granitoids.

507 Our preferred interpretation is that subduction remains the most probable process for
508 producing continental crust, in the Archean as on the modern Earth. But how do we respond
509 to the argument that subduction could not have operated in a hot Archean Earth? Weller et al.
510 (2019) found that thick Archean oceanic crust was less, not more, buoyant than modern crust
511 and concluded that subduction was indeed possible early in Earth history. Furthermore, recent
512 estimates of Archean mantle temperatures are less extreme than those made in the past. Ganne
513 and Feng (2017) re-evaluated the potential temperatures of Herzberg et al.'s (2010) non-arc
514 basalts and concluded that the difference between Archean and modern ambient mantle was
515 no more than 150K. Not considered in this appraisal is the likelihood that many greenstone-
516 belt basalts erupted in plume-generated oceanic plateaus (e.g. Abouchami et al., 1990; Storey
517 et al., 1991) and do not represent normal oceanic crust; if so, the potential temperatures
518 defined by Ganne and Feng (2017) are likely to define maximum values for the convecting
519 upper mantle. An independent estimate from eclogite xenoliths (Aulbach and Arndt, 2019)
520 also indicates a small difference of 100-150K.

521 Archean subduction would nonetheless have differed from modern subduction even if
522 temperatures in the Archean mantle were only 100-150K higher than today's. Although not
523 reaching the dimensions envisaged for a hot Archean mantle, ocean crust generated from
524 Archean mantle was thicker than modern ocean crust. In addition, temperatures in the lower
525 part of the crust and in the subduction zone were higher, and both crust and mantle probably
526 were less oxidised. The reported absence of features like those in modern subduction settings
527 (e.g. Korenaga, 2013; Bédard, 2017) may result from these differences, not from an absence
528 of subduction.

529 Many authors accept that subduction may have been possible in the Archean but propose that
530 this subduction was intermittent (e.g., O'Neill et al. 2007; Moyen and Van Hunen, 2012) or
531 that the slabs broke off before penetrating deep into the mantle (van Hunen & van den Berg,
532 2008; Sizova et al., 2010; Gerya, 2010). Other authors propose intermittent subduction during
533 slab roll-back, perhaps triggered by the ascent of plumes (Gerya et al., 2015). In modern
534 subduction zones, the aqueous fluids that flux melting in subduction zones, are released at the
535 early stages of subduction, before the plate reaches a depth of about 100 km. The operation of
536 multiple, short-lived, shallow-penetrating subducting slabs would be totally capable of
537 generating the granitoids that make up Archean continental crust.

538 Whether the onset of subduction signals the start of plate tectonics remains a subject of debate.
539 If it is accepted that subduction is one of the key features of plate tectonics, and if the oldest
540 preserved granitoids represent portions of continental crust, then it could be argued that the
541 absence of a credible alternative to subduction for generating continental crust indicates that
542 plate tectonics operated from the early Archean. Alternatively, a transitional stage could be
543
544 envisaged in which small segments of oceanic slabs subduct, without being parts of rigid
545 plates created at spreading centers.

546
547 **5. Conclusion**
548 In this study we tested the sagduction model for the production of Archean granitic crust.
549 Petrological arguments and observations of natural analogues indicate that ocean crust,
550 whether generated at spreading centers or as oceanic plateaus, was differentiated, with
551 hydrous basalt confined to the uppermost parts, and lower portions consisting of mafic-
552 ultramafic cumulates. When heat from the intrusions is taken into account, our thermal models
553 predict that the lower part of the crust was dry. Given our proposed thermal and buoyancy
554 structure of the crust, we showed that gravitational instabilities of the lower crust are unable
555 to bring water and basaltic material to pressures greater than 1 GPa, which is the basic
556 requirement for the generation of felsic magma matching the chemical features of the
557 preserved Archean continental crust. We conclude that the sagduction model, as presented in
558 previous studies, is not consistent with the generation of granitic crust, and propose that such
559 crust was generated in analogs of modern subduction zones.

560
561 **Acknowledgements**
562 AR acknowledges support from Institut de Physique du Globe de Paris and Jet Propulsion
563 Laboratory, NASA Postdoctoral Program fellowship. NTA acknowledges support from
564 Institut des Sciences de la Terre of University Grenoble Alpes.

565
566
567 **References**
568 Abbott, D. H., Mooney, W. D. & VanTongeren, J. A. (2013). The character of the Moho and
569 lower crust within Archean cratons and the tectonic implications. Tectonophysics 609,
570 690-705.
571 Abouchami, W., Boher, M., Michard, A. & Albarède, F. (1990). A major 2.1 Ga old event of
572 mafic magmatism in West Africa. J. Geophys. Res. 95, 17605-17629.
573 Aranovich, L.Y., Makhluf, A.R., Manning, C.E., Newton, R.C. (2014) Dehydration melting
574 and the relationship between granites and granulites. Precambrian Research 253, 26-37.
575 Annen, C. & Sparks, R. S. J. (2002). Effects of Repetitive Emplacement of Basaltic Intrusions
576 on Thermal Evolution and Melt Generation in the Deep Crust. Earth and Planetary
577 Science Letters 203, 937-955.
578 Annen, C. (2017). Factors affecting the thickness of thermal aureoles. Frontiers in Earth
579 Science, 5, 82.
580 Arndt, N. T., Barnes, S. J. & Lesher, C. M. (2008). Komatiite. Cambridge: Cambridge
581 University Press.
582 Artemieva, I. & Thybo, H. (2013). EUNAseis: A seismic model for Moho and crustal structure
583 in Europe, Greenland, and the North Atlantic region. Tectonophysics 609, 97-153.
584 Aulbach, S. and Arndt, N.T. (2019) Eclogites as palaeodynamic archives: Evidence for warm
585 (not hot) and depleted (but heterogeneous) Archaean ambient mantle. Earth Planetary
586 Science Letters, 505 162-172.
587 Beard, J. S. & Lofgren, G. E. (1991). Dehydration melting and water-saturated melting of
588 basaltic and andesitic greenstones and amphibolites at 1, 3, and 6.9 kb. J. Petrol. 32, 365-
589 401.
590 Bédard, J. H. (2006). A catalytic delamination-driven model for coupled genesis of Archaean crust and sub-continental lithospheric mantle. *Geochimica et Cosmochimica Acta* 70, 1188-1214.

591
592 Bédard, J. H. (2017). Stagnant lids and mantle overturns: Implications for Archaean tectonics, magmagenesis, crustal growth, mantle evolution, and the start of plate tectonics. *Geoscience Frontiers*.

593
594
595 Bohrson, W. A., & Spera, F. J. (2001). Energy-constrained open-system magmatic processes II: application of energy-constrained assimilation–fractional crystallization (EC-AFC) model to magmatic systems. *Journal of Petrology*, 42(5), 1019-1041.

599
600
601 Bouhallier, H., Choukroune, P. & Ballèvre, M. (1983). Diapirism, bulk homogeneous shortening and transcurrent shearing in the Archean Dharwar craton : the Holenarsipur area, South India. *Precambrian Research* 63, 43-58.

602
603 Bowring, S., A, & Williams, I. S. (1999). Priscoan (4.00–4.03 Ga) orthogneisses from northwestern Canada. *Contributions to Mineralogy and Petrology* 134, 3-16.

604
605 Burov, E., & Gerya, (2014) T. Asymmetric three-dimensional topography over mantle plumes. *Nature* 513(7516), 85-89.

606 Cawthorn, R. G. (ed.) (1996). Layered Intrusions. Amsterdam: Elsevier

607
608 Chapman, D. S. (1992). Thermal state of the continental lower crust. *Continental lower crust, 179-199*.

609
610
611 Choukroune, P., Bouhallier, H. & Arndt, N. T. (1995). Soft lithosphere during periods of Archaean crustal growth or crustal reworking. In: Coward, M. P. & Ries, A. C. (eds.) *Early Precambrian Processes: Geol. Soc Lond, Spec Publ.*, 67-86.

612
613
614 Condie, K. C. & Kröner, A. (2008). When did plate tectonics begin? Evidence from the geologic record. In: Condie, K. C. & Pease, V. (eds.) When did plate tectonics begin on planet Earth? *Boulder: Geological Society of America Special Paper*.

615
616
617 Conrad, C. P. & Molnar, P. (1997). The growth of Rayleigh–Taylor-type instabilities in the lithosphere for various rheological and density structures. *Geophysical Journal International* 129, 95-112.

618 Cox, K. G. (1980). A model for continental flood vulcanism. *J. Petrol.* 21, 629-650.

619
620 Crisp, J. A. (1984). Rates of magma emplacement and volcanic output. *J. Volcanol. Geotherm. Res.* 20, 177-211.

621
622 Dhuime, B., Hawkesworth, C. J., Cawood, P. A. & Storey, C. D. (2012). A change in the geodynamics of continental growth 3 billion years ago. *Science* 335, 1334-1336.

623
624 Dhuime, B., Wuestefeld, A. & Hawkesworth, C. J. (2015). Emergence of modern continental crust about 3 billion years ago. *Nature Geoscience* 8, 552-555.

625 Faccenda, M. (2014). Water in the slab: A trilogy. *Tectonophysics* 614, 1-30.

626
627 Fitton, J. G., Mahoney, J. J., Wallace, P. J. & Saunders, A. D. (2004). Origin and Evolution of the Ontong Java Plateau. *Geological Society, Special Publication* 229, 384pp

628
629
630 Fontaine, F. J., Rabinowicz, M., & Cannat, M. (2017). Can high-temperature, high-heat flux hydrothermal vent fields be explained by thermal convection in the lower crust along fast-spreading Mid-Ocean Ridges? *Geochemistry, Geophysics, Geosystems*, 18(5), 1907-1925.

631
632 Furlong, K. P., Hanson, R. B., & Bowers, J. R. (1991). Modeling thermal regimes. *Reviews in Mineralogy and Geochemistry*, 26(1), 437-505

633
634 Ganne, J. & Feng, X. (2017). Primary magmas and mantle temperatures through time. *Geochem. Geophys. Geosyst.* 18, DOI: 10.1002/2016GC006787.

635
636 Gerya, T.V. (2014) Precambrian geodynamics: Concepts and models. *Gondwana Research* 25, 442-463
637 Hacker, B.R., Abers, G.A., Peacock, S.M., (2003). Subduction factory 1. Theoretical
638 mineralogy, densities, seismic wave speeds, and H₂O contents. J. Geophys. Res. 108.

639 Herzberg, C., & Asimow, P. D. (2008). Petrology of some oceanic island basalts: PRIMELT2.
640 XLS software for primary magma calculation. Geochemistry, Geophysics,
641 Geosystems, 9(9).

642 Herzberg, C., Condie, K. & Korenaga, J. (2010). Thermal history of the Earth and its
643 petrological expression. Earth and Planetary Science Letters 292, 79-88.

644 Houseman, G. A. & Molnar, P. (1997). Gravitational (Rayleigh-Taylor) instability of a layer
645 with non-linear viscosity and convective thinning of continental lithosphere. Geophysical
646 Journal International 128, 125-150.

647 Ito, K. & C. Kennedy, G. (1971). An Experimental Study of the Basalt-Garnet Granulite-
648 Eclogite Transition.

649 Jahn, B. M., Glikson, A. Y., Peucat, J.-J. & Hickman, A. H. (1981). REE geochemistry and
650 isotopic data of Archaean silicic volcanics and granitoids from the Pilbara Block, western
651 Australia: implications for the early crustal evolution. Geochemica et Cosmochimica
652 Acta 45, 1633-1652.

653 Jain, C., Rozel, A. B., Tackley, P. J., Sanan, P., & Gerya, T. V. (2019). Growing primordial
654 continental crust self-consistently in global mantle convection models. Gondwana
655 Research. 73, 96-122

656 Jarrard, R.D. (2003). Subduction fluxes of water, carbon dioxide, chlorine, and potassium.
657 Geochem. Geophys. Geosyst. 4.

658 Jaupart, C., Mareschal, J. C., & Roman, A. (2018). The Formation of Continental Crust from
659 a Physics Perspective. Geochemistry International, 56(13), 1289-1321.

660 Johnson, T. E., Brown, M., Kaus, B. J. P. & VanTongeren, J. A. (2013). Delamination and
661 recycling of Archaean crust caused by gravitational instabilities. Nature Geoscience 7,
662 47-52.

663 Jull, M. & Kelemen, P. A. (2001). On the conditions for lower crustal convective instability.
664 Journal of Geophysical Research 106, 6423-6446.

665 Keller, C. B. & Schoene, B. (2012). Statistical geochemistry reveals disruption in secular
666 lithospheric evolution about 2.5 Gyr ago. Nature 485, 490-493.

667 Korenaga, J. (2006). Archean geodynamics and the thermal evolution of earth. In: Benn, K.,
668 Mareschal, J.-C. & Condie, K. C. (eds.) Archean geodynamics and environments.
669 Washington: Americal Geophysical Union Monograph, 7-32.

670 Korenaga, J. (2013). Initiation and evolution of plate tectonics on Earth: Theories and
671 observations. Annual Review of Earth and Planetary Sciences 41, 117-151.

672 McKenzie, D. & Bickle, M. J. (1988). The volume and composition of melt generated by
673 extension of the lithosphere. J. Petrol. 29, 625-679.

674 Michaut, C., & Jaupart, C. (2006). Ultra-rapid formation of large volumes of evolved magma.
675 Earth and Planetary Science Letters, 250(1-2), 38-52.

676 Moore, W. B. & Webb, A. A. G. (2013). Heat-pipe Earth. Nature 501, 501-505.

677 Morgan, J. P. & Chen, Y. J. (1993). The genesis of oceanic crust: Magma injection,
678 hydrothermal circulation, and crustal flow. Journal of Geophysical Research: Solid Earth
679 98, 6283-6297.

680 Moyen, J.-F. & Martin, H. (2012). Forty years of TTG research. Lithos 148, 312-336.

681 Moyen, J.-F. & Stevens, G. (2006). Experimental studies on TTG petrogenesis: implications
682 for Archean geodynamics. In: Benn, K., Mareschal, J.-C. & Condie, K. C. (eds.) Archean
683 geodynamics and environments. Washington: Americal Geophysical Union Monograph,
149-176.

685 Moyen, J.-F. & Van Hunen, J. (2012) Short-term episodicity of Archaean plate tectonics.
686 Geology 40, 451-454.

687 Müntener, O., Kelemen, P. B., & Grove, T. L. (2001). The role of H 2 O during crystallization
688 of primitive arc magmas under uppermost mantle conditions and genesis of igneous
689 pyroxenites: an experimental study. *Contributions to Mineralogy and Petrology*, *141*(6),
690 643-658.

691 Neal, C.R., Sager, W.W., Sano, T., Erba, E. (2015). The origin, evolution, and environmental
692 impact of oceanic large igneous provinces. *Geological Society of America, Special Paper,
693 511*.

694 O'Neill, C, Lenardic, A., Moresi, L., Torsvik, T.H., Lee, C.-T. (2007) Episodic Precambrian
695 subduction. Earth and Planetary Science Letters 262(3-4):552-562

696 Piccolo, A., Palin, R.M. Kaus, B.J. and White, R.W. (2019) Generation of Earth's early
697 continents from a relatively cool Archean mantle. Geochemistry, Geophysics,
698 Geosystems, 20, doi.org/10.1029/2018GC008079

699 Robin, C. M. I. & Bailey, R. C. (2009). Simultaneous generation of Archean crust and
700 subcratonic roots by vertical tectonics. Geology 37, 523-526.

701 Rozel, A. B., Golabek, G. J., Jain, C., Tackley, P. J., & Gerya, T. (2017). Continental crust
702 formation on early Earth controlled by intrusive magmatism. Nature, 545(7654), 332.

703 Roman, A., & Jaupart, C. (2016). The fate of mafic and ultramafic intrusions in the continental
704 crust. Earth and Planetary Science Letters, 453, 131-140.

705 Shirey, S. B. & Richardson, S. H. (2011). Start of the Wilson cycle at 3 Ga shown by diamonds
706 from subcontinental mantle. Science 333, 434-436.

707 Shirley, D. N. (1987). Differentiation and Compaction of the Palisades Sill, New Jersey.
708 Journal of Petrology 28, 835-865.

709 Sizova, E., Gerya, T., Brown, M. & Perchuk, L. (2010). Subduction styles in the Precambrian:
710 Insight from numerical experiments. Lithos 116, 209-229.

711 Sizova, E., Gerya, T., Stüwe, K. & Brown, M. (2015). Generation of felsic crust in the
712 Archean: A geodynamic modeling perspective. Precambrian Research 271, 198-224.

713 Sleep, N. H. & Windley, B. F. (1982). Archaean plate tectonics: constraints and inferences. J.
714 Geol. 90, 363-379.

715 Gerya, T.V., Stern, B., Baes, M., Sobolev, S.V., Wjatten, S.A. (2015) Plate tectonics on the
716 Earth triggered by plume-induced subduction initiation. Nature 527(7577):221-225

717 Storey, M., Mahoney, J. J., Kroenke, L. W. & Saunders, A. D. (1991). Are oceanic plateaus
718 sites of komatiite formation? Geology 19, 376-379.

719 Tang, M., Chen, K. & Rudnick, R. L. (2016). Archean upper crust transition from mafic to
720 felsic marks the onset of plate tectonics. Science 351, 372-375.

721 Thébaud, N. & Rey, P. F. (2013). Archean gravity-driven tectonics on hot and flooded
722 continents: Controls on long-lived mineralised hydrothermal systems away from
723 continental margins. Precambrian Research 229, 93-104.

724 Theissen-Krah, S., Rüpke, L. H. & Hasenclever, J. (2016). Modes of crustal accretion and
725 their implications for hydrothermal circulation. Geophysical Research Letters 43, 3.

726 van Hunen, J. & van den Berg, A. P. (2008). Plate tectonics on the early Earth limitations
727 imposed by strength and buoyancy of subducted lithosphere. Lithos 103, 217-235.

728 Van Kranendonk, M. J., Collins, W. J., Hickman, A. H. & Pawley, M. J. (2004). Critical tests
729 of vertical vs. horizontal tectonic models for the Archaean East Pilbara Granite-
730 Greenstone Terrane, Pilbara Craton, Western Australia. Precambrian Research 131, 173-211.

731 Weller, O.M. Copley, A., Miller, W.R.G., Palin, R.M., Dyck, B. (2019) The relationship between mantle potential temperature and oceanic lithosphere buoyancy, Earth Planetary Science Letters, 518, 86-99

735 White, S. M., Crisp, J. A. & Spera, F. R. (2006). Long-Term Volumetric Eruption Rates and Magma Budgets. Geochem. Geophys. Geosyst. 7, DOI: 10.1029/2005GC001002.

737 Wilson, A. (2012). A chill sequence to the Bushveld Complex: Insight into the first stage of emplacement and implications for the parental magmas. Journal of Petrology, 53, 1123-1168.

740 Zhang, C., Holtz, F., Koepke, J., Wolff, P. E., Ma, C. & Bédard, J. H. (2013). Constraints from experimental melting of amphibolite on the depth of formation of garnet-rich restites, and implications for models of Early Archean crustal growth. Precambrian Research 231, 206-217.

745 Appendix
746 A. Dimensionless parameters

747 If we assume a linear thermal gradient near the Moho, the variation of viscosity with depth can be written

749 $\eta(z) = \eta_0 \exp(-1/L z)$ (A1)

751 where $\eta_0$ is the viscosity at the Moho defined as

753 $\eta_0=1/2 A^{(-1/n)}(\epsilon_{xx})^{(1-n)/n} \exp(Q\beta/ nRT_0)$

755 L is the length-scale over which the viscosity decreases by a factor e, which is given by:

757 $L = nRT_0^2/Q\beta$ (A2)

759 n is the stress exponent, R is the gas constant, T₀ is the temperature at the Moho, Q is the activation energy and β is the thermal gradient dT/dz, A is a constant and $\epsilon_{xx}$ is the background strain-rate. Taking n=3.5, R=8.31 J mol⁻¹ K⁻¹, Q =500 kJ mol⁻¹, varying Moho temperature between 800 and 1200°C and temperature gradients around the Moho between 10-35 °C km⁻¹, we find that, in most of the cases, the characteristic length-scale is less than 10 km.

765 We define two dimensionless parameters describing the dynamics of the problem: the rheological parameter: $\gamma = L / h_{lc}$

768 and the buoyancy number: $B = (\rho_{lc} - \rho_{uc})/(\rho_{lc} - \rho_{m})$ where the subscript lc, uc, m, denotes the lower crust, upper crust and mantle respectively.

770 B. Numerical Code and Setup

772 The numerical calculations presented in this paper were performed using the thermo-mechanical Flamar-Paravoz code, which has been extensively described elsewhere and used in the past for a wide range of geodynamical problems (Burov and Gerya 2014, Roman and Jaupart, 2016). The code solves the Newtonian equations of motion coupled with the energy
776 conservation. These equations are solved using an explicit, time-marching scheme in a
777 Lagrangian grid, which allows us to precisely follow different interfaces, and implement a
778 remeshing procedure when the grid becomes too distorted.

779
780 In our calculations we set a zero displacement condition on both the upper and lower sides of
781 the computational domain. A comparison with some test runs revealed only minor differences
782 by imposing free surface condition on the top. The left and right boundaries are reflecting with
783 zero normal displacement and zero shear stress. The grid has a uniform and equal resolution
784 in the vertical and lateral direction, with the vertical dimension of the domain equal to ten
785 times the thickness of the cumulate layer. The lateral extent was varied between 0.5-10 times
786 the thickness of the cumulate thickness in order to study modes with wavelengths between 1-
787 20 times the cumulate thickness. An initial sinusoidal perturbation with an amplitude equal to
788 1% of the cumulate layer thickness was imposed at the crust-mantle interface in order to start
789 the instability.

790
791 Although the code allows us to simulate the flow in visco-elasto-plastic materials, we have
792 used only viscoelastic, Maxwell type material in order to isolate the interplay between
793 buoyancy and viscous forces. While the mantle may undergo to both diffusion creep (n = 1)
794 and dislocation creep (n>1) depending on the stress magnitudes, for the sake of conciseness
795 we present here only results for n = 3.5, for which the late phase non-linear effects are expected
796 to be the strongest. Growth rates are exponential in the case of n=1 and super-exponential for
797 n>1, but the effects in the entrainment of the upper layer are only minor, as will be shown in
798 a companion paper.

799
800
801 **Figure Captions**

802
803 Figure 1: a) Schematic section of an internally differentiated Archean crust b) Amount of
804 fractionation as function of the MgO content of primary liquids. Data are Herzberg et al.'s
805 (2010) Precambrian non-arc basalts. The amount of fractionation required to match primary
806 liquid composition is calculated with the method of Herzberg and Asimow, (2008). c) Purple
807 dots: density of the olivines fractionated from the primary liquids. Green dots: density of the
808 residual mantle left behind primary melt extraction. Densities are calculated from the
809 compositions obtained with the method of Herzberg and Asimow (2008) and are expressed at
810 ambient temperature and atmospheric pressure. Dots correspond to the average density of the
811 compositional range of crystallized olivines, while error-bars correspond to their standard
812 deviation.

813
814 Figure 2: Final temperature as function of depth of a crust emplaced in 1 Myr (blue line) and
815 10 Myr (green line). In this calculation, the amount of intrusive rock is set to 0. The initial
816 lithospheric thickness is 40 km at the temperature at the bottom 1400 °C. The final structure
817 is shown on the right.

818
819 Figure 3: Temperature evolution in a reference calculation in which intrusion is taken into
820 account. Temperatures are plotted as function of depth and time for 1 Myr. The intrusive to
821 extrusive ratio is set to 1, so that the final crust is composed of equal proportions of lavas and
822 intrusions (the final structure at 1 Myr is shown in the right column). The dehydration
823 temperature window is marked by red solid lines.

824
825 Figure 4: a) Final temperature profiles of crust formed in 1 Myr (blue), 5 Myr (red) and 10
826 Myr (green). Intrusions are emplaced at the Moho, and the intrusive to extrusive ratio is equal
827 to 1. Initial lithospheric thickness is set to 40 km. b) Temperatures of the Moho as function of
828 depth for the calculations shown on the left. In each calculation, the Moho depth increases
829 with time as it is buried by intrusions and flows.

830
831 Figure 5: Comparison of the final thermal profiles for intrusive to extrusive ratios equal to 1
832 (blue line) and two (green line) The emplacement time is 5 Myr.

833
834 Figure 6: a) Snapshots depicting the evolution of an experiment with B = 3 and the same
835 viscosity in the lower/upper crust and mantle. The upper crust is shown as alternating dark
836 grey and light grey layers (layering is only used to track the deformation and does not
837 correspond to any change in the crustal rock properties). b) Dimensionless amplitude of the
838 perturbation of the lower crust (solid violet line), of the thickened part of the upper crust (solid
839 black line) and of the thinned part of the upper crust (dash dot line) as function of
840 dimensionless time. The locations at which measurements are made for each line are shown
841 by black arrows. All the quantities have been scaled by the thickness of the cumulate layer h<sub>lc</sub>.

842
843 Figure 7: Maximum downward displacement of the upper crust as function of the buoyancy number
844 B. The maximum downward deflection corresponds to the positive peak observed in Figure 6b.

845
846 Figure 8: Snapshot of an experiment with an exponential decrease in viscosity. L/h<sub>lc</sub> = 0.65 and
847 B=1.0.
<figure><img src="image_2.png" /></figure>

b)

<figure><img src="image_3.png" /></figure>

c)

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

Extrusive crust

Lithospheric
mantle
Time [Myr]

<figure><img src="image_6.png" /></figure>

Extrusive crust

Intrusive crust

Lithospheric
mantle
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

Int2ext=1

<figure><img src="image_10.png" /></figure>

Int2ext=2

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>