Welcome!!  Here you will find information about the latest versions of the "LitMod" codes and other software I have developed and/or contributed to.

If you are looking for the global model LithoRef18, please go to the bottom of the page.

LitMod_2D and LitMod_T

LitMod_2D is an interactive software to perform multi-observable modelling of geophysical data for the whole lithospheric structure of the Earth and/or other terrestrial planets (in Windows and Linux). It was developed by J. C. Afonso but there are a few versions available from different groups. Currently, it is maintained mainly by the Group of Dynamics of the Lithosphere at the Institute of Earth Sciences J. Almera in Barcelona. You should get in touch with them if you are interested in using the latest versions of the code.

LitMod_T is the "transient" equivalent of LitMod_2D (i.e. it solves the transient heat conduction equation rather than the steady-state equation). Contact me if you are interested in using this code.

Here are some relevant publications using LitMod_2D:

- Afonso, J.C., Fernandez, M., Ranalli, G., Griffin, W. L., Connolly, J. (2008). Combined geophysical-petrological modelling of the lithospheric-sublithospheric upper mantle: methodology and applications. Geochem. Geophys. Geosyst., 9, Q05008, doi:10.1029/2007GC001834. PDF

- Tunini, L., Jimenez-Munt, I., Fernandez, M., Verges, Villasenor, A. Melchiorre, M.,  Afonso, J.C., Geophysical-petrological model of the crust and upperm mantle in the India-Eurasia collision zone. Tectonics , 35, 1642-1669, doi:10.1002/2016TC004161 PDF

- Carballo, A., Fernandez, M., Jimenez-Munt, I., Torne, M., Verges, J., Melchiorre, M., Pedreira, D.,  Afonso, J.C., Garcia-Castellanos, D., Diaz, J., Villasenor, A., Pulgar, J.A., Quintana, L. (2015), From the Norht-Iberian Margin to the Alboran Basin: a lithosphere geo-transect crossing the Iberian Plate. Tectonophysics. , 663, 399-418. PDF

- Pedreira, D., Afonso, J.C. et al., (2015), Geophysical-petrological modeling of the lithosphere beneath the Cantabrian Mountains and North-Iberian margin: geodynamic implications. Lithos, 230, 46-68. PDF

- Fernandez, M., Afonso, J.C., Ranalli, G. (2010). The deep lithospheric structure of the Namibian volcanic margin Tectonophysics, 481, 68-81. PDF

- Alasonati-Tasarova, Z., Afonso, J.C., Bielik, M., Gotze, H-J. Hak, J. (2009). The lithospheric structure of the Western Carpathian-Pannonian Basin region based on the CELEBRATION 2000 seismic experiment and gravity modelling, Tectonophysics, 475, 454-469. PDF 

Thermochemical model of western Tibet obtained with LitMod_2D by fitting surface heat flow, bouguer anomalies, geoid anomalies, elevation and seismic anomalies. f) anf g) show a comparison of predicted and observed P-wave anomalies. Results from Tunini et al. (2015). 


LitMod_3D is the big brother of LitMod_2D and it's used in industry and academia. It is also fully interactive and it has been developed in collaboration with J. Fullea. The currently supported version run in Linux only. The latest version of the code and the user guide can be downloaded from the button below. LitMod_3D is distributed as a compressed tar file that contains all the required executables, sources, and shell script files. All sources are freely available. If you are interested in LitMod_3D, here are some relevant publications you may want to check:

Fullea, J., Afonso, J.C., Connolly, J.A.D., Fernandez, M., Garcia-Castellanos, D., Zeyen, H. (2009). LitMod3D: an interactive 3D software to model the thermal, compositional, density, rheological, and seismological structure of the lithosphere and sublithospheric upper mantle. Geochem. Geophys. Geosyst, 10, Q08019, doi:10.1029/2009GC002391. PDF

- Fullea, J., Fernandez, M., Afonso, J.C., Verges, J., Zeyen, H. (2010). The structure and evolution of the lithosphere-asthenosphere boundary beneath the Atlanitc-Mediterranean Transition Region. Lithos, 120, 74-95. PDF 

- Alasonati Tašárová, Z., J. Fullea, M. Bielik, and P. Środa (2016), Lithospheric structure of Central Europe: Puzzle pieces from Pannonian Basin to Trans-Europe an Suture Zone resolved by geophysical-petrological modeling, Tectonics, 35, doi:10.1002/2015TC003935. PDF 

- Gradmann, S. Ebbing, J., Fullea, J. (2013) Integrated geophysical modelling of a lateral transition zone in the lithospheric mantle under Norway and Sweden, Geophys. J. Int., 194,1358–1373,

A 3D rendering of the thermal structure of the Atlantic-Mediterranean region from the study of Fullea et al. (2010) using LitMod_3D.


The code is fully functional, but under constant testing and development by J.C. Afonso's group at Macquarie University. Besides academics, LitMod_4INV is being used (and funded) by a number of industry partners and government agencies. Send me an email if you are interested in using it. Some relevant publications are:

- Jones, A.G.,  Afonso, J.C., Fullea, J., (2017), Geochemical and geophysical constrains on the dynamic topography of the Southern African Plateau., Geochem. Geophys. Geosyst., 18, 3556–3575, doi:10.1002/2017GC006908. PDF

- Qashqai, T.M.,  Afonso, J.C., Yang. Y. (2016), The crustal structure of the Arizona Transition Zone and southern Colorado Plateau from multi-observable probabilistic inversion..Geophys.Geochem.Geosys. , 17, doi:10.1002/2016GC006463. PDF

- Afonso, J.C., Rawlinson N., Yang, Y., Schutt, D., Jones, A.G., Fullea, J., Griffin, W.L. (2016), 3D multi-observable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle III: Thermochemical Tomography in the Western-Central US.J. Geophys. Res., 121, doi:10.1002/2016JB013049 PDF

This is the "big daddy" of the LitMod suite. It is a fully nonlinear probabilistic inversion code(s) in spherical coordinates for the compositional and thermal structure of the lithosphere and upper mantle. It is capable of simultaneously inverting

gravity gradients, gravity anomalies, geoid height, surface heat flow, magnetotelluric data, receiver functions, surface-wave data, absolute elevation (including both static and dynamic effects) and body-wave data, together with 

petrological information. It requires a decent cluster to run, although it is faster/cheaper when not all the above data sets are inverted together. 

Some results obtained with LitMod_4INV from Afonso et al. (2016) JGR. 

- Guo, Z.,  Afonso, J.C., Qashqai, M., Yang, Y., Chen, J. (2016), Thermochemical structure of the North China Craton from multi-observable probabilistic inversion: extent and causes of cratonic lithosphere modification. Gondwana Res. ,37, 252-265. PDF

- Shan, B.,  Afonso, J.C., Yang, Y., Grose, C.J., Zheng, Y., Xiong, X. (2014), The compositional and thermal structure of the lithosphere and upper mantle beneath South China: Results from multi-observable probabilistic inversion. J. Geophys. Res., 119, 8417-8441, doi:10.1002/2014JB011412 PDF

​- Afonso, J.C., Fullea J., Yang, Y., Connolly, J.A.D., Jones, A.G. (2013), 3D multi-observable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle II: General methodology and resolution analysis. J. Geophys. Res., 

doi:10.1002/jgrb.50123. PDF 

- Afonso, J.C., Fullea J., Griffin, W.L., Yang, Y., Jones, A.G., Connolly, J.A.D., O'Reilly, S.Y. (2013), 3D multi-observable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle I: a priori information and geophysical observables.   J. Geophys. Res., doi:10.1002/jgrb.50124. PDF

Efficient and general algorithm for implementing thermodynamic information into geodynamic and geophysical studies

Our approach is based on tensor rank decomposition methods, which transform the original multi-dimensional discrete information into a separated representation that contains significantly fewer terms, thus drastically reducing the amount of information to be stored in memory during a numerical simulation or geophysical inversion. Accordingly, the amount and resolution of the
thermodynamic information that can be used in a simulation or inversion increases substantially. Any thermodynamic software or database can be used to obtain the primary thermodynamic information.


Please refer to: Afonso, J.C., Zlotnik, S., Diez, P. (2015) “An efficient and general approach for implementing thermodynamic phase-equilibria information in geophysical and geodynamic studies”, Geochem, Geophys., Geosys., 16, doi:10.1002/2015GC006031

Density fields obtained from the energy minimization problem for water contents from 0 to 4% (from left to right) are shown in the first row. Same fields provided by the decomposition procedure with 100 terms are shown in the second row. Their absolute difference (full solution - separated solution) is shown in the third row. 

ADR FEM code

This algorithm/code was presented in the paper: 

- Oliveira, B., Afonso, J.C., Zlotnik, S. (2016), A Lagrangian-Eulerian finite element algorithm for advection-diffusion-reaction problems with phase change, Comput. Methods Appl. Mech. Engrg., 300, 375-401.

It is a matlab code that implements the algorithm described in the paper above, which is a very powerful and accurate algorithm to model a large number of problems that involve ADR with phase changes (e.g. metal casting, groundwater transport problems, magmatic systems, etc). The algorithm combines a Lagrangian formulation for the advection + reaction problem with the Eularian-based heat source method for the diffusion + phase change problem. 

This algorithm is the basis for solving the energy equation of the multi-phase multi-component reactive transport problem tackled in Oliveira et al. (2017) and Oliveira et. al. (2018)... see below.

Contour plots of the three-body rotation problem. Advection (3 revolutions) and advection+diffusion (1 revolution) problems (by rows) for different interpolation schemes (by columns). The reference solutions are shown in the rightmost column.

Code for multi-phase multi-component reactive transport 

This multi-algorithm method and code was presented in the paper: 

- Oliveira, B., Afonso, J.C., Zlotnik, S., Diez, P.  (2018), Numerical Modelling of Multi-Phase Multi-Component Reactive Transport in the Earth`s interior, Geophys. J. Int., 212, 345-388.

It is a (so far matlab) code implementing the general strategies described in the paper above. The approach is truly multiphase in the sense that each dynamic phase is explicitly modelled with an individual set of mass, momentum, energy and chemical mass balance equations coupled via interfacial interaction terms. It is also truly multicomponent in the sense that the compositions of the system and its constituent phases are expressed by a full set of fundamental chemical components (e.g. SiO2, Al2O3, MgO, trace elements and isotopes.) rather than proxies. These chemical components evolve, react with and partition into different phases according to an internally consistent thermodynamic model. We combine concepts from Ensemble Averaging and Classical Irreversible Thermodynamics to obtain sets of macroscopic balance equations that describe the evolution of systems governed by multiphase multicomponent reactive transport (MPMCRT). Equilibrium mineral assemblages, their compositions and physical properties, and closure relations for the balance equations are obtained via a ‘dynamic’ Gibbs free-energy minimization procedure (i.e. minimizations are performed on-the-fly as needed by the simulation).

We are currently using this code to simulate disequilibrium processes during the generation and segregation of melt in the mantle, metasomatism and isotopic fractionation among others. This code is constantly being benchmarked and developed by B. Oliveira and J.C. Afonso at Macquarie University. Please contact us if you are interested in using it for doing research... 

On the left: a schematic of the philosophy of the approach... where different "gears" are coupled to make the entire machine/code works

On the right: an example from Oliveira et al. (2018) where a thermal plume develops a second "melt" plume on the top due to the melt percolating through the solid matrix at a faster speed than the actual thermal plume.


This is a joint inversion code for inverting Vs and anisotropy data (Rayleigh and Love phase and/or group velocities, ellipticity (Z/H ratio) and receiver functions) with a Markov Chain Monte Carlo method. This code is a collaboration with Z. Guo and will be released soon

LithoRef18: A global reference model of the lithosphere and upper mantle from joint inversion and analysis of multiple data sets

This model was presented in the paper: 

- Afonso, J.C., Salajeghegh, F., Szwillus, W., Ebbing, J. Gaina, C. (2019), A global reference model of the lithosphere and upper mantle from joint inversion and analysis of multiple datasets. Geophys. J. Int., ggz094,

LithoRef18 is a new global model for the Earth’s lithosphere and upper mantle obtained through a formal joint inversion of 3D gravity anomalies, geoid height, satellite-derived gravity gradients and absolute elevation complemented with seismic, thermal and petrological prior information . The model includes crustal thickness, average crustal density, lithospheric thickness, depth-dependent density of the lithospheric mantle, lithospheric geotherms, and average density of the sublithospheric mantle down to 410 km depth with a surface discretization of 2°×2°. It is simple, but good...


We also provide some matlab codes to compute far-field effects and interrogate LithoRef18 (see below) 

"LithoRef model" contains the global model as a text file

"LithoRef code" is a matlab code to compute far-field effects.

This code (cube.m) will allow you to extract profiles along any direction at any location.

Some teaching/toy codes 
Free-energy minimization using the simplex/Perple_x approach

This simple matlab code implements the main idea of doing Gibbs free energy minimization using linear programming, an approach used by the famous code PerPle_X (e.g. Connolly, 2005). The example allows the user to use different equations of state and mineral phases. The default thermodynamic data is that of Fabrychnaya (1995) for the "olivine" system.

Download code
Parallel Tempering, Multiple-try Metropolis sampling

This matlab code illustrates the implementation of the Parallel Tempering (PT) and Mutliple-try Metropolis (MTM) algorithms in the context of MCMC. The user can use PT or MTM separately or together.  

Delayed Rejection Adaptive Metropolis (DRAM) in Fortran

At some point I decided to write a simple teaching code in Fortran (but it's very useful for real inversions too!) of the DRAM method of Haario et al. (2006) based on the matlab implementation by Marko Laine. 

Download code
Download code
Approximate Bayesian Computation (ABC) 

This is a matlab code that implements ABC to sample from a target distribution. It was written by T. Connell during his MRes project and it is based on the DRAM matlab code by M. Laine.

Download code
Gravity modelling and inversion in 2D

A very simple matlab code to illustrate the basics of gravity modelling and inversion using the original algorithm of Talwani et al. (1977). The inversion method is the so-called Compact Gravity Inversion (REF?).

Download code
More coming soon...

Make sure you visit this website in a desktop as the mobile version does not contain all software.

I also recommend visiting the websites of the following books, which contain lots of simple matlab examples for Earth sciences:

- Garya, T. (2012) Introduction to Numerical Geodynamic Modelling, CUP.

- Simpson, G. (2016) Practical Finite Element Modelling in earth Sciences using Matlab, Wiley.