current research projects
We are interested in understanding the physics of geological and geophysical processes. As these processes typically occur on a million-year timescale and we only observe the end-result, they are best simulated on a computer. We actively develop new numerical approaches for this purpose and apply them to a wide variety of topics, some of which are listed below. See also our publications page.
How did subduction initiate and how does plate tectonics work on Earth?
Earth is the only terrestrial planet in our solar system that currently has plate tectonics. Why this is the case remains incompletely understood, partly because we have not been able to reproduce plate tectonics in the computer with realistic parameters. Recently, we demonstrated that a free surface boundary condition, which was ignored in most previous models of mantle convection, crucially changes the physics of convection models with plates and results in asymmetric subduction. Yet, the yield stresses that have to be employed in such models remain significantly below constraints from rock experiments and from in-situ stress measurements in boreholes.
A different line of research in our group therefore focusses on better understanding how a subduction zone is created in the first place. In a series of papers we have demonstrated that shear-heating is a viable mechanism to create a large-scale, lithospheric-cutting, shear zone. Under some conditions, this also results in a well-developed subduction zone. These lithospheric-scale models involve visco-elasto-plastic rheologies and thermo-mechanical feedbacks, whereas the large-scale convection models use simpler rheologies. Ongoing research focusses on understanding the effects of more realistic rheologies in convection setups (with the aim to test whether shear-heating is a viable localization mechanism in a convecting system as well).
Development of state-of-the-art geodynamic modelling codes
Numerical modelling of geodynamic processes is challenging as the rheology of rocks is nonlinear and changes from viscous (in the mantle where temperatures are higher) to elasto-plastic under surface-near conditions. Moreover, geological processes result in very large deformations and involve processes such as phase transitions and (partial) melting. As existing modelling packages are not suitable or slow, we have developed our own packages over the last few years.
In two dimensions and for moderate resolutions we mainly use MILAMIN_VEP these days, which is a marker-in-cell finite element package that is written in almost native MATLAB, yet is relatively efficient (see www.milamin.org for a viscous version of the solver). Ongoing research in 2D is in coupling two-phase flow with lithospheric deformation, and in developing faster nonlinear solvers.
In three dimensions, the problem is considerably more challenging and we are actively developing software that runs on massively parallel machines and employs efficient multigrid techniques. The main package we are working on, LaMEM (Lithosphere and Mantle Evolution Model) is a finite element/finite difference package build on top of PETSC and we have used it to perform simulations with over 2 billion degrees of freedom on >16'000 cores on the Bluegene/Q at FZ Jülich. Current development is focussed on adding nonlinearities, coupling it to surface processes and improving the robustness of the solvers.
What is the effective rheology of melt?
Understanding the effective rheology of melts is important for developing volcanic eruption models. Experimental results on deformation of partially molten rocks indicate that the rheology is mostly linear for a fully molten rock but becomes extremely nonlinear if crystals are present. Understanding the reasons for this is challenging, which is why we reproduced the experiments in the computer which has the advantage that we can precisely control the rheology which helps to understand the cause of the observed nonlinearities in lab experiments.
How are fold and thrust belts formed?
Fold and thrust belts are of considerable economic importance but how they form mechanically is incompletely understood. Some fold and thrust belts are dominated by folds (e.g. the Zagros), whereas others are dominated by thrusts (or faults). With numerical simulations we could demonstrate that the Zagros is likely folded because it has several weak detachment horizons within the sedimentary overburden. In additional papers we looked at this in more detail, and ongoing work involves better understanding the interaction between folds in three dimensions.
Can we constrain the rheology of the lithosphere from surface observables?
Understanding the rheology of the crust and lithosphere on geological timescales remains one of the big unresolved questions in geodynamics. We are coupling Monte Carlo inverse models with geodynamic forward models in order to constrain the present-day rheology of the lithosphere by combining geophysical observations (such as gravity and GPS velocity) of the lithosphere.
How does melt migrate through the lithospheric mantle and crust?
The main interaction between the solid earth and the surface is through melt (and gases) that are erupted in volcanoes. After two centuries of observations, we know quite a bit about the field relationships between batholiths and country-rocks. Fairly sophisticated two-phase flow models have been developed to study the migration of melt towards mid-oceanic ridges. Yet, much less is known about how melt migrates through the continental lithosphere.
Endmember models include diking, diapirism, and stoping. Since field observations can be found to support each of these models, the debate about emplacement mechanisms is unlikely to be resolved anywhere soon. An additional problem is that models about dike propagation include elastic deformation of the host rocks, whereas most diapiric models include viscous deformation.
In order to better understand under which conditions each of these endmember models are most likely to occur, we need a numerical approach which can have all of them in a self-consistent manner. We have recently developed a two-phase flow approach in a visco-elasto-plastic lithosphere that fulfils these requirements and that is currently being used to understand the reletive importance of each of these observed endmember models.
In other projects, we are incorporating diking in a parameterized manner in numerical models of magma emplacement and we develop semi-analytical solutions that help to understand the physics of melt migration.
What is the interaction between tectonic and erosion?
Under which conditions does erosion affect large-scale tectonic processes is question that is being studied by a number of groups world-wide. We have looked at this with two-dimensional coupled mantle/lithosphere/surface processes models in Taiwan, where we could show that it is the high erosion rates that is probably responsible for creating Taiwan in the first place. We are currently coupling a realistic erosion model with a three-dimensional crustal scale folding model to study the interaction of erosion with tectonics in the Zagros mountains.
How does mantle flow interact with lithospheric deformation?
We are also interested in better understanding the coupling between lithospheric deformation and mantle flow. In one project we used numerical models to study the potential causes for net rotation of the lithosphere on Earth. In another, more regional, study we combine numerical models with semi-analytical approaches to understand the physics of eduction (which is the return motion of an oceanic plate after slab breakoff occurred).
How do continents collide?
Continental collision and mountain belt formation are among the most spectacular surface expressions of tectonics on Earth. Numerous models have been developed to study for example the India-Asia collision, but all models are simplified as they are two-dimensional in one way or the other. We have developed the first fully three dimensional models which simultaneously take lithospheric processes, mantle flow and surface topography into account. In the first of these studies, we compared the fully 3D models with a widely used simplification (the thin-viscous sheet approximation, which is a 2D model that assumes that the lithosphere is homogeneous in the vertical direction). In subsequent work we use more sophisticated models to study India-Asia collision, and ongoing studies involve modelling the collision of India-Asia and the formation of the European Alps on a million year timescale with fully three dimensional models.
What causes shear localization?
The underlying physics of shear localization in the brittle and ductile domain remains incompletely understood. We have performed studies in which we addressed localization in the brittle domain (what controls the angle of shear zones in numerical models of brittle crustal deformation?) and in the ductile domain (mainly focussed on shear heating). We have also performed the first fully 3D models of the formation of the San Andreas fault. In ongoing studies we are also looking at he interactions between multiple shear zones, and are looking at the relative importance of grain size evolution and shear heating on shear localization.