Geodynamics 101

Subduction through the mantle transition zone: sink or stall?

Subduction through the mantle transition zone: sink or stall?

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 latest ‘Geodynamics 101’ post, Saskia Goes, Reader at Imperial College London, UK, discusses the fate of subducting slabs at the mantle transition zone.

Saskia Goes

Subducting plates can follow quite different paths in their life times. While some sink straight through the upper into the lower mantle, others appear to stall in the mantle transition zone above 660 km depth. Geodynamicists have long puzzled about what controls these different styles of behaviour, especially because there appear to be correlations between sinking or stalling with faster or slower plate motions and mountain building or ocean basin formation, respectively. In the long run, how easily slabs sink through the transition zone controls how efficiently material and heat are circulated in the mantle.

The word subduction derives from the Latin verb subducere, which means pulled away from below, but metaphorically can mean to lose footing or remove secretly. Definitely, when Wegener first proposed continental drift, people were unaware that subduction is removing plates from the Earth’s surface. We now know this process is not quite so secret. The plates creak in earthquakes as they sink into the mantle, in some cases all the way through the mantle transition zone to about 700 km depth. Furthermore, where the subducting plate bends below the overriding plate, it creates deep-sea trenches with prominent gravity and geoid signals. This bending is a very important part of subduction dynamics, as I’ll explain below.

The seismic Wadati-Benioff zones and gravity expressions were sufficient clues of the location of the downwelling limbs of a mantle convection system to help acceptance of plate tectonics in the 1960s. However, it took another twenty odd years until seismology yielded images of cold plates sinking into the mantle, and it turned out that the plates extend beyond the seismic Wadati-Benioff zones [Van der Hilst et al., 1991; Zhou and Clayton, 1990]. These images showed that some subducting plates flatten in the mantle transition zone (e.g. below Japan and Izu-Bonin), while others continue with little to no deflection into the lower mantle (e.g., below the Northern Kuriles and Marianas) (Fig. 1). Soon after, it was realised that many of the places where the slabs are flat in the transition zone have a history of trench retreat [Van der Hilst and Seno, 1993]. Furthermore, mapping of seafloor ages revealed that flat slabs tend to form where plates older than about 50 Myr are subducted [Karato et al., 2001; King et al., 2015].

Figure 1: Variable modes of slab-transition zone interaction

Many mechanisms have been proposed for the variable slab transition-zone interaction. We recently reviewed the geodynamic and observational literature and combined these insights with those from our own set of mechanical and thermo-mechanical subduction models [Goes et al., 2017]. This effort shows that not one single mechanism, but an interplay of several mechanisms is the likely cause of the observed variable subduction behaviour.

It has long been realised that viscosity increases with depth into the mantle, quite possibly including jumps at the major phase transitions in the mantle transition zone. The ringwoodite-postspinel transition that is responsible for the global 660 km seismic discontinuity, usually taken as the base of the upper mantle, is an endothermic transition under most of the conditions prevailing in the mantle today. This means that the transition will take place at a higher pressure and thus depth in the subducting plate than the surrounding mantle, rendering the plate locally buoyant with respect to the mantle. Both these factors hamper the descent of the subducting plate through the transition zone. However, a viscosity increase within acceptable bounds (as derived from geoid and postglacial rebound modelling) can slow sinking, but does not lead to stalling material. By contrast, the phase transition can lead to stalling, as well as an alternation of periods of accumulation of material in the transition zone and periods where this material flushes rapidly into the lower mantle, at least in convection models without strong plates. But does this work with strong plates?

Making dynamic models of subduction with strong plates is challenging because the models need to capture strong strength gradients between the core of the plate and the underlying mantle, allow for some form of plate yielding, maintain a weak zone between the two plates and adequately represent the effect of plate bending (a free-surface effect). Most models prescribe at least part of the system by imposing velocities and/or plate geometries. This however needs to be done with great care and consideration for what forcing such imposed conditions imply.

“Pulled away from below” is a good description of the dynamics of subduction. Subduction is primarily driven by slab pull, the gravitational force on the dense subducting plate [Forsyth and Uyeda, 1975]. And to “lose footing” reminds us that gravity is the main driving force. Gravity tries to pull the plate straight down (Fig. 2), so the easiest way for a plate to subduct is to fall into the mantle, a process that leads to trench retreat [Garfunkel et al., 1986; Kincaid and Olson, 1987]. Besides letting the plate follow the path of gravity, subduction by trench retreat has the other advantage that the plate does not need to bend too much. Bending a high-strength plate takes significant energy. Some studies have shown that if plates are assigned laboratory-based rheologies, such bending can easily take up all of the gravitational potential energy of the subducting plate [Conrad and Hager, 1999], so if plates are to sink into the mantle, they have to do this by minimising the amount of energy used for bending into the trench. As a consequence, strong and dense plates prefer to subduct at smaller dip angles while weaker and lighter plates can be bent to subduct more vertically [Capitanio et al., 2007].

Figure 2: If subduction occurs freely, i.e., driven by the pull of gravity on the dense slab with sinking resisted by the viscous mantle, it is usually energetically most favourable to subduct by trench retreat.

The angle at which plates subduct strongly affects how they subsequently interact with viscosity or phase interfaces (Fig. 3). Steeply dipping plates will buckle and thicken when they encounter resistance to sinking. This deformation facilitates further sinking, as a bigger mass. But plates that reach the interface at a lower dip may be deflected. Such deflected plates have a harder time sinking onwards, both because the high viscosity resistance is now distributed over a wider section of the plate and due to the spread-out additional buoyancy from the depressed endothermic phase boundary.

Figure 3: The subduction angle largely determines how the slab interacts with viscosity and phase changes.

So, variable plate density and strength can lead to variable behaviour of subduction in the transition zone. And we know plates have variable density and strength. Older plates are denser and if strength is thermally controlled, as most lab experiments predict, also stronger than younger plates. This implies that older plates can drive trench retreat more easily than young plates. And indeed this matches observations that significant trench retreat has only taken places where old plates subduct. Furthermore, significant trench retreat will facilitate plate flattening in the transition zone, consistent with the observation that flat plates tends to underlie regions with a history of trench retreat (even if that does not always mean trench motions are high at the present day). This mechanism can also explain why flat slabs tend to be associated with old plate subduction.

So what about the role of other proposed mechanisms? Our models with strong slabs show that only when slabs encounter both an increase in viscosity (which forces the slabs to deform or flatten) and an endothermic phase transition (which can lead to stalling of material in the transition zone) do we find the different modes of slab dynamics. Neither a viscosity increase alone, nor an endothermic phase transition alone leads to mixed slab dynamics.

Other factors likely contribute to the regional variability. In the cold cores of the slabs, some phases may persist metastably, thus delaying the transformations to higher density phases to a larger depth. Metastability will be more pervasive in colder old plates thus making older plates more buoyant and hence resistant to sinking than young ones. In combination with trench retreat facilitated by a strong slab at the trench, this can further encourage slab flattening [Agrusta et al., 2014; King et al., 2015]. Phase transformations may also lead to slab weakening in the transition zone because they can cause grain size reduction. Such weakening can aid slab deflection [Čížková et al., 2002; Karato et al., 2001]. However, several studies have shown that transition zone slab strength is less important than slab strength at the trench, which governs how a slab starts sinking through the transition zone.

The Earth is clearly more complex than the models discussed. For example, present-day plate dip angles display various trends with plate age at the trench. Lateral variations in plate strength and buoyancy can complicate subduction behaviour. Furthermore, forces on the upper plate and large-scale mantle flow may also impede or assist trench motions and may thus affect or trigger changes in how slabs interact with the transition zone [Agrusta et al., 2017]. All these factors remain to be fully investigated. However, the first order trends of subduction-transition zone interaction can be understood as a consequence of plates of various ages interacting with a viscosity increase and endothermic phase change.

 Agrusta, R., J. van Hunen, and S. Goes (2014), The effect of metastable pyroxene on the slab dynamics, Geophys. Res. Lett., 41, 8800-8808.
 Agrusta, R., S. Goes, and J. van Hunen (2017), Subducting-slab transition-zone interaction: stagnation, penetration and mode switches, Earth Planet. Sci. Let., 464, 10-23.
 Capitanio, F. A., G. Morra, and S. Goes (2007), Dynamic models of downgoing plate buoyancy driven subduction: subduction motions and energy dissipation, Earth Planet. Sci. Lett., 262, 284-297.
 Čížková, H., J. van Hunen, A. P. van der Berg, and N. J. Vlaar (2002), The influence of rheological weakening and yield stress on the interaction of slabs with the 670 km discontinuity, Earth Plan. Sci. Let., 199(3-4), 447-457.
 Conrad, C. P., and B. H. Hager (1999), Effects of plate bending and fault strength at subduction zones on plate dynamics, J. Geophys. Res., 104(B8), 17551-17571.
 Forsyth, D. W., and S. Uyeda (1975), On the relative importance of driving forces of plate motion. , Geophys. J. R. Astron. Soc. , 43, 163-200.
 Garfunkel, Z., C. A. Anderson, and G. Schubert (1986), Mantle circulation and the lateral migration of subducted slab, J. Geophys. Res., 91(B7), 7205-7223.
 Goes, S., R. Agrusta, J. van Hunen, and F. Garel (2017), Subduction-transition zone interaction: a review, Geosphere, 13(3. Subduction Top to Bottom 2), 1-21.
 Karato, S. I., M. R. Riedel, and D. A. Yuen (2001), Rheological structure and deformation of subducted slabs in the mantle transition zone: implications for mantle circulation and deep earthquakes, Phys. Earth Plan. Int., 127, 83-108.
 Kincaid, C., and P. Olson (1987), An experimental study of subduction and slab migration, J. Geophys. Res., 92(B13), 13,832-813,840.
 King, S. D., D. J. Frost, and D. C. Rubie (2015), Why cold slabs stagnate in the transition zone, Geology, 43, 231-234.
 Van der Hilst, R. D., and T. Seno (1993), Effects of relative plate motion on the deep structure and penetration depth of slabs below the Izu-Bonin and Mariana island arcs, Earth Plan. Sci. Let., 120, 395-407.
 Van der Hilst, R. D., E. R. Engdahl, W. Spakman, and G. Nolet (1991), Tomographic imaging of subducted lithosphere below northwest Pacific island arcs, Nature, 353, 37-43.
 Zhou, H.-w., and R. W. Clayton (1990), P and S Wave Travel Time Inversions for Subducting Slab Under the Island Arcs of the Northwest Pacific, J. Geophys. Res., 95(B5), 6829-6851.

Being both strong and weak

Being both strong and weak

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 latest ‘Geodynamics 101’ post, Postdoc Anthony Osei Tutu of GFZ Potsdam shares the outcomes of his PhD work, showing us that, like the lithosphere, it is OK to be weak sometimes!

Strength is not everything in achieving one’s goal. The lithospheric plate acts both strong and weak at times. This dual characteristic of the outermost part of the Earth, the crustal-lithospheric shell, is thought to have sustained plate tectonics throughout Earth’s history, in the presence of other controlling mechanisms such as the weak asthenospheric layer (Bercovici et al. 2000; Karato 2012). In the world of the lithospheric plates there is the saying “I might be strong and unbreakable, but sometimes and somewhere, I am very weak, soft and brittle” and this allows the plates to accommodate each other in their relative movements.

We all sometimes need to bring out the soft part in us to accommodate others such as friends, family or colleagues. For example, my graduate school, the Helmholtz-Kolleg GEOSIM, an experiment by the Helmholtz Association, GFZ-Potsdam, University of Potsdam and Free University of Berlin, brought together two or more experts in mathematics and geosciences to collaborate on and serve as PhD supervisors for answering some of Earth Sciences’ pressing questions. The many, many benefits of this multidisciplinary PhD supervising approach also came with challenges. Sometimes, the different supervisors would make opposing/contrasting suggestions to investigate a particular problem according to the experience of some students and myself. Then it falls on you as the student to stand firm (i.e. be strong) on what you believe works for your experiments and at the same time to be receptive (i.e. flexible or soft) to the different suggestions, while keeping in mind the limited time you have as a PhD student.

Figure 1: Schematic plot of the conditions in a subduction system (left) aiding or (right) hindering global plate motions.

The both strong and weak behavior of the lithospheric plates was one of the conclusions of my PhD study. Besides the strong plate interiors (Zhong and Watts 2013), weak regions along the plate boundaries, aided by sediment and water (see Fig. 1), are required to give the low friction between the subducting and overriding plates (Moresi and Solomatov 1998; Sobolev and Babeyko 2005), combined with a less viscous sublithospheric mantle. This combination was key to match the magnitude and direction of present-day global plate motions in the numerical modeling study (Osei Tutu et al. 2018). I used the global 3D lithosphere-asthenosphere numerical code SLIM3D (Popov and Sobolev 2008) with visco-elasto-plastic rheology coupled to a mantle flow code (Hager and O’Connell 1981) for the investigation. To understand the influence of intra-plate friction (brittle/plastic yielding) and asthenospheric viscosity on present-day plate motions, I tested a range of strengths of the plate boundary. Past numerical modeling studies (Moresi and Solomatov 1998; Crameri and Tackley 2015) have suggested that small friction coefficients (μ < 0.1, yield stress ~100 MPa) can lead to plate tectonics in models of mantle convection. This study shows that in order to match present-day plate motions and net rotation, the static frictional parameter must be less than 0.05 (15 MPa yield stress). I am able to obtain a good fit with the magnitude and orientation of observed plate velocities (NUVEL-1A) in a no-net-rotation reference frame with μ < 0.04 and a minimum asthenosphere viscosity of 5•1019 Pas to 1•1020 Pas (Fig. 2). The estimates of net-rotation (NR) of the lithosphere suggest that amplitudes of ~0.1– 0.2 °/My, similar to most observation-based estimates, can be obtained with asthenosphere viscosity cutoff values of ~1•1019 Pas to 5•1019 Pas and a friction coefficient μ < 0.05.

Figure 2: Set of predicted global plate motions for varying asthenosphere viscosity and plate boundary frictions, modified after Osei Tutu et al. (2018). Rectangular boxes show calculations with RMS velocities comparable to the observed RMS velocity of NUVEL-1A (DeMets et al. 2010).

The second part of my PhD study focused on the responses of the strong plate interiors to the convecting mantle below by evaluating the influence of shallow and deep mantle heterogeneities on the lithospheric stress field and topography. I explored the sensitivity of the considered surface observables to model parameters providing insights into the influence of the asthenosphere and plate boundary rheology on plate motion by testing various thermal-density structures to predict stresses and topography. Lithospheric stresses and dynamic topography were computed using the model setup and rheological parameters that gave the best fit to the observed plate motions (see rectangular boxes in Fig. 2). The modeled lithosphere stress field was compared the World Stress Map 2016 (Heidbach et al. 2016) and the modeled dynamic topography to models of observed residual topography (Hoggard et al. 2016; Steinberger 2016). I tested a number of upper mantle thermal-density structures. The thermal structure used to calculate the plate motions before is considered the reference thermal-density structure, see also Osei Tutu et al. (2017). This reference thermal-density structure is derived from a heat flow model combined with a sea floor age model. In addition I used three different thermal-density structures derived from global S-wave velocity models to show the influence of lateral density heterogeneities in the upper 300 km on model predictions. These different structures showed that a large portion of the total dynamic force generating stresses in the crust/lithosphere has its origin in the deep mantle, while topography is largely influenced by shallow heterogeneities. For example, there is hardly any difference between the stress orientation patterns predicted with and without consideration of the heterogeneities in the upper mantle density structure across North America, Australia and North Africa. However, inclusion of crustal thickness variations in the stress field simulations (as shown in Fig. 3a) resulted in crustal dominance in areas of high altitude in terms of stress orientation, for example in the Andes and Tibet, compared to the only-deep mantle contributions (as shown in Fig. 3b).

Figure 3: Modeled lithosphere stress field in the Andes considering (a) crustal thickness variations from the CRUST 1.0 model as well as lithospheric variations and (b) uniform crustal and lithospheric thicknesses.

The outer shell of the solid Earth is complex, exhibiting different behaviors on different scales. In our quest to understand its dynamics, we can learn from the lithospheric plate’s life cycle how to live our lives and preserve our existence as scientist-humans by accommodating one another. After all, they have existed for billions of years.



Bercovici, David, Yanick Ricard, and Mark A. Richards. 2000. “The Relation Between Mantle Dynamics and Plate Tectonics: A Primer.” 5–46.

Crameri, Fabio and Paul J. Tackley. 2015. “Parameters Controlling Dynamically Self-Consistent Plate Tectonics and Single-Sided Subduction in Global Models of Mantle Convection.” Journal of Geophysical Research: Solid Earth 120(5):3680–3706, 10.1002/2014JB011664.

DeMets, Charles, Richard G. Gordon, and Donald F. Argus. 2010. “Geologically Current Plate Motions.” Geophys. J. Int 181:1–80.

Hager, BH and RJ O’Connell. 1981. “A Simple Global Model of Plate Dynamics and Mantle Convection.” Journal of Geophysical Research: Solid Earth, 86(B6):4843–4867, 10.1029/JB086iB06p04843.

Heidbach, Oliver, Mojtaba Rajabi, Moritz Ziegler, Karsten Reiter, and Wsm Team. 2016. “The World Stress Map Database Release 2016 -Global Crustal Stress Pattern vs. Absolute Plate Motion.” Geophysical Research Abstracts EGU General Assembly 18:2016–4861.

Hoggard, M. J., N. White, and D. Al-Attar. 2016. “Global Dynamic Topography Observations Reveal Limited Influence of Large-Scale Mantle Flow.” Nature Geoscience 9(6):456–63, 10.1038/ngeo2709.

Karato, Shun-Ichiro. 2012. “On the Origin of the Asthenosphere.” Earth and Planetary Science Letters 321–322:95–103.

Moresi, Louis and Viatcheslav Solomatov. 1998. “Mantle Convection with a Brittle Lithosphere: Thoughts on the Global Tectonic Styles of the Earth and Venus.” Geophysical Journal International 133(3):669–82, 10.1046/j.1365-246X.1998.00521.x.

Osei Tutu, A., S. V Sobolev, B. Steinberger, A. A. Popov, and I. Rogozhina. 2018. “Evaluating the Influence of Plate Boundary Friction and Mantle Viscosity on Plate Velocities.” Geochemistry, Geophysics, Geosystems n/a-n/a, 10.1002/2017GC007112.

Popov, A. A. and S. V. Sobolev. 2008. “SLIM3D: A Tool for Three-Dimensional Thermomechanical Modeling of Lithospheric Deformation with Elasto-Visco-Plastic Rheology.” Physics of the Earth and Planetary Interiors 171(1–4):55–75.

Sobolev, S. V. and A. Y. Babeyko. 2005. “What Drives Orogeny in the Andes?” Geology 33(8).

Steinberger, Bernhard. 2016. “Topography Caused by Mantle Density Variations: Observation-Based Estimates and Models Derived from Tomography and Lithosphere Thickness.” Geophysical Journal International 205(1):604–21, 10.1093/gji/ggw040.

Osei Tutu, A., B. Steinberger, S. V Sobolev, I. Rogozhina, and A. Popov. 2017. "Effects of Upper Mantle Heterogeneities on Lithospheric Stress Field and Dynamic Topography." Solid Earth Discuss.,, in review, 2017

Zhong, Shijie and A. B. Watts. 2013. “Lithospheric Deformation Induced by Loading of the Hawaiian Islands and Its Implications for Mantle Rheology.” Journal of Geophysical Research: Solid Earth 118(11):6025–48, 10.1002/2013JB010408.

Finding the forces in continental rifting

Finding the forces in continental rifting

Luke Mondy

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 latest ‘Geodynamics 101’ post, PhD candidate Luke Mondy from the EarthByte Group at the University of Sydney blogs about some impressively high-resolution numerical models of ‘rotational rifting,’ and the role of gravity. Luke also shares a bit about the journey behind this work, which recently appeared in Geology.


In geodynamic modelling, we’re always thinking about forces. It’s a balancing act of plate driving forces potentially interacting with the upwelling mantle, or maybe sediment loading, or thermal relaxation… the list goes on.

Figure 1: A summary of the forces interacting during continental rifting, from Brune, 2018.

But the thing that underpins all of these forces, fundamentally, is our favourite but oft forgotten force: gravity. Here, I’ll tell the story of investigating a numerical model of continental rifting and discovering – or rather, rediscovering – the importance of gravity as a fundamental force in driving Earth dynamics.

How it started – a side project!

A few years ago, my colleagues and I were granted access to not just one, but two, big supercomputers in Australia: Raijin, and Magnus. Both were brand new and raring to go – but we needed something big to test them out on. At the time, 3D geodynamic models were typically limited to quite low resolution, since they can be so computationally demanding, but since we had access to this new power, we decided to see how far we could push the computers to address a fundamentally 3D problem.

2D vs 3D

Historically, subduction and rifting have been ideal settings to model as they can be constrained to two dimensions while still retaining most of their characteristic properties.

Figure 2. A 2D subduction model. Despite being ‘only’ two dimensions, the fundamental and interesting aspects of the problem are still captured by the model. Figure from Rey et al., 2014.

However, as tremendously useful as these models have been, many interesting problems in geodynamics are fundamentally three dimensional. The obvious example is global mantle convection, but we are starting to see more and more papers addressing both rifting and subduction problems that require 3D contexts, for example: continental accretion (Moresi et al., 2014), metamorphic core complex formation (Rey et al., 2017), or oblique rifting (Brune et al., 2012).

Typically when we model a rift in 2D, the dimensionality implies that we are looking at orthogonal rifting – that the plates move away from each other perpendicular to the rift axis. Since 2D models cannot account for forces in the third dimension, they are only suitable for when the applied tectonic forces pull within the plane of the model – that is, when the 2D model lies along a small circle of an Euler pole.

Euler poles have another interesting geometric property – the velocity of extension between two plates changes as we move closer or further away from the Euler pole: zero velocity at the pole itself, and fastest at the equator to the pole (Lundin et al., 2014).

Figure 3. Left: From Lundin et. al. (2014), the figure shows the geometric relationship of increasing rifting velocity as the distance from the pole increases. Right: the same relationship graphed out, showing the cosine curve (Kearey 2009).

This leads to differing extension velocities along the length of the rift axis. Extension velocities are a huge control on the resulting geodynamics (e.g., Buck et al., 1999). Employing a series of 2D models along a rift axis (Brune et al., 2014) has been used to show how these dynamics change, but misses out on the three-dimensionality of the problem – how do these differing and diachronous dynamics interact with each along other the rift margin as it forms?

Rotational Rifting

We decided to attempt to model this sort of rifting, as we termed it “rotational rifting”. Essentially we linked up the 2D slices along the rift axis into one big 3D model – so that we have slow extension towards the Euler pole, and fast extension away from it.

To do this, we ended up using the code Underworld (at the time version 1.8 – but their 2.0 version is the best place to start!), and a framework developed inside the EarthByte group at the University of Sydney called the ‘Lithospheric Modelling Recipe’, or LMR.


Figure 4. Map view of the two experiments. Arrows show the velocity boundary conditions applied. Note they are perpendicular to the model domain – we thought long and hard about this choice, and explain it fully in the Data Repository.

Using the LMR, we set up two 3D experiments: both are 1000 km by 500 km along the surface, and 180 km deep. The ‘orthogonal’ experiment is modelled at the equator to the pole – so the velocities along the walls are the same all the way along the rift axis. The ‘rotational’ experiment is very close to the Euler pole (where the rate of extension velocity change is greatest), from 89 degrees to 79 degrees (90 degrees being the Euler pole), which gives an imposed velocity at the slow end (89 degrees) of 0.5 cm/yr and at the fast end (79 degrees) 5.0 cm/yr.


Since we wanted to stress test the supercomputers, we ran these experiments at just under 2 km grid resolution (256 x 512 x 96). This meant each experiment ended up using about 2.5 billion particles to track the materials! The 2 km grid size is an important milestone – to properly resolve faulting, sub-2 km grid sizes are required (Gerya, 2009).

The results!

So we ran the experiments, and compared the results! To give a broad overview of what we found, here’s a nice animation:

Figure 5. Top: Animation showing the orthogonal experiment from a south-west perspective (with the Euler pole being the ‘north’ pole). The light grey layers show the upper crust, dark grey the lower crust. Half of the crust has been removed to show the lithospheric mantle topography. The blue to the white colours show the lithospheric mantle temperature, and from white to red shows the asthenospheric temperature. Bottom: As above but for the rotational experiment. Notice that the asthenospheric dome migrates along the rift towards the Euler pole.

What to do now?

Cool looking experiments, of course! The supercomputers had been able to handle the serious load we put on them (it took about 2 weeks per experiment, on ~800 cpus), so that part of the project was a success. But what about the experiments themselves – did switching to 3D actually tell us anything useful?

What we expected…

The things we expected were there. The orthogonal experiment behaved identically to a 2D model. For the rotational experiment, we found the style of faulting changed and evolved along the rift axis, and seemed to match up nicely with the 2D work about differing extension rates. We were able to identify phases of rifting via strain patterns, which were similar to those described by Lavier and Manatschal (2006), and seemed to match the outputs of the series of 2D models along a rift axis.

Figure 6. Map view of strain-rate of the rotational experiment through time. The phases (1 through 4, representing different modes of deformation) migrate along the rift towards the Euler pole.

What we didn’t expect…

Almost on a whim, we decided to start looking into the tectonic regime. Using the visualization program Paraview, we calculated the eigenvectors of the deviatoric stress and assigned a tectonic regime (blue for extension, red for compression, green for strike-slip, and white for undetermined), following a similar scheme to the World stress map (Zoback, 1992). Apologies to colour blind folk!

Here’s what a selected section of the orthogonal experiment surface looks like through time:

Figure 7. The stress regimes at the surface of the orthogonal experiment (clipped to y = ~400 to ~600 km).

Not really that surprising – we found mostly extension everywhere, with a bit of compression when the central graben sinks down and gets squeezed. However, it was a little bit surprising to see the compression come back on the rift flanks.

But when we applied the same technique to the rotational experiment, we found this on the surface:

Figure 8. The stress regimes at the surface of the rotational experiment. The three numbers at the top represent the total extension at y = 0 km, 500 km, and 1000 km respectively.

Now all of a sudden we’re seeing strike-slip stress regimes in different areas of the experiment!

The above figures displaying the stress in the experiments so far have been of the surface – where one of the principal stress axes must be vertical – but our colouring technique does not limit us to just the surface. We noticed when looking at cross-sections that the lithospheric mantle was also showing unexpected stress regimes!

Figure 9. Slices at y = 500 km across the rift axis (right in the middle). Coloured areas show where the plunge of one principal stress axis is >60 degrees. Both experiments have the same applied extension velocity at y = 500 km, and so total extension is equivalent between experiments.

In most of the lithosphere, the strain rate is still very small, not enough to notice much deformation (1e-16 to 1e-18 1/s). But a few puzzling questions were raised: why do we see compressional tectonic regimes in the orthogonal experiment; and why do we also see strike-slip regimes in the rotational experiment?

Gravitational Potential Energy (GPE)

It quickly became apparent that these stress changes were related to the upwelling asthenosphere, as the switch between regimes was well timed to when the asthenosphere would approach the Moho – about 40 km depth. This gave us the hint that perhaps buoyancy forces were at play. We used Paraview again to calculate the gravitational potential energy at each point on the surface (taking into account all the temperature dependent densities, detailed topography, and so on), and produced these maps:

Figure 10. A time series showing the gravitational potential energy (GPE) at each point on the surface of the rotational experiment. Only half the surface is shown because it is symmetrical. The small triangle notch is where we determined the rift tip to be located (where 1/(beta factor) < 0.2).

What we saw confirmed our suspicions – the rise of the asthenospheric dome induces a gravitational force that radiates outwards. The juxtaposition of the hot, yet still quite heavy, asthenospheric material, next to practically unthinned crust on both the rift flanks and ahead of the rift tip, produces a significant force.

But why the switch to compression or strike-slip tectonic regimes in an otherwise extensional setting? In the case of the orthogonal model, the force (aka the difference in GPE) is perpendicular to the rift axis, since the dome rises synchronously along the axis. When this force overcomes the far-field tectonic force (essentially the force required to drive our experiment boundary conditions), the stress regime changes from extension to compression.

However, in the rotational experiment, the dome is larger the further away from the Euler pole, and so instead the gravitational force radiates outwards from the dome. Now the stress in the lithospheric mantle has to deal with not only the force induced from the upwelling asthenosphere right next to it, but also from along the rift axis (have a look at the topography of the lithospheric mantle in Fig. 5). These combined forces end up rotating the principal stresses such that sigma_2 stands vertical and a strike-slip regime is generated.

We also see the gravitational force manifest in other ways. Looking at the along axis flow in the asthenosphere, the experiment initially predicts a suction force towards the rapidly opening end of the model (away from the Euler pole), similar to Koopmann et al. (2014). But once the dome is formed, we see a reversal of this flow, back towards the Euler pole, driven by gravitational collapse. This flow appears to apply a strong stress to the crust surrounding the dome, reaching upwards of 50 MPa in some places.

Figure 11. A: The direction of flow at the lithosphere-asthenosphere boundary in the centre of the rift. Early in the experiment, we see suction towards the fast end of the rift, while later in the experiment, we see a return flow. The dashed line shows the flow after the tectonic boundary conditions have been removed. B,C: cross-sections showing stress and velocity arrows from the experiment just after the tectonic boundary conditions have been removed.

How do we know it’s gravity?

To test this idea further, we ran some additional experiments. First, we let the rotational experiment run for about 3.6 Million years, and then ‘stopped’ the tectonics (changed the side velocity boundary conditions to 0 cm/yr) – leaving gravity as the only driving force. We saw that the return flow towards the Euler pole was still present (though reduced). By running some more rotational experiments with either doubled or halved Euler pole rotational rate, we saw that the initial suction magnitude correlates with the change in opening velocity, but the return flow to the Euler pole is almost identical, giving further evidence that this is gravity driven.

What about the real world?

We numerical modellers love to stay in the world of numbers – but alas sometime we must get our hands dirty and look at the real world – just to make sure our models actually tell us something useful!

Despite our slightly backwards methodology (model first, check nature second), it did give us an advantage: our experiments were producing predictions for us to go and test. We had our hypothesis – now to see if it could be validated.

So we went out and looked for examples of rifting near an Euler pole, and the two most notable we found were in the Woodlark Basin, Papua New Guinea, and the Galapagos Rise in the Pacific. Despite the ‘complications’ of the natural world (things like sediment loading, pre-existing weakness in the crust, etc. – things that get your hands dirty), we found a striking first order relationship between the earthquake focal mechanisms present in both areas, and what our experiments predicted:

Figure 12. Top: the Woodlark Basin, PNG. Bottom: the Galapagos Rise. Both show earthquake focal mechanisms, coloured the same way as our experiments: blue for extension, red for compression, and green for strike-slip.

Furthermore, much work has been done investigating the Hess deep, a depression that sits ahead of the rift tip in the Galapagos. We found in our rotational experiment a similar ‘deep’ that moves ahead of the rift tip through time, giving us greater confidence in our experimental predictions.


There are a few things I’ve taken away from this experience. The first is that it’s important to remember the fundamentals. I’ve found that, generally, geodynamicists initially think about the force-balances going on in a particular setting, but gravity was staring me in the face for a while before I understood its critical role.

The second take-away was that exploratory modelling – playing around with experiments just for fun – is a great thing to do. Probably most of us do this anyway as part of the day-to-day activities, but putting aside some time to think about what sort of things to try out allowed us to find something really interesting. Furthermore, we then had a whole host of predictions we could go out and look for, rather than trying to tweak out experiment parameters to match something we already had found.

Finally, the 3D revolution we’re going through at the moment is exciting! Now that there are computers available to us that are able to run these enormous calculations, it gives us a chance to explore these fundamental problems in a new way and hopefully learn something about the world!

If you would like to checkout our paper, you can see it here. We made all of our input files open-source (and the code Underworld is already open-source), so please check them out too!


Brune, S., Popov, A. A., & Sobolev, S. V. (2012). Modeling suggests that oblique extension facilitates rifting and continental break‐up. Journal of Geophysical Research: Solid Earth, 117(B8).

Brune, S., Heine, C., Pérez-Gussinyé, M., & Sobolev, S. V. (2014). Rift migration explains continental margin asymmetry and crustal hyper-extension. Nature Communications, 5, 4014.

Brune, S. (2018). Forces within continental and oceanic rifts: Numerical modeling elucidates the impact of asthenospheric flow on surface stress. Geology, 46(2), 191-192.

Buck, W. R., Lavier, L. L., & Poliakov, A. N. (1999). How to make a rift wide. Philosophical Transactions - Royal Society of London Series A, Mathematical Physical and Engineering Sciences, 671-689.

Gerya, T. (2009). Introduction to numerical geodynamic modelling. Cambridge University Press.

Kearey, P. (Ed.). (2009). The Encyclopedia of the solid earth sciences. John Wiley & Sons.

Lundin, E. R., Redfield, T. F., Péron-Pindivic, G., & Pindell, J. (2014, January). Rifted continental margins: geometric influence on crustal architecture and melting. In Sedimentary Basins: Origin, Depositional Histories, and Petroleum Systems. 33rd Annual GCSSEPM Foundation Bob F. Perkins Conference. Gulf Coast Section SEPM (GCSSEPM), Houston, TX (pp. 18-53).

Koopmann, H., Brune, S., Franke, D., & Breuer, S. (2014). Linking rift propagation barriers to excess magmatism at volcanic rifted margins. Geology, 42(12), 1071-1074.

Lavier, L. L., & Manatschal, G. (2006). A mechanism to thin the continental lithosphere at magma-poor margins. Nature, 440(7082), 324.

Mondy, L. S., Rey, P. F., Duclaux, G., & Moresi, L. (2018). The role of asthenospheric flow during rift propagation and breakup. Geology.

Moresi, L., Betts, P. G., Miller, M. S., & Cayley, R. A. (2014). Dynamics of continental accretion. Nature, 508(7495), 245.

Rey, P. F., Coltice, N., & Flament, N. (2014). Spreading continents kick-started plate tectonics. Nature, 513(7518), 405.

Rey, P. F., Mondy, L., Duclaux, G., Teyssier, C., Whitney, D. L., Bocher, M., & Prigent, C. (2017). The origin of contractional structures in extensional gneiss domes. Geology, 45(3), 263-266.

Zoback, M. L. (1992). First‐and second‐order patterns of stress in the lithosphere: The World Stress Map Project. Journal of Geophysical Research: Solid Earth, 97(B8), 11703-11728.

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.