Thickening of cratonic lithosphere: Implications for craton growth and kimberlite eruption trends

2026-03-14

Cratons are thought to be the stable cores of continental lithosphere that have survived for 3000 Myr. Such long term survival is often attributed to the excess thickness and elevated viscosity of the cratonic lithosphere. Yet, the evolution of craton thickness during these 3000 Myr has remained highly debated. Several studies have explored three possible scenarios. First, cratons may have been thicker in the past and thinned to ~200 km in the present day. Second, cratons may have thickened slowly or third, they have maintained their current thickness since their origin. In this study we explore the evolution of craton thickness in the past 3000 Myr using 2-D thermo-mechanical numerical models. We initiate each model with a thick and compositionally lighter (1.5% less dense) craton of 200 km in a hot convecting mantle and let it run for 3000 Myr. We impose periodic compression and extension on the craton to mimic supercontinental cycles. We run a total of 24 models exploring a range of initial thicknesses, density contrasts, radioactive heating, and mantle cooling parameters, in order to test multiple evolutionary scenarios. The main results suggest that due to its lower density, the craton is initially flattened. As the craton cools, thermal density overcomes the compositional density, and the craton thickness increases. Viscosity increases concurrently and the mantle flow is diverted along the cratonic edges to self-compress the craton gradually. Due to periodic compression and extension in the model, craton topography varies within a few hundred meters, consistent with observations suggesting basin opening and erosion during and after the assembly and break up of supercontinents. However, the continental lithosphere remains stable. After 1500 Myr, the craton becomes thicker than 160 km depth, a crucial depth for generating kimberlites. Kimberlites are volatile-rich ultramafic rocks that are generated within a depth range of 160-250 km, and are only found above thick continental cratons. Importantly, most kimberlite ages cluster within the last 300 Myr, and available databases suggest that kimberlites were scarce between 3000 and 2000 Myr. Eruptions began occurring more continuously after ~1500 Ma, and accelerated after ~1100 Ma. This pattern is consistent with our models of a slowly growing craton thickness. We find that before 1500 Ma cratons were mostly thinner than the critical depth for kimberlite generation. After 1500 Myr, their thickness increased, allowing them to host more kimberlites. Although previous hypotheses emphasize mantle temperature and carbon availability as primary controls on kimberlite eruptions in the later part of Earth’s history, our results suggest that craton thickness also exerts a strong control on the eruption of kimberlite magmas.

Widespread felsic volcanism as a possible step towards Archean subaerial landmass: Insights from combined oxygen and hafnium isotopes in zircon

2026-03-14

Subaerial land today is mainly formed by continental crust, but before the stabilization of the first cratons, at ca. 3 Ga, volcanic structures (e.g., oceanic islands) may have been the first subaerial regions of the early Earth. Understanding the onset of felsic magmatism is crucial for constraining the formation of both continental crust and hypothetical early volcanic islands. Studies of ancient zircons suggest that subaerial land likely emerged at least by 3.5 Ga, but how long before that it began remains unknown. Our work examines circa 3.5 Ga old felsic volcanic rocks, the oldest known in the Kaapvaal (South Africa) and Singhbhum (India) cratons. We analyzed oxygen and Lu-Hf isotopes in zircon as they are effective proxies for distinguishing the melt source between mantle-derived and crustal (remelting of altered rocks and sediments). Oxygen isotopes ratios (δ18O) were measured by Secondary Ion Mass Spectrometry (SIMS) in coeval felsic units of the Kaapvaal Craton (i.e. Theespruit, Sandspruit, and Toggekry formations), and of the Singhbhum Craton (Daitari and Gorumahisani greenstone belts). This new data was compared with a newly compiled global Archean δ18O dataset (ca. 13,000 data points). Our felsic volcanic rocks display the averaged δ18O values ranging between 5.1 and 5.8 ± 0.24 ‰ (2 sd), which are purely mantle-like values. The only exception is a Toggekry formation sample (δ18O 3.9 ± 0.24 ‰), which reflects remelting of hydrothermally altered rocks. Published εHf values for the same rocks fall between CHUR and Depleted Mantle trends, implying juvenile melt signatures. In this context, we highlight the significance of the early Earth's first felsic rocks, whose formation is usually attributed to partial melting of a hydrated basaltic oceanic crust. In contrast, our data emphasizes the importance of purely mantle-derived felsic melts in the Archean. These felsic melts can be a result of extensive fractional crystallization (ca. 80%) of a stalled basaltic melt. Such relatively dry melting (possessing only juvenile water) requires elevated heat flow, and thick lithosphere. During the Archean, these conditions may have prevailed in a thick basaltic oceanic plateau setting. Reworking (i.e., melting) of such ancient oceanic plateaus could have led to the renewed generation of felsic melts producing buoyant silicic rocks and ultimately result in the consolidation and emergence of the earliest continental crust. The global Archean δ18O values compilation suggests that the mantle and seawater-altered rocks are both important sources of felsic melts during the Archean. This highlights the significance of global Archaean tectonic regimes that may have led to the formation of the first subaerial landmass in brief stints.

Modelled mantle texture evolution during periodic geodynamic cycles

2026-03-14

Many geodynamic processes on Earth occur with a certain periodicity, which can range from decades to centuries for earthquake cycles to hundreds of millennia for glacial cycles and up to hundreds of millions of years for plate tectonic cycles. From the simplified view of a rock, all of these geodynamic cycles induce a deformation during a loading/opening phase followed by deformation in the opposite direction during the unloading/closing phase. These periodic cycles thus produce deformation without any net strain on the rock.In this work, we use simple models to determine the types of rock texture that can develop within mantle rocks after multiple cycles of dynamic processes, and to understand how such textures can influence the effective viscosity of the mantle.Our simplified setup consists of an olivine polycrystal aggregate ( = our mantle rock) that has an initial (either isotropic or anisotropic) texture at the start of the model. We impose a velocity gradient representing either simple or pure shear in a given direction. The aggregate is sheared with the given velocity gradient for a prescribed amount of strain and then the deformation is reversed. To be impartial, we test the same setting with multiple texture evolution models, including the MDM, the D-REX, the SpecFab and the VPSC models.Our results show that the frequency of deformation cycling and the magnitude of the deformation (in the measure of strain) can dramatically impact both the stability and the type of texture that forms after a few or many deformation cycles. Because these textures are viscously anisotropic, the strain achieved in a deformation cycle thus greatly influences how the mechanical anisotropy of the mantle evolves, and in turn, influences different geodynamic processes. As an example, during glacial cycles one expects small amounts of strain in the mantle ranging from 0.0001 to a maximum of 0.01 units of strain during loading and unloading of ice on the surface. Our results suggest that i) a given piece of mantle needs to experience the same glacial-interglacial cycles hundreds to thousands of times to experience enough strain cycles to develop a significant texture, and ii) when this happens the developed texture is very different than what one would get for continuously deforming the mantle in the same direction. Instead of developing a point maximum in the shear direction, we observe a girdle-type texture with a small maxima normal to the shear plane that remains stable once developed. At lower frequencies, for which shear direction reversals occur less frequently and with larger amounts of strain, the texture does not stabilize. Instead, the texture initially develops toward a point maxima that becomes partially destroyed by the subsequent reverse deformation. Given these trends, we conclude that periodic geodynamic processes may significantly influence the formation of upper mantle rock textures, and that the deformation frequency exerts a particularly important control on the eventual rock texture.

Dataset for the paper titled 'CPO2Hill: A new model linking olivine texture parameters to anisotropic viscous behavior' by Király et al.

Open MIND · 2026-01-15

This dataset includes input files, MATLAB scripts and functions, and results associated with the paper CPO2Hill: A new model linking olivine texture parameters to anisotropic viscous behavior'’ by Király et al., submitted to Geodynamica

Depth rotations of azimuthal seismic anisotropy associated with relative importance of Couette/Poiseuille flow in the asthenosphere

2026-03-13

Azimuthal seismic anisotropy in the upper mantle is crucial for understanding the spatial patterns of past and present upper mantle deformation. Traditional interpretation of such anisotropy attributes to relative shear between surface plates and mantle. This requires the orientation of anisotropy azimuths to remain constant with depth. However, inferences of azimuthal anisotropy based on surface wave tomographic models often reveals depth-dependent azimuths. To this end, the existence of mechanically weak, thin asthenosphere beneath the lithosphere facilitates the channelization of plate-driven Couette flow and pressure-driven Poiseuille flow. The combination of two flows, especially when misaligned, yields depth rotations of asthenospheric shear. This provides a geodynamically plausible link between asthenospheric flow properties and depth rotations of azimuthal seismic anisotropy. In this submission, we utilize publicly available azimuthal seismic anisotropy models together with predictions from a global mantle flow model that incorporates Couette/Poiseuille flow. We find that Poiseuille flow profoundly affects depth rotations of seismically inferred azimuthal anisotropy. Prominent depth rotations are under the Atlantic basin and the Nazca plate, where Poiseuille flow dominates the modeled asthenospheric flow regime. Significant Poiseuille flow may exist beneath the Indian basin, yet with small depth rotation, probably because of its directional alignment with Couette flow. Our results indicate that interpretation of azimuthal seismic anisotropy cannot be simply tied to relative shearing between plates and mantle. Instead, the relative importance of Couette and Poiseuille flows must be taken into account.

Mantle Flow and Anisotropy in Subduction Zones: Modeling and Clustering of Olivine Textures

Geochemistry Geophysics Geosystems · 2025-07-01 · 1 citations

Abstract The mantle near Earth's subduction zones endures intense deformation that generates anisotropic rock textures. These textures can be observed seismically and modeled geodynamically, but the complexity of this deformation makes analyses of these textures difficult. In this study, we apply time‐series clustering analysis to tracers within subduction models, allowing for the identification of regions in the subduction zone with common deformation histories and olivine crystallographic‐preferred orientation development. We compare olivine texture evolution predicted using different methods in both retreating and stationary‐trench settings. Our results reveal distinct variations in olivine texture, indicating that both seismic and viscous anisotropy can exhibit substantial heterogeneity within the mantle wedge, sub‐slab, and subducting plate regions. For retreating trenches, olivine textures are strongest in the mid‐depth mantle wedge region about 200 km away from the trench between 100 and 300 km depth. Our study shows that trench‐normal olivine a ‐axis orientations dominate in the center of subduction zones. Toroidal flow around slab edges generates a mix of trench‐normal, trench‐parallel, and oblique fast seismic directions. Textures and anisotropy are stronger for the retreating trench model than for the stationary trench model since more deformation has been accumulated due to trench motion. These findings provide insights for interpreting seismic anisotropy in subduction zones and highlight the importance of considering texture heterogeneity, as characterized by clustering algorithms, when analyzing both geodynamic models and seismic observations of subduction zones.

Linking Plate Kinematics and True Polar Wander over the last 250 Myrs

2025-03-14

The flux of subducting slabs into the mantle is an essential component of the Earth’s mantle convection. However, the slab flux remains poorly known for pre-Jurassic times because of the absence of preserved oceanic seafloor. Sinking of subducted slabs within the mantle perturbs Earth’s moment of inertia, which, in addition to perturbations related to upwellings, results in long-term motion of the solid Earth relative to the rotation axis, resulting in so-called True Polar Wander (TPW). This motion, which can be inferred using paleomagnetic data, should therefore yield crucial information about the large-scale subduction kinematics back in time. However, it is not yet clear how to separate the numerous contributions to TPW, since these result from the superimposition of a complex distribution of mantle mass heterogeneities that are advected through time. In this study, we developed a new approach to assess the impact of subducting slabs on TPW based on the harmonic decomposition of plate kinematics into large-scale patterns. We constructed simple plate models that yielded pure dipole and pure quadrupole and net stretching kinematics, which represent the spherical harmonic degree 1 and degree 2 components of relative plate motions, respectively. We then implemented these three patterns of large-scale plate motions, and their subduction zones, into three simple mechanistic models and computed mantle mass heterogeneities through time. We then calculated changes to Earth’s moment of inertia tensor to predict the resulting TPW. In this contribution, we will first show the results of these sensitivity experiments highlighting the evolution of inertia perturbations associated to each of these three large-scale patterns. We will then show the calculated TPW using the harmonic decomposition of full-plate models over the last 250 Myrs and discuss the influence of each of these three plate kinematic components on the observed TPW. Finally, we will discuss how the observed TPW can help better constrain the evolution of mantle mass heterogeneities and rates of subduction flux for past times.

Subduction dynamics and mantle anisotropy: modeling and clustering of olivine textures

2025-03-15

The mantle near Earth's subduction zones experiences significant deformation, forming anisotropic rock textures. These textures can be detected using seismic methods and simulated in geodynamic models. This study employs time-series clustering to examine tracers in subduction models, identifying regions with similar deformation histories, olivine crystallographic-preferred orientation (CPO) development, and CPO-induced anisotropic viscosity. We compare the evolution of olivine textures predicted by various numerical methods (e.g. D-Rex, MDM, and MDM+AV) for both retreating and stationary trench subduction settings.Our modeling shows notable variations in olivine texture around the slab and as a function of subduction dynamics. These variations, which are illuminated by the clustering analysis, show that texture, seismic, and viscous anisotropy can vary greatly within the mantle wedge, sub-slab, and subducting plate regions of the upper mantle. In the retreating-trench model, the strongest textures are observed in the mid-depth mantle wedge region and beneath the slab at the 660 km transition zone. Trench-normal olivine a-axis orientations are predominant in the center of subduction zones, while toroidal flow around slab edges produces a mix of trench-normal, trench-parallel, and oblique fast seismic directions. On the other hand, in the stationary-trench model, the trench-normal signal in front of the slab is weaker while there are stronger trench-normal signals behind the slab at shallow depths between 100 and 300 km. At the edge of the slab, weak toroidal flow produces trench-oblique orientations while trench-parallel and trench-normal orientations are missing. In general, the retreating trench model exhibits stronger textures and anisotropy due to increased deformation from trench motion.These results provide valuable insights into seismic anisotropy in subduction zones and underscore the importance of considering texture heterogeneity when interpreting geodynamic models and seismic data. The use of time-series clustering algorithms highlights the intricate pattern of evolution and the relationship between deformation history, CPO, and CPO-induced viscous anisotropy occurring within subduction zones.

Benchmarking dynamic topography across geodynamical codes

2025-03-14

During the ESA funded 4D Dynamic Earth project, different sensitivity studies are performed to understand the applicability of current ground and satellite datasets available to study the dynamical behavior of the solid Earth, in particular the complete mantle. This project is a joint effort between ESA and many European universities and is lead by Delft University of Technology (https://4ddynamicearth.tudelft.nl/). The project consists of ten work packages, many of them relying on some form of forward geodynamical modelling. Given the diversity of participants multiple codes are used in the project: a 2D axisymmetric Python code developed by C.T. at the Utrecht University, a 3D Matlab code developed by O.O-G. and the 3D massively parallel C++ community code ASPECT.One recurring quantity that is of paramount importance for some work packages is dynamic topography, i.e. the outer surface expression to dynamic mantle flow. We have therefore designed a simple isothermal experiment of an anomalous sphere present in the mantle of a planet (the core is ignored as is customary in whole-Earth geodynamic modelling). The sphere itself can be positively or negatively buoyant, and the mantle can be isoviscous or characterized by a radial viscosity profile. Boundary conditions at the core-mantle boundary and at the surface are either no-slip or free-slip. Dynamic topography calculations involve the radial stress which is derived from the primitive variables velocity (actually, its gradient) and pressure which are found to be sensitive to mesh size in both radial and lateral directions. We therefore report on the root mean square velocity, the surface strain rate, stress and dynamic topography and the gravity anomaly for a range of experiments. Our objective is two-fold: characterize the accuracy of our codes and provide the community with a benchmark. All three codes are Finite Element codes and all rely on the Taylor-Hood element but they are also quite different with respect to meshing and solver architecture. Nevertheless we find that all measured quantities converge within approx. 1% for radial resolutions of at least 30km.

ANIMA the journey: how we model olivine CPO-related anisotropic viscosity

2025-03-15

The long-term fluid-like movements in the Earth’s mantle largely depend on the rheological behaviour of olivine, the main rock-forming mineral in the upper mantle. Although the average viscosity of the mantle can be estimated from post-glacial rebound or geoid anomalies, the micromechanical mechanisms that facilitate the deformation of the solid mantle have been identified from rock mechanics experiments. Dislocation creep emerges as the predominant deformation mechanism in the uppermost mantle, aligning olivine crystals into a crystallographic preferred orientation (CPO) parallel to the flow, while this alignment of crystals also results in anisotropic viscous behaviour. Thus, anisotropic viscosity and CPO evolve hand in hand, and this interaction may impact many geodynamic processes. For example, beneath tectonic plates CPO evolves parallel to the plate motion direction, weakening the asthenosphere in that direction. However, if the plate motion direction changes, the asthenosphere will resist this change, leading to smaller velocities, less deformation and therefore a slow evolution of the CPO towards the new plate motion direction. In the ANIMA project, we aimed to find an efficient way of modelling CPO evolution and the related anisotropic viscosity in a fully coupled way within a geodynamic simulation. We developed a method that tracks CPO evolution on advected particles based on the D-REX method and utilizes the eigenvalues of the mean CPO orientation matrices to predict the anisotropic viscous parameters. These parameters allow us to calculate a tensor form of the viscosity, which we then feed back into our model solution. This method can be applied in combination with other rheologies, although with a cost of having to represent the viscosity as a tensor in the entire model domain, regardless of the dominant deformation mechanism. Despite an estimated increase in computational cost by up to an order of magnitude, incorporating anisotropic viscosity coupled to CPO evolution stands feasible for regional geodynamic models. This development will facilitate the study of a broad new range of geodynamics problems that involve olivine texture and anisotropic viscosity.