3D geodynamic modelling of the present-day and long-term deformation of the Alps and Adria

Linking geophysical data with geological constraints to understand the dynamics of Alpine Mountain building was one of the main goals of the 4D-MB SPP project. Whereas seismic tomography inversions give a snapshot of how the seismic velocity structure may look like today, geological constraints give (incomplete) pieces of information of how it may have evolved over time. Linking such information with the physics of the lithosphere requires geodynamic numerical models and due to the geometric complexity of the region, 3D is a must. A problem with geodynamic models is that the driving forces are usually the density differences of subducting plates with the surrounding mantle (or far-field forces), whereas the (uncertain) rheology of mantle and crustal rocks plays a crucial role as well. As such, there are quite a few uncertainties in the input parameters, even if the geometry is well-constrained. During the first phase of the 4D-MB project, we focused on present-day models of the Alpine system. Whereas it is possible to invert for the rheology of the lithosphere in mountain belts using probabilistic Bayesian methods (Baumann and Kaus, 2015), this method requires a large number of forward simulations and thus remains infeasible in 3D. An alternative approach is to employ gradient-based inversion methods, in which the adjoint method is employed as a particularly efficient method to compute the gradient of the misfit of the model and data (usually GPS data) versus model parameters (Reuber et al., 2020; Reuber, 2021). Since the adjoint gradient method is computationally cheap (compared to a forward simulation), it can also be used to quickly determine the key model parameters of a particular simulation (Reuber et al., 2018b) or can be combined with gravity and seismic inversions (Reuber and Simons, 2020). We had initially applied this to a case where the starting forward model setup was already giving a reasonably good fit to the uplift data, in which case the method rapidly converged (Reuber et al., 2018a). It thus seemed straightforward to do the same for the Alps. Yet, several issues were encountered in the process: a) we need to consider a much wider region than just the Alps to avoid issues with the lateral boundaries; b) even with high-resolution P-wave tomographic models at hand, one still needs to interpret the seismic velocity anomalies to create an initial model setup, which is a highly non-unique step and results in various possible interpretations; c) coming up with an initial forward model that gives a reasonably good fit to the GPS data turned out to be a significant challenge. Despite running well over 350 forward simulations, we failed to obtain forward simulations that provided a well-enough fit of the velocity in Adria, and without a good starting model, gradient based geodynamic inversions do not converge to a meaningful solution (Reuber, 2020). More recently, we made another attempt in which seismic velocity was directly translated to density and viscosity anomalies using a simple, linear, scaling law, while also prescribing the far-field velocities at the model boundaries (such as that of the N. Anatolian plate). Results give a better fit in Adria (Fig. 1A), but also show that the details of the slab geometry underneath the Alps do not have much impact on the model results, while the model fit within the Alps remains unsatisfying (perhaps because of the small velocities there). Instead of just focusing at the present day structure of the Alps, it is also interesting to see how the system evolved over the last 20-30 million years, which was the focus of our project during the 2nd phase of 4D-MB. The idea was to start with a plate tectonic reconstruction (Le Breton et al., 2021) and let the model evolve forward in time. As for the present-day models, there are many uncertainties in the plate tectonic reconstructions, such as: What was the slab dip? What were the lengths of the slabs? Were they laterally broken or not? What was the thermal and rheological structure of the plates? Given the difficulties with the present-day models, and the increased computational demands of time-dependent simulations, it is unreasonable to expect model results that magically fit all available constraints. Yet, after performing many hundreds of forward simulations, we do get some consistent results and in some of the simulations Adria moves northward and rotates anticlockwise relative to Europe by about the correct amount. The Gibraltar slab arrives at the correct place (Fig. 1B), and the models clearly show that the northward motion of Africa has little impact on the dynamics of Adria, which is instead mostly driven by the interaction of the Hellenic and Calabrian slabs while being pulled northwards by the retreating Western Alpine slab. The size and thermal structure of the Ionian oceanic lithosphere is important as well. Figure 1: A) Example of present-day geodynamic models, B) Snapshot of a forward geodynamic simulation that started at 30 Ma. We also made various technical advances, which includes the Julia package GeophysicalModelGenerator.jl to create complicated 3D geodynamic model setups from geophysical/geological data, DataPicker.jl which provides a GUI for GMG, and LaMEM.jl which is the Julia interface to LaMEM and allows installing and running LaMEM in parallel (either directly from Julia or via Jupyter or Pluto notebooks). We also extended LaMEM to include a continuous integration, adjoint inversion (Reuber et al., 2020; Reuber, 2021) and sensitivity testing (Reuber et al., 2018b)

Unraveling the physics of the Yellowstone magmatic system using geodynamic simulations

The Yellowstone magmatic system is one of the largest magmatic systems on Earth, and thus an ideal location to study magmatic processes. Whereas previous seismic tomography results could only image a shallow magma reservoir, a recent study using more seismometers showed that a second and massive partially molten mush reservoir exists above the Moho (Huang et al., 2015). To understand the measurable surface response of this system to visco-elasto-plastic deformation, it is thus important to take the whole system from the mantle plume up to the shallow magma reservoirs into account. Here, we employ lithospheric-scale 3D visco-elasto-plastic geodynamic models to test the influence of parameters such as the connectivity of the reservoirs and rheology of the lithosphere on the dynamics of the system. A gravity inversion is used to constrain the effective density of the magma reservoirs, and an adjoint modeling approach reveals the key model parameters affecting the surface velocity. Model results show that a combination of connected reservoirs with plastic rheology can explain the recorded slow vertical surface uplift rates of around 1.2 cm/year, as representing a long term background signal. A geodynamic inversion to fit the model to observed GPS surface velocities reveals that the magnitude of surface uplift varies strongly with the viscosity difference between the reservoirs and the crust. Even though stress directions have not been used as inversion parameters, modeled stress orientations are consistent with observations. However, phases of larger uplift velocities can also result from magma reservoir inflation which is a short term effect. We consider two approaches: (1) overpressure in the magma reservoir in the asthenosphere and (2) inflation of the uppermost reservoir prescribed by an internal kinematic boundary condition. We demonstrate that the asthenosphere inflation has a smaller effect on the surface velocities in comparison with the uppermost reservoir inflation. We show that the pure buoyant uplift of magma bodies in combination with magma reservoir inflation can explain (varying) observed uplift rates at the example of the Yellowstone volcanic system

Neutral Atom Quantum Computing Hardware: Performance and End-User Perspective

We present an industrial end-user perspective on the current state of quantum computing hardware for one specific technological approach, the neutral atom platform. Our aim is to assist developers in understanding the impact of the specific properties of these devices on the effectiveness of algorithm execution. Based on discussions with different vendors and recent literature, we discuss the performance data of the neutral atom platform. Specifically, we focus on the physical qubit architecture, which affects state preparation, qubit-to-qubit connectivity, gate fidelities, native gate instruction set, and individual qubit stability. These factors determine both the quantum-part execution time and the end-to-end wall clock time relevant for end-users, but also the ability to perform fault-tolerant quantum computation in the future. We end with an overview of which applications have been shown to be well suited for the peculiar properties of neutral atom-based quantum computers

Multi-physics adjoint modeling of Earth structure: combining gravimetric, seismic, and geodynamic inversions

We discuss the resolving power of three geophysical imaging and inversion techniques, and their combination, for the reconstruction of material parameters in the Earth’s subsurface. The governing equations are those of Newton and Poisson for gravitational problems, the acoustic wave equation under Hookean elasticity for seismology, and the geodynamics equations of Stokes for incompressible steady-state flow in the mantle. The observables are the gravitational potential, the seismic displacement, and the surface velocity, all measured at the surface. The inversion parameters of interest are the mass density, the acoustic wave speed, and the viscosity. These systems of partial differential equations and their adjoints were implemented in a single Python code using the finite-element library FeNICS. To investigate the shape of the cost functions, we present a grid search in the parameter space for three end-member geological settings: a falling block, a subduction zone, and a mantle plume. The performance of a gradient-based inversion for each single observable separately, and in combination, is presented. We furthermore investigate the performance of a shape-optimizing inverse method, when the material is known, and an inversion that inverts for the material parameters of an anomaly with known shape

Subduction metamorphism in the Himalayan ultrahigh-pressure Tso Moran massif: An integrated geodynamic and petrological modelling approach

The Tso Morari massif is one of only two regions where ultrahigh-pressure (UHP) metamorphism of subducted crust has been documented in the Himalayan Range. The tectonic evolution of the massif is enigmatic, as reported pressure estimates for peak metamorphism vary from similar to 2.4 GPa to similar to 4.8 GPa. This uncertainty is problematic for constructing large-scale numerical models of the early stages of India-Asia collision. To address this, we provide new constraints on the tectonothermal evolution of the massif via a combined geodynamic and petrological forward-modelling approach. A prograde-to-peak pressure-temperature-time (P-T-t)path has been derived from thermomechanical simulations tailored for Eocene subduction in the northwestern Himalaya. Phase equilibrium modelling performed along this P-T path has described the petrological evolution of felsic and mafic components of the massif crust, and shows that differences in their fluid contents would have controlled the degree of metamorphic phase transformation in each during subduction. Our model predicts that peak P-T conditions of similar to 2.6-2.8 GPa and -600-620 °C, representative of 90-100 km depth (assuming lithostatic pressure), could have been reached just similar to 3 Myr after the onset of subduction of continental crust. This P-T path and subduction duration correlate well with constraints reported for similar UHP eclogite in the Kaghan Valley, Pakistan Himalaya, suggesting that the northwest Himalaya contains dismembered remnants of what may have been a similar to 400-km-long UHP terrane comparable in size to the Western Gneiss Region, Norway, and the Dabie-Sulu belt, China. A maximum overpressure of similar to 0.5 GPa was calculated in our simulations for a homogeneous crust, although small-scale mechanical heterogeneities may produce overpressures that are larger in magnitude. Nonetheless, the extremely high pressures for peak metamorphism reported by some workers (up to 4.8 GPa) are unreliable owing to conventional thermobarometry having been performed on minerals that were likely not in equilibrium. Furthermore, diagnostic high-P mineral assemblages predicted to form in Tso Morari orthogneiss at peak metamorphism are absent from natural samples, which may reflect the widespread metastable preservation of lower-pressure assemblages in the felsic component of the crust during subduction. If common in such subducted continental terranes, this metastability calls into question the reliability of geodynamic simulations of orogenesis that are predicated on equilibrium metamorphism operating continuously throughout tectonic cycles

Plume — Lid interactions during the Archean and implications for the generation of early continental terranes

Many Archean terranes are interpreted to have a tectonic and metamorphic evolution that indicates intra-crustal reorganization driven by lithospheric-scale gravitational instabilities. These processes are associated with the production of a significant amount of felsic and mafic crust, and are widely regarded to be a consequence of plume-lithosphere interactions. The juvenile Archean felsic crust is made predominantly of rocks of the tonalite–trondhjemite–granodiorite (TTG) suite, which are the result of partial melting of hydrous metabasalts. The geodynamic processes that have assisted the production of juvenile felsic crust, are still not well understood. Here, we perform 2D and 3D numerical simulations coupled with the state-of-the-art of petrological thermodynamical modelling to study the tectonic evolution of a primitive Archean oceanic plateau with particular regard on the condition of extraction of felsic melts. In our numerical simulations, the continuous emplacement of new, dry mafic intrusions and the extraction of the felsic melts, generate an unstable lower crust which drips into the mantle soon after the plume arrival. The subsequent tectonic evolution depends on the asthenosphere TP. If the TP is high enough (≥ 1500 ∘C) the entire oceanic crust is recycled within 2 Myrs. By contrast at low TP, the thin oceanic plateau slowly propagates generating plate-boundary like features