Numerical modelling of thermochemical convection
Thermochemical compressible convection is studied using the code StagYY23 in two-dimensional spherical annulus geometry24. Vergelijkingen van behoud van momentum, massa en energie worden opgelost met behulp van een parallelle MUMPS solver beschikbaar in het PETSc pakket25. De eindige-verschilbenadering wordt gebruikt op een verspringend rooster26 met een radiaal variërende resolutie (hoger in de grenslagen en bij de 660 faseovergang). Het domein wordt gediscretiseerd door 64 (radiaal) maal 512 (lateraal) cellen. Voor zowel de buiten- als de binnengrenzen worden vrije sliprandvoorwaarden gebruikt. De oppervlaktetemperatuur is vastgesteld op 300 K en de begintemperatuur van de kern is 6.000 K.
Een tracerveld wordt door de maas geleid met behulp van tweede-orde divergentievrije ruimtelijke interpolatie van het snelheidsveld en een vierde-orde Runge-Kutta schema door de tijd. Elke tracer bevat verschillende grootheden, zoals watergehalte, samenstelling, temperatuur en radiogene verhittingsgraad. Het gebruik van deze velden op tracers beperkt de numerieke diffusie. Radiogene elementen worden gepartitioneerd tijdens het smelten: basaltsmelten zijn 100 keer meer verrijkt dan hun vaste residu. Eenvoudigheidshalve wordt ervan uitgegaan dat water volledig doordringt tot de bovenste 10 km en door de hele mantel wordt gevecteerd. Het samenstellingsmodel beschouwt vaste en gesmolten gesteenten als een lineaire combinatie van basalt (korstmateriaal) en harzburgiet (verarmde mantel). De oermantel begint met een petrologische samenstelling van 20% basalt en 80% harzburgiet, een mengsel van 60% olivijn en 40% pyroxeen-garnet. Voor numerieke efficiëntie wordt de samenstelling opgeslagen op tracers, terwijl de smeltberekeningen worden uitgevoerd op het gaas. Daartoe wordt een gemaasd samenstellingsveld berekend in de middelpunten van de cellen door het gemiddelde te nemen van de tracers binnen de cel. In het netwerk wordt bij elke tijdstap een smeltfractie berekend door de druk-temperatuuromstandigheden te vergelijken met een samenstellingafhankelijke solidusfunctie, waarbij een latente warmte van 600 kJ kg-1 in aanmerking wordt genomen. Zodra de hoeveelheid smelt in het netwerk is berekend, wordt de overeenkomstige massa gesmolten materiaal gecreëerd in het tracerveld. Als het smelten op minder dan 300 km diepte plaatsvindt, worden gesmolten tracers van volledig basaltische samenstelling naar de top van het domein (eruptie) of naar de bodem van de korst (intrusie) getransporteerd. De temperatuur van de geërodeerde tracers wordt op de oppervlaktetemperatuur gebracht, waardoor snel een koude lithosfeer wordt gevormd. Bij de temperatuur van de binnengedrongen tracers wordt alleen rekening gehouden met adiabatische decompressie, wat resulteert in een warme lithosfeer. De materiaalkolom tussen het smeltbrongebied en de intrusie- of eruptielocatie wordt onmiddellijk naar beneden getransporteerd om massa te behouden. Het vaste residu dat door de eruptie-intrusieprocedure wordt achtergelaten, is zuiver harzburgiet. Voor meer details over de uitvoering zie refs 27 en 28. De viscositeit is zowel diepte- als temperatuurafhankelijk volgens de Arrheniusformule (zie ref. 18 voor details). De dichtheid wordt berekend als de som van druk, thermische en samenstellingseffecten, met inbegrip van vaste stof/vaste faseovergangen. Olivijn- en pyroxeen-garnetfasen worden afzonderlijk behandeld, zodat rekening kan worden gehouden met de dichtheidstoename die samenhangt met de basalt-eclogiet faseovergang. De eclogitisatie drijft het druppelproces aan dat wordt waargenomen gedurende de eerste paar honderd miljoen jaar na de kristallisatie van de magma-oceaan (zie ref. 18 voor details).
Wij nemen aan dat de viskeuze vervorming werkt via het diffusiekrimp-proces, dat een temperatuur- en drukafhankelijke Arrhenius-wet volgt:
waarbij η0 de referentieviscositeit is (1020 Pa s) bij nuldruk en referentietemperatuur T0 (1600 K), Δη een factor is die wordt gebruikt om viscositeitssprongen tussen de verschillende lagen i op te leggen (tot 5 of 6), Ei de activeringsenergie in laag i is, P de druk is, Vi het activeringsvolume is, R de gasconstante is (8.314 J mol-1 K-1), Ei de activeringsenergie in laag i is, P de druk is, Vi het activeringsvolume is, R de gasconstante is (8.314 J mol-1 K-1).314 J mol-1 K-1) en T is de absolute temperatuur. Verschillende waarden voor Ei en Vi worden gebruikt voor de boven- en ondermantel en voor de post-perovskietlaag dicht bij de kern-mantelgrens29,30,31 (zie uitgebreide gegevenstabel 3). In de onderste mantel neemt het activeringsvolume af met toenemende druk volgens:
Pi wordt gegeven in de uitgebreide gegevenstabel 3 (of is oneindig indien niet gegeven). Een viscositeitssprong van 30 wordt opgelegd bij de overgang tussen de boven- en ondermantel32.
Plastische rek in de lithosfeer wordt berekend met een vloeispanning τy volgens de Byerlee-benadering20:
waarbij f de wrijvingscoëfficiënt is. De effectieve viscositeit wordt berekend als:
waar de tweede invariant van de rekstrekspannings tensor is. De viscositeit is niet afhankelijk van de smeltfractie of samenstelling.
Omdat we oplossen voor samendrukbare convectie, zijn de adiabatische temperatuur, dichtheid, thermische geleidbaarheid, thermische expansiviteit en warmtecapaciteit druk-afhankelijk volgens een derde-orde Birch-Murnaghan toestandsvergelijking, die de bulkmodulus relateert aan de druk (zie ref. 31 voor een gedetailleerde uitleg). Het belangrijkste voor de huidige studie is dat de mantelmaterialen worden behandeld als een mengsel van olivijn-, pyroxeen-garnet- en smeltfasesystemen, volgens een parametrisatie die is gebaseerd op mineraalfysische gegevens33,34. Elk fasesysteem heeft zijn eigen faseovergangsdiepten en fysische eigenschappen. De smeltdichtheid neemt geleidelijk toe met de druk. De dichtheden van olivijn en pyroxeen-garnet in vaste vorm nemen ook toe met de diepte, maar vertonen ook sprongen bij faseovergangen, zoals te zien is in uitgebreide gegevenstabel 4. Met name de basalt-eclogiet faseovergang op een diepte van ongeveer 60 km heeft een eerste-orde effect op de dynamica van de lithosfeer18. Boven de eclogiet faseovergang is basalt 160 kg m-3 lichter dan olivijn. Onder de faseovergang wordt eclogiet gevormd uit basalt met 190 kg m-3 dichter dan olivijn, wat gemakkelijk mechanische instabiliteiten in de lithosfeer veroorzaakt. In de diepere mantel vertonen zowel pyroxeen-garnet als olivijn ook dichtheidstoenames met de diepte aan de grenzen van de overgangszone en bovenop de post-perovskietlaag (zie uitgebreide gegevenstabel 4). De druk PPT waarbij de faseovergang optreedt is temperatuurafhankelijk: PPT = P0 + γT, waarbij P0 een referentiedruk is. De Clapeyron-hellingen γ van elke overgang (zie uitgebreide gegevenstabel 4) kwantificeren dit effect.
Raming van de geproduceerde volumes TTG
De volumes TTG-rotsen die in fig. 2c worden gepresenteerd, worden verkregen via een tweestapsproces. Eerst berekenen we het volume van gehydrateerd basalt dat overeenkomt met de druk-temperatuurcondities voor TTG-vorming (zie uitgebreide gegevenstabel 1 en ref. 11). Om het volume van TTG-rotsen in driedimensionale bolgeometrie te verkrijgen, berekenen we de gemiddelde mondiale dikte die wordt geproduceerd en vermenigvuldigen dit met het aardoppervlak voor elke tijdstap t. Vervolgens sommeren we alle cellen i van het berekeningsdomein voor een bepaalde tijd t en verkrijgen we een momentaan volume TTG:
waarbij S het oppervlak van de aarde is, δzi de dikte is die de cel i zou hebben als deze over het gehele oppervlak van het berekeningsdomein zou zijn uitgespreid. is een dimensieloze waterconcentratie (1 voor volledig gehydrateerd gesteente) en is een dimensieloze ‘basaltconcentratie’ (1 voor basalt, 0 voor harzburgiet). Het volume wordt door 2 gedeeld in de veronderstelling dat de gemiddelde graad van partieel smelten van gehydrateerde basaltkorst die TTG-smelt produceert, overeenkomt met 50%35. De celsommatie wordt alleen uitgevoerd voor cellen met zowel water- als basaltconcentraties hoger dan de helft van het opgelegde oppervlaktewatergehalte (dat wil zeggen, die een grote hoeveelheid gehydrateerde basaltkorst bevatten).
Het cumulatieve volume V(t) dat in Fig. 2c wordt verkregen door gewogen tijdintegratie van het momentane volume v(t):
waarbij VTTG-P het volume is van gesteenten bij drukcondities die relevant zijn voor TTG-vorming (een constante waarde) en Vhbs het volume is van gehydrateerd basalt in het hele domein. We verwaarlozen zeldzame situaties waarin ∂Vhbs/∂t negatief is. Vergelijking (6) zorgt ervoor dat TTG alleen wordt geproduceerd als aan de voorwaarden voor de vorming ervan is voldaan (v is niet nul) en ook als er daadwerkelijk nieuw gehydrateerd basalt wordt gevormd (∂Vhbs/∂t is positief). Deze strategie is veel robuuster dan het eenvoudigweg integreren van het volume van het gesteente bij condities voor TTG-vorming zonder enige weging. Een evenwichtsconvectiesituatie waarin geen basalt meer wordt gevormd, maar de geotherm in de korst overeenkomt met de omstandigheden voor TTG-vorming, zou een dramatische overschatting geven van het volume TTG’s dat daadwerkelijk wordt geproduceerd.
Het extraheren van kwantitatieve gegevens uit onze numerieke simulaties moet echter met de nodige voorzichtigheid worden overwogen. Ons vereenvoudigde algoritme voor het berekenen van het TTG-volume in de korst (Methodes) maakt gebruik van een gemiddelde graad van gehydrateerd smelten van de basaltkorst en creëert en traceert niet direct TTG-rotsregio’s binnen de korst. Regionale modellen met directe tracering van TTG-gesteente laten over het algemeen hogere groeisnelheden van de korst zien10 , maar een dergelijke tracering is nog niet haalbaar op wereldschaal vanwege het gebrek aan numerieke resolutie. Bovendien kan ons vereenvoudigd model de duur overschatten van episodes van stilstaand deksel (bijna 50% van de simulatietijd, Fig. 2), die verwaarloosbare hoeveelheden TTG-gesteente produceren. In regionale modellen met een hogere resolutie is de duur van dergelijke perioden beperkt tot ongeveer 80 miljoen jaar9,10 als gevolg van de plaatselijke verzwakking van de lithosfeer door magmatische processen, die niet in ons vereenvoudigde mondiale model is geïmplementeerd. Hieronder bespreken we de mogelijke implicaties van deze studie voor de dynamica van Venus.
Beschikbaarheid van code
Vergelijking van onze numerieke simulaties met veldwaarnemingen
Vergelijking van onze numerieke resultaten met de verscheidenheid en complexiteit van veldgegevens van Archeaanse kratons is ambitieus. Hoewel we kunnen evalueren of bepaalde reologische condities en magmatische gedragingen kunnen leiden tot de vorming van felsische gesteenten op de vroege Aarde, kunnen we zeker niet de complexiteit verklaren van de petrologische en geochemische gegevens die in de laatste decennia zijn gerapporteerd. Wij zullen ons daarom richten op het eerste-ordebeeld waarover in de Archeaanse veldgeologie overeenstemming bestaat. Wij zijn ervan overtuigd dat onze numerieke modellen een robuust antwoord geven op de vraag hoe groot de intrusie-efficiëntie moet zijn om een redelijke verdeling van TTG-gesteenten te verkrijgen. Toch moeten we ook bedenken dat in onze modellen nog belangrijke processen ontbreken, waardoor een gedetailleerde vergelijking met veldgegevens gevaarlijk en waarschijnlijk misleidend zou zijn, hoewel de eenvoudige boodschap van deze studie nog steeds geldig zou zijn.
Ten eerste is onze reologie niet compositie-afhankelijk, waardoor we geen zeer stijve subcontinentale lithosferische mantel kunnen produceren. Onze viscositeit is echter sterk temperatuurafhankelijk (zoals door experimentele gegevens wordt gesuggereerd), wat op zich al een zeer visceuze subkristallische lithosfeer oplevert. Deze niet-triviale functionaliteit is alleen mogelijk omdat wij zeer robuuste parallelle directe oplossers gebruiken18,25. Ten tweede kan, aangezien wij de evolutie van de aardmantel over een miljard jaar bestuderen, de resolutie van onze modellen niet hoog genoeg zijn om de details van de processen in de aardkorst te bestuderen, die zeer waardevol zijn voor veldgeologen die de top-down-benadering gebruiken (waarbij gegevens van het aardoppervlak worden gekoppeld om een globaal geodynamisch beeld te schetsen). Deze belangrijke beperking stelt ons niet in staat om regionale processen in de korst (zoals koepel- en kielstructuren) te modelleren die het temperatuurprofiel in de korst zouden wijzigen. Dergelijke simulaties zijn uitgevoerd met behulp van hoge-resolutie driedimensionale simulaties9, en suggereren een beeld dat vergelijkbaar is met het onze, maar dan op regionale schaal. Bovendien is onze benadering zeker valide, aangezien de omwenteling van de korst alleen zal optreden als de korst daadwerkelijk wordt geproduceerd, hetgeen wij kunnen modelleren. Ten derde is ons petrologisch model relatief eenvoudig omdat we willen bepalen welke globale geodynamische omgeving zal leiden tot gunstige omstandigheden voor TTG-vorming. Tenslotte zijn onze simulaties uitgevoerd in twee dimensies, omdat driedimensionale simulaties computationeel erg duur zijn en ons niet in staat zouden hebben gesteld om onze parameterruimte systematisch te bestuderen. Met deze modelinherente beperkingen in het achterhoofd, vergelijken we kenmerken die in onze numerieke modellen zijn waargenomen met geodynamische scenario’s gebaseerd op veldgegevens.
Hoewel er maar zeer weinig Archeens gesteente op aarde bewaard is gebleven, lijkt het algemeen aanvaard dat de tektonische modus sterk veranderde kort na het begin van het Paleoarcheïsch tijdperk, rond 3,55 miljard jaar geleden36. Al het materiaal dat ouder is dan deze leeftijd (in het Eoarchean) kan worden gevonden in gestapelde gneisterranen van hoge kwaliteit, terwijl felsische gesteenten die later (tijdens het grootste deel van het Paleoarchean) zijn geproduceerd, lijken te zijn ontstaan in een minder catastrofale tektonische omgeving, die stabiel genoeg was om korstdiafirisme te genereren16. Er bestaat momenteel geen consensus over het geodynamische/tectonische proces dat de verschillende structuren heeft gevormd die door veldstudies zijn gerapporteerd7,16. Eoarcheaanse structuren vertonen intense vervorming en worden soms geïnterpreteerd als subductie-accretie-terranen37 of het resultaat van hete orogenen38. De enige bekende kratons uit deze periode zijn gevonden in het Itsaq gneiscomplex in West-Groenland en in de Superior Province in Canada7. Paleoarcheens felsisch materiaal is te vinden in het Kaapvaalkraton (voornamelijk in Zuid-Afrika) en in het Pilbara-kraton (West-Australië). Op beide plaatsen heeft de korst blijkbaar een interne kanteling ondergaan, waarbij TTG-koepels zijn gevormd, omgeven door komatiitische en basaltische ‘kielen’, gevuld met sedimenten. Het ontstaan van deze koepel- en kielgeometrieën wordt alleen mogelijk geacht door de vorming van een zeer visceuze subcontinentale lithosferische mantel, die sterk genoeg is om de felsische korst op zijn plaats te houden ondanks zijn diapirisme39. Bovendien lijkt het felsische plutonisme dat de koepel- en kielstructuren genereert, verband te houden met verschillende magmatiseringspulsen die worden toegeschreven aan een stationaire pluim onder het kraton40. Zo melden veldgeologen twee hoofdtypen van Archeens materiaal: Eoarcheïsche gestapelde terranen die zeer intense vervorming vertonen en Paleoarcheïsche koepel- en kielstructuren die blijkbaar meer gematigde convectie vereisen.
De schijnbaar rustige fase van het Paleoarcheïsche tijdperk lijkt te eindigen bij het begin van het Mesoarcheïsche tijdperk (3,2 miljard jaar geleden), op welk moment meer vervorming wordt waargenomen. Dit is geïnterpreteerd als het begin van platentektoniek7 , maar deze vraag valt buiten het bestek van deze brief en we merken hier alleen op dat belangrijke deformatie naar verluidt begint bij 3,2 miljard jaar geleden.
Onze numerieke simulaties komen goed overeen met het belangrijkste resultaat van deze veldwaarnemingen. Uitgebreide gegevens Fig. 1 toont een vergelijking van onze simulaties (onderste paneel) met de eerder genoemde veldgegevens (bovenste paneel) door de tijd heen. In overeenstemming met de gemilonitiseerde gestapelde terranen uit het Eoarchean in West-Groenland en de Superior Province, zien we in onze simulaties steeds een beginfase van intense deformatie. Deze heftige periode is te zien in Fig. 2a (links, de ‘druipende’ fase) en duurt honderd miljoen jaar tot wel zevenhonderd miljoen jaar, afhankelijk van de reologie van de lithosfeer (zie de rood- en bruinachtige gebieden onderin de Extended Data Fig. 1). De druppelfase treedt op als gevolg van de volgende modelaannames: de zeer hoge begintemperatuur van de kern41 produceert intense pluimen en de bovenmantel is volledig verrijkt. Gedurende deze honderden miljoenen jaren wordt een grote hoeveelheid mafische korst geproduceerd. Hierna zien we zowel het uitdruipen van eclogitische brokken van honderden kilometers groot als proto-subductie-achtige oppervlaktetektoniek, zoals te zien is in Fig. 2. Dergelijk gedrag suggereert dat de Eoarcheaanse geodynamica een combinatie zou kunnen zijn geweest van verticale tectonomagmatische processen en kortstondige plaattektoniek-achtige gebeurtenissen die moeilijk te interpreteren zouden zijn als slechts één paradigma zou worden toegepast. De meeste TTG-gesteenten met lage en middelhoge druk worden in deze periode gevormd.
De druipfase wordt altijd gevolgd door een ‘stilstaand deksel’-periode die de stabiliteit zou bieden die nodig is voor de vorming van de koepel- en kielstructuren (gele regio’s in het onderste paneel van Extended Data Fig. 1). Bovendien zien we een verrassende ruimtelijke stabiliteit van de pluimen, als gevolg van de afzetting en langzame opwarming van de enorme eclogitische blokken op de kern-mantel grens. Voornamelijk hoge-druk TTG’s worden gevormd tijdens deze fase, hoewel een kleine hoeveelheid midden-druk TTG’s ook worden geproduceerd. Het is hier te verwachten dat regionale aardkorst kantelingen de vorming van meer midden- en lagedruk TTG gesteenten zouden veroorzaken.
Interessant is dat zeer intense opduikingen (paarse gebieden in Extended Data Fig. 1) soms volgen op de rustige fase, waarbij het grootste deel van de mafische korst nog steeds aanwezig is. Uitgebreide data Fig. 1 laat zien dat grote opduikingen plaatsvinden na 800 miljoen jaar evolutie van de Aarde (dat is precies 3,2 miljard jaar geleden) voor een intrusie efficiëntie van 90%. Deze secundaire opduikingen als gevolg van de heractivering van de convectie (zodra de eclogietblokken op de grens tussen kern en mantel warm genoeg zijn) kunnen echter op verschillende tijdstippen van de evolutie van de Aarde plaatsvinden. Over het geheel genomen genereren de simulaties met een intrusie-efficiëntie van 90% geodynamische overgangen die consistent zijn met de algemene geologische gegevens en produceren zij de verwachte hoeveelheden TTG gesteenten.
Samenvattend laten we zien dat onze simulaties niet alleen het belang van intrusief magmatisme voor de productie van proto-continentale kratons bevestigen, maar ook in overeenstemming zijn met geologische waarnemingen voor Archeaanse kratons. Verdere modelverbeteringen zijn echter nodig (en in voorbereiding) voor een meer gedetailleerde vergelijking van onze modellen met veldgegevens. We verwachten dat weinig van de TTG’s met lage en middelhoge druk die tijdens de druipfase zijn geproduceerd, bewaard zouden blijven (wat we in onze modellen nog niet kunnen zien) en dat meer van deze TTG’s zouden kunnen worden geproduceerd tijdens het koepel- en kielproces in de rustfase. Wij hebben aangetoond dat de initiële druppelfase een mix van verticale en kortstondige horizontale tektoniek genereert gedurende enkele honderden miljoenen jaren. De blokken mafisch materiaal die tijdens deze eerste episode worden geproduceerd en naar de kern-mantelinggrens zinken, kunnen verantwoordelijk zijn voor een stagnerende fase van enkele honderden miljoenen jaren. En tenslotte kan een meer globale convectiefase worden verwacht wanneer het eclogiet bij de kern-mantelingrens warm genoeg is. Over het geheel genomen komen deze kwalitatieve waarnemingen goed overeen met kwantitatieve metingen van volumes TTG gesteenten die geproduceerd worden wanneer (en alleen wanneer) een hoge intrusie-efficiëntie in aanmerking wordt genomen.
Implicaties van de huidige studie voor de dynamica van Venus
Onlangs is gesuggereerd dat Venus – hoewel het geen continentale korst heeft – een betere analogie zou kunnen zijn voor de Archeaanse Aarde dan Jupiters maan Io7,42,43. Daartoe geven we een beknopt overzicht van de huidige kennis over de tektoniek op Venus en vergelijken die met de resultaten van numerieke modellen, als volgt. Het kleine aantal kraters op Venus leidde tot de conclusie dat de planeet in de voorbije miljarden jaren een globale heropflakkering onderging21,42,44,45,46. Venus-achtige gebeurtenissen kunnen waargenomen worden in globale modellen die zowel extrusieve22 als intrusieve korstvormingsmodellen47 gebruiken (zie ook Fig. 2), dus zijn bijkomende geologische beperkingen nodig om een onderscheid te maken tussen korstvormingsscenario’s. Spectrometrische metingen op de landingsplaatsen Venera en Vega toonden aan dat het oppervlak van Venus gedomineerd wordt door basaltisch materiaal48 , dus een bepaalde fractie van de gegenereerde smelt moet geëxtrudeerd zijn op Venus. Er zijn echter ook duidelijke tekenen van lithosfeervervorming, zoals novae49,50, coronae51,52,53,54 en tessera terreinen55, terwijl het bestaan van subductiezones onzeker blijft56. Er is gesuggereerd dat dijkintrusies belangrijk zijn voor de vorming van novae57 en structurele analyse geeft aan dat novae in de loop van de tijd kunnen evolueren tot coronae-achtige structuren50. Modellen voor coronavorming waarbij een dikke en sterke lithosfeer wordt gebogen58,59 hebben een deel van de negen waargenomen topografische groepen van coronae gereproduceerd; deze modellen hebben echter de neiging de grootte van coronae te overschatten. Het model van een dikke en koude lithosfeer is in vraag gesteld door recente studies die een lage elastische dikte van minder dan 20 km suggereren voor grote delen van Venus (wat een dunne en warme lithosfeer impliceert60,61). Een recente studie die rekening hield met deze nieuwe gegevens door smeltintrusie in de lagere korst van Venus te modelleren, toonde aan dat coronae op zelfconsistente wijze kunnen evolueren uit novastructuren62, in overeenstemming met geologische gegevens. Samenvattend lijkt een korstvormingsmodel dat zowel extrusieve als intrusieve korstvorming integreert noodzakelijk om de waargenomen tektoniek op het huidige Venus te verklaren.