Finite Element method
Applied to




FANTOM is a finite element code for the solution of the Stokes and heat transport equations. It has been purposely designed to address geological flow problems in two and three dimensions. A variety of rheologies has been implemented including linear and non-linear thermally activated creep and brittle (or plastic) frictional models. A cloud of particles is used to track materials in the simulation domain which allows to record the integrated history of deformation; its density is variable and dynamically adapted. The solution to the large system of algebraic equations resulting from the finite element discretization and linearization of the set of coupled partial differential equations to be solved is obtained by means of an efficient (sequential or massively parallel) direct solver. The code can be run in sequential or parallel mode.

Three typical simulations illustratring a different version of the code are presented and briefly discussed: the 2D numerical sandbox, a 2D lithospheric extension experiment, and a 3D crustal scale orogeny experiment.


In order to test the large deformation viscous-plastic behaviour of the code, the extensional numerical sandbox experiment [1] was carried out with the two-dimensional sequential version of the code. The central idea is to compare analogue experiments observations with numerical results, the former being conducted in several laboratories, the latter obtained with different numerical codes.

We focus here on the extensional setup: this experiment consists of a sand layer which overlies a viscous central layer. A thin rigid sheet covering the base of the model extends from the middle to the right lateral wall and is attached to it. The outward mouvement of the right wall leads to an extensional system. The sand is frictional plastic and characterised by ρ=1560kg/m3, c=10Pa, φ=36o, φsw=31o, ε1=0.5, ε2=1.0. The silicone is a viscous material characterised by ρ=965kg/m3 and μ=5x104Pa.s.

[1] The numerical sandbox: comparison of model results for a shortening and an extension experiment S. Buiter and A.Y. Babeyko and S. Ellis and T.V. Gerya and B.J.P. Kaus and A. Kellner and G. Schreurs and Y. Yamada, Analogue and Numerical Modelling of Crustal-Scale Processes. Geological Society, London. Special Publications, p29-64, 253, 2006.


This experiment is nearly identical to the one discussed in [2] in which plane strain thermo-mechanical finite-element model experiments have been used to investigate the effects of frictional plastic strain softening and inherited weakness on the style of lithospheric extension.

The setup, shown hereunder consists of three layers:

The size of the numerical domain is Lx=1200km, Ly=600km and the boundary conditions are as follows: the temperature is set to T=1330o at the base of the model and to 0o at the top. At startup, a constant geotherm T=550o is placed at the base of the crust

The extensional velocity applied to the sides of the crust is vext=0.5cm/yr and a re-entrant velocity field is applied on the rest of the boundary so as to lead to a zero net-flux through the vertical sides of the box. A weak seed is placed in the upper part of the lithosphere and represents the weakest zone of inherited damage, therefore controlling the strain localisation process.

All materials see their density depend on the temperature field with a thermal expansion coefficient α=3.1x 10-5 oC-1.

[2] Roles of lithospheric strain softening and heterogeneity in determining the geometry of rifts and continental margins R.S. Huismans and C. Beaumont, Geological Society, London, Special Publications, 282, 2007.


Let us consider a three-dimensional domain of dimensions Lx x Ly x Lz, filled with the same wet quartz which constituted the crust of the previous experiment. It is submitted to compressive velocity boundary conditions: u=+0.5cm/yr on face x=0 and u=-0.5cm/yr on face x=L_x, while v=0 is imposed on both faces y=0 and y=Ly and w=0 is imposed at the bottom of the box.

A random weak seed is introduced as shown hereunder: cloud points are given a random strain value between 0 and 1.75 ponderated by the function (1-\cos(2π x/Lx)) (1+\cos(2π y/Ly)) (1-\cos(2 π z/Lz))/8. Having set ε1=0.25 and ε2=1.25, this insures that some elements will be in the strain weakening regime, while not being fully strain weakened, i.e. ε1 < ε < ε2. Analogously to the previous experiment, T=0oC is imposed on the free surface while T=550oC is imposed at the bottom.