GD
Geodynamics

Geodynamics 101

The past is the key

The past is the key

Lorenzo Colli

“The present is the key to the past” is a oft-used phrase in the context of understanding our planet’s complex evolution. But this perspective can also be flipped, reflected, and reframed. In this Geodynamics 101 post, Lorenzo Colli, Research Assistant Professor at the University of Houston, USA, showcases some of the recent advances in modelling mantle convection.  

 

Mantle convection is the fundamental process that drives a large part of the geologic activity at the Earth’s surface. Indeed, mantle convection can be framed as a dynamical theory that complements and expands the kinematic theory of plate tectonics: on the one hand it aims to describe and quantify the forces that cause tectonic processes; on the other, it provides an explanation for features – such as hotspot volcanism, chains of seamounts, large igneous provinces and anomalous non-isostatic topography – that aren’t accounted for by plate tectonics.

Mantle convection is both very simple and very complicated. In its essence, it is simply thermal convection: hot (and lighter) material goes up, cold (and denser) material goes down. We can describe thermal convection using classical equations of fluid dynamics, which are based on well-founded physical principles: the continuity equation enforces conservation of mass; the Navier-Stokes equation deals with conservation of momentum; and the heat equation embodies conservation of energy. Moreover, given the extremely large viscosity of the Earth’s mantle and the low rates of deformation, inertia and turbulence are utterly negligible and the Navier-Stokes equation can be simplified accordingly. One incredible consequence is that the flow field only depends on an instantaneous force balance, not on its past states, and it is thus time reversible. And when I say incredible, I really mean it: it looks like a magic trick. Check it out yourself.

With four parameters I can fit an elephant, and with five I can make him wiggle his trunk

This is as simple as it gets, in the sense that from here onward every additional aspect of mantle convection results in a more complex system: 3D variations in rheology and composition; phase transitions, melting and, more generally, the thermodynamics of mantle minerals; the feedbacks between deep Earth dynamics and surface processes. Each of these additional aspects results in a system that is harder and costlier to solve numerically, so much so that numerical models need to compromise, including some but excluding others, or giving up dimensionality, domain size or the ability to advance in time. More importantly, most of these aspects are so-called subgrid-scale processes: they deal with the macroscopic effect of some microscopic process that cannot be modelled at the same scale as the macroscopic flow and is too costly to model at the appropriate scale. Consequently, it needs to be parametrized. To make matters worse, some of these microscopic processes are not understood sufficiently well to begin with: the parametrizations are not formally derived from first-principle physics but are long-range extrapolations of semi-empirical laws. The end result is that it is possible to generate more complex – thus, in this regard, more Earth-like – models of mantle convection at the cost of an increase in tunable parameters. But what parameters give a truly better model? How can we test it?

Figure 1: The mantle convection model on the left runs in ten minutes on your laptop. It is not the Earth. The one on the right takes two days on a supercomputer. It is fancier, but it is still not the real Earth.

Meteorologists face similar issues with their models of atmospheric circulation. For example, processes related to turbulence, clouds and rainfall need to be parametrized. Early weather forecast models were… less than ideal. But meteorologists can compare every day their model predictions with what actually occurs, thus objectively and quantitatively assessing what works and what doesn’t. As a result, during the last 40 years weather predictions have improved steadily (Bauer et al., 2015). Current models are better at using available information (what is technically called data assimilation; more on this later) and have parametrizations that better represent the physics of the underlying processes.

If time travel is possible, where are the geophysicists from the future?

We could do the same, in theory. We can initialize a mantle convection model with some best estimate for the present-day state of the Earth’s mantle and let it run forward into the future, with the explicit aim of forecasting its future evolution. But mantle convection evolves over millions of years instead of days, thus making future predictions impractical. Another option would be to initialize a mantle convection model in the distant past and run it forward, thus making predictions-in-the-past. But in this case we really don’t know the state of the mantle in the past. And as mantle convection is a chaotic process, even a small error in the initial condition quickly grows into a completely different model trajectory (Bello et al., 2014). One can mitigate this chaotic divergence by using data assimilation and imposing surface velocities as reconstructed by a kinematic model of past plate motions (Bunge et al., 1998), which indeed tends to bring the modelled evolution closer to the true one (Colli et al., 2015). But it would take hundreds of millions of years of error-free plate motions to eliminate the influence of the unknown initial condition.

As I mentioned before, the flow field is time reversible, so one can try to start from the present-day state and integrate the governing equations backward in time. But while the flow field is time reversible, the temperature field is not. Heat diffusion is physically irreversible and mathematically unstable when solved back in time. Plainly said, the temperature field blows up. Heat diffusion needs to be turned off [1], thus keeping only heat advection. This approach, aptly called backward advection (Steinberger and O’Connell, 1997), is limited to only a few tens of millions of years in the past (Conrad and Gurnis, 2003; Moucha and Forte, 2011): the errors induced by neglecting heat diffusion add up and the recovered “initial condition”, when integrated forward in time (or should I say, back to the future), doesn’t land back at the desired present-day state, following instead a divergent trajectory.

Per aspera ad astra

As all the simple approaches turn out to be either unfeasible or unsatisfactory, we need to turn our attention to more sophisticated ones. One option is to be more clever about data assimilation, for example using a Kalman filter (Bocher et al., 2016; 2018). This methodology allow for the combining of the physics of the system, as embodied by the numerical model, with observational data, while at the same time taking into account their relative uncertainties. A different approach is given by posing a formal inverse problem aimed at finding the “optimal” initial condition that evolves into the known (best-estimate) present-day state of the mantle. This inverse problem can be solved using the adjoint method (Bunge et al., 2003; Liu and Gurnis, 2008), a rather elegant mathematical technique that exploits the physics of the system to compute the sensitivity of the final condition to variations in the initial condition. Both methodologies are computationally very expensive. Like, many millions of CPU-hours expensive. But they allow for explicit predictions of the past history of mantle flow (Spasojevic & Gurnis, 2012; Colli et al., 2018), which can then be compared with evidence of past flow states as preserved by the geologic record, for example in the form of regional- and continental-scale unconformities (Friedrich et al., 2018) and planation surfaces (Guillocheau et al., 2018). The past history of the Earth thus holds the key to significantly advance our understanding of mantle dynamics by allowing us to test and improve our models of mantle convection.

Figure 2: A schematic illustration of a reconstruction of past mantle flow obtained via the adjoint method. Symbols represent model states at discrete times. They are connected by lines representing model evolution over time. The procedure starts from a first guess of the state of the mantle in the distant past (orange circle). When evolved in time (red triangles) it will not reproduce the present-day state of the real Earth (purple cross). The adjoint method tells you in which direction the initial condition needs to be shifted in order to move the modeled present-day state closer to the real Earth. By iteratively correcting the first guess an optimized evolution (green stars) can be obtained, which matches the present-day state of the Earth.

1.Or even to be reversed in sign, to make the time-reversed heat equation unconditionally stable.

How to make a subduction zone on Earth

How to make a subduction zone on Earth

Subduction zones are ubiquitous features on Earth, and an integral part of plate tectonics. They are known to have a very important role in modulating climate on Earth, and are believed to have played an essential part in making the Earth’s surface habitable, a role that extends to present-day. This week, Antoniette Greta Grima writes about the ongoing debate on how subduction zones form and persist for millions of years, consuming oceanic lithosphere and transporting water and other volatiles to the Earth’s mantle.

Antoniette Greta Grima. PhD Student at Dept. of Earth Sciences, University College London, UK.

Before we can start thinking about how subduction zones form, we need to be clear on what we mean by the term subduction zone. In the most generic sense, this term has been described by White et al. (1970) as “an abruptly descending or formerly descended elongate body of lithosphere, together with an existing envelope of plate deformation”. In simple words this defines subduction zones as places where pieces of the Earth’s lithosphere bend downwards into the Earth’s interior. This definition however, does not take into account the spatio-temporal aspect of subduction zone formation. It also, does not differentiate between temporary, episodic lithosphere ‘peeling’ or ‘drips’, thought to precede the modern-day ocean plate-tectonic regime (see van Hunen and Moyen, 2012; Crameri et al., 2018; Foley, 2018, and references therein) and the rigid self-sustaining subduction, which we see on the present-day Earth (Gurnis et al., 2004).

A self-sustaining subduction zone is one where the total buried, rigid slab length extends deep into the upper mantle and is accompanied at the surface by back-arc spreading (Gurnis et al., 2004). The latter is an important surface observable indicating that the slab has overcome the resistive forces impeding its subduction and is falling quasi-vertically through the mantle. Gurnis et al. (2004) go on to say that if one or the other of these defining criteria is missing then subduction is forced rather than self-sustaining. Forced or induced subduction (Stern, 2004; Leng and Gurnis, 2011; Stern and Gerya, 2017; Baes et al., 2018) is described by Gurnis et al. (2004) as a juvenile, early stage system, that cannot be described as a fully fledged subduction zone. These forced subduction zones are characterised by incipient margins, short trench-arc distance, narrow trenches and a volcanically inactive island arc and/or trench. Furthermore, although these juvenile systems might be seismically active they will lack a well defined Benioff Zone. Examples of forced subduction include the Puysegur-Fiordland subduction along the Macquarie Ridge Complex, the Mussau Trench on the eastern boundary of the Caroline plate and the Yap Trench south of the Marianas, amongst others. On the other hand, Cenozoic (<66 Ma) subductions, shown in figure 2, are by this definition self-sustaining and mature subduction zones. These subduction zones including the Izu-Bonin-Mariana, Tonga-Kermadec and Aleutians subduction zones, are characterised by their extensive and well defined trenches (see figure 2) (Gurnis et al., 2004). However, despite their common categorization subduction zones can originate through various mechanisms and from very different tectonic settings.

Figure 1: Map the Earth’s subduction zones and tectonic plates from W.K and E.H. (2003). See how subduction zones dominate the figure.

1. How are subduction zones formed?

We know from the geological record that the formation of subduction zones is an ongoing process, with nearly half of the present day active subduction zones initiating during the Cenozoic (<66 Ma) (see Gurnis et al., 2004; Dymkova and Gerya, 2013; Crameri et al., 2018, and references therein). However, it is less clear how subduction zones originate, nucleate and propagate to pristine oceanic basins.

 

Figure 2: The oldest, still active subduction zones on Earth can be dated back to 66 million years ago. These are self- sustaining mature subduction zones with well defined trenches and trench lengths (modified from Gurnis et al., 2004).

Crameri et al. (2018, and references therein) list a number of mechanisms, some which are shown in figure 3, that may work together to weaken and break the lithosphere including:

  • Meteorite impact
  • Sediment loading
  • Major episode of delamination
  • Small scale convection in the sub-lithospheric mantle
  • Interaction of thermo-chemical plume with the overlying lithosphere
  • Plate bending via surface topographic variations
  • Addition of water or melt to the lithosphere
  • Pre-existing transform fault or oceanic plateau
  • Shear heating
  • Grain size reduction

Some of these mechanisms, particularly those listed at the beginning of the list are more appropriate to early Earth conditions while others, such as inherited weaknesses or fracture zones, transform faults and extinct spreading ridges are considered to be prime tectonic settings for subduction zone formation in the Cenozoic (<66 Ma) (Gurnis et al., 2004). As the oceanic lithosphere grows denser with age, it develops heterogeneity which facilitates its sinking into the mantle to form new subduction zones. However, it is important to keep in mind that without inherited, pre-existing weaknesses, it is extremely difficult to form subduction zones at passive margins. This is because as the oceanic lithosphere cools and becomes denser, it also becomes stronger and therefore harder to bend into the mantle that underlies it (Gurnis et al., 2004; Duarte et al., 2016). Gurnis et al. (2004) note that the formation of new subduction zones alters the force balance on the plate and suggest that the strength of the lithosphere during bending is potentially the largest resisting component in the development of new subduction zones. Once that resistance to bending is overcome, either through the negative buoyancy of the subducting plate and/or through the tectonic forces acting on it, a shear zone extending through the plate develops (Gurnis et al., 2004; Leng and Gurnis, 2011). This eventually leads to plate failure and subduction zone formation.

Figure 3: Different ways to form a new subduction zone (from Stern and Gerya, 2017).

2. Where can new subduction zones form?

From our knowledge of the geological record, observations of on-going subduction, and numerical modelling (Baes et al., 2011; Leng and Gurnis, 2011; Baes et al., 2018; Beaussier et al., 2018) we think that subduction zone initiation primarily occurs through the following:

In an intra-oceanic setting through surface weakening processes

An intra-oceanic setting refers to a subduction zone forming right within the oceanic plate itself. Proposed weakening mechanisms include weakening of the lithosphere due to melt and/or hydration (e.g., Crameri et al., 2018; Foley, 2018, and references therein), localised lithospheric shear heating (Thielmann and Kaus, 2012) and density variations within oceanic plate due to age heterogeneities, where its older and denser portions flounder and sink (Duarte et al., 2016). Another mechanism proposed by Baes et al. (2018) suggests that intra-oceanic subduction can also be induced by mantle suction flow. These authors suggest that mantle suction flow stemming from either slab remnants and/or from slabs of active subduction zones can act on pre-existing zones of weakness, such as STEP (subduction-transfer edge propagate) faults to trigger a new subduction zone, thus facilitating spontaneous subduction initiation (e.g. figure 3) (Stern, 2004). The Sandwich and the Tonga-Kermadec subduction zones are often cited as prime examples of intra-oceanic subduction zone formation due to mantle suction forces (Baes et al., 2018). Ueda et al. (2008) and Gerya et al. (2015) also suggest that thermochemical plumes can break the lithosphere and initiate self-sustaining subduction, provided that the overlying lithosphere is weakened through the presence of volatiles and melt (e.g. figure 3). This mechanism can explain the Venusian corona and could have facilitated the initation of plate tectonics on Earth (Ueda et al., 2008; Gerya et al., 2015). Similarly Burov and Cloetingh (2010) suggest that in the absence of plate tectonics, mantle lithospheric interaction through plume-like instabilities, can induce spontaneous downwelling of both continental and oceanic lithosphere.

Through subduction infection or invasion

Subduction invasion/infection (Waldron et al., 2014; Duarte et al., 2016) occurs when subduction migrates from an older system into a pristine oceanic basin. Waldron et al. (2014) suggest that the closure of the Iapetus Ocean is due to the encroachment of old lithosphere into a young ocean. These authors suggests that subduction initiated at the boundary between old and new oceanic lithosphere and was introduced to the area through trench rollback. This process is thought be similar to the modern day Caribbean, Scotia and Gibraltar Arcs (Duarte et al., 2016). This suggests that the older subductions of the Pacific are invading the younger Atlantic basin, which might potentially lead to collision, orogeny and closure of the Atlantic ocean (Duarte et al., 2016).

Following a subduction polarity reversal

Subduction polarity reversal describes a process where the trench jumps from the subducting plate to the overriding one, flipping its polarity in the process (see figure 3). This can result from the arrival at the trench of continental lithosphere (McKenzie, 1969) or young positively buoyant lithosphere (Crameri and Tackley, 2015). Subduction polarity reversal is often invoked to explain and justify the two juxtaposed Wadati-Benioff zones and their opposite polarities, in the Solomon Island Region (Cooper and Taylor, 1985). Indications of a polarity reversal are also exhibited below the Alpine and Apennine Belts (Vignaroli et al., 2008). Furthermore, Crameri and Tackley (2014) also suggest that the continental connection between South America and the Antarctic peninsula has been severed through a subduction polarity reversal, resulting in the lateral detachment of the South Sandwich subduction zone.

Subduction initiation at ancient/ inherited zones of lithospheric weakness

Subduction zones can also initiate at ancient, inherited zones of weakness such as old fracture zones, transform faults, extinct subduction boundaries and extinct spreading ridges (Gurnis et al., 2004). Gurnis et al. (2004) suggest that the Izu-Bonin-Mariana subduction zone initiated at a fracture zone, while the Tonga-Kermadec subduction initiated at an extinct subduction boundary. The same study also proposes that the incipient Puysegur-Fiordland subduction zone nucleated at an extinct spreading centre.

In conclusion, we can say that subduction zone formation is a complex and multi layered process that can stem from a variety of tectonic settings. However, it is clear that our planet’s current convection style, mode of surface recycling and its ability to sustain life are interlinked with subduction zone formation. Therefore, to understand better how subduction zones form is to better understand what makes the Earth the planet it is today.

 

References:

Baes, M., Govers, R., and Wortel, R. (2011). Subduction initiation along the inherited weakness zone at the edge of a slab: Insights from numerical models. Geophysical Journal International, 184(3):991–1008.

Baes, M., Sobolev, S. V., and Quinteros, J. (2018). Subduction initiation in mid-ocean induced by mantle suction flow. Geophysical Journal International, 215(3):1515–1522.

Beaussier, S. J., Gerya, T. V., and Burg, J.-p. (2018). 3D numerical modelling of the Wilson cycle: structural inheritance of alternating subduction polarity. Fifty years of the Wilson Cycle concept in plate tectonics, page First published online.

Burov, E. and Cloetingh, S. (2010). Plume-like upper mantle instabilities drive subduction initiation. Geophys. Res. Lett., 37(3).

Cooper, P. and Taylor, B. (1985). Polarity reversal in the Solomon Island Arc. Nature, 313(6003):47–48.

Crameri, F., Conrad, C. P., Mont ́esi, L., and Lithgow-Bertelloni, C. R. (2018). The dynamic life of an oceanic plate. Tectonophysics.

Crameri, F. and Tackley, P. J. (2014). Spontaneous development of arcuate single-sided subduction in global 3-D mantle convection models with a free surface. Journal of Geophysical Research: Solid Earth, 119(7):5921–5942.

Crameri, F. and Tackley, P. J. (2015). Parameters controlling dynamically self-consistent plate tectonics and single-sided subduction in global models of mantle convection. Journal of Geophysical Research: Solid Earth, 3(55):1–27.

Duarte, J. C., Schellart, W. P., and Rosas, F. M. (2016). The future of Earth’s oceans: Consequences of subduction initiation in the Atlantic and implications for supercontinent formation. Geological Magazine, 155(1):45–58.

Dymkova, D. and Gerya, T. (2013). Porous fluid flow enables oceanic subduction initiation on Earth. Geophysical Research Letters, 40(21):5671–5676.

Foley, B. J. (2018). The dependence of planetary tectonics on mantle thermal state : applications to early Earth evolution. 376.

Gerya, T. V., Stern, R. J., Baes, M., Sobolev, S. V., and Whattam, S. A. (2015). Plate tectonics on the Earth triggered by plume-induced subduction initiation. Nature, 527(7577):221–225.

Gurnis, M., Hall, C., and Lavier, L. (2004). Evolving force balance during incipient subduction. Geochemistry, Geophysics, Geosystems, 5(7).

Leng, W. and Gurnis, M. (2011). Dynamics of subduction initiation with different evolutionary pathways. Geochemistry, Geophysics, Geosystems, 12(12).

McKenzie, D. P. (1969). Speculations on the Consequences and Causes of Plate Motions. Geophys. J. R. Astron. Soc., 18(1):1–32.

Stern, R. J. (2004). Subduction initiation: spontaneous and induced. Earth and Planetary Science Letters, 226(3-4):275–292.

Stern, R. J. and Gerya, T. (2017). Subduction initiation in nature and models: A review. Tectonophysics.

Thielmann, M. and Kaus, B. J. (2012). Shear heating induced lithospheric-scale localization: Does it result in subduction? Earth Planet. Sci. Lett., 359-360:1–13.

Ueda, K., Gerya, T., and Sobolev, S. V. (2008). Subduction initiation by thermal-chemical plumes: Numerical studies. Phys. Earth Planet. Inter., 171(1-4):296–312.

van Hunen, J. and Moyen, J.-F. (2012). Archean Subduction: Fact or Fiction? Annual Review of Earth and Planetary Sciences, 40(1):195–219.

Vignaroli, G., Faccenna, C., Jolivet, L., Piromallo, C., and Rossetti, F. (2008). Subduction polarity reversal at the junction between the Western Alps and the Northern Apennines, Italy. Tectonophysics, 450(1-4):34–50.

Waldron, J. W., Schofield, D. I., Brendan Murphy, J., and Thomas, C. W. (2014). How was the iapetus ocean infected with subduction? Geology, 42(12):1095–1098.

White, D. A., Roeder, D. H., Nelson, T. H., and Crowell, J. C. (1970). Subduction. Geological Society of America Bulletin, 81(October):3431–3432.

W.K, H. and E.H., C. (2003). Earth’s Dynamic Systems. Prentice Hall; 10 edition.

Tomography and plate tectonics

Tomography and plate tectonics

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 ‘Geodynamics 101’ post for 2019, Assistant Prof. Jonny Wu from the University of Houston explains how to delve into the subduction record via seismic tomography and presents some fascinating 3D workflow images with which to test an identified oceanic slab. 

Jonny Wu, U. Houston

Tomography… wait, isn’t that what happens in your CAT scan? Although the general public might associate tomography with medical imaging, Earth scientists are well aware that ‘seismic tomography’ has enabled us to peer deeper, and with more clarity, into the Earth’s interior (Fig. 1). What are some of the ways we can download and display tomography to inform our scientific discoveries? Why has seismic tomography been a valuable tool for plate reconstructions? And what are some new approaches for incorporating seismic tomography within plate tectonic models?

Figure 1: Tomographic transect across the East Asian mantle under the Eurasian-South China Sea margin, the Philippine Sea and the western Pacific from Wu and Suppe (2018). The displayed tomography is the MITP08 global P-wave model (Li et al., 2008).

Downloading and displaying seismic tomography

Seismic tomography is a technique for imaging the Earth’s interior in 3-D using seismic waves. For complete beginners, IRIS (Incorporated Research Institutions for Seismology) has an excellent introduction that compares seismic tomography to medical CT scans.

A dizzying number of new, high quality seismic tomographic models are being published every year. For example, the IRIS EMC-EarthModels catalogue  currently contains 64 diverse tomographic models that cover most of the Earth, from global to regional scales. From my personal count, at least seven of these models have been added in the past half year – about one new model a month. Aside from the IRIS catalog, a plethora of other tomographic models are also publicly-available from journal data suppositories, personal webpages, or by an e-mail request to the author.

Downloading a tomographic model is just the first step. If one does not have access to custom workflows and scripts to display tomography, consider visiting an online tomography viewer. I have listed a few of these websites at the end of this blog post. Of these websites, a personal favourite of mine is the Hades Underworld Explorer built by Douwe van Hinsbergen and colleagues at Utrecht University, which uses a familiar Google Maps user interface. By simply dragging a left and right pin on the map, a user can display a global tomographic section in real time. The displayed tomographic section can be displayed in either a polar or Cartesian view and exported to a .svg file. Another tool I have found useful are tomographic ‘vote maps’, which provide indications of lower mantle slab imaging robustness by comparison of multiple tomographic models (Shephard et al., 2017). Vote maps can be downloaded from the original paper above or from the SubMachine website (Hosseini et al. (2018); see more in the website list below).

Using tomography for plate tectonic reconstructions

Tomography has played an increasing role in plate tectonic studies over the past decades. A major reason is because classical plate tectonic inputs (e.g. seafloor magnetic anomalies, palaeomagnetism, magmatism, geology) are independent from the seismological inputs for tomographic images. This means that tomography can be used to augment or test classic plate reconstructions in a relatively independent fashion. For example, classical plate tectonic models can be tested by searching tomography for slab-like anomalies below or near predicted subduction zone locations. These ‘conventional’ plate modelling workflows have challenges at convergent margins, however, when the geological record has been significantly destroyed from subduction. In these cases, the plate modeller is forced to describe details of past plate kinematics using an overly sparse geological record.

Figure 2: Tomographic plate modeling workflow proposed by Wu et al. (2016). The final plate model in c) is fully-kinematic and makes testable geological predictions for magmatic histories, terrane paleolatitudes and other geology (e.g. collisions) that can be compared against the remnant geology in d), which are relatively independent.

A ‘tomographic plate modelling’ workflow (Fig. 2) was proposed by Wu et al. (2016) that essentially reversed the conventional plate modelling workflow. In this method, slabs are mapped from tomography and unfolded (i.e. retro-deformed) (Fig. 2a). The unfolded slabs are then populated into a seafloor spreading-based global plate model. Plate motions are assigned in a hierarchical fashion depending on available kinematic constraints (Fig. 2b). The plate modelling will result in either a single unique plate reconstruction, or several families of possible plate models (Fig. 2c). The final plate models (Fig. 2c) are fully-kinematic and make testable geological predictions for magmatic histories, palaeolatitudes and other geological events (e.g. collisions). These predictions can then be systematically compared against remnant geology (Fig. 2d), which are independent from the tomographic inputs (Fig. 2a).

The proposed 3D slab mapping workflow of Wu et al. (2016) assumed that the most robust feature of tomographic slabs is likely the slab center. The slab mapping workflow involved manual picking of a mid-slab ‘curve’ along hundreds (and sometimes thousands!) of variably oriented 2D cross-sections using software GOCAD (Figs. 3a, b). A 3-D triangulated mid-slab surface is then constructed from the mid-slab curves (Fig. 3c). Inspired by 3D seismic interpretation techniques from petroleum geoscience, the tomographic velocities can be extracted along the mid-slab surface for further tectonic analysis (Fig. 3d).


Figure 3: Slab unfolding workflow proposed by Wu et al. (2016) shown for the subducted Ryukyu slab along the northern Philippine Sea plate. The displayed tomography in a), d) and e) is from the MITP08 global P-wave model (Li et al., 2008).

For relatively undeformed upper mantle slabs, a pre-subduction slab size and shape can be estimated by unfolding the mid-slab surface to a spherical Earth model, minimizing distortions and changes to surface area (Fig. 3e). Interestingly, the slab unfolding algorithm can also be applied to shoe design, where there is a need to flatten shoe materials to build cut patterns (Bennis et al., 1991).  The three-dimensional slab mapping within GOCAD allows a self-consistent 3-D Earth model of the mapped slabs to be developed and maintained. This had advantages for East Asia (Wu et al., 2016), where many slabs have apparently subducted in close proximity to each other (Fig. 1).

Web resources for displaying tomography

Hades Underworld Explorer : http://www.atlas-of-the-underworld.org/hades-underworld-explorer/

Seismic Tomography Globe : http://dagik.org/misc/gst/user-guide/index.html

SubMachine : https://www.earth.ox.ac.uk/~smachine/cgi/index.php

 

References

Bennis, C., Vezien, J.-M., Iglesias, G., 1991. Piecewise surface flattening for non-distorted texture mapping. Proceedings of the 18th annual conference on Computer graphics and interactive techniques 25, 237-246.

Hosseini, K. , Matthews, K. J., Sigloch, K. , Shephard, G. E., Domeier, M. and Tsekhmistrenko, M., 2018. SubMachine: Web-Based tools for exploring seismic tomography and other models of Earth's deep interior. Geochemistry, Geophysics, Geosystems, 19. 

Li, C., van der Hilst, R.D., Engdahl, E.R., Burdick, S., 2008. A new global model for P wave speed variations in Earth's mantle. Geochemistry, Geophysics, Geosystems 9, Q05018.

Shephard, G.E., Matthews, K.J., Hosseini, K., Domeier, M., 2017. On the consistency of seismically imaged lower mantle slabs. Scientific Reports 7, 10976.

Wu, J., Suppe, J., 2018. Proto-South China Sea Plate Tectonics Using Subducted Slab Constraints from Tomography. Journal of Earth Science 29, 1304-1318.

Wu, J., Suppe, J., Lu, R., Kanda, R., 2016. Philippine Sea and East Asian plate tectonics since 52 Ma constrained by new subducted slab reconstruction methods. Journal of Geophysical Research: Solid Earth 121, 4670-4741

Inversion 101 to 201 – Part 3: Accounting for uncertainty – Bayes and friends

Inversion 101 to 201 – Part 3: Accounting for uncertainty – Bayes and friends

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 time, Lars Gebraad, PhD student in seismology at ETH Zürich, Switzerland, talks about inversion and how it can be used in geodynamics! Due to the ambitious scope of this topic, we have a 3-part series on this! You can find the first post here. This week, we have our final post in the series, where Lars discusses probabilistic inversion.

One integral part of doing estimations on parameters is an uncertainty analysis. The aim of a general inverse problem is to find the value of a parameter, but it is often very helpful to indicate the measure of certainty. For example in the last figure of my previous post, the measurement values at the surface are more strongly correlated to the upper most blocks. Therefore, the result of an inversion in this set up will most likely be more accurate for these parameters, compared to the bottom blocks.

In linear deterministic inversion, the eigenvalues of the matrix system provide an indication of the resolvability of parameters (as discussed in the aforementioned work by Andrew Curtis). There are classes of methods to compute exact parameter uncertainty in the solution.

From what I know, for non-linear models, uncertainty analysis is limited to the computation of second derivatives of the misfit functional in parameter space. The second derivatives of X (the misfit function) are directly related to the standard deviations of the parameters. Thus, by computing all the second derivatives of X, a non-linear inverse problem can still be interrogated for its uncertainty. However, the problem with this is its linearisation; linearising the model and computing derivatives may not be truly how the model reacts in model space. Also, for strongly non-linear models many trade-offs (correlations) exist which influence the final solution, and these correlations may very strongly depending on the model to be inverted.

Bayes’ rule

Enter reverend Thomas Bayes

This part-time mathematician (he only ever published one mathematical work) from the 18th century formulated the Bayes’ Theorem for the first time, which combines knowledge on parameters. The mathematics behind it can be easily retrieved from our most beloved/hated Wikipedia, so I can avoid getting to caught up in it. What is important is that it allows us to combine two misfit functions or probabilities. Misfits and probabilities are directly interchangeable; a high probability of a model fitting our observations corresponds to a low misfit (and there are actual formulas linking the two). Combining two misfits allows us to accurately combine our pre-existing (or commonly: prior) knowledge on the Earth with the results of an experiment. The benefits of this are two-fold: we can use arbitrarily complex prior knowledge and by using prior knowledge that is bounded (in parameter space) we can still invert underdetermined problems without extra regularisation. In fact, the prior knowledge acts as regularisation.

One probabilistic regularised bread, please

Let’s give our friend Bayes a shot at our non-linear 1D bread. We have to come up with our prior knowledge of the bread, and because we did not need that before I’m just going to conjure something up! We suddenly find the remains of a packaging of 500 grams of flour

This is turning in quite the detective story!

However, the kitchen counter that has been worked on is also royally covered in flour. Therefore, we estimate that probably this pack was used; about 400 grams of it, with an uncertainty (standard deviation) of 25 grams. Mathematically we can formulate our prior knowledge as a Gaussian distribution with the aforementioned standard deviation and combine this with our misfit of the inverse problem (often called the likelihood). The result is given here:

Prior and original misfit

Combined misfit

One success and one failure!

First, we successfully combined the two pieces of information to make an inverse problem that is no longer non-unique (which was a happy coincidence of the prior: it is not guaranteed). However, we failed to make the problem more tractable in terms of computational requirements. To get the result of our combined misfit, we still have to do a systematic grid search, or at least arrive at a (local) minimum using gradient descent methods.

We can do the same in 2D. We combine our likelihood (original inverse problem) with rather exotic prior information, an annulus in model space, to illustrate the power of Bayes’ theorem. The used misfit functions and results are shown here:

Original misfit for a bread of 500 grams

Prior knowledge misfit

Combined misfit

This might also illustrate the need for non-linear uncertainty analysis. Trade-offs at the maxima in model space (last figure, at the intersection of the circle and lines) distinctly show two correlation directions, which might not be fully accounted for by using only second derivative approximations.

Despite this ‘non-progress’ of still requiring a grid search even after applying probability theory, we can go one step further by combining the application of Bayesian inference with the expertise of other fields in appraising inference problems…

Jump around!

Up to now, using a probabilistic (Bayesian) approach has only (apparently) made our life more difficult! Instead of one function, we now have to perform a grid search over the prior and the original problem. That doesn’t seem like a good deal. However, a much used technique in statistics deals with exactly the kind of problems we are facing here: given a very irregular and high dimensional function

How do we extract interesting information (preferably without completely blowing our budget on supercomputers)?

Let’s first say that with interesting information I mean minima (not necessarily restricted to global minima), correlations, and possibly other statistical properties (for our uncertainty analysis). One answer to this question was first applied in Los Alamos around 1950. The researches at the famous institute developed a method to simulate equations of state, which has become known as the Metropolis-Hastings algorithm. The algorithm is able to draw samples from a complicated probability distribution. It became part of a class of methods called Markov Chain Monte Carlo (MCMC) methods, which are often referred to as samplers (technically they would be a subset of all available samplers).
The reason that the Metropolis-Hastings algorithm (and MCMC algorithms in general) is useful, is that a complicated distribution (e.g. the annulus such as in our last figure) does not easily allow us to generate points proportional to its misfit. These methods overcome this difficulty by starting at a certain point in model space and traversing a random path through it – jumping around – but visiting regions only proportional to the misfit. So far, we have only considered directly finding optima to misfit functions, but by generating samples from a probability distribution proportional to the misfit functions, we can readily compute these minima by calculating statistical modes. Uncertainty analysis subsequently comes virtually for free, as we can calculate any statistical property from the sample set.

I won’t try to illustrate any particular MCMC sampler in detail. Nowadays many great tools for visualising MCMC samplers exist. This blog by Alex Rogozhnikov does a beautiful job of both introducing MCMC methods (in general, not just for inversion) and illustrating the Metropolis Hastings Random Walk algorithm as well as the Hamiltonian Monte Carlo algorithm. Hamiltonian Monte Carlo also incorporates gradients of the misfit function, thereby even accelerating the MCMC sampling. Another great tool is this applet by Chi Feng. Different target distributions (misfit functions) can be sampled here by different algorithms.

The field of geophysics has been using these methods for quite some time (Malcom Sambridge writes in 2002 in a very interesting read that the first uses were 30 years ago), but they are becoming increasingly popular. However, strongly non-linear inversions and big numerical simulations are still very expensive to treat probabilistically, and success in inverting such a problem is strongly dependent on the appropriate choice of MCMC sampler.


Summary of probabilistic inversions

In the third part of this blog we saw how to combine any non-linear statistical model, and how to sample these complex functions using MCMC samplers. The resulting sample sets can be used to do an inversion and compute statistical moments of the inverse problem.


Some final thoughts

If you reached this point while reading most of the text, you have very impressively worked your way through a huge amount of inverse theory! Inverse theory is a very diverse and large field, with many ways of approaching a problem. What’s discussed here is, to my knowledge, only a subset of what’s being used ‘in the wild’. These ramblings of a aspiring seismologist might sound uninteresting to the geodynamiscists at the other side of the geophysics field. Inverse methods seem to be not nearly discussed as much in geodynamics as they are in seismology. Maybe it’s the terminology that differs, and that all these concepts are well known and studied under different names and you recognise some of the methods. Otherwise, I hope I have given an insight in the wonderful and sometimes ludicrous mathematical world of (some) seismologists.

Interested in playing around with inversion yourself? You can find a toy code about baking bread here.