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).
