Last Summer while I was interning at BGI (Bavarian Research Institute of Experimental Geochemistry and Geophysics), I had the opportunity to attend lectures from Professor Golabek on the subject of numerical modelling of geophysical processes. This covered the theory of finite differences, and included examples of explicit and implicit finite differences relating to geological scenarios such as diffusion of heat from igneous intrusions.

For the final project, I attempted to reproduce some of the results from a paper by Pascal and Olesen, [2009]. This paper involved modelling the heat flow from an asthenospheric diapir that interacts with the base of a cold crustal lithosphere. A model setup was created (see Figure 1), and results were gathered using 2D finite element thermal modelling.

I attempted to recreate this in MATLAB, using the finite differences skills I learnt during lectures. Here are some of my results, in gif form!

This was created using a moderate mesh of nodes and a large timestep size, to observe the large-scale thermal diffusion created by the diapir.

This simulation was created with a fine mesh of nodes and a small timestep, to see the fine features of the heat diffusion. This would have been very computationally expensive, had it not been for my changes to the efficiency in my code (see Figure 5).

I also computed the temperature evolution with depth for the centre of the diapir (see Figure 4).

Here it is possible to observe the predicted increase in shallow geothermal gradient due to the heat diffusion, as well as how the initial sharp jump in temperature is smoothed out over time.

A major part of these calculations was creating efficient scripts. Figure 5 shows an example of these results.

Each time I am creating the same sparse matrix for use in the finite differences code. The first run I use an non-vectorised script, for the second run I vectorise part of it, and for the third run I have completely vectorised the code. Using this, I was able to run simulations at fine grid spacing without dealing with inordinate computation time.

In the future, I plan to try to use my knowledge if finite differences methods to recreate results of deformation of a thin viscous sheet, from England and McKenzie [1982].

References:

Pascal, C., & Olesen, O. (2009). Are the Norwegian mountains compensated by a mantle thermal anomaly at depth?. Tectonophysics, 475(1), 160-168.

England, P., & McKenzie, D. (1982). A thin viscous sheet model for continental deformation. Geophysical Journal International, 70(2), 295-321.