Modélisation numérique de la convection thermochimique
La convection thermochimique compressible est étudiée à l’aide du code StagYY23 en géométrie bidimensionnelle d’anneau sphérique24. Les équations de conservation de la quantité de mouvement, de la masse et de l’énergie sont résolues à l’aide d’un solveur parallèle MUMPS disponible dans le paquet PETSc25. L’approximation des différences finies est utilisée sur une grille décalée26 avec une résolution variant radialement (plus élevée dans les couches limites et à la transition de phase 660). Le domaine est discrétisé par 64 (radial) fois 512 (latéral) cellules. Des conditions limites de glissement libre sont utilisées pour les frontières extérieures et intérieures. La température de surface est fixée à 300 K et la température initiale du noyau est à 6 000 K.
Un champ de traceurs est advecté à travers le maillage en utilisant une interpolation spatiale sans divergence de second ordre du champ de vitesse et un schéma Runge-Kutta de quatrième ordre à travers le temps. Chaque traceur transporte plusieurs quantités, telles que la teneur en eau, la composition, la température et le taux de chauffage radiogénique. L’utilisation de ces champs sur les traceurs limite la diffusion numérique. Les éléments radiogéniques sont partitionnés lors de la fusion : les fontes basaltiques sont 100 fois plus enrichies que leur résidu solide. Pour simplifier, on considère que l’eau pénètre entièrement dans les 10 km supérieurs et qu’elle est advectée dans tout le manteau. Le modèle de composition traite les roches solides et fondues comme étant une combinaison linéaire de basalte (matériau crustal) et de harzburgite (manteau appauvri). Le manteau primordial a une composition pétrologique initiale de 20 % de basalte et de 80 % de harzburgite, soit un mélange de 60 % d’olivine et de 40 % de pyroxène-grenat. Pour des raisons d’efficacité numérique, la composition est stockée sur des traceurs, tandis que les calculs de fusion sont effectués sur la maille. Dans ce but, un champ de composition maillé est calculé au centre des cellules en faisant la moyenne des traceurs dans la cellule. Sur la maille, une fraction fondue est calculée à chaque pas de temps, en comparant les conditions de pression-température à une fonction de solidus dépendant de la composition, en considérant une chaleur latente de 600 kJ kg-1. Une fois que la quantité de matière fondue est calculée sur le maillage, la masse correspondante de matière fondue est créée sur le champ de traçage. Si la fusion se produit à moins de 300 km de profondeur, les traceurs fondus de composition entièrement basaltique sont transportés soit vers le sommet du domaine (éruption), soit vers le bas de la croûte (intrusion). La température des traceurs éruptifs est fixée à la température de surface, formant rapidement une lithosphère froide. La température des traceurs intrudés ne prend en compte que la décompression adiabatique, ce qui donne une lithosphère chaude. La colonne de matière entre la région de la source de fusion et le lieu de l’intrusion ou de l’éruption est instantanément advectée vers le bas pour conserver la masse. Le résidu solide laissé par la procédure d’éruption-intrusion est de la harzburgite pure. Pour plus de détails concernant l’implémentation, voir les réf. 27 et 28. La viscosité dépend à la fois de la profondeur et de la température selon la formulation d’Arrhenius (voir réf. 18 pour plus de détails). La densité est calculée comme la somme des effets de pression, thermiques et compositionnels, y compris les transitions de phase solide-solide. Les phases d’olivine et de pyroxène-grenat sont traitées séparément, de sorte que l’augmentation de la densité associée à la transition de phase basalte-eclogite peut être prise en compte. L’éclogitisation conduit le processus d’égouttage observé pendant les premières centaines de millions d’années après la cristallisation de l’océan magmatique (voir réf. 18 pour plus de détails).
Nous supposons que la déformation visqueuse opère par le processus de fluage par diffusion, qui suit une loi d’Arrhenius dépendant de la température et de la pression :
où η0 est la viscosité de référence (1020 Pa s) à pression nulle et température de référence T0 (1600 K), Δη est un facteur utilisé pour imposer des sauts de viscosité entre les différentes couches i (jusqu’à 5 ou 6), Ei est l’énergie d’activation dans la couche i, P est la pression, Vi est le volume d’activation, R est la constante des gaz (8.314 J mol-1 K-1) et T est la température absolue. Diverses valeurs pour Ei et Vi sont utilisées pour le manteau supérieur et inférieur et pour la couche post-perovskite proche de la limite noyau-manteau29,30,31 (voir le tableau 3 des données étendues). Dans le manteau inférieur, le volume d’activation diminue avec l’augmentation de la pression selon:
Pi est donné dans le tableau de données étendues 3 (ou est infini s’il n’est pas donné). Un saut de viscosité de 30 est imposé à la transition entre le manteau supérieur et le manteau inférieur32.
La déformation plastique dans la lithosphère est calculée en utilisant une contrainte d’écoulement τy suivant l’approche de Byerlee20:
où f est le coefficient de friction. La viscosité effective est calculée comme:
où 
est le second invariant du tenseur des taux de déformation. La viscosité ne dépend pas de la fraction fondue ou de la composition.
Puisque nous résolvons la convection compressible, la température adiabatique, la densité, la conductivité thermique, l’expansivité thermique et la capacité thermique dépendent de la pression suivant une équation d’état de Birch-Murnaghan du troisième ordre, qui relie le module apparent à la pression (voir réf. 31 pour une explication détaillée). Plus important encore pour la présente étude, les matériaux du manteau sont traités comme un mélange de systèmes de phases olivine, pyroxène-grenat et fonte, suivant une paramétrisation basée sur des données de physique minérale33,34. Chaque système de phase a ses propres profondeurs de transition de phase et ses propriétés physiques. La densité de la masse fondue augmente régulièrement avec la pression. Les densités de l’olivine solide et du pyroxène-grenat augmentent également de façon régulière avec la profondeur mais contiennent également des sauts aux transitions de phase, comme le montre le tableau 4 des données étendues. En particulier, la transition de phase basalte-éclogite située à une profondeur d’environ 60 km a un effet de premier ordre sur la dynamique de la lithosphère18. Au-dessus de la transition de phase éclogite, le basalte est plus léger que l’olivine de 160 kg m-3. En dessous de la transition de phase, l’éclogite formée à partir du basalte devient plus dense que l’olivine de 190 kg m-3, ce qui entraîne facilement des instabilités mécaniques dans la lithosphère. Dans le manteau plus profond, le pyroxène-grenat et l’olivine présentent également des augmentations de densité avec la profondeur aux limites de la zone de transition et au sommet de la couche post-pérovskite (voir le tableau de données étendues 4). La pression PPT à laquelle la transition de phase se produit dépend de la température : PPT = P0 + γT, où P0 est une pression de référence. Les pentes de Clapeyron γ de chaque transition (voir le tableau des données étendues 4) quantifient cet effet.
Estimation des volumes de TTG produits
Les volumes de roches TTG présentés sur la figure 2c sont obtenus via un processus en deux étapes. Nous calculons d’abord le volume de basalte hydraté qui correspond aux conditions de pression-température pour la génération de TTG (voir le tableau 1 des données étendues et la réf. 11). Pour obtenir ce que serait le volume des roches TTG dans une géométrie sphérique tridimensionnelle, nous calculons l’épaisseur moyenne globale produite et la multiplions par la surface de la Terre pour chaque pas de temps t. Ensuite, nous faisons la somme de toutes les cellules i du domaine de calcul pour un certain temps t et obtenons un volume instantané de TTG :
où S est la surface de la Terre, δzi est l’épaisseur qu’aurait la cellule i si elle était étalée sur toute la surface du domaine de calcul. 
est une concentration d’eau sans dimension (1 pour une roche complètement hydratée) et 
est une ‘concentration de basalte’ sans dimension (1 pour le basalte, 0 pour la harzburgite). Le volume est divisé par 2 en supposant que le degré moyen de fusion partielle de la croûte basaltique hydratée produisant le TTG-melt correspond à 50%35. La sommation des cellules n’est effectuée que pour les cellules ayant à la fois des concentrations en eau et en basalte supérieures à la moitié de la teneur en eau de surface imposée (c’est-à-dire contenant une grande quantité de croûte basaltique hydratée).
Le volume cumulé V(t) représenté sur la Fig. 2c est obtenu par intégration temporelle pondérée du volume instantané v(t):
où VTTG-P est le volume des roches aux conditions de pression pertinentes pour la formation de TTG (une valeur constante) et Vhbs est le volume de basalte hydraté dans l’ensemble du domaine. Nous négligeons les rares situations dans lesquelles ∂Vhbs/∂t est négatif. L’équation (6) garantit que les TTG ne sont produits que lorsque les conditions pertinentes pour leur formation sont réunies (v est non nul) et également lorsque du nouveau basalte hydraté est effectivement formé (∂Vhbs/∂t est positif). Cette stratégie est beaucoup plus robuste que la simple intégration du volume de roche aux conditions de formation du TTG sans aucune pondération. Une situation de convection à l’équilibre dans laquelle il ne se forme plus de basalte mais où la géothermie dans la croûte correspond aux conditions de formation de TTG surestimerait considérablement le volume de TTG réellement produit.
Cependant, l’extraction de données quantitatives de nos simulations numériques doit être considérée avec prudence. Notre algorithme simplifié de calcul du volume de la croûte TTG (Méthodes) utilise un degré moyen de fusion de la croûte basaltique hydratée et ne crée pas directement et ne trace pas de régions de roches TTG au sein de la croûte. Les modèles régionaux avec traçage direct des roches TTG montrent généralement des taux de croissance crustale plus élevés10, mais un tel traçage n’est pas encore réalisable à l’échelle mondiale en raison du manque de résolution numérique. En outre, notre modèle simplifié peut surestimer la durée des épisodes de couvercle stagnant (près de 50 % du temps de simulation, figure 2), qui produisent des quantités négligeables de roche TTG. Dans les modèles régionaux à plus haute résolution, la durée de tels épisodes est limitée à environ 80 millions d’années9,10 en raison de l’affaiblissement localisé de la lithosphère par des processus magmatiques, ce qui n’est pas mis en œuvre dans notre modèle global simplifié. Nous discutons ci-dessous des implications potentielles de la présente étude pour la dynamique de Vénus.
Disponibilité du code
Comparaison de nos simulations numériques avec les observations de terrain
Comparer nos résultats numériques avec la variété et la complexité des données de terrain des cratons archéens est ambitieux. Si nous pouvons évaluer si certaines conditions rhéologiques et certains comportements magmatiques peuvent conduire à la formation de roches felsiques sur la Terre primitive, nous ne pouvons certainement pas expliquer la complexité des données pétrologiques et géochimiques rapportées au cours des dernières décennies. Nous nous concentrerons donc sur l’image de premier ordre convenue dans la communauté de la géologie de terrain de l’Archéen. Nous sommes convaincus que nos modèles numériques donnent une réponse robuste concernant l’efficacité élevée nécessaire de l’intrusion pour produire une distribution raisonnable des roches TTG. Pourtant, nous devons également garder à l’esprit que nos modèles manquent encore d’importants processus, ce qui rendrait une comparaison détaillée avec les données de terrain hasardeuse et probablement trompeuse, bien que le message simple de cette étude soit toujours valable.
Premièrement, notre rhéologie ne dépend pas de la composition, ce qui nous empêche de produire un manteau lithosphérique sub-continental très rigide. En revanche, notre viscosité est fortement dépendante de la température (comme le suggèrent les données expérimentales), ce qui génère en soi une lithosphère sub-crustale très visqueuse. Cette fonctionnalité non triviale n’est possible que parce que nous utilisons des solveurs directs parallèles très robustes18,25. Deuxièmement, étant donné que nous considérons un milliard d’années d’évolution globale du manteau, la résolution de nos modèles ne peut pas être suffisamment élevée pour étudier les détails des processus crustaux, qui sont très précieux pour les géologues de terrain qui utilisent l’approche descendante (reliant les données de « surface » pour suggérer une image géodynamique globale). Cette limitation importante ne nous permet pas de modéliser les processus crustaux régionaux (tels que les structures de dôme et de quille) qui modifieraient le profil de température dans la croûte. De telles simulations ont été réalisées à l’aide de simulations tridimensionnelles à haute résolution9, et suggèrent une image similaire à la nôtre, mais à une échelle régionale. De plus, notre approche est certainement valable puisque le renversement crustal ne se produira que si la croûte est réellement produite en premier lieu, ce que nous avons la capacité de modéliser. Troisièmement, notre modèle pétrologique est relativement simple car nous voulons déterminer quel environnement géodynamique global conduira à des conditions favorables pour la formation de TTG. Enfin, nos simulations ont été réalisées en deux dimensions, car les simulations tridimensionnelles sont très coûteuses en termes de calcul et ne nous auraient pas permis d’étudier systématiquement notre espace de paramètres. En gardant à l’esprit ces limites inhérentes à la modélisation, nous comparons les caractéristiques observées dans nos modèles numériques à des scénarios géodynamiques basés sur des données de terrain.
Bien que très peu de roches archéennes aient été préservées sur Terre, il semble largement admis que le mode tectonique a grandement changé peu après le début du Paléo-archéen, à environ 3,55 milliards d’années36. En effet, tout le matériel plus ancien que cet âge (à l’Eoarchéen) se trouve dans des terranes de gneiss à haute teneur empilés, alors que les roches felsiques produites plus tard (pendant la majeure partie du Paléo-archéen) semblent provenir d’un environnement tectonique moins catastrophique, suffisamment stable pour générer du diapirisme crustal16. Il n’y a actuellement aucun consensus sur le processus géodynamique/tectonique formant les différentes structures rapportées par les études de terrain7,16. Les structures éoarchéennes montrent une déformation intense et sont parfois interprétées comme des terranes de subduction-accrétion37 ou le résultat d’orogènes chauds38. Les seuls cratons connus de cette période ont été trouvés dans le complexe gneissique d’Itsaq à l’ouest du Groenland et dans la Province du Supérieur au Canada7. On trouve du matériel felsique paléo-archéen dans le craton de Kaapvaal (principalement en Afrique du Sud) et dans le craton de Pilbara (Australie occidentale). Dans ces deux endroits, la croûte a apparemment subi un renversement interne, formant des dômes TTG entourés de » quilles » komatiitiques et basaltiques remplies de sédiments. On pense que l’apparition de ces géométries de dômes et de quilles n’est possible que grâce à la formation d’un manteau lithosphérique sous-continental très visqueux, suffisamment fort pour maintenir la croûte felsique en place malgré son diapirisme39. De plus, il semble que le plutonisme felsique générant les structures de dôme et de quille soit lié à plusieurs pulsations de magmatisme que l’on pense être dues à un panache stationnaire sous le craton40. Ainsi, les géologues de terrain signalent deux types principaux de matériel archéen : Les terranes empilés de l’Éoarchéen montrant une déformation très intense et les structures en dôme et en quille du Paléo-archéen nécessitant apparemment une convection plus modérée.
La phase apparemment calme du Paléo-archéen semble se terminer au début du Mésoarchéen (il y a 3,2 milliards d’années), moment où l’on observe davantage de déformation. Cela a été interprété comme le début de la tectonique des plaques7, mais cette question dépasse le cadre de la présente Lettre et nous nous contenterons de noter ici que l’on rapporte qu’une déformation importante commence à 3,2 milliards d’années.
Nos simulations numériques sont en bon accord avec le résultat principal de ces observations de terrain. Données étendues La figure 1 représente une comparaison de nos simulations (panneau inférieur) avec les données de terrain précédemment mentionnées (panneau supérieur) à travers le temps. En accord avec les terranes empilés milonitisés de l’Eoarchéen dans l’ouest du Groenland et la Province du Supérieur, nous observons toujours une phase initiale de déformation intense dans nos simulations. Cette période vigoureuse est visible sur la Fig. 2a (à gauche, la phase d’égouttage) et dure de cent millions d’années à sept cents millions d’années selon la rhéologie de la lithosphère (voir les régions rouge et brunâtre en bas de la Fig. 1 des données étendues). La phase d’égouttement se produit à la suite des hypothèses suivantes du modèle : la température initiale très élevée du noyau41 produit des panaches intenses et le manteau supérieur est entièrement enrichi. Une grande quantité de croûte mafique est produite pendant ces centaines de millions d’années. Par la suite, nous observons à la fois l’égouttement de morceaux d’éclogite de plusieurs centaines de kilomètres et une tectonique de surface de type proto-subduction, comme on peut le voir sur la figure 2. Ce comportement suggère que la géodynamique de l’Éoarchéen pourrait avoir été une combinaison de processus tectonomagmatiques verticaux et d’événements de courte durée de type tectonique des plaques, qu’il serait difficile d’interpréter en appliquant un seul paradigme. La plupart des roches TTG de basse et moyenne pression sont formées pendant cette période.
La phase d’égouttement est toujours suivie d’une période de » couvercle stagnant » qui offrirait la stabilité nécessaire à la formation des structures de dôme et de quille (régions jaunes dans le panneau inférieur de la Fig. 1 des données étendues). De plus, nous observons une stabilité spatiale surprenante des panaches, due au tassement et au réchauffement lent des énormes blocs éclogitiques à la limite noyau-manteau. Principalement des TTG de haute pression sont formés pendant cette phase, bien qu’une petite quantité de TTG de moyenne pression soit également produite. On peut s’attendre ici à ce que les événements régionaux de renversement de la croûte provoquent la formation d’un plus grand nombre de roches TTG de moyenne et basse pression.
Il est intéressant de noter que des événements de resurfaçage très intenses (régions violettes dans les Données étendues Fig. 1) suivent parfois la phase quiescente, recyclant la majeure partie de la croûte mafique encore présente. La figure 1 des données étendues montre que les événements majeurs de resurfaçage se produisent après 800 millions d’années d’évolution de la Terre (c’est-à-dire, il y a exactement 3,2 milliards d’années) pour une efficacité d’intrusion de 90%. Pourtant, cet épisode de resurfaçage secondaire dû à la réactivation de la convection globale (une fois que les blocs d’éclogite à la limite noyau-manteau sont suffisamment chauds) peut se produire à différents moments de l’évolution de la Terre. Dans l’ensemble, les simulations avec une efficacité d’intrusion de 90% génèrent des transitions géodynamiques cohérentes avec le registre géologique général et produisent les quantités attendues de roches TTG.
En résumé, nous montrons que nos simulations confirment non seulement l’importance du magmatisme intrusif pour la production de la croûte proto-continentale, mais qu’elles sont également en accord avec les observations géologiques pour les cratons archéens. Cependant, d’autres améliorations du modèle sont nécessaires (et en cours) pour une comparaison plus détaillée de nos modèles avec les données de terrain. Nous nous attendons à ce que peu de TTG de basse et moyenne pression produits pendant la phase d’égouttement soient préservés (ce que nous ne pouvons pas encore voir dans nos modèles) et que davantage de ces TTG puissent être produits pendant le processus de dôme et de quille dans la phase de quiescence. Nous avons montré que la phase initiale d’égouttement génère un mélange de tectonique verticale et horizontale de courte durée pendant quelques centaines de millions d’années. Les blocs de matériel mafique produits pendant cet épisode initial et descendant jusqu’à la limite noyau-manteau pourraient être responsables d’une phase stagnante de quelques centaines de millions d’années. Enfin, une phase de convection plus globale peut être attendue lorsque l’éclogite à la limite noyau-manteau est suffisamment chaude. Dans l’ensemble, ces observations qualitatives sont en bon accord avec les mesures quantitatives des volumes de roches TTG produites lorsque (et seulement lorsque) une efficacité d’intrusion élevée est considérée.
Implications de la présente étude pour la dynamique de Vénus
Il a été suggéré récemment que Vénus – bien que dépourvue de croûte continentale – pourrait être un meilleur analogue pour la Terre archéenne que la lune de Jupiter Io7,42,43. Dans ce but, nous compilons brièvement les connaissances actuelles sur la tectonique de Vénus et les comparons aux résultats des modèles numériques, comme suit. Le petit nombre de cratères sur Vénus a conduit à la conclusion que la planète a connu un événement de resurfaçage global au cours du dernier milliard d’années21,42,44,45,46. Des événements de resurfaçage global semblables à ceux de Vénus peuvent être observés dans des modèles globaux prenant en compte à la fois des modèles de formation de croûte extrusive22 et intrusive47 (voir aussi Fig. 2), de sorte que des contraintes géologiques supplémentaires sont nécessaires pour distinguer les scénarios de formation de croûte. Les mesures spectrométriques effectuées sur les sites d’atterrissage Venera et Vega ont montré que la surface de Vénus est dominée par du matériel basaltique48, de sorte qu’une certaine fraction de la fonte générée doit être extrudée sur Vénus. Cependant, il existe également des signes clairs de déformation de la lithosphère, comme les novae49,50, les coronae51,52,53,54 et les terrains en tesselles55, alors que l’existence de zones de subduction reste incertaine56. On a suggéré que les intrusions de type dike sont importantes pour la formation des novae57 et l’analyse structurelle indique que les novae peuvent évoluer vers des structures de type corona au fil du temps50. Les modèles de formation de coronées impliquant la flexion d’une lithosphère épaisse et solide58,59 ont reproduit une partie des neuf groupes topographiques observés de coronées ; cependant, ces modèles ont tendance à surestimer la taille des coronées. Le modèle d’une lithosphère épaisse et froide a été remis en question par des études récentes suggérant une faible épaisseur élastique de moins de 20 km pour de grandes parties de Vénus (ce qui implique une lithosphère de Vénus mince et chaude60,61). Une étude récente prenant en compte ces nouvelles données en modélisant l’intrusion de matière fondue dans la croûte inférieure de Vénus a montré que les coronae peuvent évoluer de manière auto-consistante à partir des structures de novae62, en accord avec les données géologiques. En conclusion, un modèle de formation de croûte incorporant à la fois la formation de croûte extrusive et intrusive semble être nécessaire pour expliquer la tectonique observée sur la Vénus actuelle.