Mantle Convection and Earth’s Thermal Evolution

My PhD research, supervised by W.R. Peltier at the University of Toronto, involved simulating convection in Earth’s mantle and, in particular, investigating the effects of the endothermic phase transition at 660 km depth in the mantle and its efficacy for producing layered or partially layered convection. I then became interested in in the effect of partial layering on Earth’s surface heat flow which led to an interest in understanding Earth’s thermal history.

I was also briefly a post-doc at the York University in Toronto working with Prof. Garry Jarvis at York University, where we looked at the magnitude of stresses imparted on continents caused by mantle flow.

My former PhD student, Simona Costin, developed a parameterized thermal and magnetic model for Earth’s core that we coupled to the mantle convection model to create a thermal evolution model for Earth. A copy of Simona's thesis can be found here .

fig 1

fig2

fig 3

fig 4

Here are some old movies of mantle convection in the axisymmetric spherical shell (fig 1) that we used to calculate thermal evolutions. The first has no phase transitions while in the next movie weak phase transitions (fig 2)are acting and in this one earth-like phase transitions (fig 3). There is also a movie with a very strong phase transition effect (fig 4). Note the mantle "avalanches" that are taking place in the last two movies, particularly the last one.

Another former PhD student, Gunjan Sinha, used mantle convection models to investigate why the mean mantle temperature gradient might be less than adiabatic and to investigate the combined effects of continents and phase transitions on the planform of mantle convection. A copy of Gunjan's thesis can be obtained here .

As an undergraduate, Christine Shiels used mantle convection models to study the creation of long wavelength flow in models with stiff mantle plates and a low viscosity asthenosphere. In the image below, the results of a simulation with a surface plate and a low viscosity asthenosphere is shown in an aspect ratio 5 box. Temperatures show color while arrows show velocity. It can be seen that the velocity in the asthenosphere exceeds that in the plate.

I have recently been using Comsol Multiphysics to carry out some 3D mantle convection simulations in spherical shells. The image below shows a temperature isosurface from one of these simulations.

Post-doctoral fellow Sean Trim worked with me in 2018 and 2019. A paper was published, in collaboration with Julian Lowman, on improving mass conservation in thermo-chemical convection using tracer repositioning.

Papers on Mantle Convection and Earth's Thermal Evolution

Flows in Porous Media

In 2002, I spent some time at Cambridge University working at the Institute for Theoretical Geophysics with Herbert Huppert and Grae Worster on modeling mushy layers. Mushy layers are porous media where solidification and melting can take place. A great deal of interesting dynamical behaviour can occur because of the change in density that can occur because of solidification and melting and because of the effects on permeability.

A porous layer melting due to heating from above. The colours show the porosity and the top layer becomes completely molten.

An interesting issue in porous media is that a dissolved solute and heat both advect and diffuse at very different rates. This can lead to interesting behaviour if a liquid with a different solute concentration and temperature is injected into a porous layer. Julia Milne, a summer student, and I carried out a number of simulations to better understand this problem. Take a look at some simulationsthat Julia made in the summer of 2003.

In a mushy layer, the composition and thermal fields are coupled by the liquidus condition and so the effective advection and diffusion rates will be affected. In one research project, I calculated these, as well as transport induced melting rates in a mushy layer.

In 2007, I took a sabbatical leave to Brown University in Providence Rhode Island where I became interested in issues related to compaction in porous layers. Theory and experiments have shown that if the viscosity of the solid matrix decreases with porosity, that if the matrix is strained, the matrix will spontaneously separate into regions of high and low porosity. When buoyancy is also present, porosity waves can also occur. I created a numerical model to investigate the combined effects of these two phenomena and showed that porosity localization still occurs in bands, but there may be more than one set of porosity bands.

fig 4

fig 5

fig 6
(fig 4) a movie of a compacting porous layer undergoing pure shear with horizontal compression and buoyancy where the viscosity depends only on porosity. (fig 5) is a movie of the same thing with vertical compression. Note the change in the oscillations because of the different orientation of the bands relative to gravity. (fig 6) is a model with vertical compression and gravity where the viscosity of the matrix depends on the strain-rate of the matrix as well. Notice how there are now two sets of bands at roughly 20 degrees to vertical.

The deformation of partially molten samples also results in anisotropic arrangement of melt pockets at the grain scale. It has been proposed that this arrangement may cause the matrix to have an anisotropic viscosity. I have carried out simulations of melt bands with anisotropic viscosity and these can lead to low angle bands.

Current graduate student, David Gebhardt, is working on modeling melt bands where the corner flow of a mid-ocean ridge is used as the background flow. The aim is to evaluate how effective melt bands might be for channeling melt towards mid-ocean ridges.

For his MSc research, Michael Bird imported synchrotron X-ray tomography derived images into Comsol and created meshes of the pore-space. He then solved the Navier-Stokes equations for fluid flow and Laplace's equation for electrical conduction in order to calculate the permeability and the electrical formation factor. An image of the streamlines, color coded with the magnitude of the velocity is shown in the image below.

Here is a link to a lecture that I gave as part of a Melt in Earth's Mantle workshop at the Isaac Newton Institute at Cambridge University.

Papers on Flows in Porous Media

Tektites

Splashform tektites are glassy rocks that were created from the splash of molten rock that occurred after certain large Earth impacts.

Splashform tektites are found in a number of interesting shapes including near spheres, disks and dumbbells. The shapes of tektites are thought to arise from an interaction of centrifugal and surface tension forces on rotating molten drops. My colleague, Mel Stauffer , has compiled data on the dimensions of thousands of tektites from southeast Asia while my colleague, Ray Spiteri , and I have created a model of rotating liquid drops with surface tension. We have shown that the time evolution of this model agrees well with the shape distribution of real splash-form tektites.

A movie from one of our simulations of a drop going from a near sphere to a flattened ellipsoid to a dumbbell can be seen here.

Claire Samson's group at Carleton University has imaged tektites using a laser scanning camera. These images allow us to calculate the volume and moment of inertia of tektites which allows us to calculate their density and constrain their rotation rates. We have also used the resulting image in Comsol to flow air over the tektites in order to calculate their drag coefficients which help us to determine whether aerodynamic effects are important in shaping tektites. The image below shows the streamlines and magnitude of velocity (colors) from a simulation of a tektite moving to the left through air.

I have also provided numerical model results for comparison laboratory-produced synthetic tektites. Here is a link to a recent physicsworld article regarding artificial tektites.

Papers on Splashform Tektites

Applied Geophysical Techniques

I teach a course in gravity, magnetics, electromagnetics and radiation methods. As part of the course, I have created virtual labs, using Comsol Multiphysics, that go along with the real physical labs for the course. If you are interested in using any of the labs, please contact me. I have published a couple of papers regarding simulating various techniques of applied geophysics using Comsol. MSc student, Michael Zhang, worked on simulating VLF waves and investigated the effects of topography and lakes.

I have also spent some time delving into the theory of the resistivity method. As a result of an enquiry from a former student, we have looked into the effects of an leakage current or of an unintentional additional potential electrode.

It has long been held that there is no simple, general, formula to estimate the depth of investigation for the resistivity method that takes into account the positions of all four electrodes. I have found that if you consider the mean of the sensitivity function of the resistivity method for a homogeneous infinite half space, that it does have a simple formula which is shown below. Here, kg is the geometrical factor and rAN is the distance from the A electrode to the N electrode. The other distances are defined similarly. This estimate works for all array types except a pole-pole.

Papers on Applied Geophysics

Click here for Comsol files in support of our 2015 Computers and Geosciences paper "Forward Modeling of Geophysical Electromagnetic Methods Using Comsol"