One month to AGU!

One month to AGU!
As the leaves are falling; the sun is going down before you leave the office; and the stores are stacking up on Christmas decorations, it’s time to face the facts: it’s almost AGU! It shouldn’t come as a surprise, but just in case. Don’t worry, there is still time to reread your abstract to see what you’re supposed to be presenting, figure out how to do that in the several weeks that are left and wrap it all up in a convincing poster or talk.
But first, check out these tips on presenting your work: our EGU GD blogs on creating prize-winning posters and selecting proper color schemes, EGU resources on making posters and these Science magazine posts on creating and giving great oral presentations. And while you’re at it, research New Orleans  highlights, such as the French quarter and the excellent cuisine. There is still time to register for some of the AGU field trips too! And wait, did you book your hotel and flights? Then, get to work!

The jelly sandwich lithosphere: elastic bread, the jelly, and gummy bears

The jelly sandwich lithosphere: elastic bread, the jelly, and gummy bears

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 Vojtěch Patočka from the Charles University in Prague, Czech Republic, discusses the rheology of the lithosphere and its food analogies. Do you want to talk about your research? Contact us!

It is becoming increasingly obvious that geodynamics and cooking are closely related, especially to the participants of a recent workshop in the Netherlands. Food analogies are helping students to get a physical grasp of continuum mechanics and to forget their lunches at university canteens. Here we take a basic look at a field where talking about fine cuisine has long been established: the rheology of the lithosphere. But first we step back a little.

Elastic solids deform when a force is applied and return to their original shape when the force is removed. Viscous fluids eventually take the shape of a container that is applied to them: they spontaneously flow to a state of zero shear stress. One can hardly imagine more different materials than these two, and yet the Earth’s mantle is sometimes modelled as an elastic shell and sometimes considered to be a viscous fluid.

You may be thinking:

Sure, it is a matter of what time scale is relevant for the process at hand

– just as J.C. Maxwell already thought: “Hence, a block of pitch may be so hard that you cannot make a dent in it by striking it with your knuckles; and yet it will in the course of time flatten itself by its weight, and glide downhill like a stream of water.” To visualise Maxwell’s dreams, the University of Queensland has been continuously running a pitch drop experiment for the past 90 years.

It is for the reason above that seismologists vibrate a Hookean Earth in their computers and geodynamicists play with viscous fluids. Are both of them right? Surprisingly, only the seismologists are. There is direct evidence that the outer parts of planetary mantles have an important elastic component even on geological time scales, which implies that treating the entire mantle as viscous is wrong. A textbook example can be found near the deepest part of the world’s oceans. Fig. 1, adopted from the bible by Turcotte & Schubert, shows how one can fit the bathymetric profile across the Mariana trench to the shape of a bent elastic plate. Note that the slab is subducting at the rate of a few centimeters per year, meaning that each segment of the slab is loaded for tens of millions of years before it disappears into the mantle and it still retains the elastic strain.

Figure 1: Comparison of a bathymetric profile across the Mariana trench (solid line) with the deflection of a thin elastic plate subject to end loading (dashed line). Distance xbx0 is the half-width of the forebulge from which the thickness of the plate can be inferred. Adopted from section 3 of Turcotte & Schubert, 2002.

Flexure studies and effective elastic thickness

Other similar examples, referred to as flexure studies, are summarised nicely by Watts et al., 2013. They include deflections under seamounts that litter oceanic floors and also some much longer lasting structures, such as continental foreland basins with the Ganges basin representing a particularly popular case. The free parameter that plays first fiddle in flexure studies is the thickness of elastic plate that matches observation. What does the resulting value, known as the effective elastic thickness, actually tell us?

I will use a dirty trick of bad journalism and quote an innocent geology blog slightly out of context:

The Indian crust is cold and rigid. Clever folk can do the maths on the shape of the crust as it bends down. This confirms that the pattern matches the model for rigid, elastic deformation. It also allows us to calculate the plate’s flexural rigidity, which is a measure of its strength. This means quantifying the rheology of real bits of the earth, which is a very useful trick

The math referred to involves purely elastic deformation and the flexural rigidity in its definition depends only on Young’s modulus, thickness, and Poisson ratio of the plate. Are the clever folk trying to fool us into believing that the lithosphere is an elastic plate? From lab measurements and geology we know for sure that it is not. Brittle failures are an abundant feature in both the crust and upper mantle, and various solid state creep mechanisms must be active in the deeper parts of the lithosphere.

The best fit for the Indo-Australian plate subducted below the Ganges basin is obtained with an elastic plate of circa 90 km thickness and for the Pacific plate at the Mariana trench it is circa 30 km. In both cases the plates are actually much thicker. What the computed values of 90 and 30 km tell us, is how much of elastic energy is present within the plate over the process of subduction, despite the brittle failures and despite the ductile creep activated in response to shear stresses within the plate.

To move away from sinking slabs, think about the story of a seamount that popped out on an ocean floor. Imagine a fresh, unstressed segment of oceanic plate that gets suddenly loaded by the uninvited underwater volcano. Small intra-plate fractures may immediately form, depending on the size of the load, and ductile creep will begin to continuously deform the plate’s deeper parts. The elastic energy present in the material upon loading gets partially released, either immediately through cracking or gradually via ductile creep. Measuring effective (or equivalent) elastic thickness merely tells us the amount of elastic energy left at the time of measurement.

This is all well known to the authors of flexure studies. In fact, they often re-draw the Christmas trees (see Fig. 2) and provide constraints on rheological zonation of the Earth. Pioneers in the field are E.B. Burov and A.B. Watts, whose papers are often accidentaly googled by chefs searching for latest trends in the dessert industry due to the extensive use of the words jelly sandwich and crème brûlée ([4] and [5], for breakfast see also [6]). The main point of the discussion is to determine the brittle/ductile transition and the active creep mechanism in a realistic, compositionally stratified lithosphere. This is usually complicated by the fact that lithologies measured in lab experiments strongly depend on composition, water content and temperature – and these are not well constrained in the real Earth.

Figure 2: The total force per unit width necessary to break or viscously deform a lithospheric section at a given strain rate. Plots like these are known as the Christmas tree plots, here adopted from Basin Analysis by P.A. Allen and J.R. Allen (without the food).

A geodynamical paradox

It is a paradox that in geodynamical modeling we often use the constraints from flexure studies and at the same time forget about elasticity, without which there would be no such studies at all. Recall that the primary result of shape fitting, for instance of the one depicted in Fig. 1, is that some elastic support is present. As I warned above, one can only hardly distinguish between a purely elastic plate of a given thickness and some other, thicker elasto-brittle or visco-elastic plate when looking at its surface flexure. However, there must be elasticity involved in some form: purely viscous or visco-plastic plates would not form the observed forebulge. Forebulges are related to the way elastic rods and plates transfer bending moment throughout the medium.

The tendency to disregard elasticity may be related to the way we use the word ‘rigid’ in the context of plate tectonics. In physics, rigid means ‘not deforming’. For plate tectonics to work the plates do not have to be rigid. They may not flow apart on geological time scales (they must remain plates), but some relatively small reversible deformation is well compliant with the concept. So the main point is: Be the lithosphere a sandwich or a crème brûlée, it is also flexible, even on geological time scales.

Batman and gummy bears

If you are in numerical modelling then there is good news for you. In the past two decades, several Robins, including myself, have enhanced the Batman codes to account for elasticity. Maybe all you have to do is to switch on the right button. The implementation of visco-elasticity is usually based on a method developed by L. Moresi, who also, according to a previous Geodynamics 101 blog post, coined the hero terminology I just borrowed. Visco-elasticity works in quite a simple way, just like gummy bears. I am currently running an experiment with them. They quickly bend and squeeze when tortured, their resistance governed by their shear modulus, and recover when let go. At the same time they can flow, the resistance being controlled by their viscosity. Let’s put a book on top of one and see if we can find its nose after a week…

The elasticity button

Don’t be afraid to push the elasticity button, if your code has it. Usually it won’t do anything dramatic to your simulations, but exceptions exist [7]. In thermal convection models without plate tectonics there is not much feedback between the lid and the underlying mantle, and so only the build-up of dynamic topography is affected [8]. In subduction modeling your slabs may obtain different dipping angles. In continental extensions the total amount of extension will become more important than the divergence rate [9]. And in the simulation of continental shortening mentioned above [7], the elastic energy accumulated in the entire model gets partially released upon the onset of a shear zone. In such cases, i.e. when large scale elastic strains suddenly influence a much smaller region, one can expect some earthquakes to shake the conventional view of elasticity in geodynamical modeling. And if you still do not care about elasticity but yet you made it all the way here, then you deserve a bonus: the convergence of Stokes solvers is way better for visco-elastic rheologies than for the viscous ones – if you are numerically troubled by large viscosity contrasts of your model, elasticity is the way to go for you.

(1) D.L. Turcotte and G. Schubert (2002), Geodynamics
(2) A.B. Watts, S.J. Zhong, and J. Hunter (2013), The Behavior of the Lithosphere on Seismic to Geologic Timescales, doi: 10.1146/annurev-earth-042711-105457
(3) P.A. Allen and J.R. Allen (2005), Basin analysis : principles and applications
(4) E.B. Burov and A.B. Watts (2006), The long-term strength of continental lithosphere: 'jelly sandwich' or 'creme brulée'?, doi: 10.1130/1052-5173(2006)016
(5) E.B. Burov (2009), Time to burn out creme brulee?, doi: 10.1016/j.tecto.2009.06.013
(6) E.H. Hartz, Y.Y. Podladchikov (2008), Toasting the jelly sandwich: The effect of shear heating on lithospheric geotherms and strength, doi: 10.1130/G24424A.1
(7) Y. Jaquet, T. Duretz, and S.M. Schmalholz (2016), Dramatic effect of elasticity on thermal softening and strain localization during lithospheric shortening, doi: 10.1093/gji/ggv464
(8) V. Patocka, O. Cadek, P.J. Tackley, and H. Cizkova (2017), Stress memory effect in viscoelastic stagnant lid convection, doi: 10.1093/gji/ggx102
(9) J.A. Olive, M.D. Behn, E. Mittelstaedt, G. Ito, and B.Z. Klein, The role of elasticity in simulating long-term tectonic extension, doi: 10.1093/gji/ggw044

Going with the toroidal mantle flow

Going with the toroidal mantle flow

Subduction zones host one of the most complex and fascinating tectonic systems on the planet. Numerical models by Király and colleagues recently published in Earth and Planetary Science Letters reveal that the strength of the toroidal flow depends on the mantle viscosity and the magnitude of the slab pull force while the characteristic size of the toroidal cells mainly depends on the size of the convecting mantle.


The motion of Earth’s tectonic plates—the lithosphere—is driven by the subduction of relatively cold and dense oceanic plates into the mantle. Subduction zones are some of the most striking features on Earth. They represent one of the two types of convergent plate boundaries, in which one tectonic plate sinks underneath another one into the Earth’s mantle. The resulting forces induce mantle flow around the subducting plate, but the manner in which this happens is still a matter of debate.

Below the tectonic plates, the mantle moves in a slow and unseen manner—as cold slabs sink, hot upwellings rise, and convection slowly rids the Earth’s interior of its primordial heat. Much of our knowledge of mantle flow patterns is indirect, inferred from geodynamical models or seismic tomography images.

Seismological and geochemical investigations suggest 3D subduction-induced mantle flow around lateral slab edges from the sub-slab zone towards the mantle wedge[1]. Such flow has also been observed in laboratory experiments[2] and numerical models of subduction[3]. In particular, geodynamic numerical models of subduction demonstrate that back-arc extension at narrow subduction zones is driven by rollback-induced mantle flow[4].

Observations of seismic anisotropy—the directional dependence of seismic wave speed—provide us with tantalizingly direct information about mantle flow direction. For example, crystals of olivine, which is the most common upper-mantle mineral, tend to become aligned by the mantle flow, and seismic anisotropy is an indicator of this alignment.

To apply the laboratory-derived viscosity laws to nature, they must be extrapolated over 10 orders of magnitude, which introduces uncertainty. Moreover, it is unclear to what extent experiments on centimeter-scale samples are representative of the crust and lithosphere.

The viscosity distribution of the lithosphere and the surrounding mantle is therefore one of the least certain parameters in geodynamics. Alternative ways to determine viscosity on geological time scales are thus needed.

Schematic diagram of a subduction zone, showing the dominance of 3D flow beneath the slab and the competing influence of 2D and 3D flow fields in the mantle wedge. Credit: Long and Silver, 2008, Science.


Writing in Earth and Planetary Science Letters, Király and colleagues[5] present a cutting-edge numerical model that shows how the strength and length scale of the toroidal flow vary with the mantle viscosity and the magnitude of the slab pull force. Király and co-workers highlight some remarkable implications of these effects: around subducting plates, the characteristic length, axis, and shape of the toroidal cell are almost independent of the slab’s properties and mainly depend on the thickness of the convecting mantle. The independence of the shape of the toroidal cell on the slab width can appear controversial with respect to previous studies[6–7] that showed that the overall mantle flow is dependent on the slab width. However, as Király et al. point out, these differences can be explained with the model setup adopted and the range of slab widths investigated.

In order to characterize the flow—in terms of its geometry and strength—Király and co-workers analysed the vertical component of the mantle vorticity as well as the ratio between the vertical and the trench parallel component of the vorticity vector. With a series of numerical experiments, they find that subduction-induced mantle flow is highly three-dimensional, and that the toroidal component is sub-horizontal with some vertical flow components. According to the authors, this vertical flow around the slab edges has significant implications, as it can be responsible for the presence of the off-arc volcanism around several laterally confined subduction zones.

The numerical experiments by Király et al. represent another step forward in our understanding of how mantle circulation plays a relevant role in the shaping of tectonic features around subduction zones, and they provide evidence that the slab properties only impact the vigour of the flow around subducting slabs. The extension of these numerical experiments to a range of different parameters will allow for a fuller characterization of the olivine fabric, which in turn will allow seismologists to relate their measurements of seismic anisotropy to flow directions and, ultimately, to mantle processes.



1. Long, M. D., and P. G. Silver (2008), The subduction zone flow field from seismic anisotropy: A global view, Science, 319, 315-318.
2. Funiciello, F., C. Faccenna, and D. Giardini (2004), Role of lateral mantle flow in the evolution of subduction systems: insights from laboratory experiments, Geophy. J. Int., 157, 1393-1406.
3. Jadamec, M. A., and M. I. Billen (2010), Reconciling surface plate motions with rapid three-dimensional mantle flow around a slab edge, Nature, 465, 338-342.
4. Sternai, P., L. Jolivet, A. Menant, and T. Gerya (2014), Driving the upper plate surface deformation by slab rollback and mantle flow, Earth Planet. Sci. Lett., 405, 110-118.
5. Király, A. Capitanio, F.A., Funiciello, F., Faccenna, C. (2017). Subduction induced mantle flow: Length-scales and orientation of the toroidal cell. Earth Planet. Sci. Lett., 479, 284–297.
6. Funiciello, F., Moroni, M., Piromallo, C., Faccenna, C., Cenedese, A., Bui, H.A., 2006. Mapping mantle flow during retreating subduction: laboratory models analyzed by feature tracking. J. Geophys. Res., Solid Earth 111, 1–16.
7. Stegman, D.R., Freeman, J., Schellart, W.P., Moresi, L., May, D., 2006. Influence of trench width on subduction hinge retreat rates in 3-D models of slab rollback. Geochem. Geophys. Geosyst. 7.

Planting seeds of deformation in numerical models

Planting seeds of deformation in numerical models

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 we continue the conversation that was started at NetherMod 2017 by discussing how we initiate deformation in numerical models of the lithosphere and more importantly, does it matter how we start such models? Do you want to talk about your research? Contact us!

Geodynamic modelling often concerns itself with the study of localized deformation in the crust and lithosphere. The models generally start with an initial geometry, boundary conditions and a prescribed set of initial conditions that can be either thermal or mechanical. This initial setup is usually a standard representation of the lithosphere and asthenosphere, as inferred from geological and geophysical observations and laboratory rock experiments. The boundary conditions usually drive the deformation in the system. However, as a synthetic (computer) model is pristine at the start, the localization of deformation can take a long model time of up to millions of years, as numerical disturbances need to accumulate. Besides that, the deformation will likely localise near or at the boundary of the system because of the boundary conditions. In order to avoid this long starting phase, and to exert some control on the location of the initial deformation, several approaches are widely used to initiate and localize deformation. Examples of these are the S-point velocity discontinuity at the bottom of the system (e.g. Braun and Beaumont, 1995; Ellis et al., 1995; Willett, 1999; Beaumont et al., 2000; Buiter et al., 2006; Thieulot et al., 2008; Braun and Yamato, 2009) and the use of a weak seed. The latter are usually small zones that are weaker than the surrounding crust or lithosphere. As the lithosphere and crust are never homogeneous, using weak seeds in a model can be easily justified. They could represent regions of different material properties (e.g., heat production), inherited faults, inherited crustal thickness changes and/or plumes impacting the lithosphere. The popular NetherMod 2017 potatoes (yep, still referring back to that one. If you don’t know what this is about, check out this post) are an extreme case in which multiple weak seeds are used to reflect local geology. Traditionally, simpler, single weak seeds are used.

Here, I will give a brief (and by no means comprehensive) overview of the different weak seed methods used to initiate deformation in models of continental extension, and I will conclude with a discussion on how these different initial conditions affect the model evolution.

Modified figure from Burg & Podladchikov (1999): a thermal perturbation is used in the middle of their model to localise deformation there.

Seeding through thermal effects

The weak zone could be implemented as a temperature anomaly which could reflect a region in the crust of higher radiogenic heat production or a region in the lithosphere of locally reduced viscosity. They can be implemented by elevating the temperature or basal heat flux at the crust-lithosphere boundary or lithosphere-asthenosphere boundary with a certain amplitude over a finite region. Examples of studies using thermal weak seeds are Burg and Podladchikov (1999), Frederiksen and Braun (2001), Hansen and Nielsen (2003), and Brune and Autin (2013).

Seeding by mechanical inhomogeneity

Modified figure from Kaus (2009): Model setup with a rectangular, viscous inclusion in the middle of the domain to initiate and localise deformation. Dimensions of the weak zone were varied, but maintained the aspect ratio 2:1.

A weak seed could be composed of a material with a lower rheological strength than its surroundings. There are multiple ways of achieving this, but often used methods include:
• A weak seed with a lower viscosity (e.g., Gray and Pysklywec, 2010; Gray and Pysklywec, 2012b; Gray and Pysklywec, 2012a; Kaus, 2009, and Mishin, 2011)
• A weak seed with a different angle of internal friction (e.g., Pysklywec et al., 2002; Kaus and Podlachikov, 2006; Thieulot, 2011; Gray and Pysklywec, 2013; Chenin and Beaumont, 2013)
• A weak seed with a different density (e.g., Tirel et al., 2008)
• A weak seed consisting of a different material (e.g., Pysklywec et al., 2000; Huismans and Beaumont, 2007)
• A weak seed with more accumulated strain than its surroundings (e.g., Lavier et al., 2000; Huismans et al., 2005; Warren et al., 2008a; Warren et al., 2008b; Petrunin and Sobolev, 2008; Beaumont et al., 2009; Allken et al., 2011; Kneller et al., 2013; Allken et al., 2013)

Apart from these different ways of making a seed weak, you can also find weak seeds of many different shapes and sizes in the literature. Most commonly, you will find
• Square weak seeds (e.g., Gray and Pysklywec, 2013)
• Rectangular weak seeds (e.g., Huismans et al., 2005) with different aspect ratios
• Fault-shaped weak seeds (e.g., Currie et al., 2007, and Currie and Beaumont, 2011)

Modified from Gray & Pysklywec, 2013: Model setup with a square, frictionally weak zone of dimensions 10×10 km.

Seeding through geometrical discontinuity
Another method to create a zone of different rheological strength is by varying the thickness of the crust or lithosphere. A thickened lithosphere could represent a remnant of a previous mountain building phase, whereas a thinned lithosphere would represent a remnant of a previous rifting phase. Studies using this method of weakening include Gac et al. (2014) and Burg and Schmalholz (2008).

Modified from Gac et al. (2014): Model setup with an inherited thin crust.

Influence of weak seeds on the model evolution

As mentioned at NetherMod 2017, you would ideally either have a generic model for which you should determine the influence of different initial weak seeds to check how robust your model is, or you would have a region-specific model for which you find the optimal initial conditions to get your desired model output. Only few studies have investigated the former in detail. Dyksterhuis et al. (2013) found that a single weak seed typically produces symmetric narrow rifts; multiple seeds produce a wide rift; and an initial fault-shaped weak zone produces an asymmetric rift.

It also raises the question about whether or not there are differences between similar codes when the same initiation method is used.

To shine a very preliminary light on this problem, I ran some models of continental extension using the SULEC and ELEFANT codes with different initial conditions. As both codes are based on the same physics, and have similar implementation, they should show a high degree of similarity when using the same deformation initiation method. The results show that different initiation methods indeed result in different results, particularly with respect to the timing of the deformation (see figure below). Besides that, SULEC tends to show more asymmetric behaviour than ELEFANT. For a more complete overview of the results and the model setup, please look here (my old (first!) poster for GeoMod 2014).

Models of continental extension after 10 Myr of extension for SULEC (top) and ELEFANT (bottom) for different initial, frictionally weak, weak seeds with aspect ratios of 6×6 elements, 12×3 elements, and 3×12 elements (i.e., the weak seed consists of the same amount of elements in each model).

In conclusion, I hope to add to the discussion of ‘how we start our models’ with this post by affirming that the model evolution is affected by our choice of weak seed (if only by the amount of waiting until deformation starts) and its effect can differ slightly between codes, even if the codes are very similar. Taking into account the vast variability of methods to initiate deformation, one really should be careful when assessing model results.

Allken, V., Huismans, R., Fossen, H., and Thieulot, C. (2013). 3D numerical modelling of graben interaction and linkage: a case study of the Canyonlands grabens, Utah. Basin Research.

Allken, V., Huismans, R., and Thieulot, C. (2011). Three-dimensional numerical modeling of upper crustal extensional systems. Journal of Geophysical Research, page doi:10.1029/2011JB008319.

Beaumont, C., Jamieson, R., Butler, J., and Warren, C. (2009). Crustal structure: A key constraint on the mechanism of ultra-high-pressure rock exhumation. EPSL, 287:116-129.

Beaumont, C., Munoz, J., Hamilton, J., and Fullsack, P. (2000). Factors controlling the alpine evolution of the central pyrenees inferred from a comparison of observations and geodynamical models. Journal of Geophysical Research, 105:8121-8145.

Braun, J. and Beaumont, C. (1995). Three-dimensional numerical experiments of strain partitioning at oblique plate boundaries: Implications for contrasting tectonic styles in the southern Coast Ranges, California, and central South Island, New Zealand. Journal of Geophysical Research, 100(B9):18,059-18,074.

Braun, J. and Yamato, P. (2009). Structural evolution of a three-dimensional, finite-width crustal wedge.Tectonophysics, 484:181-192.

Brune, S. and Autin, J. (2013). The rift to break-up evolution of the Gulf of Aden: Insights from 3D numerical lithospheric-scale modelling. Tectonophysics, 607(0):65-79. The Gulf of Aden rifted margins system: Special Issue dedicated to the YOCMAL project (Young Conjugate Margins Laboratory in the Gulf of Aden).

Buiter, S., Babeyko, A., Ellis, S., Gerya, T., Kaus, B., Kellner, A., Schreurs, G., and Yamada, Y. (2006). The numerical sandbox: comparison of model results for a shortening and an extension experiment. Analogue and Numerical Modelling of Crustal-Scale Processes. Geological Society, London. Special Publications, 253:29-64.

Burg, J.-P. and Podladchikov, Y. (1999). Lithospheric scale folding: numerical modelling and application to the Himalayan syntaxes. International Journal of Earth Sciences, 88(2):190-200.

Burg, J.-P. and Schmalholz, S. (2008). Viscous heating allows thrusting to overcome crustal-scale buckling: Numerical investigation with application to the Himalayan syntaxes. Earth and Planetary Science Letters, 274(1):189-203.

Chenin, P. and Beaumont, C. (2013). Influence of offset weak zones on the development of rift basins: Activation and abandonment during continental extension and breakup. JGR, 118:1-23.

Currie, C. and Beaumont, C. (2011). Are diamond-bearing Cretaceous kimberlites related to low-angle subduction beneath western North America. EPSL, 303:59-70.

Currie, C., Beaumont, C., and Huismans, R. (2007). The fate of subducted sediments: a case for backarc intrusion and underplating. Geology, 35(12):1111-1114.

Dyksterhuis, S., Rey, P., Mueller, R., and Moresi, L. (2013). Effects of initial weakness on rift architecture. Geological Society, London, Special Publications, 282:443-455.

Ellis, S., Fullsack, P., and Beaumont, C. (1995). Oblique convergence of the crust driven by basal forcing: implications for length-scales of deformation and strain partitioning in orogens. Geophys. J. Int., 120:24-44.

Frederiksen, S. and Braun, J. (2001). Numerical modelling of strain localisation during extension of the continental lithosphere. Earth and Planetary Science Letters, 188(1):241-251.

Gac, S., Huismans, R. S., Simon, N. S., Faleide, J. I., and Podladchikov, Y. Y. (2014). Effects of lithosphere buckling on subsidence and hydrocarbon maturation: A case-study from the ultra-deep East Barents Sea basin. Earth and Planetary Science Letters, 407:123-133.

Gray, R. and Pysklywec, R. N. (2010). Geodynamic models of archean continental collision and the formation of mantle lithosphere keels. Geophysical Research Letters, 37(19).

Gray, R. and Pysklywec, R. (2012a). Geodynamic models of mature continental collision: Evolution of an orogen from lithospheric subduction to continental retreat/delamination. JGR, 117(B03408).

Gray, R. and Pysklywec, R. N. (2012b). Influence of sediment deposition on deep lithospheric tectonics. Geophysical Research Letters, 39(11).

Gray, R. and Pysklywec, R. (2013). Influence of viscosity pressure dependence on deep lithospheric tectonics during continental collision. JGR, 118.

Hansen, D. and Nielsen, S. (2003). Why rifts invert in compression. Tectonophysics, 373(1):5-24.

Huismans, R. and Beaumont, C. (2007). Roles of lithospheric strain softening and heterogeneity in determining the geometry of rifts and continental margins. Geological Society, London, Special Publications, 282(1):111-138.

Huismans, R., Buiter, S., and Beaumont, C. (2005). Effect of plastic-viscous layering and strain softening on mode selection during lithospheric extension. Journal of Geophysical Research, 110:B02406.

Kaus, B. (2009). Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation. Tectonophysics, 484(1), 36-47.

Kaus, B. and Podlachikov, Y. (2006). Initiation of localized shear zones in viscoelastoplastic rocks. Journal of Geophysical Research: Solid Earth 111.B4.

Kneller, E. A., Albertz, M., Karner, G. D., , and Johnson, C. A. (2013). Testing inverse kinematic models of paleocrustal thickness in extensional systems with high- resolution forward thermo-mechanical models.

Lavier, L., Buck, W., and Poliakov, A. (2000). Factors controlling normal fault offset in an ideal brittle layer. 105(B10):23,431--23,442.

Mishin, Y. (2011). Adaptive multiresolution methods for problems of computational geodynamics. PhD thesis, ETH Zurich.

Petrunin, A. and Sobolev, S. (2008). Three-dimensional numerical models of the evolution of pull-apart basins. Physics of the Earth and Planetary Interiors, 171:387-399.

Pysklywec, R., Beaumont, C., and Fullsack, P. (2000). Modeling the behavior of continental mantle lithosphere during plate convergence. Geology, 28(7):655-658.

Pysklywec, R., Beaumont, C., and Fullsack, P. (2002). Lithospheric deformation during the early stages of continental collision: Numerical experiments and comparison with South Island, New Zealand. JGR, 107(B72133).

Thieulot, C., Fullsack, P., and Braun, J. (2008). Adaptive octree-based finite element analysis of two- and three-dimensional indentation problems. Journal of Geophysical Research, 113:B12207.

Thieulot, C. (2011). FANTOM: two-and three-dimensional numerical modelling of creeping
 flows for the solution of geological problems. Physics of the Earth and Planetary Interiors, 188(1):47-68.

Tirel, C., Brun, J.-P., and Burov, E. (2008). Dynamics and structural development of metamorphic core complexes. Journal of Geophysical Research: Solid Earth (1978-2012), 113(B4).

Warren, C., Beaumont, C., and Jamieson, R. (2008a). Formation and exhumation of ultra-highpressure rocks during continental collision: Role of detachment in the subduction channel. Gcubed.

Warren, C., Beaumont, C., and Jamieson, R. (2008b). Modelling tectonic styles and ultra-high pressure (UHP) rock exhumation during the transition from oceanic subduction to continental collision. EPSL, 267:129-145.

Willett, S. D. (1999). Rheological dependence of extension in wedge models of convergent orogens. Tectonophysics, 305:419--435.