We are actively developing software to simulate long-term tectonic and geodynamic process and tools to better link these models with geological and geophysical constraints.


LaMEM (Lithosphere and Mantle Evolution Model) is a highly scalable 3D code that can be used to simulate multiphase thermo-mechanical geodynamic problems which require visco-elasto-plastic rheologies. It is based on a staggered finite difference discretization/marker-and-cell method and has an internal free surface to study the of mantle flow and lithospheric deformation on surface processes. It has various linear and nonlinear solvers and uses the PETSc framework, which ensures that it runs both on our laptops and on large parallel machines (our largest simulation used over 450'000 processors), with either direct or Galerkin multigrid iterative solvers. Thanks to the use of a stable and computational efficient discretization, we can perform high-resolution 3D simulations in approximately the same time as our earlier 2D simulations (few hours to a few days per simulation).

Check out our publication page for some scientific applications of the code.

LaMEM was developed with funding from the European Research Council.
Download it here.



MVEP2 is a 2D Lagrangian thermo-mechanical finite element code, combined with a marker and cell approach to simulate geodynamic processes. It includes thermo-mechanical feedback, phase transitions, erosion and sedimentation, melt extraction and crust formation, chemical evolution of partially molten rocks, and can model both viscous flow in the mantle and localized deformation in the elasto-plastic crust. One of the main advantages of the code is that it is easy to setup models, change the rheology, visualize results and perform production runs. As a result, we have used the code and its predecessor (MILAMIN_VEP) both for teaching and research projects (see our publication list).

The main improvements in MVEP2 compared to earlier code versions is that remeshing is faster and that it is now very simple to add new rheologies (e.g., grain size evolution & feedback with deformation). Despite being written in MATLAB and being based on direct solvers, the code is relatively efficient thanks to the use of fast matrix assembly routines and MUTILS (see It can use both quadrilateral (structured) meshes (with quadratic or linear elements) or unstructured triangular meshes. We are working on adding full two-phase flow models of magma migration to the code.

Download MVEP2 here.


Performing geodynamic simulations requires both efficient/accurate solvers and realistic input models. Particularly in 3D, it can be quite challenging to create such input models from the multitude of available published data. For this purpose, we developed geomIO, which allows you to draw 2D (vertical or horizontal) cross-sections in Inkscape or Illustrator and create a 3D model from that. Changing aspects of the 3D model geometry is now a matter of a few clicks, which makes it much easier to test the effect of the model geometry on the results of geodynamic simulations (an aspect that is often ignored). It works in 2D as well, and can be used with any (geodynamic) software.

Check out for the software and documentation. The geomIO masters (Tobias Baumann & Arthur Bauville) are preparing a paper on it.


Geophysical data give important constraints on the present-day geometry of the lithosphere. Yet, in order to understand how the lithosphere deforms, we also need to have constraints on the rheology of the lithosphere. Unfortunately, experimental data have significant uncertainties. An alternative method to constrain the dynamics of the lithosphere is to couple geodynamic simulations with geophysical data. We have tested this geodynamic inverse modeling approach over the last few years and could demonstrate that it indeed works, using both analytical models and synthetic 2D and 3D simulations. For the inverse models, we use a geometrically constrained Markov-Chain Monte Carlo method, the Neighborhood Algorithm (NA) developed by Malcolm Sambridge (ANU). Yet, as geodynamic forward simulations are relatively expensive and as the speed of a forward simulation depends on the viscosity contrast and can vary significantly, it was necessary to fully rewrite the parallel part of the NA to remove any scalability bottlenecks. The resulting code, NAplus, is thus similar to the original NA method, but is highly scalable (we have tested in on over 16'000 cores), has better (parallel) random number generators, and can be combined either serial or parallel forward models written in Fortran, C or MATLAB. You only need to provide a routine that takes the input parameters and returns a misfit value. The code can be used for any Monte-Carlo based inversion problem.

NAplus was developed with funding from the European Research Council and the code was described in Baumann et al. (2014) and Baumann and Kaus (2015). Download a copy of it here.