Modelización numérica de la convección termoquímica
La convección compresible termoquímica se estudia utilizando el código StagYY23 en una geometría de anillo esférico bidimensional24. Las ecuaciones de conservación de momento, masa y energía se resuelven utilizando un solver MUMPS paralelo disponible en el paquete PETSc25. La aproximación por diferencias finitas se utiliza en una malla escalonada26 con una resolución que varía radialmente (mayor en las capas límite y en la transición de fase 660). El dominio se discretiza mediante 64 (radial) por 512 (lateral) celdas. Se utilizan condiciones de contorno de deslizamiento libre para los límites exteriores e interiores. La temperatura de la superficie se fija en 300 K y la temperatura inicial del núcleo es de 6.000 K.
Un campo de trazadores es advectado a través de la malla utilizando una interpolación espacial de segundo orden libre de divergencias del campo de velocidad y un esquema Runge-Kutta de cuarto orden a través del tiempo. Cada trazador lleva varias cantidades, como el contenido de agua, la composición, la temperatura y la tasa de calentamiento radiogénico. El uso de estos campos en los trazadores limita la difusión numérica. Los elementos radiogénicos se particionan durante la fusión: los fundidos basálticos están 100 veces más enriquecidos que su residuo sólido. Para simplificar, se considera que el agua penetra completamente en los 10 km superiores y es advectada por todo el manto. El modelo composicional trata las rocas sólidas y fundidas como una combinación lineal de basalto (material de la corteza) y harzburgita (manto agotado). El manto primordial comienza con una composición petrológica inicial de 20% de basalto y 80% de harzburgita, siendo una mezcla de 60% de olivino y 40% de piroxeno-granate. Para una mayor eficiencia numérica, la composición se almacena en los trazadores, mientras que los cálculos de fusión se realizan en la malla. Para ello, se calcula un campo de composición en la malla en los centros de las celdas, promediando los trazadores dentro de la celda. En la malla, se calcula una fracción de fusión en cada paso de tiempo, comparando las condiciones de presión-temperatura con una función de solidus dependiente de la composición, considerando un calor latente de 600 kJ kg-1. Una vez calculada la cantidad de fusión en la malla, se crea la correspondiente masa de material fundido en el campo de trazado. Si la fusión se produce a menos de 300 km de profundidad, los trazadores fundidos de composición totalmente basáltica son transportados hacia la parte superior del dominio (erupción) o hacia el fondo de la corteza (intrusión). La temperatura de los trazadores erupcionados se ajusta a la temperatura de la superficie, formando rápidamente una litosfera fría. La temperatura de los trazadores intruidos sólo tiene en cuenta la descompresión adiabática, lo que da lugar a una litosfera caliente. La columna de material entre la región de la fuente de fusión y el lugar de la intrusión o erupción se advierte instantáneamente hacia abajo para conservar la masa. El residuo sólido que deja el procedimiento de erupción-intrusión es harzburgita pura. Para más detalles sobre la aplicación, véanse las referencias 27 y 28. La viscosidad depende tanto de la profundidad como de la temperatura siguiendo la formulación de Arrhenius (ver ref. 18 para más detalles). La densidad se calcula como la suma de los efectos de presión, térmicos y de composición, incluyendo las transiciones de fase sólido-sólido. Las fases de olivino y piroxeno-granate se tratan por separado, de modo que el aumento de densidad asociado a la transición de fase basalto-eclogita puede tenerse en cuenta. La eclogitización impulsa el proceso de goteo observado durante los primeros cientos de millones de años tras la cristalización del océano magmático (véase la ref. 18 para más detalles).
Suponemos que la deformación viscosa opera a través del proceso de fluencia por difusión, que sigue una ley de Arrhenius dependiente de la temperatura y la presión:
Donde η0 es la viscosidad de referencia (1020 Pa s) a presión cero y temperatura de referencia T0 (1.600 K), Δη es un factor utilizado para imponer saltos de viscosidad entre las distintas capas i (hasta 5 o 6), Ei es la energía de activación en la capa i, P es la presión, Vi es el volumen de activación, R es la constante de los gases (8.314 J mol-1 K-1) y T es la temperatura absoluta. Se utilizan varios valores para Ei y Vi para el manto superior e inferior y para la capa postperovskita cercana al límite entre el núcleo y el manto29,30,31 (véase la Tabla 3 de datos ampliados). En el manto inferior, el volumen de activación disminuye con el aumento de la presión según:
Pi se da en la Tabla de Datos Extendidos 3 (o es infinito si no se da). Se impone un salto de viscosidad de 30 en la transición entre el manto superior y el inferior32.
La cesión plástica en la litosfera se calcula utilizando una tensión de fluencia τy siguiendo el enfoque de Byerlee20:
donde f es el coeficiente de fricción. La viscosidad efectiva se calcula como:
donde 
es la segunda invariante del tensor de velocidad de deformación. La viscosidad no depende de la fracción de fusión ni de la composición.
Dado que resolvemos para la convección compresible, la temperatura adiabática, la densidad, la conductividad térmica, la expansividad térmica y la capacidad calorífica dependen de la presión siguiendo una ecuación de estado de tercer orden de Birch-Murnaghan, que relaciona el módulo de volumen con la presión (ver ref. 31 para una explicación detallada). Lo más importante para el presente estudio es que los materiales del manto se tratan como una mezcla de sistemas de fase de olivino, piroxeno-granate y fundido, siguiendo una parametrización basada en datos de física mineral33,34. Cada sistema de fase tiene sus propias profundidades de transición de fase y propiedades físicas. La densidad del fundido aumenta suavemente con la presión. Las densidades del olivino sólido y del piroxeno-granate también aumentan suavemente con la profundidad, pero también contienen saltos en las transiciones de fase, como se muestra en la Tabla 4 de Datos Extendidos. En particular, la transición de fase basalto-eclogita situada a una profundidad de unos 60 km tiene un efecto de primer orden en la dinámica de la litosfera18. Por encima de la transición de fase de la eclogita, el basalto es más ligero que el olivino en 160 kg m-3. Por debajo de la transición de fase, la eclogita formada a partir del basalto se vuelve más densa que el olivino en 190 kg m-3, lo que impulsa fácilmente las inestabilidades mecánicas dentro de la litosfera. En el manto más profundo, tanto el piroxeno-granate como el olivino también muestran aumentos de densidad con la profundidad en los límites de la zona de transición y en la parte superior de la capa post-perovskita (véase la Tabla de Datos Extendidos 4). La presión PPT a la que se produce la transición de fase depende de la temperatura: PPT = P0 + γT, donde P0 es una presión de referencia. Las pendientes de Clapeyron γ de cada transición (ver Tabla de Datos Extendidos 4) cuantifican este efecto.
Estimación de los volúmenes de TTG producidos
Los volúmenes de rocas TTG presentados en la Fig. 2c se obtienen mediante un proceso de dos pasos. Primero calculamos el volumen de basalto hidratado que se ajusta a las condiciones de presión-temperatura para la generación de TTG (véase la Tabla 1 de Datos Extendidos y la ref. 11). Para obtener lo que sería el volumen de las rocas TTG en una geometría esférica tridimensional, calculamos el espesor global medio producido y lo multiplicamos por la superficie de la Tierra para cada paso de tiempo t. Luego sumamos sobre todas las celdas i del dominio computacional para un determinado tiempo t y obtenemos un volumen instantáneo de TTG:
Donde S es la superficie de la Tierra, δzi es el espesor que tendría la celda i si se extendiera por toda la superficie del dominio computacional. 
es una concentración de agua adimensional (1 para roca totalmente hidratada) y 
es una «concentración de basalto» adimensional (1 para basalto, 0 para harzburgita). El volumen se divide por 2 bajo el supuesto de que el grado medio de fusión parcial de la corteza basáltica hidratada que produce la fusión TTG corresponde al 50%35. La suma de celdas se realiza sólo para las celdas que tienen concentraciones tanto de agua como de basalto superiores a la mitad del contenido de agua superficial impuesto (es decir, que contienen una gran cantidad de corteza basáltica hidratada).
El volumen acumulado V(t) mostrado en la Fig. 2c se obtiene por integración temporal ponderada del volumen instantáneo v(t):
donde VTTG-P es el volumen de las rocas en las condiciones de presión relevantes para la formación del TTG (un valor constante) y Vhbs es el volumen del basalto hidratado en todo el dominio. Despreciamos las raras situaciones en las que ∂Vhbs/∂t es negativo. La ecuación (6) asegura que el TTG se produce sólo cuando se dan las condiciones relevantes para su formación (v es distinto de cero) y también cuando se forma realmente nuevo basalto hidratado (∂Vhbs/∂t es positivo). Esta estrategia es mucho más robusta que la simple integración del volumen de roca en condiciones de formación de TTG sin ninguna ponderación. Una situación de convección de equilibrio en la que no se forme más basalto pero la geotermia en la corteza se corresponda con las condiciones de formación de TTGs sobreestimaría drásticamente el volumen de TTGs realmente producido.
Sin embargo, la extracción de datos cuantitativos de nuestras simulaciones numéricas ha de considerarse con precaución. Nuestro algoritmo simplificado de cálculo del volumen de TTG en la corteza (Métodos) utiliza un grado medio de fusión de la corteza basáltica hidratada y no crea ni traza directamente regiones de roca TTG dentro de la corteza. Los modelos regionales con trazado directo de la roca TTG suelen demostrar mayores tasas de crecimiento de la corteza10 , pero dicho trazado aún no es factible a escala global debido a la falta de resolución numérica. Además, nuestro modelo simplificado puede sobreestimar la duración de los episodios de estancamiento de la tapa (casi el 50% del tiempo de simulación, Fig. 2), que producen cantidades insignificantes de roca TTG. En los modelos regionales de mayor resolución, la duración de estos episodios se limita a unos 80 millones de años9,10 debido al debilitamiento localizado de la litosfera por procesos magmáticos, que no está implementado en nuestro modelo global simplificado. A continuación, discutimos las potenciales implicaciones del presente estudio para la dinámica de Venus.
Disponibilidad de código
Comparación de nuestras simulaciones numéricas con las observaciones de campo
Comparar nuestros resultados numéricos con la variedad y complejidad de los datos de campo de los cratones arcaicos es ambicioso. Si bien podemos evaluar si ciertas condiciones reológicas y comportamientos magmáticos pueden conducir a la formación de rocas félsicas en la Tierra primitiva, ciertamente no podemos explicar la complejidad de los datos petrológicos y geoquímicos reportados en las últimas décadas. Por lo tanto, nos centraremos en la imagen de primer orden acordada en la comunidad de geología de campo del Arcaico. Confiamos en que nuestros modelos numéricos ofrezcan una respuesta robusta respecto a la alta eficiencia de intrusión necesaria para producir una distribución razonable de las rocas TTG. Sin embargo, también debemos tener en cuenta que a nuestros modelos todavía les faltan procesos importantes, lo que haría que una comparación detallada con los datos de campo fuera peligrosa y probablemente engañosa, aunque el simple mensaje de este estudio seguiría siendo válido.
En primer lugar, nuestra reología no depende de la composición, lo que nos impide producir un manto litosférico subcontinental muy rígido. Sin embargo, nuestra viscosidad es fuertemente dependiente de la temperatura (como sugieren los datos experimentales), lo que en sí mismo genera una litosfera subcustal muy viscosa. Esta funcionalidad no trivial sólo es posible porque utilizamos solvers directos paralelos muy robustos18,25. En segundo lugar, dado que estamos analizando mil millones de años de evolución global del manto, la resolución de nuestros modelos no puede ser lo suficientemente alta como para estudiar los detalles de los procesos de la corteza, que son muy valiosos para los geólogos de campo que utilizan el enfoque descendente (que vincula los datos «superficiales» para sugerir una imagen geodinámica global). Esta importante limitación no permite modelar los procesos regionales de la corteza (como las estructuras de domos y quillas) que modificarían el perfil de temperatura en la corteza. Tales simulaciones se han realizado mediante simulaciones tridimensionales de alta resolución9, y sugieren una imagen similar a la nuestra, pero a escala regional. Además, nuestro enfoque es ciertamente válido, ya que el vuelco de la corteza sólo se producirá si la corteza se produce realmente en primer lugar, lo que tenemos la capacidad de modelar. En tercer lugar, nuestro modelo petrológico es relativamente sencillo porque queremos determinar qué entorno geodinámico global dará lugar a condiciones favorables para la formación de TTG. Por último, nuestras simulaciones se realizaron en dos dimensiones, porque las simulaciones tridimensionales son muy costosas desde el punto de vista informático y no nos habrían permitido estudiar sistemáticamente nuestro espacio de parámetros. Teniendo en cuenta estas limitaciones inherentes a la modelización, comparamos las características observadas en nuestros modelos numéricos con los escenarios geodinámicos basados en datos de campo.
Aunque se han conservado muy pocas rocas del Arcaico en la Tierra, parece ampliamente aceptado que el modo tectónico cambió en gran medida poco después del comienzo del Paleoarqueano, hace unos 3.550 millones de años36. En efecto, todo el material anterior a esta edad (en el Eoarqueo) se encuentra en terranos de gneis apilados de alto grado, mientras que las rocas félsicas producidas posteriormente (durante la mayor parte del Paleoarqueo) parecen haber surgido de un entorno tectónico menos catastrófico, lo suficientemente estable como para generar diapirismo en la corteza16. En la actualidad no hay consenso sobre el proceso geodinámico/tectónico que forma las distintas estructuras señaladas por los estudios de campo7,16. Las estructuras eoarqueanas muestran una intensa deformación y a veces se interpretan como terranos de subducción-acreción37 o como el resultado de orógenos calientes38. Los únicos cratones conocidos de este periodo se han encontrado en el complejo de gneis Itsaq en el oeste de Groenlandia y en la Provincia Superior de Canadá7. Se puede encontrar material félsico paleoarqueano en el cratón Kaapvaal (principalmente en Sudáfrica) y en el cratón Pilbara (Australia Occidental). En ambos lugares, la corteza ha experimentado aparentemente un vuelco interno, formando domos TTG rodeados de «quillas» komáticas y basálticas rellenas de sedimentos. Se cree que la aparición de estas geometrías de domos y quillas sólo es posible debido a la formación de un manto litosférico subcontinental muy viscoso, lo suficientemente fuerte como para mantener la corteza félsica en su lugar a pesar de su diapirismo39. Además, parece que el plutonismo félsico que genera las estructuras de domos y quillas está vinculado a varios pulsos de magmatismo que se cree que se deben a una pluma estacionaria por debajo del cratón40. Así, los geólogos de campo informan de dos tipos principales de material arcaico: Los terranos apilados del Eoarqueo que muestran una deformación muy intensa y las estructuras de cúpula y quilla del Paleoarqueo que aparentemente requieren una convección más moderada.
La fase aparentemente quiescente del Paleoarqueo parece terminar al inicio del Mesoarqueo (hace 3.200 millones de años), momento en el que se observa más deformación. Esto se ha interpretado como el inicio de la tectónica de placas7, pero esta cuestión está fuera del alcance de la presente Carta y sólo señalaremos aquí que se informa de que la deformación importante comienza hace 3.200 millones de años.
Nuestras simulaciones numéricas están en buen acuerdo con el resultado principal de estas observaciones de campo. Datos ampliados La Fig. 1 representa una comparación de nuestras simulaciones (panel inferior) con los datos de campo anteriormente mencionados (panel superior) a lo largo del tiempo. De acuerdo con los terranos apilados milonitizados del Eoarqueano en el oeste de Groenlandia y la Provincia Superior, siempre observamos una fase inicial de intensa deformación en nuestras simulaciones. Este período vigoroso puede verse en la Fig. 2a (izquierda, la fase de «goteo») y dura desde cien millones de años hasta hasta setecientos millones de años, dependiendo de la reología de la litosfera (ver las regiones rojas y marrones en la parte inferior de la Fig. 1 de Datos Extendidos). La fase de goteo se produce como resultado de las siguientes suposiciones del modelo: la altísima temperatura inicial del núcleo41 produce plumas intensas y el manto superior está totalmente enriquecido. Durante estos cientos de millones de años se produce una gran cantidad de corteza máfica. Tras esto, se observa tanto el goteo de trozos eclogíticos de cientos de kilómetros de tamaño como una tectónica superficial de tipo proto-subducción, como puede verse en la Fig. 2. Este comportamiento sugiere que la geodinámica eoarqueana podría haber sido una combinación de procesos tectonomagmáticos verticales y eventos de tipo tectónico-placa de corta duración que serían difíciles de interpretar aplicando sólo un paradigma. La mayor parte de las rocas de baja y media presión del TTG se forman durante este periodo de tiempo.
La fase de goteo siempre va seguida de un periodo de «tapa estancada» que ofrecería la estabilidad necesaria para la formación de las estructuras de domo y quilla (regiones amarillas en el panel inferior de la Fig. 1 de Datos Extendidos). Además, observamos una sorprendente estabilidad espacial de las plumas, debido al asentamiento y lento calentamiento de los enormes bloques eclogíticos en el límite entre el núcleo y el manto. Durante esta fase se forman principalmente TTG de alta presión, aunque también se produce una pequeña cantidad de TTG de media presión. Se puede esperar aquí que los eventos regionales de vuelco de la corteza causen la formación de más rocas TTG de media y baja presión.
Interesantemente, eventos de resurgimiento muy intensos (regiones púrpuras en la Fig. 1 de Datos Extendidos) a veces siguen a la fase quiescente, reciclando la mayor parte de la corteza máfica aún presente. La Fig. 1 de Extended Data muestra que los principales eventos de resurgimiento ocurren después de 800 millones de años de evolución de la Tierra (es decir, hace exactamente 3.200 millones de años) para una eficiencia de intrusión del 90%. Sin embargo, este episodio de resurgimiento secundario debido a la reactivación de la convección global (una vez que los bloques de eclogita en el límite entre el núcleo y el manto se calientan lo suficiente) puede ocurrir en varios momentos de la evolución de la Tierra. En general, las simulaciones con una eficiencia de intrusión del 90% generan transiciones geodinámicas consistentes con el registro geológico general y producen las cantidades esperadas de rocas TTG.
En resumen, mostramos que nuestras simulaciones no sólo confirman la importancia del magmatismo intrusivo para la producción de la corteza protocontinental, sino que también están de acuerdo con las observaciones geológicas para los cratones arcaicos. Sin embargo, se requieren más mejoras en los modelos (que están en curso) para una comparación más detallada de nuestros modelos con los datos de campo. Esperamos que se conserven pocos de los TTG de baja y media presión producidos durante la fase de goteo (lo que todavía no podemos ver en nuestros modelos) y que se produzcan más de estos TTG durante el proceso de domo y quilla en la fase de reposo. Hemos demostrado que la fase de goteo inicial genera una mezcla de tectónica vertical y horizontal de corta duración durante unos cientos de millones de años. Los bloques de material máfico producidos durante este episodio inicial y que se hunden hasta el límite entre el núcleo y el manto podrían ser responsables de una fase de estancamiento de unos cientos de millones de años. Y, por último, cabe esperar una fase de convección más global cuando la eclogita en el límite entre el núcleo y el manto se caliente lo suficiente. En general, estas observaciones cualitativas están en buen acuerdo con las mediciones cuantitativas de los volúmenes de rocas TTG producidas cuando (y sólo cuando) se considera una alta eficiencia de intrusión.
Implicaciones del presente estudio para la dinámica de Venus
Recientemente se ha sugerido que Venus -aunque carece de corteza continental- podría ser un mejor análogo de la Tierra Arcaica que la luna Io de Júpiter7,42,43. Para ello, recopilamos brevemente los conocimientos actuales sobre la tectónica de Venus y los comparamos con los resultados de los modelos numéricos. El pequeño número de cráteres en Venus llevó a la conclusión de que el planeta experimentó un evento de resurgimiento global en los últimos mil millones de años21,42,44,45,46. Los eventos de rejuvenecimiento global similares a los de Venus pueden observarse en modelos globales que toman tanto modelos de formación de corteza extrusiva22 como intrusiva47 (véase también la Fig. 2), por lo que son necesarias restricciones geológicas adicionales para distinguir entre los escenarios de formación de corteza. Las mediciones espectrométricas realizadas en los lugares de aterrizaje de Venera y Vega mostraron que la superficie de Venus está dominada por material basáltico48, por lo que una cierta fracción del fundido generado debe ser extruido en Venus. Sin embargo, también hay signos claros de deformación de la litosfera, como novas49,50, coronas51,52,53,54 y terrenos de teselas55, mientras que la existencia de zonas de subducción sigue siendo incierta56. Se ha sugerido que las intrusiones de diques son importantes para la formación de novas57 y el análisis estructural indica que las novas pueden evolucionar hacia estructuras similares a las coronas con el tiempo50. Los modelos de formación de coronas que implican la flexión de una litosfera gruesa y fuerte58,59 reproducen parte de los nueve grupos topográficos observados de las coronas; sin embargo, estos modelos tienden a sobreestimar el tamaño de las coronas. El modelo de una litosfera gruesa y fría ha sido cuestionado por estudios recientes que sugieren un bajo espesor elástico de menos de 20 km para grandes partes de Venus (lo que implica una litosfera de Venus delgada y cálida60,61). Un estudio reciente que tiene en cuenta estos nuevos datos mediante la modelización de la intrusión de fundido en la corteza inferior de Venus mostró que las coronas pueden evolucionar de forma autoconsistente a partir de estructuras de novae62, de acuerdo con los datos geológicos. En conclusión, un modelo de formación de la corteza que incorpore tanto la formación de la corteza extrusiva como la intrusiva parece ser necesario para explicar la tectónica observada en el Venus actual.