Conference Agenda

Advancements in numerical methods and the growing computational power have allowed us to incorporate numerical forward models into geodynamic inverse problems. We now have several strategies to constrain the rheological properties of the crust and lithosphere. However, among these approaches, the geometric uncertainties of deep geological structures have not been a primary objective. The initial geometry of geological formations (e.g., salt bodies, magma bodies, crust, lithosphere) has either been assumed to be constant or, has been represented by simplified geometries (e.g., planes, spheres, ellipsoids). Given the number of geophysical and geological constraints, it is quite an effort to account for the uncertainties within a geological interpretation of these structures.

In this study, we present a method to describe and manipulate complex, three-dimensional bodies with as few parameters as possible, which allows us to incorporate the geometry of the geological structures into the inversion process. Our approach relies on a partly automated procedure to create an entire ensemble of possible input geometries to account for uncertain geometrical constraints.

The formation of the alpine mountain belt is an indicator of the active geological processes in the Mediterranean that causes earthquakesin Italy and other hazards. To understand this complex system the European project AlpArray, with the German contribution 4D-MB, was funded in order to investigate the structure, dynamics and geology of the alpine area in more detail and through all scales.
We focus on the large scale geodynamic processes that drive this complex system of multiple subduction zones, ranging from the surface to the mantle. The preliminary geodynamic modelling results are based on large scale seismic models of the Mediterranean. The main aim is to investigate the effect of each subduction zone on the observable surface deformation. Additionally, the gradients of these surface observables will be computed using an adjoint technique. These gradients can be used to highlight the pointwise sensitivity of the surface velocity to the material parameters at depth. Starting from a reference model we perform gradient based inversions to better fit the surface observations.

State-of-the-art lithosphere dynamics simulations take (nonlinear) visco-elasto-plastic rheologies into account, along with many different rocktypes. As a result, the simulations typically have a large number of parameters (often >50), all of which have uncertainties. By running systematic forward simulations in which the parameters are changed one obtains an understanding of the relative importance of each of them. With experience, computational geodynamicists have built up an intuition for which parameters are usually important, such as the effective viscosity or density of the lithosphere and mantle. Yet, for a different case study or setup, this may change. Performing such parameter studies thus require many forward simulations which is time-consuming, particularly in 3D.

It would thus be useful to have an additional method that gives some automated guidance into which model parameters are important for a particular model simulation. We discuss here that this can be done by computing scaling law exponents, which relies on computing the gradient of the velocity at a certain observation point with respect to all material parameters in the model. This gradient can be computed using a standard finite difference approach, but it can also be computed using the adjoint method which is computationally very efficient. We will illustrate how this works using both simple examples (falling blocks, folding, faulting), as well as complicated 3D models of magmatic systems. The results indeed show that the largest exponents play a key role in the model dynamics, similar to what we found by manually varying arameters which gives useful guidance on how to perform parametric studies.

Another application is to compute the misfit between modelled and observed velocities at certain observation points. In this case, we can also employ the adjoint method to compute the gradient of the misfit value with respect to the model parameters. By embedding this in a gradient-based algorithm we can perform an inversion to fit geodynamic models to observations, which include velocities, strain rate tensor values or principal stress directions.

All these methods have been implemented in the 3D open-source, lithosphere dynamics code LaMEM, available from bitbucket.

The Central Alps are one of the continental collision zones that cannot be explained by conventional collision tectonics. The striking isostatic disequilibrium inspired the formulation of the slab rollback orogenesis. This model revolves around the notion that vertical forces dominate over horizontal forces. It has been applied to the Alps kinematically, and to some extent in a dynamic context as well. This research uses quantitative (seismo-)thermomechanical models to analyse the physical viability of the slab rollback orogeny model for the Central Alps. Specifically, the role of the continental lower crust is analysed, along with sensitivity studies. Approximately one hundred model runs indicate several parameters to be of key importance, such as frictional strength of sediments versus upper crust. Strong sediments generate a decoupled orogen, while weak sediments promote upper crustal subduction. Additionally, a plagioclase lower crust promotes large displacement along a retroshear, such as the Insubric Line. An Alpine depth of slab detachment (≈200 km) was obtained in very few models, as no Peierls creep was used in most models. Incorporation of this mechanism would help reaching those depths. Convergence rate affects mostly the accretionary wedge and slab temperatures. 5 cm/yr gives very similar results to 10 cm/yr, while slower models generate Rayleigh-Taylor instabilities. Most rollback was achieved in models with decoupled orogens, while coupled orogens result in less, but still measurable, amounts. Finally, the initial temperature distribution is important. A moho temperature ≥ 450°C is required to decouple the upper from the lower continental crust.

A number of Mediterranean subduction systems are assumed to be the product of narrow oceanic lithosphere subductions. In fully coupled two-dimensional thermo-mechanical numerical models such oceanic domains are too narrow to generate the features that can be observed in these subduction systems.

Modelling results show that extent of the Ceahlau-Severin Ocean most commonly proposed by paleogeographic reconstruction is incapable of bringing about roll-back during subduction.

To properly explain the features of the Carpathians Mountain Range we started to work towards developing a numerical model, where the oceanic domain has an inherited component from a previous orogenic evolution.

The model setup was constrained by contemporary mantle structure and geodynamic reconstructions. Several parameters, such as the total kinematic convergence, temperature distribution, and crustal rheological parameters, were explored.

Extension of the upper plate occurs due to the roll-back of the oceanic slab in nearly all of our models, this demonstrates that the inherited component is capable of explaining the roll-back subduction of narrow oceans.

Slab detachment occurs in models, where the continental crust is weaker, however, detachment happens before continental collision, which is not in line with established knowledge of the region.

In our contribution we discuss the most important results of our efforts and report on the issues that still need to be addressed. We examine the possibility of seismo-thermo-mechanical modelling based on our optimal model, directed at unraveling the seismic cycle of the Vrancea-zone under the Eastern Carpathians.

Convergent plate boundaries, and in particular the interface between subducting and overriding plates, are responsible for a large quantity of the global seismicity. Deformation at and above the subduction interface can be categorized in different stages, which are part of the so-called subduction seismic cycle. This cycle includes a stress accumulation phase, a seismic event or earthquake, and a subsequent phase of relaxation. An ideal model thus takes tectonic loading, rate dependent friction, and viscoelastic stress relaxation into account.

This study aims at better understanding the interaction between frictional behaviour in the seismogenic zone and long term deformation structures. Hence, this study focusses on the velocity weakening section of a subduction zone megathrust and the overlying wedge. The approach includes physical analogue modelling of the upper 10 km of a subduction zone and a quantitative study of frictional behaviour of analogue granular materials at experimental conditions. Analogue models are monitored by digital cameras, a load cell to measure loading forces, and a velocity meter to measure displacement and velocity of the downgoing plate. PIVlab and Python based scripts are used to link seismic behaviour to slip on individual faults and to the large scale deformation.

Preliminary results show the formation of splay faults (linking subduction fault to the surface), frontal thrusts, and proto-backthrusts. (Re-)activation of these faults depends on slip size and therefore on the frictional properties of the seismogenic zone and loading rate.

The in-depth behaviour of magmatic systems is still poorly constrained due to their lack of accessibility and the difficulty of finding good analogue representations. Numerical models are thus useful to better understand and interpret these constraints. The ERC-funded MAGMA project aims to develop tools and software for studying a range of magmatic processes in the lithosphere. On the way to building a general framework able to model the behaviour of a fluid-solid chemically coupled magmatic system, we here give an overview of the current development of a mechanical staggered-grid finite difference code using robust analytical linear and non-linear solvers via the PETSc infrastructure, that is able to run in parallel on high-performance computers. The mesh is assembled using the recently developed DMStag framework, which is part of PETSc. This code solves the Stokes equations for elasto-visco-plastic rheologies. We will present the equations and implementations used and show initial results and benchmarks.

A recent focus of studies in geodynamic modeling and magmatic petrology is to understand the coupled behavior between deformation and magmatic processes. Here, we present a 2D numerical model of an upper crustal magma (or mush) chamber in a visco-elastic host rock, with coupled thermal, mechanical and chemical (TMC) processes. The magma chamber is isolated from deeper sources of magma and it is cooling, and thus shrinking. We quantify the mechanical interaction between the shrinking magma chamber and the surrounding host rock, using a compressible visco-elastic formulation, considering several geometries of the magma chamber.

We present a self-consistent system of the conservation equations for coupled TMC processes, under the assumptions of slow (negligible inertial forces), visco-elastic deformation and constant chemical bulk composition. The thermodynamic melting/crystallization model is based on a pelitic melting model calculated with Perple_X, assuming a granitic composition and is incorporated as a look-up table. We will discuss the numerical implementation, show the results of systematic numerical simulations, and illustrate the effect of volume changes due to crystallization on stresses in the host rocks.