Modelação numérica de convecção termoquímica
Convecção termoquímica compressiva é estudada usando o código StagYY23 em geometria de anel esférico bidimensional24. As equações de momento, massa e conservação de energia são resolvidas utilizando um solver MUMPS paralelo disponível no pacote PETSc25. A aproximação da diferença finita é utilizada numa grelha escalonada26 com uma resolução radialmente variável (mais elevada nas camadas de fronteira e na transição de fase 660). O domínio é discretizado por 64 (radial) vezes 512 (lateral) células. As condições de limite de deslizamento livre são utilizadas tanto para os limites exteriores como para os interiores. A temperatura da superfície é fixada em 300 K e a temperatura inicial do núcleo é de 6.000 K.
Um campo de traçador é advançado através da malha utilizando uma interpolação espacial de segunda ordem livre de divergências do campo de velocidade e um esquema de Runge-Kutta de quarta ordem através do tempo. Cada traçador transporta várias quantidades, tais como o conteúdo de água, composição, temperatura e taxa de aquecimento radiogénico. A utilização destes campos em traçadores limita a difusão numérica. Os elementos radiogénicos são particionados durante a fusão: os derretimentos basálticos são 100 vezes mais enriquecidos do que os seus resíduos sólidos. Por simplicidade, a água é considerada como penetrando totalmente nos 10 km superiores e é aconselhada em todo o manto. O modelo composicional trata as rochas sólidas e fundidas como sendo uma combinação linear de basalto (material crustoso) e harzburgite (manto esgotado). O manto primordial começa com uma composição petrológica inicial de 20% basalto e 80% harzburgite, sendo uma mistura de 60% olivina e 40% piroxeno-garnet. Para eficiência numérica, a composição é armazenada em traçadores, enquanto os cálculos de fusão são efectuados na malha. Para este efeito, um campo de composição em malha é calculado nos centros celulares, calculando a média dos marcadores dentro da célula. Na malha, é calculada uma fracção de fusão em cada passo, comparando as condições de pressão-temperatura com uma função de solidus dependente da composição, considerando um calor latente de 600 kJ kg-1. Uma vez calculada a quantidade de derretimento na malha, a massa correspondente de material fundido é criada no campo do traçador. Se o derretimento ocorrer a menos de 300 km de profundidade, os traçadores fundidos de composição totalmente basáltica são transportados para o topo do domínio (erupção) ou para o fundo da crosta (intrusão). A temperatura dos traçadores em erupção é ajustada à temperatura de superfície, formando rapidamente uma litosfera fria. A temperatura dos traçadores intrudidos tem em conta apenas a descompressão adiabática, resultando numa litosfera quente. A coluna de material entre a região da fonte de fusão e o local de intrusão ou erupção é instantaneamente aconselhada para baixo para conservar a massa. O resíduo sólido deixado pelo procedimento de intrusão da erupção é puro harzburgite. Para mais pormenores sobre a implementação ver as ref. 27 e 28. A viscosidade é dependente tanto da profundidade como da temperatura, seguindo a formulação de Arrhenius (ver ref. 18 para detalhes). A densidade é calculada como a soma dos efeitos de pressão, térmico e composicional, incluindo as transições de fase sólida e sólida. As fases Olivina e piroxeno-garnet são tratadas separadamente, pelo que o aumento de densidade associado à transição de fase basal-eclogite pode ser tido em conta. A eclogitização impulsiona o processo de gotejamento observado durante as primeiras centenas de milhões de anos após a cristalização do oceano magma (ver ref. 18 para detalhes).
Partimos do princípio de que a deformação viscosa opera através do processo de difusão fluente, que segue uma lei Arrhenius dependente da temperatura e da pressão:
onde η0 é a viscosidade de referência (1020 Pa s) à pressão zero e temperatura de referência T0 (1.600 K), Δη é um factor utilizado para impor saltos de viscosidade entre as várias camadas i (até 5 ou 6), Ei é a energia de activação na camada i, P é a pressão, Vi é o volume de activação, R é a constante de gás (8.314 J mol-1 K-1) e T é a temperatura absoluta. Vários valores para Ei e Vi são utilizados para o manto superior e inferior e para a camada pós-perovskite próxima do limite do core-mantle29,30,31 (ver Quadro de Dados Estendidos 3). No manto inferior, o volume de activação diminui com o aumento da pressão de acordo com:
Pi é dado na Tabela de Dados Estendida 3 (ou é infinito se não for dado). Um salto de viscosidade de 30 é imposto na transição entre o manto superior e inferior32.
Rendimento plástico na litosfera é calculado usando uma tensão de rendimento τy seguindo a abordagem Byerlee20:
onde f é o coeficiente de fricção. A viscosidade efectiva é calculada como:
>p>where 
é o segundo invariante do tensor de tensão. A viscosidade não depende da fracção de fusão ou composição.
Desde que resolvemos para convecção compressiva, a temperatura adiabática, densidade, condutividade térmica, expansibilidade térmica e capacidade térmica são dependentes da pressão seguindo uma equação de estado de terceira ordem Birch-Murnaghan, que relaciona o módulo de massa com a pressão (ver ref. 31 para uma explicação detalhada). Mais importante para o presente estudo, os materiais do manto são tratados como uma mistura de olivina, piroxeno-garnet e sistemas em fase de fusão, seguindo uma parametrização baseada em dados de física mineral33,34. Cada sistema de fase tem as suas próprias profundidades de transição de fase e propriedades físicas. A densidade de fusão aumenta suavemente com a pressão. As densidades sólidas de olivina e piroxeno-garnet também aumentam suavemente com a profundidade, mas também contêm saltos nas transições de fase, como mostra a Tabela de Dados Estendidos 4. Em particular, a transição de fase basal-eclogite localizada a uma profundidade de cerca de 60 km tem um efeito de primeira ordem sobre a dinâmica da litosfera18. Acima da transição de fase de eclogite, o basalto é mais leve que a olivina por 160 kg m-3. Abaixo da transição de fase, a eclogite formada a partir do basalto torna-se mais densa que a olivina por 190 kg m-3, o que facilmente acciona instabilidades mecânicas no interior da litosfera. No manto mais profundo, tanto a pirogena-garnet como a olivina também apresentam aumentos de densidade com a profundidade nos limites da zona de transição e no topo da camada pós-perovskite (ver Quadro de Dados Estendidos 4). A pressão PPT em que ocorre a transição de fase é dependente da temperatura: PPT = P0 + γT, onde P0 é uma pressão de referência. Os declives de Clapeyron γ de cada transição (ver Quadro de Dados Alargado 4) quantificam este efeito.
Estimativa dos volumes de TTG produzidos
Os volumes de rochas de TTG apresentados na Fig. 2c são obtidos através de um processo em duas etapas. Primeiro calculamos o volume de basalto hidratado que corresponde às condições de pressão-temperatura para a geração de TTG (ver Quadro de Dados Estendidos 1 e ref. 11). Para obter o que seria o volume de rochas TTG em geometria esférica tridimensional, calculamos a espessura média global produzida e multiplicamo-la pela superfície da Terra para cada passo de tempo t. Depois somamos sobre todas as células i do domínio computacional durante um determinado tempo t e obtemos um volume instantâneo de TTG:
onde S é a superfície da Terra, δzi é a espessura que teria a célula se estivesse espalhada por toda a superfície do domínio computacional. 
é uma concentração de água sem dimensão (1 para rocha totalmente hidratada) e 
é uma ‘concentração de basalto’ sem dimensão (1 para basalto, 0 para harzburgite). O volume é dividido por 2 sob a hipótese de que o grau médio de fusão parcial da crosta basáltica hidratada que produz a mistura TTG corresponde a 50%35. A soma celular é realizada apenas para células com concentrações de água e basalto superiores a metade do teor imposto de água superficial (ou seja, contendo uma grande quantidade de crosta basáltica hidratada).
O volume acumulado V(t) mostrado na Fig. 2c é obtido pela integração temporal ponderada do volume instantâneo v(t):
onde VTTG-P é o volume de rochas em condições de pressão relevantes para a formação de TTG (um valor constante) e Vhbs é o volume de basalto hidratado em todo o domínio. Negligenciamos situações raras em que ∂Vhbs/∂t é negativo. A equação (6) assegura que o TTG é produzido apenas quando as condições relevantes para a sua formação são combinadas (v não é zero) e também quando o novo basalto hidratado é efectivamente formado (∂Vhbs/∂t é positivo). Esta estratégia é muito mais robusta do que simplesmente integrar o volume de rocha em condições para a formação de TTG sem qualquer ponderação. Uma situação de convecção de equilíbrio em que não se forma mais basalto mas a geotermia na crosta corresponde a condições de formação de TTG sobrestimaria dramaticamente o volume de TTG efectivamente produzido.
No entanto, a extracção de dados quantitativos das nossas simulações numéricas tem de ser considerada com cautela. O nosso algoritmo simplificado de cálculo do volume de TTG-crust (Métodos) utiliza um grau médio de fusão da crosta basáltica hidratada e não cria e traça directamente regiões de rocha TTG dentro da crosta. Os modelos regionais com rastreio directo de rocha TTG demonstram geralmente taxas de crescimento de crosta mais elevadas10 mas tal rastreio ainda não é viável à escala global devido à falta de resolução numérica. Além disso, o nosso modelo simplificado pode sobrestimar a duração dos episódios de estagnação da tampa (quase 50% do tempo de simulação, Fig. 2), que produzem quantidades negligenciáveis de rocha TTG. Em modelos regionais de maior resolução, a duração de tais episódios é limitada a cerca de 80 milhões de anos9,10 devido ao enfraquecimento localizado da litosfera por processos mágicos, o que não é implementado no nosso modelo global simplificado. Abaixo, discutimos as implicações potenciais do presente estudo para a dinâmica de Vénus.
Disponibilidade de código
Comparação das nossas simulações numéricas com observações de campo
Comparar os nossos resultados numéricos com a variedade e complexidade dos dados de campo dos cratões Archean é ambicioso. Embora possamos avaliar se certas condições reológicas e comportamentos mágicos podem levar à formação de rochas geladas no início da Terra, não podemos certamente explicar a complexidade dos dados petrológicos e geoquímicos relatados nas últimas décadas. Iremos, portanto, concentrar-nos no quadro de primeira ordem acordado na comunidade geológica de campo Archean. Estamos confiantes de que os nossos modelos numéricos dão uma resposta robusta relativamente à alta eficiência de intrusão necessária para produzir uma distribuição razoável de rochas TTG. Contudo, devemos também ter presente que os nossos modelos ainda carecem de processos importantes, o que faria uma comparação detalhada com dados de campo perigosos e provavelmente enganadores, embora a simples mensagem deste estudo ainda fosse válida.
Primeiro, a nossa reologia não é dependente da composição, impedindo-nos de produzir um manto sub-continental muito rígido da litosfera. No entanto, a nossa viscosidade é fortemente dependente da temperatura (como sugerem os dados experimentais), o que por si só gera uma litosfera subcrostal muito viscosa. Esta funcionalidade não trivial só é possível porque utilizamos solventes directos paralelos muito robustos18,25. Em segundo lugar, uma vez que estamos a olhar para mil milhões de anos de evolução do manto global, a resolução nos nossos modelos não pode ser suficientemente alta para estudar os detalhes dos processos de crosta, que são muito valiosos para os geólogos de campo que utilizam a abordagem de cima para baixo (ligando dados de ‘superfície’ para sugerir uma imagem geodinâmica global). Esta importante limitação não nos permite modelar processos regionais da crosta (tais como estruturas de cúpula e quilha) que modificariam o perfil de temperatura na crosta. Tais simulações foram realizadas utilizando simulações tridimensionais de alta resolução9, e sugerem uma imagem semelhante à nossa, mas a uma escala regional. Além disso, a nossa abordagem é certamente válida, uma vez que a inversão da crosta só ocorrerá se a crosta for realmente produzida em primeiro lugar, o que temos a capacidade de modelar. Em terceiro lugar, o nosso modelo petrológico é relativamente simples porque queremos determinar que ambiente geodinâmico global conduzirá a condições favoráveis para a formação de TTG. Finalmente, as nossas simulações foram realizadas em duas dimensões, porque as simulações tridimensionais são computacionalmente muito caras e não nos teriam permitido estudar sistematicamente o nosso espaço de parâmetros. Tendo em mente estas limitações de modelação, comparamos as características observadas nos nossos modelos numéricos com cenários geodinâmicos baseados em dados de campo.
Embora muito poucas rochas arqueanas tenham sido preservadas na Terra, parece amplamente aceite que o modo tectónico mudou muito pouco depois do início do Paleoarqueano, há cerca de 3,55 mil milhões de anos36. De facto, todo o material mais antigo que esta idade (no Eoarchean) pode ser encontrado em terranos gneiss empilhados de alta qualidade, enquanto que as rochas felsicas produzidas mais tarde (durante a maior parte do Paleoarchean) parecem ter surgido de um ambiente tectónico menos catastrófico, suficientemente estável para gerar o diapirismo crustoso16. Não existe actualmente consenso sobre o processo geodinâmico/tectónico que forma as várias estruturas relatadas pelos estudos de campo7,16. As estruturas eoarquianas apresentam deformações intensas e são por vezes interpretadas como terranos de subducção-acriação37 ou o resultado de orogênios quentes38. Os únicos cratões conhecidos deste período foram encontrados no complexo Itsaq gneiss na Gronelândia Ocidental e na Província Superior no Canadá7. Material paleo-arqueano felsic pode ser encontrado no cratão Kaapvaal (principalmente na África do Sul) e no cratão Pilbara (Austrália Ocidental). Em ambos os locais, a crosta parece ter sofrido uma reviravolta interna, formando cúpulas TTG rodeadas de ‘quilhas’ komatiíticas e basálticas cheias de sedimentos. Pensa-se que o aparecimento destas geometrias de cúpula e quilha só é possível devido à formação de um manto sub-continental litosférico muito viscoso, suficientemente forte para manter a crosta félica no lugar, apesar do seu diapirismo39. Além disso, parece que o plutonismo félico gerador das estruturas da cúpula e da quilha está ligado a vários impulsos de magmatismo que se pensa serem devidos a uma pluma estacionária abaixo do cratão40. Assim, os geólogos de campo relatam dois tipos principais de material Arqueano: Os terráqueos empilhados de Eoarchean mostrando uma deformação muito intensa e as estruturas de cúpula e quilha de Paleoarquean que aparentemente requerem uma convecção mais moderada.
A fase aparentemente quiescente do período Paleoarqueano parece terminar no início da Mesoarquean (3,2 mil milhões de anos atrás), altura em que mais deformação é observada. Isto foi interpretado como o início da tectónica de placas7 , mas esta questão está para além do âmbito da presente Carta, e apenas iremos notar aqui que a deformação importante é relatada como tendo começado há 3,2 mil milhões de anos atrás.
As nossas simulações numéricas estão em bom acordo com o resultado principal destas observações de campo. Dados Ampliados Fig. 1 representa uma comparação das nossas simulações (painel inferior) com os dados de campo anteriormente mencionados (painel superior) ao longo do tempo. De acordo com os terranos empilhados milonitizados do Eoarchean na Gronelândia Ocidental e na Província Superior, observamos sempre uma fase inicial de deformação intensa nas nossas simulações. Este período vigoroso pode ser visto na Fig. 2a (esquerda, a fase de ‘gotejamento’) e dura de cem milhões de anos até setecentos milhões de anos, dependendo da reologia da litosfera (ver regiões vermelhas e acastanhadas no fundo dos dados alargados Fig. 1). A fase de gotejamento ocorre como resultado dos seguintes pressupostos do modelo: a temperatura central inicial muito elevada41 produz plumas intensas e o manto superior é totalmente enriquecido. Uma grande quantidade de crosta mafiosa é produzida durante estas centenas de milhões de anos. A seguir, observamos tanto gotejamento de pedaços eclogíticos de centenas de quilómetros de tamanho como tectónica de superfície tipo proto-subducção, como se pode ver na Fig. 2. Tal comportamento sugere que a geodinâmica eoarqueana pode ter sido uma combinação de processos tectonomagmáticos verticais e eventos de curta duração tipo placa-tectónica que seriam difíceis de interpretar aplicando apenas um paradigma. A maioria das rochas TTG de baixa e média pressão são formadas durante este período de tempo.
A fase de gotejamento é sempre seguida por um período de “tampa estagnada” que ofereceria a estabilidade necessária para a formação da cúpula e estruturas de quilha (regiões amarelas no painel inferior dos Dados Estendidos Fig. 1). Além disso, observamos uma surpreendente estabilidade espacial das plumas, devido ao assentamento e ao aquecimento lento dos enormes blocos eclogíticos no limite do core-mantle. Formam-se principalmente TTGs de alta pressão durante esta fase, embora se produza também uma pequena quantidade de TTGs de média pressão. Pode esperar-se aqui que eventos regionais de inversão de crosta causariam a formação de mais TTG de média e baixa pressão.
Interessantemente, eventos de ressurfacing muito intensos (regiões roxas em Dados Estendidos Fig. 1) seguem por vezes a fase quiescente, reciclando a maior parte da crosta mafiosa ainda presente. Dados alargados Fig. 1 mostra que grandes eventos de repavimentação ocorrem após 800 milhões de anos de evolução da Terra (ou seja, há exactamente 3,2 mil milhões de anos) para uma eficiência de intrusão de 90%. No entanto, este episódio secundário de ressurfacing devido à reactivação da convecção global (uma vez que os blocos de eclogite no limite do core-mantle são suficientemente quentes) pode acontecer em vários momentos da evolução da Terra. Globalmente, as simulações com uma eficiência de intrusão de 90% geram transições geodinâmicas consistentes com o registo geológico geral e produzem as quantidades esperadas de rochas TTG.
Em resumo, mostramos que as nossas simulações não só confirmam a importância do magmatismo intrusivo para a produção de crosta proto-continental, como também estão de acordo com as observações geológicas para as crateras de Arqueano. No entanto, são necessárias novas melhorias dos modelos (e em curso) para uma comparação mais detalhada dos nossos modelos com os dados de campo. Esperamos que poucos dos TTGs de baixa e média pressão produzidos durante a fase de gotejamento sejam preservados (o que ainda não podemos ver nos nossos modelos) e que mais destes TTGs possam ser produzidos durante o processo de cúpula e quilha na fase de quiescência. Mostrámos que a fase de gotejamento inicial gera uma mistura de tectónica horizontal vertical e de curta duração durante cerca de centenas de milhões de anos. Os blocos de material mafioso produzidos durante este episódio inicial e afundando até ao limite do core-mantle podem ser responsáveis por uma fase estagnada de algumas centenas de milhões de anos. E finalmente, uma fase de convecção mais global pode ser esperada quando a eclogite no limite do core-mantle estiver suficientemente quente. Globalmente, estas observações qualitativas estão em boa sintonia com as medições quantitativas dos volumes de rochas TTG produzidas quando (e apenas quando) é considerada uma elevada eficiência de intrusão.
Implicações do presente estudo para a dinâmica de Vénus
Tem sido sugerido recentemente que Vénus – embora carecendo de crostas continentais – pode ser um melhor análogo para a Terra Arqueana do que a lua de Júpiter Io7,42,43. Para este efeito, compilamos brevemente os conhecimentos actuais sobre a tectónica de Vénus e comparamo-los com os resultados dos modelos numéricos, como se segue. O pequeno número de crateras em Vénus levou à conclusão de que o planeta viveu um evento global de ressurgimento nos últimos mil milhões de anos21,42,44,45,46. Eventos de ressurfacing globais semelhantes aos de Vénus podem ser observados em modelos globais que tomam modelos tanto extrusivos22 como intrusivos de formação de crosta47 (ver também Fig. 2), pelo que são necessárias restrições geológicas adicionais para distinguir entre cenários de formação de crosta. As medições espectométricas nos locais de aterragem de Venera e Vega mostraram que a superfície de Vénus é dominada por material basáltico48 , pelo que uma certa fracção do derretimento gerado deve ser extrudido em Vénus. Contudo, existem também sinais claros de deformação da litosfera, tais como novae49,50, coronae51,52,53,54 e terrenos de tessera55, enquanto a existência de zonas de subducção permanece incerta56. As intrusões de diques têm sido sugeridas como importantes para a formação de novae57 e a análise estrutural indica que as novae podem evoluir para estruturas do tipo corona ao longo do tempo50. Os modelos de formação de coroas envolvendo flexão de uma litosfera espessa e forte58,59 reproduziram parte dos nove grupos topográficos de coroas observados; contudo, estes modelos tendem a sobrestimar o tamanho das coroas. O modelo de uma litosfera espessa e fria tem sido desafiado por estudos recentes que sugerem uma baixa espessura elástica de menos de 20 km para grandes partes de Vénus (o que implica uma litosfera de Vénus fina e quente60,61). Um estudo recente tendo em conta estes novos dados através da modelização da intrusão de fusão na crosta inferior de Vénus mostrou que as coroas podem evoluir de forma autoconsistente a partir de estruturas de novae62 , em concordância com os dados geológicos. Em conclusão, um modelo de formação de crosta incorporando tanto a formação de crosta extrusiva como a intrusiva parece ser necessário para explicar a tectónica observada na actual Vénus.