Modellazione numerica della convezione termochimica
La convezione termochimica comprimibile è studiata usando il codice StagYY23 nella geometria dell’anello sferico bidimensionale24. Le equazioni di conservazione di quantità di moto, massa ed energia sono risolte usando un solutore MUMPS parallelo disponibile nel pacchetto PETSc25. L’approssimazione alle differenze finite è usata su una griglia sfalsata26 con una risoluzione variabile radialmente (più alta negli strati limite e alla transizione di fase 660). Il dominio è discretizzato da 64 (radiali) volte 512 (laterali) celle. Le condizioni limite di scivolamento libero sono usate sia per i confini esterni che per quelli interni. La temperatura della superficie è fissata a 300 K e la temperatura iniziale del nucleo è di 6.000 K.
Un campo di traccianti viene fatto avanzare attraverso la mesh usando un’interpolazione spaziale del secondo ordine senza divergenza del campo di velocità e uno schema Runge-Kutta del quarto ordine nel tempo. Ogni tracciante porta diverse quantità, come il contenuto d’acqua, la composizione, la temperatura e il tasso di riscaldamento radiogenico. L’utilizzo di questi campi sui traccianti limita la diffusione numerica. Gli elementi radiogeni sono partizionati durante la fusione: le fusioni basaltiche sono 100 volte più arricchite del loro residuo solido. Per semplicità si considera che l’acqua penetri completamente nei primi 10 km e che venga avvelenata in tutto il mantello. Il modello composizionale tratta le rocce solide e fuse come una combinazione lineare di basalto (materiale crostale) e harzburgite (mantello impoverito). Il mantello primordiale inizia con una composizione petrologica iniziale del 20% di basalto e dell’80% di harzburgite, essendo una miscela del 60% di olivina e del 40% di pirosseno-granato. Per l’efficienza numerica la composizione è memorizzata sui traccianti, mentre i calcoli di fusione sono eseguiti sulla mesh. A questo scopo, un campo di composizione meshato viene calcolato al centro della cella facendo la media dei traccianti all’interno della cella. Sulla mesh, una frazione di fusione viene calcolata ad ogni passo temporale, confrontando le condizioni di pressione-temperatura con una funzione solidale dipendente dalla composizione, considerando un calore latente di 600 kJ kg-1. Una volta calcolata la quantità di fusione sulla mesh, la massa corrispondente di materiale fuso viene creata sul campo del tracciante. Se la fusione avviene a meno di 300 km di profondità, i traccianti fusi di composizione completamente basaltica sono trasportati o verso la cima del dominio (eruzione) o verso il fondo della crosta (intrusione). La temperatura dei traccianti eruttati è impostata sulla temperatura superficiale, formando rapidamente una litosfera fredda. La temperatura dei traccianti intrusi tiene conto solo della decompressione adiabatica, risultando in una litosfera calda. La colonna di materiale tra la regione della sorgente di fusione e il luogo dell’intrusione o dell’eruzione viene istantaneamente spinta verso il basso per conservare la massa. Il residuo solido lasciato dalla procedura di eruzione-intrusione è harzburgite pura. Per maggiori dettagli sull’implementazione vedi i rif. 27 e 28. La viscosità è sia dipendente dalla profondità che dalla temperatura seguendo la formulazione di Arrhenius (vedi rif. 18 per dettagli). La densità è calcolata come somma di pressione, effetti termici e composizionali, incluse le transizioni di fase solido-solido. Le fasi olivina e pirosseno-granato sono trattate separatamente, in modo che l’aumento di densità associato alla transizione di fase basalto-eclogite possa essere preso in considerazione. L’eclogitizzazione guida il processo di gocciolamento osservato durante le prime centinaia di milioni di anni dopo la cristallizzazione dell’oceano magmatico (vedi rif. 18 per i dettagli).
Assumiamo che la deformazione viscosa operi attraverso il processo di creep per diffusione, che segue una legge di Arrhenius dipendente dalla temperatura e dalla pressione:
dove η0 è la viscosità di riferimento (1020 Pa s) a pressione zero e temperatura di riferimento T0 (1600 K), Δη è un fattore usato per imporre salti di viscosità tra i vari strati i (fino a 5 o 6), Ei è l’energia di attivazione nello strato i, P è la pressione, Vi è il volume di attivazione, R è la costante dei gas (8.314 J mol-1 K-1) e T è la temperatura assoluta. Diversi valori per Ei e Vi sono utilizzati per il mantello superiore e inferiore e per lo strato post-perovskite vicino al confine tra nucleo e mantello29,30,31 (vedi Tabella 3 dei dati estesi). Nel mantello inferiore, il volume di attivazione diminuisce all’aumentare della pressione secondo:
Pi è dato nella Extended Data Table 3 (o è infinito se non è dato). Un salto di viscosità di 30 è imposto alla transizione tra il mantello superiore e inferiore32.
Lo snervamento plastico nella litosfera è calcolato usando una tensione di snervamento τy seguendo l’approccio di Byerlee20:
dove f è il coefficiente di attrito. La viscosità effettiva è calcolata come:
dove 
è il secondo invariante del tensore della velocità di deformazione. La viscosità non dipende dalla frazione di fusione o dalla composizione.
Siccome risolviamo per la convezione comprimibile, la temperatura adiabatica, la densità, la conduttività termica, l’espansività termica e la capacità termica sono dipendenti dalla pressione secondo un’equazione di stato di Birch-Murnaghan del terzo ordine, che mette in relazione il modulo di massa con la pressione (vedi rif. 31 per una spiegazione dettagliata). La cosa più importante per il presente studio è che i materiali del mantello sono trattati come una miscela di olivina, pirosseno-granato e sistemi di fase di fusione, secondo una parametrizzazione basata su dati di fisica minerale33,34 . Ogni sistema di fase ha la propria profondità di transizione di fase e le proprie proprietà fisiche. La densità del fuso aumenta dolcemente con la pressione. Anche le densità dell’olivina solida e del pirosseno-granato aumentano dolcemente con la profondità, ma contengono anche salti alle transizioni di fase, come mostrato nella tabella 4 dei dati estesi. In particolare, la transizione di fase basalto-eclogite situata ad una profondità di circa 60 km ha un effetto di primo ordine sulla dinamica della litosfera18. Sopra la transizione di fase dell’eclogite, il basalto è più leggero dell’olivina di 160 kg m-3. Al di sotto della transizione di fase, l’eclogite formata dal basalto diventa più densa dell’olivina di 190 kg m-3 , il che determina facilmente instabilità meccaniche all’interno della litosfera. Nel mantello più profondo, sia il pirosseno-granato che l’olivina mostrano anch’essi aumenti di densità con la profondità ai confini della zona di transizione e in cima allo strato di post-perovskite (vedi Tabella 4 dei dati estesi). La pressione PPT a cui avviene la transizione di fase è dipendente dalla temperatura: PPT = P0 + γT, dove P0 è una pressione di riferimento. Le pendenze di Clapeyron γ di ogni transizione (vedi Tabella 4 dei dati estesi) quantificano questo effetto.
Stima dei volumi di TTG prodotti
I volumi delle rocce TTG presentati in Fig. 2c sono ottenuti attraverso un processo in due fasi. Prima calcoliamo il volume del basalto idratato che corrisponde alle condizioni di pressione-temperatura per la generazione di TTG (vedi Tabella 1 dei dati estesi e rif. 11). Per ottenere quello che sarebbe il volume delle rocce TTG nella geometria sferica tridimensionale, calcoliamo lo spessore globale medio prodotto e lo moltiplichiamo per la superficie della Terra per ogni passo temporale t. Poi sommiamo su tutte le celle i del dominio di calcolo per un certo tempo t e otteniamo un volume istantaneo di TTG:
dove S è la superficie della Terra, δzi è lo spessore che la cella i avrebbe se fosse distribuita sull’intera superficie del dominio di calcolo. 
è una concentrazione di acqua adimensionale (1 per la roccia completamente idratata) e 
è una “concentrazione di basalto” adimensionale (1 per il basalto, 0 per l’harzburgite). Il volume è diviso per 2 assumendo che il grado medio di fusione parziale della crosta basaltica idratata che produce TTG-fusa corrisponda al 50%35. La somma delle celle viene eseguita solo per le celle che hanno sia acqua che concentrazioni di basalto superiori alla metà del contenuto di acqua superficiale imposto (cioè, contenenti una grande quantità di crosta basaltica idratata).
Il volume cumulativo V(t) mostrato in Fig. 2c è ottenuto dall’integrazione temporale ponderata del volume istantaneo v(t):
dove VTTG-P è il volume delle rocce alle condizioni di pressione rilevanti per la formazione di TTG (un valore costante) e Vhbs è il volume di basalto idratato nell’intero dominio. Trascuriamo le rare situazioni in cui ∂Vhbs/∂t è negativo. L’equazione (6) assicura che le TTG siano prodotte solo quando le condizioni rilevanti per la loro formazione sono soddisfatte (v non è zero) e anche quando si forma effettivamente nuovo basalto idratato (∂Vhbs/∂t è positivo). Questa strategia è molto più robusta della semplice integrazione del volume di roccia alle condizioni di formazione del TTG senza alcuna ponderazione. Una situazione di convezione in equilibrio in cui non si forma più basalto ma la geotermia nella crosta corrisponde alle condizioni di formazione delle TTG sovrastimerebbe drammaticamente il volume delle TTG effettivamente prodotte.
Tuttavia, l’estrazione di dati quantitativi dalle nostre simulazioni numeriche deve essere considerata con cautela. Il nostro algoritmo semplificato di calcolo del volume delle TTG-crosta (Metodi) utilizza un grado medio di fusione della crosta basaltica idratata e non crea e traccia direttamente le regioni di TTG-rock all’interno della crosta. I modelli regionali con il tracciamento diretto della roccia TTG dimostrano generalmente tassi di crescita crostale più elevati10 , ma tale tracciamento non è ancora fattibile su scala globale a causa della mancanza di risoluzione numerica. Inoltre, il nostro modello semplificato può sovrastimare la durata degli episodi di ristagno del coperchio (quasi il 50% del tempo di simulazione, Fig. 2), che producono quantità trascurabili di roccia TTG. Nei modelli regionali ad alta risoluzione la durata di tali episodi è limitata a circa 80 milioni di anni9,10 a causa dell’indebolimento localizzato della litosfera da processi magmatici, che non è implementato nel nostro modello globale semplificato. Di seguito, discutiamo le potenziali implicazioni del presente studio per la dinamica di Venere.
Disponibilità del codice
Confronto delle nostre simulazioni numeriche con le osservazioni di campo
Confrontare i nostri risultati numerici con la varietà e la complessità dei dati di campo dei cratoni Archeani è ambizioso. Mentre possiamo valutare se certe condizioni reologiche e comportamenti magmatici possono portare alla formazione di rocce felsiche sulla Terra primitiva, non possiamo certamente spiegare la complessità dei dati petrologici e geochimici riportati negli ultimi decenni. Ci concentreremo quindi sul quadro di primo ordine concordato nella comunità della geologia archeana sul campo. Siamo fiduciosi che i nostri modelli numerici diano una risposta robusta per quanto riguarda l’alta efficienza di intrusione necessaria per produrre una distribuzione ragionevole di rocce TTG. Tuttavia, dovremmo anche tenere a mente che i nostri modelli mancano ancora di processi importanti, che renderebbero un confronto dettagliato con i dati di campo pericoloso e probabilmente fuorviante, anche se il semplice messaggio di questo studio sarebbe ancora valido.
In primo luogo, la nostra reologia non è dipendente dalla composizione, impedendoci di produrre un mantello litosferico sub-continentale molto rigido. Tuttavia, la nostra viscosità è fortemente dipendente dalla temperatura (come suggerito dai dati sperimentali), che di per sé genera una litosfera sub-crustale molto viscosa. Questa funzionalità non banale è possibile solo perché usiamo solutori diretti paralleli molto robusti18,25. In secondo luogo, poiché stiamo guardando un miliardo di anni di evoluzione globale del mantello, la risoluzione dei nostri modelli non può essere abbastanza alta per studiare i dettagli dei processi crostali, che sono molto preziosi per i geologi del campo che usano l’approccio top-down (collegando i dati “di superficie” per suggerire un quadro geodinamico globale). Questa importante limitazione non ci permette di modellare i processi crostali regionali (come le strutture a cupola e a chiglia) che modificherebbero il profilo della temperatura nella crosta. Tali simulazioni sono state effettuate utilizzando simulazioni tridimensionali ad alta risoluzione9, e suggeriscono un quadro simile al nostro, ma su scala regionale. Inoltre, il nostro approccio è certamente valido, poiché il rovesciamento della crosta si verificherà solo se la crosta viene effettivamente prodotta in primo luogo, cosa che abbiamo la capacità di modellare. In terzo luogo, il nostro modello petrologico è relativamente semplice perché vogliamo determinare quale ambiente geodinamico globale porterà a condizioni favorevoli per la formazione delle TTG. Infine, le nostre simulazioni sono state effettuate in due dimensioni, perché le simulazioni tridimensionali sono computazionalmente molto costose e non ci avrebbero permesso di studiare sistematicamente il nostro spazio dei parametri. Con queste limitazioni inerenti alla modellazione, confrontiamo le caratteristiche osservate nei nostri modelli numerici con gli scenari geodinamici basati sui dati di campo.
Anche se sulla Terra si sono conservate pochissime rocce archeane, sembra ampiamente accettato che la modalità tettonica sia cambiata notevolmente poco dopo l’inizio del Paleoarcheo, a circa 3,55 miliardi di anni fa36 . Infatti, tutto il materiale più antico di questa età (nell’Eoarcheano) si trova in terranei impilati di gneiss di alto grado, mentre le rocce felsiche prodotte più tardi (durante la maggior parte del Paleoarcheano) sembrano essere nate da un ambiente tettonico meno catastrofico, abbastanza stabile da generare il diapirismo crostale16. Al momento non c’è consenso sul processo geodinamico/tettonico che forma le varie strutture riportate dagli studi sul campo7,16. Le strutture eoarcheane mostrano un’intensa deformazione e sono a volte interpretate come terranes di subduzione-accrezione37 o il risultato di orogeni caldi38. Gli unici cratoni conosciuti di questo periodo sono stati trovati nel complesso Itsaq gneiss nella Groenlandia occidentale e nella Provincia Superiore in Canada7. Materiale felsico paleoarcheo può essere trovato nel cratone del Kaapvaal (principalmente in Sudafrica) e nel cratone del Pilbara (Australia occidentale). In entrambe queste località, la crosta ha apparentemente sperimentato un rovesciamento interno, formando cupole TTG circondate da “chiglie” komatiitiche e basaltiche riempite di sedimenti. Si pensa che la comparsa di queste geometrie di cupole e chiglie sia possibile solo a causa della formazione di un mantello litosferico subcontinentale molto viscoso, abbastanza forte da tenere la crosta felsica in posizione nonostante il suo diapirismo39. Inoltre, sembra che il plutonismo felsico che genera le strutture a cupola e a chiglia sia legato a diversi impulsi di magmatismo che si pensa siano dovuti a un pennacchio stazionario sotto il cratone40. Così, i geologi del campo riportano due tipi principali di materiale archeano: Terranes impilati dell’Eoarcheano che mostrano una deformazione molto intensa e strutture a cupola e a chiglia del Paleoearcheano che apparentemente richiedono una convezione più moderata.
La fase apparentemente quiescente del Paleoarcheano sembra finire all’inizio del Mesoarcheano (3,2 miliardi di anni fa), a questo punto si osserva una maggiore deformazione. Questo è stato interpretato come l’inizio della tettonica a placche7, ma questa questione va oltre lo scopo della presente lettera e noteremo solo qui che è stato riportato che la deformazione importante inizia a 3,2 miliardi di anni fa.
Le nostre simulazioni numeriche sono in buon accordo con il risultato principale di queste osservazioni sul campo. Dati estesi La Fig. 1 rappresenta un confronto delle nostre simulazioni (pannello inferiore) con i dati di campo precedentemente menzionati (pannello superiore) nel tempo. In accordo con i terranei impilati milonizzati dell’Eoarcheano nella Groenlandia occidentale e nella Provincia Superiore, osserviamo sempre una fase iniziale di intensa deformazione nelle nostre simulazioni. Questo periodo vigoroso può essere visto in Fig. 2a (a sinistra, la fase di “gocciolamento”) e dura da cento milioni di anni fino a settecento milioni di anni a seconda della reologia della litosfera (vedi le regioni rosse e marroni nella parte inferiore della Extended Data Fig. 1). La fase di dripping si verifica come risultato delle seguenti ipotesi del modello: l’altissima temperatura iniziale del nucleo41 produce intensi pennacchi e il mantello superiore è completamente arricchito. Una grande quantità di crosta mafica viene prodotta durante queste centinaia di milioni di anni. In seguito, si osserva sia lo sgocciolamento di pezzi eclogitici di centinaia di chilometri di dimensione, sia una tettonica superficiale simile alla proto-subduzione, come si può vedere nella Fig. 2. Tale comportamento suggerisce che la geodinamica dell’Eoarcheo potrebbe essere stata una combinazione di processi tettonomagmatici verticali e di eventi di breve durata simili alla tettonica a placche, che sarebbe difficile da interpretare applicando un solo paradigma. La maggior parte delle rocce TTG a bassa e media pressione si formano durante questo periodo di tempo.
La fase di gocciolamento è sempre seguita da un periodo di ‘coperchio stagnante’ che offrirebbe la stabilità necessaria per la formazione delle strutture a cupola e a chiglia (regioni gialle nel pannello inferiore della Extended Data Fig. 1). Inoltre, osserviamo una sorprendente stabilità spaziale dei pennacchi, dovuta all’assestamento e al lento riscaldamento degli enormi blocchi eclogitici al confine nucleo-mantello. Durante questa fase si formano principalmente TTG ad alta pressione, anche se viene prodotta anche una piccola quantità di TTG a media pressione. Ci si può aspettare che gli eventi regionali di ribaltamento crostale causino la formazione di più rocce TTG a media e bassa pressione.
Interessante, eventi di resurfacing molto intensi (regioni viola in Extended Data Fig. 1) a volte seguono la fase di quiescenza, riciclando la maggior parte della crosta mafica ancora presente. Extended Data Fig. 1 mostra che i maggiori eventi di resurfacing si verificano dopo 800 milioni di anni di evoluzione della Terra (cioè esattamente a 3,2 miliardi di anni fa) per un’efficienza di intrusione del 90%. Tuttavia, questo episodio secondario di riemersione dovuto alla riattivazione della convezione globale (una volta che i blocchi di eclogite al confine nucleo-mantello sono abbastanza caldi) può avvenire in vari momenti dell’evoluzione della Terra. Nel complesso, le simulazioni con un’efficienza di intrusione del 90% generano transizioni geodinamiche coerenti con il record geologico generale e producono le quantità previste di rocce TTG.
In sintesi, dimostriamo che le nostre simulazioni non solo confermano l’importanza del magmatismo intrusivo per la produzione della crosta proto-continentale, ma sono anche in accordo con le osservazioni geologiche per i cratoni Archeani. Tuttavia, ulteriori miglioramenti del modello sono necessari (e in corso) per un confronto più dettagliato dei nostri modelli con i dati di campo. Ci aspettiamo che pochi dei TTG a bassa e media pressione prodotti durante la fase di dripping si conservino (cosa che non possiamo ancora vedere nei nostri modelli) e che un numero maggiore di questi TTG possa essere prodotto durante il processo di cupola e chiglia nella fase di quiescenza. Abbiamo dimostrato che la fase iniziale di dripping genera un mix di tettonica verticale e orizzontale di breve durata per alcune centinaia di milioni di anni. I blocchi di materiale mafico prodotti durante questo episodio iniziale e che affondano verso il confine nucleo-mantello potrebbero essere responsabili di una fase stagnante lunga alcune centinaia di milioni di anni. Infine, ci si può aspettare una fase di convezione più globale quando l’eclogite al confine nucleo-mantello è abbastanza calda. Nel complesso, queste osservazioni qualitative sono in buon accordo con le misure quantitative dei volumi di rocce TTG prodotte quando (e solo quando) si considera un’alta efficienza di intrusione.
Implicazioni del presente studio per la dinamica di Venere
E’ stato suggerito recentemente che Venere – sebbene privo di crosta continentale – potrebbe essere un migliore analogo della Terra archeana rispetto alla luna di Giove Io7,42,43. A questo scopo compiliamo brevemente le conoscenze attuali sulla tettonica di Venere e le confrontiamo con i risultati dei modelli numerici, come segue. Il piccolo numero di crateri su Venere ha portato alla conclusione che il pianeta ha sperimentato un evento di resurfacing globale nell’ultimo miliardo di anni21,42,44,45,46. Eventi di resurfacing globale simili a quelli di Venere possono essere osservati in modelli globali che prendono sia modelli di formazione della crosta estrusiva22 che intrusiva47 (vedi anche Fig. 2), quindi sono necessari ulteriori vincoli geologici per distinguere tra gli scenari di formazione della crosta. Misure spettrometriche presso i siti di atterraggio Venera e Vega hanno mostrato che la superficie di Venere è dominata da materiale basaltico48 , quindi una certa frazione del fuso generato deve essere estruso su Venere. Tuttavia, ci sono anche chiari segni di deformazione della litosfera, come novae49,50, coronae51,52,53,54 e terreni a tessera55 , mentre l’esistenza di zone di subduzione rimane incerta56. È stato suggerito che le intrusioni a dike siano importanti per la formazione delle novae57 e l’analisi strutturale indica che le novae possono evolvere in strutture simili alle corone nel tempo50. I modelli di formazione delle corone che coinvolgono la flessione di una litosfera spessa e forte58,59 hanno riprodotto parte dei nove gruppi topografici osservati delle corone; tuttavia, questi modelli tendono a sovrastimare le dimensioni delle corone. Il modello di una litosfera spessa e fredda è stato messo in discussione da studi recenti che suggeriscono un basso spessore elastico di meno di 20 km per ampie parti di Venere (il che implica una litosfera di Venere sottile e calda60,61). Uno studio recente che tiene conto di questi nuovi dati modellando l’intrusione di fusione nella crosta inferiore di Venere ha mostrato che le corone possono evolvere in modo autoconsistente dalle strutture delle novae62 , in accordo con i dati geologici. In conclusione, un modello di formazione della crosta che incorpori la formazione della crosta sia estrusiva che intrusiva sembra essere necessario per spiegare la tettonica osservata su Venere attuale.