Grace Shephard / Tobias Meier


Find out more about the blog team here.

Rheological Laws: Atoms on the Move

Rheological Laws: Atoms on the Move

The Geodynamics 101 series serves to showcase the diversity of research topics and methods in the geodynamics community in an understandable manner. We welcome all researchers – PhD students to Professors – to introduce their area of expertise in a lighthearted, entertaining manner and touch upon some of the outstanding questions and problems related to their fields. For our first ‘101’ for 2018, we have an entry by postdoctoral researcher Elvira Mulyukova from Yale University about rheology and deformation occurring on atomic scales … it’s a fun and informative read indeed! Do you want to talk about your research? Contact us!

Elvira Mulyukova, Yale University

Most of us have an intuitive understanding that different materials resist being moved, or deformed, to different degrees. Splashing around in the mud is more energy-consuming (and fun, but never mind that) than in the water, and splashing around in the block of concrete is energy-intensive bordering on deadly. What are the physical reasons for these differences?

For Earth materials (rocks), the answer lies in the restless nature of their atoms: the little buggers constantly try to sneak out of their crystal lattice sites and relocate. Some are more successful at it than others, making those materials more easily deformable. A lattice site is really just other atoms surrounding the one that is trying to escape. You see, atoms are like a bunch of introverts: each is trying to escape from its neighbours, but doesn’t want to get near them. The ones that do escape have to overcome a temporary discomfort (or an increase in their potential energy, for those physically inclined) of getting close to their neighbours. This requires energy. The closer the neighbours – the more energy it takes to get past them. When you exert force on a material, you force some of the neighbours to be further away from our potential atomic fugitive, making it more likely for the atom to sneak in the direction of those neighbours. The fun part (well, fun for nerds like me) is that it doesn’t happen to just one atom, but to a whole bunch of them, wherever the stress field induced by the applied force is felt. A bunch of atoms escaping in some preferential direction is what we observe as material deformation. The more energy you need to supply to induce the mass migration of atoms – the stronger the material. But it’s really a question of how much energy the atom has to begin with, and how much energy is overall needed to barge through its detested neighbuors. For example, when you crank up the temperature, atoms wiggle more energetically and don’t need as much energy supplied from external forcing in order to escape – thus the material gets weaker.

“A lattice site is really just other atoms surrounding the one that is trying to escape. You see, atoms are like a bunch of introverts: each is trying to escape from its neighbours, but doesn’t want to get near them.” cartoon by Elvira Mulyukova

One more thing. Where are the atoms escaping to? Well, there happen to be sanctuaries within the crystal lattice – namely, crystalline defects such as vacancies (aka point defects, where an atom is missing from the otherwise ordered lattice), dislocations (where a whole row of atoms is missing), grain boundaries (where one crystal lattice borders another, which is tilted relative to it) and other crystalline imperfections. These regions are sanctuaries because the lattice is more disordered there, which allows for larger distances in-between the neighbors. When occupying a regular lattice site – the atom is sort of trapped by the crystalline order. Think of the lattice as an oppressive regime, and the crystalline defects as liberal countries that are welcoming refugees. I don’t know, is this not the place for political metaphors? *…whistling and looking away…*

Ok, enough anthropomorphisms, let’s get to the physics. If this is the last sentence you’ll read in this blog entry, let it be this: rocks are made up of atoms that are arranged into crystal lattices (i.e., ordered rows and columns of atoms), which are further organized into crystal grains (adjacent crystals tilted relative to each other); applying force to a material encourages atoms to move in a preferential direction of the largest atomic spacing, as determined by the direction of the applied force; the ability of the lattice sites to keep their atoms in place (call it a potential energy barrier) determines how easily a material deforms. Ok, so it was more like three sentences, but now you know why we need to get to the atomic intricacies of the matter to understand materials macroscopic behaviour.

Alright, so we’re applying a force (or stress, which is simply force per area) to a material and watch it deform (a zen-inducing activity in its own right). We say that a material behaves like a fluid when its response to the applied stress (and not just any stress, but differential stress) is to acquire a strain rate (i.e., to progressively shorten or elongate in one direction or the other at some rate). On geological time scales, rocks behave like fluids, and their continuous deformation (mass migration of atoms within a crystal lattice) under stress is called creep.

The resistance to deformation is termed viscosity (let’s denote it µ), which basically tells you how much strain rate (˙e) you get for a given applied differential stress (τ). Buckle up, here comes the math. For a given dimension (say x, and for the record – I’ll only be dealing with one dimension here to keep the math symbols simple, but bear in mind that µ, ˙e and τ are all tensors, so you’d normally either have a separate set of equations for each dimension, or some cleverly indexed symbols in a single set of equations), you have:

So if I’m holding a chunk of peridotite with a viscosity of 1020 Pa s (that’s units of Pascal-seconds, and that’s a typical upper mantle viscosity) and squeezing it in horizontal direction with a stress of 108 Pa (typical tectonic stress), it’ll shorten at a rate of 5 · 10-13 s-1 (typical tectonic rates). A lower viscosity would give me a higher strain rate, or, equivalently, with a lower viscosity I could obtain the same strain rate by applying a smaller stress. If at this point you’re not thinking “oh, cool, so what determines the viscosity then?,” I failed massively at motivating the subject of this blog entry. So I’m just gonna go ahead and assume that you are thinking that. Right, so what controls the viscosity? We already mentioned temperature (let’s call it T), and this one is a beast of an effect. Viscosity depends on temperature exponentially, which is another way of saying that viscosity depends on temperature hellavulot. To throw more math at you, here is what this dependence looks like:

where R = 8.3144598 J/K/mol (that’s Joule per Kelvin per mol) is the gas constant and E is the activation energy. Activation energy is the amount of energy that an atom needs to have in order to even start thinking about escaping from its lattice site, which of course depends on the potential energy barrier set up by its neighbours. Let’s say your activation energy is E = 530·103 J mol-1 . If I raised your temperature from 900 to 1000 K (that’s Kelvin, and those are typical mid-lithospheric temperatures), your viscosity would drop by a factor of ∼ 1000. That’s a three orders of magnitude drop.

Like I said, helluvalot. If instead you had a lower activation energy, say E = 300 · 103 J mol-1 , the same temperature experiment would bring your viscosity down by a factor of ∼ 50, which is less dramatic, but still significant. It’s like running through peanut butter versus running through chocolate syrup (running through peanut butter is a little harder… I clearly need to work on my intuition-enhancing examples). Notice, however, that while the temperature dependence is stronger for materials with higher activation energies, it is more energy-consuming to get the creep going in those materials in the first place, since atoms have to overcome higher energy barriers. There’s more to the story. Viscosity also depends on pressure (call it P), which has a say in both the energy barrier the atoms have to overcome in order to escape their neighbours, as well as how many lattice defects (called sanctuaries earlier) the atoms have available to escape to. The higher the pressure, the higher the energy barrier and the fewer lattice sanctuaries to resort to, thus the higher the viscosity. Throwing in the pressure effect, viscosity goes as:

The exact dependence of viscosity on pressure is determined by V – the activation volume.

Alright, we’re finally getting to my favourite part – the atoms’ choice of sanctuary sites. If the atomic mass migration happens mainly via point defects, i.e., by atoms hopping from one single lattice vacancy to another, the deformation regime is called diffusion creep. As atoms hop away, vacancies accumulate in regions of compressive stress, and fewer vacancies remain in regions of tensional stress. Such redistribution of vacancies can come about by atoms migrating through the bulk of a crystal (i.e., the interior of a grain, which is really just a crystal that is tilted relative to its surrounding crystals), or atoms migrating along the boundary of a crystal (i.e., a grain boundary). In both cases, the rate at which atoms and vacancies get redistributed depends on grain size (let’s denote it r). The larger the grains – the more distance an atom has to cover to get from the part of the grain that is being compressed to the part that is under tension. More math is due. Here is what the viscosity of a material deforming by diffusion creep looks like:

Exponent m depends on whether the atoms are barging through the bulk of a grain (m = 2), or along the grain boundaries (m = 3). What’s that new symbol B in the denominator, you ask? That’s creep compliance (in this case – diffusion creep compliance), and you two have already met, sort of. Creep compliance specifies how a given creep mechanism depends on pressure and temperature:

For diffusion creep of upper mantle rocks, I typically use m = 3, B0 ∼ 13 µmm MPa-1 s-1 (which is just a material-specific prefactor), Ediff = 300 · 103 J mol-1 and Vdiff = 5 cm3 mol-1 from Karato and Wu (1993). Sometimes I go bananas and set Vdiff = 0 cm3 mol-1 , blatantly ignoring pressure dependence of viscosity, which is ok as long as I’m looking at relatively modest depth-ranges, like a few tens of kilometers.

At sufficiently high stresses, a whole row of atoms can become mobilized and move through the crystal, instead of the meagre one-by-one atomic hopping between the vacancies. This mode of deformation is called dislocation creep. Dislocations are really just a larger scale glitch in the structure of atoms (compared to vacancies). They are linear lattice defects, where a whole row of atoms can be out of order, displaced or missing. It requires more energy to displace a dislocation, because you are displacing more than one atom, but once it’s on the move, it accommodates strain much more efficiently than in the each-atom-for-itself diffusion creep regime. As the material creeps, dislocations get born (nucleated), get displaced and get dead (annihilated). Dislocations don’t care about grain size. What they do care about is stress. Stress determines the rate at which dislocations appear, move and disappear. I know you saw it coming, more math, here is the dislocation creep viscosity:

Exponent n dictates the stress dependence of viscosity. Stress dependence of dislocation creep viscosity is a real pain, making the whole thing nonlinear and difficult to use in a geodynamical model. Not impossible, but rage-inducingly difficult. Say you’re trying to increase the strain rate by some amount, so you increase the stress, but then the viscosity drops, and suddenly you have a monster of a strain rate you never asked for. Ok, maybe it’s not quite this bad, but it’s not as good as if the viscosity just stayed constant. You wouldn’t be able to have strain localization, form tectonic plate boundaries and develop life on Earth then, but maaaan would you be cracking geodynamic problems like they were peanuts! I’m derailing. Just like all the other creeps, dislocation creep has its own compliance, A, that governs its dependence on pressure and temperature:

For dislocation creep of upper mantle rocks, I typically use n = 3, A0 = 1.1 · 105 MPa-n s -1 (which is just a material-specific prefactor), Edisl = 530 · 103 J mol -1 and Vdiff = 20 cm3 mol-1 , similar to Karato and Wu (1993). Just like for diffusion creep, I sometimes just set Vdisl = 0 cm3 mol-1 to keep things simple.

We’re almost done. Allow me one last remark. A rock has an insane amount of atoms, crystal grains and defects, all subject to local and far-field conditions (stress, temperature, pressure, deformation history, etc). A typical rock is therefore heterogeneous on the atomic (nano), granular (micro) and outcrop (meter) scales. Thus, within one and the same rock, deformation will likely be accommodated by more than just one mechanism. With that in mind, and sticking to just two deformation mechanisms described above, we can mix it all together to get:

This is known as composite rheology, where we assumed that the strain rates accommodated by diffusion and dislocation creep can be simply summed up, like so:

Alright. If you got down to here, I salute you! Next time you’re squeezing a peridotite, or splashing in the mud, or running through peanut butter – give a shout out to those little atoms that enable you to do such madness. And if you want to get to the physics of it all, you can find some good introductory texts in Karato (2008); Turcotte and Schubert (2002).


S. Karato. Deformation of earth materials: an introduction to the rheology of solid earth. Cambridge Univ Pr, 2008. 

S. Karato and P. Wu. Rheology of the upper mantle: a synthesis. Science, 260(5109):771–778, 1993. 

D.L. Turcotte and G. Schubert. Geodynamics. Cambridge Univ Pr, 2002.

From hot to cold – 7 peculiar planets around the star TRAPPIST-1

From hot to cold – 7 peculiar planets around the star TRAPPIST-1

Apart from Earth, there are a lot of Peculiar Planets out there! Every 8 weeks, give or take, we look at a planetary body or system worthy of our geodynamic attention. When the discovery of additional Earth-sized planets within the TRAPPIST-1 system was revealed last year, bringing the total to 7 planets, it captured the minds of audiences far and wide. This week, two of the authors from a 2017 Nature Astronomy study on the TRAPPIST-1 planets, Lena Noack from the Department of Earth Sciences at the Free University of Berlin and Kristina Kislyakova from the Department of Astrophysics at the University of Vienna, explain more about this fascinating system. 

Blog authors Lena Noack and Kristina Kislyakova

For Earth scientists, it often seems like a huge endeavour to talk about the geodynamics and other interior processes of the other planets in our Solar System like Mars or Venus. But what about exoplanets? It’s very daring! We have almost no information about the thousands of planets that have been discovered so far in other places of our galaxy. These planets orbit other stars, of which some are quite similar to our Sun whereas other stars behave very differently. But how much do we actually know about planets around these stars?

Exoplanet hunting missions like Kepler have shown that the majority of exoplanets are actually small-mass planets – not huge gas giants like Jupiter – and are often smaller than Neptune, with some being even smaller than Earth. We have a pretty good idea of what some of these planets could look like, for example we know their mass, their radius, we might even have some spectral information on their atmospheres, we know how much energy they get from their star, and we might even know something about the star’s composition. This information hints at the composition of the planetary disk from which planets are made, and how much radioactive heating they may experience during their later evolution. Putting all these pieces together gives us several clues on how the planets may have evolved over time, and is comparable to the wealth of information we had of our neighbouring planets before the age of space exploration.

However, in contrast to our Solar System, we cannot (at least not with our technological standard of today) travel to these planets. The only way we can learn more about exoplanets is if we combine geophysical, thermodynamical and astrophysical models – derived and tested for Earth and the Solar System – and apply them to exoplanet systems.


Artist’s impression of TRAPPIST-1e, ©NASA

One exoplanet system that is quite intriguing is the TRAPPIST-1 system, which has been observed by several different space and ground-based telescopes including TRAPPIST (short for TRAnsiting Planets and PlanetesImals Small Telescope, or otherwise known as a European monastery-brewed beer) and the Spitzer Space Telescope.

The system contains at least 7 small, densely-packed planets around an 8 Gyr old M dwarf. All planets have masses and radii close to Earth – from TRAPPIST-1c and -1h, which are both ¾ the radius of Earth, to TRAPPIST-1g, which is 13% larger than Earth. For comparison, Venus, our sister planet, has a radius 5% smaller than Earth, and Mars, our small brother planet, is only half the size of Earth. And the greatest news: TRAPPIST-1 is actually in our direct neighbourhood, only 39 light years away. This is literally around the corner! For comparison, our closest neighbour planet outside the Solar System is Proxima Centauri b with a distance of 4.2 light years. Its star belongs to a system of three stars, the most well-known of which is Alpha Centauri, the closest star outside the Solar System. Some day, it may actually be in our reach to travel to both the Centauri system as well as TRAPPIST-1. So we should learn now as much about these planets as possible.

What makes the TRAPPIST-1 planets truly peculiar are their tight orbits around the star – the closest planet orbits at a distance of 0.0011 AU – so only 0.1 percent of Earth’s orbit. Even the furthest planet in the system discovered so far– TRAPPIST-1h – has an orbit of only 0.0063 AU. In our Solar System, Mercury, the closest-in planet, orbits at a distance of 0.39 AU. Does this mean that the planets are boiling up due to their close orbit? Not necessarily, since TRAPPIST-1 is a very dim red M dwarf, which emits much less light than the Sun in our system. If we would place Earth around this red M dwarf star, it would actually need to orbit at a distance of about 0.0022 AU to receive the same incident flux from the star. Actually, if we look at the possible distances from the star, where (depending on the atmosphere greenhouse effect) liquid water at the surface could theoretically exist for a somewhat Earth-like atmosphere (that is, composed of gases such as CO2 and N2), TRAPPIST-1d, -1e, -1f and -1g could potentially contain liquid water at the surface and would thus be habitable places where Earth-like life could, in principle, form. Of course, for that to occur several other factors have to be just right, as well. This zone, where liquid water at the surface could exist, is called the Habitable Zone or Temperate Zone, and is indicated in green in the illustration of the TRAPPIST-1 system compared to the inner Solar System below.

TRAPPIST-1 system compared to the inner Solar System below showing the green region of a Habitable Zone. © Caltech/NASA

So, should we already book our trip to TRAPPIST-1? Well, there are other factors that may endanger the possible habitability of these otherwise fascinating planets. First of all, TRAPPIST-1 is really different from the Sun. Although it is much dimmer and redder, it still emits almost the same amount of harsh X-ray and extreme ultra violet radiation as our Sun, and in addition, produces powerful flares. For the TRAPPIST-1 planets, which are so close to their star, it means that their atmospheres are exposed to much higher levels of short wavelength radiation, which is known to lead to very strong atmospheric escape. A nitrogen-dominated atmosphere, like the one Earth has, would likely not be stable on the TRAPPIST-1 planets in the habitable zone due to exposure to this short wavelength radiation for gigayears, so carbon dioxide Venus-like atmospheres are more probable. Besides that, stellar wind of TRAPPIST-1 may be very dense at planetary orbits, powering strong non-thermal escape from planetary atmospheres and leading to further erosion of the atmosphere.

Another interesting feature of the M dwarfs, especially such low-mass ones as TRAPPIST-1, is their extremely slow evolution. On the one hand, this means very long main-sequence life times of such stars, with stable radiation levels for many gigayears. Could this maybe allow very sophisticated life forms to evolve? On the other hand, when these stars are young, they go through a contraction phase before entering main sequence, which is much longer than the contraction phase of G-dwarfs such as the Sun. During this phase, the stars are much brighter and hotter than later in their history. For TRAPPIST-1 planets this would mean they have been grilled by hot temperatures for about a billion years! Can they still retain some water after such a violent past? Can life form under such conditions? We don’t really know. In any event, it seems that water retaining and delivery might be a critical factor for TRAPPIST-1 planets’ habitability.

Since the planets are so densely-packed in the system, the masses of neighbouring planets as well as the mass of the star have a gravitational effect on each other – just as the Moon leads to high and low tides of Earth’s oceans. Only, the tidal forces acting on the TRAPPIST-1 planets would be much stronger, and could lead to immense energy being released in the interior of the planets due to tidal dissipation. Furthermore, the star itself appears to have a strong magnetic field. An electrical current is produced if a conductive material is embedded in a changing magnetic field, which is used, for example, to melt iron in induction furnaces. Similarly, the mantle of rocky planets are conductive and can experience enhanced energy release deep in the upper mantle due to induction heating.

Both induction heating and tidal heating can have a negative effect on the potential habitability of a rocky planet, since strong heating in the interior can be reflected by equally strong volcanic activity at the surface. This would lead to a hellish surface to live on! The interior may even be partly molten, leading to subsurface oceans of magma, which actually may be the case for TRAPPIST-1b and -1c. Even TRAPPIST-1d may be affected by strong volcanic events due to both induction and tidal heating of the interior. TRAPPIST-1f, -1g and -1h might be too cold at the surface to have liquid water, and might rather resemble our water-rich icy moons orbiting around Saturn and Jupiter. Hence TRAPPIST-1e, which receives only a little less stellar flux compared to Earth, may be the most interesting planet to visit in the system.

But what would life look like on such a planet?

The tidal forces described above also lead to a different effect: the planets would always face the star with only one side (this is called a tidal lock). Therefore, planets would have a day side that was always facing the star, and a night side immersed into eternal darkness and where no light ray is ever received from the star. Such a tidally-locked orbit is similar to the Moon-Earth system, as the Moon shows us always the same face – the “near-side” of the Moon. The other side, the “far-side” of the Moon, is only known to us due to lunar space missions. Can you imagine living at a place where it never gets dark? On the other hand, the luminosity from the star is very weak. Life on the TRAPPIST-1 planets might therefore actually look different than on Earth. To obtain the needed photons used in photosynthsesis (if this process would also evolve on these planets), life might evolve to favour a large variety of pigments that would enable it to make use of the full range of visible and infrared light – in other words, plants on these planets would appear black to us.

TRAPPIST-1 planets certainly still harbour many mysteries. They are a very good example how diverse the planets in the universe can be. If we set our imagination free… Black trees under the red star in the sky, which never sees a sunrise or sunset, powerful volcanoes filling the air with the ash and shaking the ground.
Very different from our Earth, isn’t it?

Further reading:

Kislyakova, K., Noack, L. et al. Magma oceans and enhanced volcanism on TRAPPIST-1 planets due to induction heating. Nature Astronomy 1, 878-885 (2017).

Gillon, M. et al. Seven temperate terrestrial planets around the nearby ultracool dwarf star TRAPPIST-1. Nature 542, 456-460 (2017).

Kiang, N. et al. The Color of Plants on Other Worlds. Scientific American, April 2008, 48-55 (2008).

Barr, A.C. et al. Interior Structures and Tidal Heating in the TRAPPIST-1 Planets. Astronomy and Astrophysics, in press.

Luger, R. et al. A seven-planet resonant chain in TRAPPIST-1. Nature Astronomy 1, 0129 (2017).

Scalo, J. et al. M stars as targets for terrestrial exoplanet searches and biosignature detection. Astrobiology 7(1), 85-166 (2007).

Ramirez, R.M. and Kaltenegger, L. The habitable zones of pre-main-sequence stars. The Astrophysical Journal Letters 797(2), L25 (2014).

On the influence of grain size in numerical modelling

On the influence of grain size in numerical modelling

The Geodynamics 101 series serves to showcase the diversity of research topics and methods in the geodynamics community in an understandable manner. We welcome all researchers – PhD students to Professors – to introduce their area of expertise in a lighthearted, entertaining manner and touch upon some of the outstanding questions and problems related to their fields. This month Juliane Dannberg from Colorado State University, discusses the influence of grain size and why it is important to consider it in numerical models. Do you want to talk about your research? Contact us!

Juliane Dannberg

When I started my PhD on geodynamic modelling, I was not aware that the size of mineral grains was something I might need to consider in my simulations. To me, grain size was something that sedimentologists need to describe rocks, and not something I had to deal with in my computations. In all the modelling papers I had read, if the mineral grain size was even mentioned, it was always assumed to be constant. However, it turns out that these tiny grains can have huge effects.

I first heard about the importance of grain size in a series of lectures given by Uli Faul when I participated in the CIDER summer program in 2014 (in case you’re interested, the lectures were recorded, and are available here and here). Primarily, I learned that for diffusion creep – the deformation mechanism people predominantly use in convection models – the viscosity does not only depend on temperature, but is also strongly controlled by the grain size, and that this grain size varies both in space and in time.

This made me think. If grain size in the mantle changes by several orders of magnitude, and the viscosity scales with the grain size to the power of 3, didn’t that mean that grain size variations could cause huge variations in viscosity that we do not account for in our models? Shouldn’t that have a major effect on mantle dynamics, and on the evolution of mantle plumes and subduction zones? How large are the errors we make by not including this effect? I was perplexed that there was such a major control on viscosity I had not thought about before, and wanted to look into this further myself.

Luckily, the multidisciplinary nature of CIDER meant that there were a number of people who could help me answer my questions. I teamed up with other participants1 interested in the topic, and from them I learned a lot about deformation in the Earth’s mantle.

How does the mantle deform?
For most of the mantle, the two important deformation mechanisms are diffusion and dislocation creep. In diffusion creep, single defects (or vacancies) – where an atom is missing in the lattice of the crystal – move through this lattice and the crystal deforms. For this type of deformation, the strain rate is generally proportional to the stress, which means that the viscosity does not depend on the strain rate. Many global mantle convection models use this kind of rheology, because it is thought to be dominant in the lower mantle, and also because it means that the problems described using this rheology are usually linear, which makes them easier to solve numerically. Furthermore, this is also the deformation mechanism that depends on the grain size. Usually, it is assumed that the viscosity scales with the grain size to the power of 3:

where ndiff is the diffusion creep viscosity, d is the grain size, m=3, T is the temperature, P is the pressure, and Adiff, E*diff, V*diff and R are constants.

However, if the stress is high (or the grain size is large, Figure 1), dislocation creep is the dominant deformation mechanism. In dislocation creep, linear defects – so-called dislocations – move through the crystal and cause deformation. In this regime, the viscosity depends on the strain rate, but not on grain size. Dislocation creep is generally assumed to be the dominant deformation mechanism in the upper mantle.

Figure 1: Deformation mechanisms in olivine

Why do grains grow or shrink?

In general, from an energy standpoint, larger crystals are more stable than smaller crystals, and so crystals tend to grow over time in a process called Ostwald ripening. The smaller the grains are, and the higher the temperature, the faster the grains grow. One the other hand, the propagation of dislocations through the grains causes so-called dynamic recrystallization, which reduces the grain size if the rock deforms due to dislocation creep. This means that there are always the competing mechanisms of grain growth and grain size reduction, and their balance depends on the dominance of either of the two deformation mechanisms described above – diffusion or dislocation creep:

The left-hand side of the equation describes the change of the grain size d over time, the first term on the right-hand side is grain growth (depending on grain size and temperature), the second term describes grain size reduction (depending the strain rate, the stress, and also on the grain size itself). The parameters Pg, kg, Eg, Vg, R, λ, c and γ are all constants.

If the flow field does not change, grains will evolve towards an equilibrium grain size, balancing grain growth and grain size reduction. In addition, the grain size may change when minerals cross a phase transition. If the mineral composition does not change upon crossing a phase transition (a polymorphic phase transition such as olivine–wadsleyite), there is almost no effect on grain size. But if the composition of the mineral that is stable after crossing the transition is different from the one before, the mineral breaks down, and the grain size is reduced, probably to the micrometer-scale [Solomatov and Reese, 2008].

And what does that mean for the dynamics of the Earth?

As there is a complex interaction between grain size evolution, mantle rheology and the deformation in the mantle, it is not straightforward to predict how an evolving grain size changes mantle dynamics. But it turned out that there had been a number of modelling studies investigating this effect. And they indeed found that grain size evolution may substantially influence the onset and dynamics of convection [Hall and Parmentier, 2003], the shape of mantle plumes [Korenaga, 2005], mixing of chemical heterogeneities [Solomatov and Reese, 2008], the seismic structure of the mantle [Behn et al., 2009], and the convection regime and thermal history of terrestrial planets [Rozel, 2012].

The long way to a working model…

But even knowing all of these things, it was still a long way to implement these mechanisms in a geodynamic modelling code, testing and debugging the implementation, and applying it to convection in the Earth. There were several reasons for that:

Firstly, large viscosity contrasts are already a problem for most solvers we use in our codes, and the strong dependence of viscosity on grain size means that viscosity varies by several orders of magnitude over a very small length scale in the model.

Secondly, considering an evolving grain size makes the problem we want to solve strongly nonlinear: Already in models with a diffusion–dislocation composite rheology and a constant grain size, the viscosity – which is needed to calculate the solution for the velocity – depends on the strain rate, making the momentum conservation equation nonlinear. But an evolving grain size introduces an additional nonlinearity: The viscosity now also depends (nonlinearly) on the grain size, whose evolution in turn depends on the velocity field. In terms of dynamics, this means that there is now another mechanism that can localize deformation. If the strain rate is large, the grain size is reduced due to dynamic recrystallization (as described above). A smaller grain size means a lower viscosity, which again enables a larger strain rate. Due to this feedback loop, velocities can become very high, up to several meters per year, which severely limits the time step size of a numerical model.

Finally, the equation (2) that describes grain size evolution is an ordinary differential equation in itself, and the time scales of grain growth and grain size reduction can be much smaller than changes in the flow field in the mantle. So, in order to model grain size evolution and mantle convection together, one has to come up with a way to separate these scales, and use a different (and probably much smaller) time step to compute how the grain size evolves. I remember that at one point, our models generated mineral grains the size of kilometers (whereas the grain sizes we expect in the mantle are on the order of millimeters), because we had not chosen the time step size properly. And on countless occasions, the code would just crash, because the problem was so nonlinear that a small change in just one parameter or a solution variable had such a large impact that material properties, velocities and pressures went outside of the range of what was physically reasonable.

However, after a lot of debugging, we could finally investigate how an evolving grain size would influence mantle dynamics. But see for yourself below. In an example from our models, plumes become much thinner when reaching the upper mantle, and cause much more vigorous small-scale convection when they interact with the lithosphere. Slabs bend rather than thicken, and accumulate as piles at the core-mantle boundary.

Figure 2: Comparison of plumes and slabs in models with and without grain size evolution. Modified from Dannberg et al., 2017

Of course, there are also many other areas where grain size evolution is important, and many recent studies are concerned with the influence of grain size on the Earth’s dynamic evolution. Dave Bercovici and his collaborators found that grain evolution and damage mechanisms may be a key factor for the onset of plate tectonics [e.g. Bercovici and Ricard, 2014, 2016]: Grain size reduction in shear zones could make them weak enough to for subduction initiation. The evolution of grain size may also be a major factor for focusing of melt to mid-ocean ridges [Turner et al., 2017], as it influences how fast the solid matrix can dilate and compact to let melt flow in and out. And if the Large Low Shear Velocity Provinces at the core-mantle boundary are indeed piles of hot material that are stable on long time scales, mineral grains would have a long time to grow and may play a crucial role for pile stability [Schierjott et al., 2017].

So if you do not include grain size evolution in your geodynamic models – which in many cases is just not feasible to do – I hope that you now have a better feeling for how that may affect your model results.

1The other researchers in my CIDER group were Zach Eilon, Ulrich Faul, Rene Gassmöller, Raj Moulik and Bob Myhill. I learned a lot about grain size in the mantle in particular from Bob Myhill and Ulrich Faul; I developed the geodynamic models together with Rene Gassmöller, and Zach Eilon and Raj Moulik investigated how the evolving grain size predicted by our models would influence seismic observations.


Solomatov, V. S., & Reese, C. C. (2008). Grain size variations in the Earth's mantle and the evolution of primordial chemical heterogeneities. Journal of Geophysical Research: Solid Earth, 113(B7).

Hall, C. E., & Parmentier, E. M. (2003). Influence of grain size evolution on convective instability. Geochemistry, Geophysics, Geosystems, 4(3).

Korenaga, J. (2005). Firm mantle plumes and the nature of the core–mantle boundary region. Earth and Planetary Science Letters, 232(1), 29-37.

Behn, M. D., Hirth, G., & Elsenbeck, J. R. (2009). Implications of grain size evolution on the seismic structure of the oceanic upper mantle. Earth and Planetary Science Letters, 282(1), 178-189.

Rozel, A. (2012). Impact of grain size on the convection of terrestrial planets. Geochemistry, Geophysics, Geosystems, 13(10).

Dannberg, J., Eilon, Z., Faul, U., Gassmöller, R., Moulik, P., & Myhill, R. (2017). The importance of grain size to mantle dynamics and seismological observations. Geochemistry, Geophysics, Geosystems.

Bercovici, D., & Ricard, Y. (2014). Plate tectonics, damage and inheritance. Nature, 508(7497), 513-516.

Bercovici, D., & Ricard, Y. (2016). Grain-damage hysteresis and plate tectonic states. Physics of the Earth and Planetary Interiors, 253, 31-47.

Turner, A. J., Katz, R. F., Behn, M. D., & Keller, T. (2017). Magmatic focusing to mid-ocean ridges: the role of grain size variability and non-Newtonian viscosity. arXiv preprint arXiv:1706.00609.

Schierjott, J., Rozel, A., & Tackley, P. (2017, April). Toward unraveling a secret of the lower mantle: Detecting and characterizing piles using a grain size-dependent, composite rheology. In EGU General Assembly Conference Abstracts (Vol. 19, p. 17433).

The lost Tethyan seaways: A deep-Earth and deep-time perspective on eastern Tethyan tectonics

The lost Tethyan seaways: A deep-Earth and deep-time perspective on eastern Tethyan tectonics

Every 8 weeks we turn our attention to a Remarkable Region that deserves a spot in the scientific limelight. Following from the first entry which showcased the Eastern Mediterranean, we move further east, and back in time, to the realm of the Tethys. The post is by postdoctoral researcher Sabin Zahirovic of the EarthByte Group and Basin GENESIS Hub, The University of Sydney.

Sabin Zahirovic

The southern and southeastern region of Eurasia (Fig. 1) represents one of the most tectonically complex areas in the world, with active deformation that appears as flurries of sometimes-deadly earthquakes, as well as the slower processes of mountain building. This region truly represents geophysical extremes, with the world’s tallest mountain ranges (namely, the Himalayas), and the planet’s deepest surface depressions (namely, the Mariana Trench), which collectively demonstrate the ongoing tectonic processes that are shaping, and re-shaping, our planet’s surface. Embedded in every rock, grain, mineral, and fossil in the region is the record of ancient tectonic activity, with a mosaic of exotic terranes sutured onto the continental margin. The Tethys was a vast oceanic domain at the southern margin of Eurasia and parts of its eastern history are now preserved in sutures onshore. The sutures are often peppered with ophiolites (Fig. 2), which are fragments of ancient ocean basins and gateways, the rest of which have been lost to subduction and now reside deep in the mantle. The opening and closure of these ocean basins has also been implicated in fundamental shifts of oceanic circulation and long-term sea level, which highlights the importance of better understanding the tectonic evolution of our planetary surface and deep churning interior.

Figure 1. Shuttle Radar Topography Mission data of the eastern Tethyan region displayed in the GPlates web portal, with 50X vertical exaggeration of elevation. (EarthByte Group and Scripps Institution of Oceanography)

Figure 2. Regional map highlighting the sutures (pink and blue lines) demarcating the terranes of the eastern Tethyan tectonic domain, as well as the ancient ocean basins that have since been consumed by subduction. Figure from Zahirovic et al. (2016b).

The first step to rewinding the tectonic clock requires undoing seafloor spreading in preserved ocean basins. Coincidentally, the geoscience community is celebrating the 50th anniversary of plate tectonics this week, with the understanding of seafloor spreading being a critical component of the plate tectonics paradigm. Through the work of numerous scientists such as Fred Vine, Harry Hess, Drummond Matthews, Dan McKenzie, amongst others, the ideas of seafloor spreading began to take shape. Another visionary who has received less acknowledgement is Marie Tharp, who painstakingly collected and manually plotted sonar measurements from voyages criss-crossing the Atlantic. Marie pieced together the morphology and topology of the mid-oceanic ridge systems (Fig. 3), which were mysterious bathymetric features at the time. When Marie superimposed earthquake epicentres on her maps, and suggested that this could be the missing piece of the continental drift puzzle, her significant insights were initially dismissed. However, Marie’s perseverance and determination laid much of the foundations in our understanding of these bathymetric features representing regions of seafloor spreading, from which a rush of subsequent work established the principles of plate tectonics.

Figure 3. Marie Tharp’s ground-breaking bathymetric map highlighting the elevated mid-oceanic ridge systems, which were later proved to be regions of seafloor spreading (© Marie Tharp 1977). Photograph inset of Marie in 2001 (Columbia University).

Following decades of data collection, including magnetic polarity reversal mapping of the oceans, the age of the oceanic crust and other features of the seafloor fabrics were compiled into global digital community models (e.g., Müller et al., 1997; Matthews et al., 2011). With these two ingredients, one can generate digital models of plate tectonic reconstructions, where the seafloor spreading of the ocean basins can be reversed to reveal the ancient configuration of continents (e.g., the reconstruction of Seton et al., 2012 using open-source community plate reconstruction software, GPlates, (Boyden et al., 2011)). However, due to the conservation of surface area (i.e., the Earth is not expanding), large gaps emerge that represent ancient subducted ocean basins.

As the Atlantic, Southern and Indian oceans opened, it occurred at the expense of the youngest (the Neo-Tethys), which was subducted along the southern margin of Eurasia. This process left a vast chain of arc volcanoes and orogens, stretching from the Alps in Europe to Sundaland in Southeast Asia (Fig. 1). However, only patchy evidence remains of how the continental terranes traversed the Tethyan ocean basins (Fig. 2), with existing methods largely relying on sedimentary affinities, fossil co-occurrences, and paleo-latitudinal estimates from paleomagnetic data (e.g., Metcalfe, 1994; Torsvik and Cocks, 2004). With the advent of modern seismic techniques (e.g., Li et al., 2008), it was possible to “image” the Earth’s mantle using seismic tomography (akin to a medical scan, Animation 1). This led teams of researchers to interpret the mantle structure (Fig. 4), and infer the history of subduction in the absence of preserved data on the surface (e.g., Hafkenscheid et al., 2006; van der Voo et al., 1999). In order to test interpretations of the India-Eurasia convergence that consumed the eastern Neo-Tethyan ocean basin, efforts were then launched to try and reproduce the India-Eurasia mantle structure using numerical models of mantle convection (e.g., Jarvis and Lowman, 2005).

Figure 4. Schematic interpretations of the India-Eurasia subduction history (A) using inferences from seismic tomography (van der Voo et al., 1999) (B), which were then tested using simple numerical experiments of mantle convection (Jarvis and Lowman, 2005) (C). Figures modified from Zahirovic et al. (2016b).

Animation 1. East-west sweep through the eastern Tethyan mantle, highlighting cold subducted slabs (blue) representing ancient Tethyan oceanic lithosphere (tomography model from Li et al., 2008).

With advances in community software platforms to model the evolution of entire plates and their boundaries through time, as well as a near-exponential increase in supercomputing resources, it became easier to model the coupled plate-mantle system in a time-evolving 3D spherical shell (e.g., using tools such as Aspect, CitcomS, and others). These advances enabled us, our colleagues, and the wider geo-community, to track subducting slabs from a complex and evolving global network of plate boundaries, including those associated with the eastern Tethyan tectonic domain. The new models provided an additional insight into the interaction of the deep Earth and surface processes, especially in such a tectonically-complex region like Southeast Asia.

For more than a decade, research has highlighted that Australia’s northern continental margin and the Sundaland continental promontory were low-lying and partially-flooded regions, likely due to the influence of sinking lithospheric plates from the eastern Tethyan slab graveyard (DiCaprio et al., 2009; DiCaprio et al., 2011; Spasojevic and Gurnis, 2012). What the models of mantle flow have revealed is that this effect, known as “dynamic topography”, is a transient signal through space and time, depending on the position of the continents in relation to large-scale mantle upwellings or downwellings. In the case of the eastern Tethys, these mantle downwellings were ephemeral, affecting the relative contribution of global sea level change and regional uplift or subsidence that drove shoreline retreat and advance over the continents. For Southeast Asia, we now think that the whole region was elevated and entirely emergent from about 80 to 40 million years ago, largely due to a hiatus in subduction along southern Eurasia between ~80 and 60 million years ago (Fig. 5). These notions complement traditional ideas of tectonic topography (where the surface experiences uplift or subsidence due to collisional or rifting processes), and allow us to consider the additional role of mantle flow (upwellings and downwellings contributing to several hundred meters in elevation change) over geological time in shaping the surface evolution of our planet.

Figure 5. Palaeogeographic reconstructions highlight that Southeast Asian topography was very different in the past, one dominated by a contiguous and emergent landmass about 80 to 40 million years ago, despite higher sea levels than today. Numerical models of the eastern Tethyan region suggests that the absence of subduction between 80 and 60 million years ago led to a “dynamic rebound” of the Sundaland continental promontory. However, the region became flooded again from about 40 million years ago because of a resumption in subduction, which led to regional subsidence, even with falling long-term sea levels. Figure adapted from Zahirovic et al. (2016a).

The plate tectonic reconstructions (Animation 2), along with the mantle flow models, have added significant insights into the evolution of the plate-mantle system. However, one emerging area of study is the exploration of how the shifting tectonic regimes have influenced other Earth systems – whether it be climate, ocean circulation, or biological evolution. In terms of the eastern Tethys, further research will help us uncover the role of tectonics and mantle flow in understanding the deep carbon cycle (at the plate-mantle scale), which has the potential to reveal the flux of carbon between shallow and deep Earth reservoirs, and their influence on deep-time climate (e.g., Jagoutz et al., 2016; Kent and Muttoni, 2013). For example, as the Neo-Tethys ocean basin was consumed, huge amounts of CO2 were emitted through a vast network of convergent plate boundaries that hosted arc volcanoes, leading to warmer greenhouse conditions. However, once India crashed into these subduction systems, the buoyant continental crust jammed and shut down subduction, turning off a major input of CO2 into the atmosphere. In addition, the mountain building processes in the Alpine-Himalayan mountain belt sequestered enormous amounts of CO2 from the atmosphere through chemical weathering of silicate rocks, influencing (and perhaps controlling) long-term cooling of the planet in the last 45 million years. These components collectively highlight the need for more integrative work that unpacks the complexity of regions like the eastern Tethys in order to understand the interaction of physical, chemical, and biological processes that have shaped planetary evolution.

Images from Animation 2 – please click the link below to download the movie file (10MB)


Animation 2 (please click the link above to download and view). The latest plate tectonic reconstructions of the eastern Tethys consider both the evolution of the continents, as well as the oceanic plate that carried it, much like a conveyor belt on the Earth’s surface. The reconstructions are embedded in global models with evolving plate boundaries (GPlates software), and are easier linked to numerical models of mantle convection (Zahirovic et al., 2016b), which provide an additional avenue to test our interpretations of past continental and oceanic arrangements

Note from author: This summary is by no means comprehensive, and represents a personal perspective of plate tectonics and eastern Tethyan geodynamics.


Boyden, J., Müller, R., Gurnis, M., Torsvik, T., Clark, J., Turner, M., Ivey-Law, H., Watson, R., and Cannon, J., 2011, Next-generation plate-tectonic reconstructions using GPlates, in Keller, G., and Baru, C., eds., Geoinformatics: Cyberinfrastructure for the Solid Earth Sciences: Cambridge, UK, Cambridge University Press, p. 95-114.

DiCaprio, L., Gurnis, M., and Müller, R. D., 2009, Long-wavelength tilting of the Australian continent since the Late Cretaceous: Earth and Planetary Science Letters, v. 278, no. 3, p. 175-185.

DiCaprio, L., Gurnis, M., Müller, R. D., and Tan, E., 2011, Mantle dynamics of continentwide Cenozoic subsidence and tilting of Australia: Lithosphere, v. 3, no. 5, p. 311-316.

Gurnis, M., Turner, M., Zahirovic, S., DiCaprio, L., Spasojevic, S., Müller, R., Boyden, J., Seton, M., Manea, V., and Bower, D., 2012, Plate Tectonic Reconstructions with Continuously Closing Plates: Computers & Geosciences, v. 38, no. 1, p. 35-42.

Hafkenscheid, E., Wortel, M., and Spakman, W., 2006, Subduction history of the Tethyan region derived from seismic tomography and tectonic reconstructions: Journal of Geophysical Research-Solid Earth, v. 111, no. B8, p. B08401.

Jagoutz, O., Macdonald, F. A., and Royden, L., 2016, Low-latitude arc–continent collision as a driver for global cooling: Proceedings of the National Academy of Sciences, v. 113, no. 18, p. 4935-4940.

Jarvis, G., and Lowman, J., 2005, Sinking slabs below fossil subduction zones: Physics of The Earth and Planetary Interiors, v. 152, no. 1-2, p. 103-115.

Kent, D. V., and Muttoni, G., 2013, Modulation of Late Cretaceous and Cenozoic climate by variable drawdown of atmospheric pCO 2 from weathering of basaltic provinces on continents drifting through the equatorial humid belt: Climate of the Past, v. 9, no. 2, p. 525-546.

Li, C., van der Hilst, R., Engdahl, E., and Burdick, S., 2008, A new global model for P wave speed variations in Earth's mantle: Geochemistry, Geophysics, Geosystems, v. 9, no. 5, p. 21.

Matthews, K. J., Müller, R. D., Wessel, P., and Whittaker, J. M., 2011, The tectonic fabric of the ocean basins: Journal of Geophysical Research, v. 116, no. B12, p. B12109.

Metcalfe, I., 1994, Gondwanaland origin, dispersion, and accretion of East and Southeast Asian continental terranes: Journal of South American Earth Sciences, v. 7, no. 3, p. 333-347.

Müller, R., Roest, W., Royer, J., Gahagan, L., and Sclater, J., 1997, Digital isochrons of the worldís ocean floor: JGR, v. 102, no. B2, p. 3211-3214.

Seton, M., Müller, R., Zahirovic, S., Gaina, C., Torsvik, T., Shephard, G., Talsma, A., Gurnis, M., Turner, M., Maus, S., and Chandler, M., 2012, Global continental and ocean basin reconstructions since 200 Ma: Earth-Science Reviews, v. 113, no. 3-4, p. 212-270.

Spasojevic, S., and Gurnis, M., 2012, Sea level and vertical motion of continents from dynamic earth models since the Late Cretaceous: AAPG bulletin, v. 96, no. 11, p. 2037-2064.

Torsvik, T., and Cocks, L., 2004, Earth geography from 400 to 250 Ma: a palaeomagnetic, faunal and facies review: Journal of the Geological Society, v. 161, no. 4, p. 555-572.

van der Voo, R., Spakman, W., and Bijwaard, H., 1999, Tethyan subducted slabs under India: Earth and Planetary Science Letters, v. 171, no. 1, p. 7-20.

Zahirovic, S., Flament, N., Müller, R. D., Seton, M., and Gurnis, M., 2016a, Large fluctuations of shallow seas in low-lying Southeast Asia driven by mantle flow: Geochem. Geophys. Geosys, v. FRONTIERS IN GEOSYSTEMS: Deep Earth - surface interactions, no. 17.

Zahirovic, S., Matthews, K. J., Flament, N., Müller, R. D., Hill, K. C., Seton, M., and Gurnis, M., 2016b, Tectonic evolution and deep mantle structure of the eastern Tethys since the latest Jurassic: Earth Science Reviews, v. 162, p. 293-337.