GD
Geodynamics

numerical modelling

What controlled the evolution of Plate Tectonics on Earth?

Great Unconformity - Immensity River, Grand Canyon
Stephan Sobolev

Prof. Dr. Stephan Sobolev. Head of the Geodynamic Modelling section of GFZ Potsdam.

Plate tectonics is a key geological process on Earth, shaping its surface, and making it unique among the planets in the Solar System. Yet, how plate tectonics emerged and which factors controlled its evolution remain controversial. The recently published paper in Nature by Sobolev and Brown suggests new ideas to solve this problem….

What makes plate tectonics possible on contemporary Earth?

It is widely accepted that plate tectonics is driven by mantle convection, but is the presence of said convection sufficient for plate tectonics? The answer is no, otherwise plate tectonics would be present on Mars and Venus and not only on Earth. The geodynamic community recognized that another necessary condition for plate tectonics is low strength at plate boundaries and particularly along the plate interfaces in subduction zones (e.g. Zhong and Gurnis 1992, Tackley 1998, Moresi and Solomatov 1998, and Bercovici 2003). To quantify the required strength at subduction interfaces, we have used global models of plate tectonics (Fig. 1A) that combine a finite element numerical technique employing visco-elasto-plastic rheology to model deformation in the upper 300 km of the Earth (Popov and Sobolev 2008) with a spectral code to model convection in the deeper mantle (Steinberger and Calderwood 2006). The model reproduces well present-day plate velocities if the effective friction at convergent plate boundaries is about 0.03 (Fig.1B). Low strength corresponds to subduction interfaces that are well lubricated by continental sediments (low friction; Lamb and Davis 2003, Sobolev and Babeyko 2005, or low viscosity; Behr and Becker 2018). In case of sediment shortages in the trenches (corresponding to a friction coefficient of 0.07-0.1), plate velocities would first decrease about two times (Fig. 1C) and then even more because of less negatively buoyant material having subducted into the mantle, leading to less convection driving force.

Effects of sediments on contemporary subduction according to global numerical models.

Figure 1. Global numerical model showing the effect of sediments on contemporary subduction. (A) The global model combines two computational domains coupled through continuity of velocities and tractions at 300 km depth. (B) NUVEL 1A plate velocities in a no-net-rotation reference frame (black arrows) versus computed velocities (blue arrows) for the global model with a friction of 0.03 at convergent plate boundaries. (C) Root mean square of computed plate velocities in the global model versus friction coefficient at convergent plate boundaries.

Hypothesis and its testing

Based on the previous discussion, we infer that continental sediments in subduction channels act as a lubricant for subduction. In addition, the presence of these sediments in trenches is a necessary condition for the stable operation of plate tectonics, particularly earlier in Earth’s evolution when the mantle was warmer and slabs were relatively weak. With this hypothesis we challenge the popular view that secular cooling of the Earth was the only major control on the evolution of plate tectonics on Earth since about 3 Ga. The hypothesis predicts that periods of stable plate tectonics should follow widespread surface erosion events, whereas times of diminished surface erosion should be associated with reduced subduction and possibly intermittent plate tectonics.

We test this prediction using geological proxies believed to identify plate tectonics activity (supercontinental cycles) and geochemical proxies that trace the influence of the continental crust on the composition of seawater (Sr isotopes in ocean sediments; Shields 2007) and continental sediments in the source of subduction-related magmas (oxygen and Hf isotopes in zircons; Cawood et al. 2013, Spencer et al. 2017). All three geochemical markers indeed show that just before or in the beginning of supercontinental cycles the influence of sediments is increasing, while it decreases before periods of diminished plate tectonic activity, like the boring billion period between 1.7 and 0.7 Ga (Cawood and Hawkesworth 2014; Fig. 2). The largest surface erosion and subduction lubrication events were likely associated with the global glaciation evens identified in the beginning (2.5-2.2 Ga) and at the end (0.7-0.6 Ga) of the Proterozoic Era (Hoffman and Schrag 2002). The latter snowball Earth glaciation event terminated the boring billion period and kick-started the modern phase of active plate tectonics.

Another prediction of our hypothesis is that in order to start plate tectonics, continents had to rise above sea level and provide sediments to the oceans. This prediction is again consistent with observations: there are many arguments for the beginning of plate tectonics between 3 and 2.5 Ga (see the review of Condie 2018) and, at the same time, this period is most likely when the continents rose above sea level (Korenaga et al. 2017).

Cartoon summarizing the factors that control the emergence and evolution of plate tectonics on Earth.

Figure 2. Cartoon summarizing the factors that control the emergence and evolution of plate tectonics on Earth. Enhanced surface erosion due to the rise of the continents and major glaciations stabilized subduction and plate tectonics for some periods after 3 Ga and particularly after 0.7 Ga after the cooling of the mantle. Blue boxes mark major glaciations; transparent green rectangles show the time intervals when all three geochemical proxies consistently indicate increasing sediment influence (major lubrication events); and, a thick black dashed curve separates hypothetical domains of stable and unstable plate tectonics. The reddish domain shows the number of passive margins (Bradley 2008), here used as a proxy for plate tectonic intensity.

What was before plate tectonics?

The earlier geodynamic regime could have involved episodic lid overturn and resurfacing due to retreating large-scale subduction triggered by mantle plumes (Gerya et al. 2015) or meteoritic impacts (O’Neill et al. 2017). Retreating slabs would bring water into the upwelling hot asthenospheric mantle, generating a large volume of magma that formed protocontinents. Extension of the protocontinental crust could have produced nascent subduction channels (Rey et al. 2014) along the edges of the protocontinents lubricated by the sediments. In this way, a global plate tectonics regime could have evolved from a retreating subduction regime.

What is next?

Despite of the support from existing data, more geochemical information is required to conclusively test our hypothesis about the role of sediments in the evolution of plate tectonics. Additionally, this hypothesis must be fully quantified, which in turn will require coupled modeling of mantle convection and plate tectonics, surface processes and climate.

References
Behr, W. M. and Becker, T. W. Sediment control on subduction plate speeds. Earth Planet. Sci. Lett. 502, 166-173 (2018).

Bercovici, D. The generation of plate tectonics from mantle convection. Earth Planet. Sci. Lett. 205, 107–121 (2003).

Bradley, D. C. Passive margins through earth history. Earth Sci. Rev. 91, 1-26 (2008).

Cawood, P. A., Hawkesworth, C. J. and Dhuime, B. The continental record and the generation of continental crust. Geol. Soc. Amer. Bull. 125, 14-32 (2013).

Cawood, P. A. and Hawkesworth, C. J. Earth's middle age. Geology 42, 503-506 (2014).

Condie, K. C. A planet in transition: The onset of plate tectonics on Earth between 3 and 2 Ga? Geosci. Front. 9, 51-60 (2018).

Gerya, T.V. et al. Plate tectonics on the Earth triggered by plume-induced subduction initiation, Nature 527, 221-225 (2015).

Hoffman, P. F. and Schrag, D. P. The snowball Earth hypothesis: testing the limits of global change. Terra Nova 14, 129–155 (2002).

Korenaga, J., Planavsky, N. J. and Evans, D. A. D. Global water cycle and the coevolution of the Earth's interior and surface environment. Phil. Trans. R. Soc. Am. 375, 20150393 (2017).

Lamb, S. and Davis, P. Cenozoic climate change as a possible cause for the rise of the Andes. Nature 425, 792-797 (2003).

Moresi, L. and Solomatov, V. Mantle convection with a brittle lithosphere: Thoughts on the global tectonic style of the Earth and Venus. Geophys. J. Int. 133, 669-682 (1998).

O’Neill, C. et al. Impact-driven subduction on the Hadean Earth. Nature Geosci. 10, 793-797 (2017).

Popov, A.A. and Sobolev, S. V. SLIM3D: A tool for three-dimensional thermo mechanical modeling of lithospheric deformation with elasto-visco-plastic rheology, Phys. Earth Planet. Inter. 171, 55-75 (2008).

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

Shields, G. A. A normalised seawater strontium isotope curve: possible implications for Neoproterozoic-Cambrian weathering rates and the further oxygenation of the Earth. eEarth 2, 35-42 (2007).

Sobolev, S. V. and Babeyko, A. Y. What drives orogeny in the Andes? Geology 33, 617-620 (2005).

Spencer, C. J., Roberts, N. M. W. and Santosh, M. Growth, destruction, and preservation of Earth's continental crust. Earth. Sci. Rev. 172, 87-106 (2017).

Steinberger, B. and Calderwood, A. Models of large-scale viscous flow in the Earth’s mantle with constraints from mineral physics and surface observations. Geophys. J. Intern., 167 1461–1481 (2006).

Tackley, P. J. Self-consistent generation of tectonic plates in three-dimensional mantle convection. Earth Planet. Sci. Lett. 157, 9-22, (1998).

Zhong, S. and Gurnis, M. Viscous flow model of a subduction zone with a faulted lithosphere: long and short wavelength topography, gravity and geoid. Geophys. Res. Lett. 19, 1891–1894 (1992).

 

Remarkable Regions – The Réunion Hotspot

Remarkable Regions – The Réunion Hotspot
Eva Bredow at Réunion caldera.

Eva Bredow in front of the caldera at Réunion Island. Credit: Simon Stähler.

This week we again turn our attention to a Remarkable Region that deserves a spot in the scientific limelight. Postdoctoral researcher Eva Bredow of Kiel University shares with us her long history with Réunion Island.

At first glance, Réunion is a relatively small tropical island, located between Madagascar and Mauritius, and from my personal experience, most Germans have never even heard of it. To be fair, it is much better known in France, because Réunion is officially a French overseas department, meaning that the eleven-hour flight from Paris is technically a domestic flight and that you can pay there with Euros (and I bet you did not know that a millimetre-sized outline of the island appears on every Euro banknote!). Besides, Réunion hosts one of the most active volcanoes in the world with one eruption per year on average. However, it rarely hits the headlines because the inhabitants live far enough away not to be overly threatened. And yet, for people interested in geodynamics, the name Réunion might actually have a familiar sound, since it regularly appears in hotspot catalogues and hotspot reference frames – a sure indication that there is more to discover.

For me, Réunion has been a very special place ever since I was a high school student lucky enough to visit the island in order to learn French. And who would have thought back then that hiking in this surreal volcanic landscape would be one of the first steps towards my decision to study geophysics? And what were the odds to stumble upon a PhD project years later, centred around the Réunion hotspot? Well, that is exactly what happened and in this article, it is my pleasure to give you at least a brief overview of why Réunion deserves to be called a remarkable spot indeed and how numerical modelling can help us to explore its geodynamic history.

NW Indian Ocean crustal thickness map.

Crustal thickness map of the north-western Indian Ocean with the entire hotspot track from Réunion Island to the Deccan Traps in India. Figure from Torsvik et al. (2013).

A deep root

The hypothesis that Réunion is an intraplate hotspot possibly fed by a hot, buoyant upwelling rooted deep in the mantle was already put forward by Jason Morgan (1971, 1972) in his famous papers outlining the classical mantle plume hypothesis. And as it happens, the Réunion plume has left a number of traces that fit the plume hypothesis extremely well and make it a kind of prototype for a deep plume and its surface manifestations. A brief look at a topographic map of the north-western Indian Ocean reveals not only the currently active hotspot at Réunion and the slightly older island of Mauritius, but also a clearly continuous (and age-progressive) hotspot track on the African and Indian plates, only split due to subsequent seafloor-spreading.

According to numerous laboratory and numerical studies that describe the mushroom-like geometry of a plume, the hotspot track is considered to be caused by the long-lived plume tail, whereas the voluminous plume head is supposed to create a huge flood basalt province in a relatively short geological time (Richards et al., 1989). In the case of the Réunion plume, the hotspot track starts at the Deccan Traps, a gigantic continental Large Igneous Province (LIP) in India. The LIP was created around 65 million years ago and the environmental changes triggered by the volcanic activities might have led to the extinction of the dinosaurs (an alternative theory to the Chicxulub impact in Mexico; Courtillot and Renne, 2003).

Further indications for a deep plume beneath Réunion include the broad topographic hotspot swell around the island, a geochemical signature of the volcanic rocks that clearly deviates from mid-ocean ridge basalts, and the present-day hotspot location above the plume generation zone at the margin of the African Large Low Shear Velocity Province (LLSVP).

Plume-ridge interaction

A more puzzling observation is the geochemical anomaly at the closest segments of the Central Indian Ridge, about 1000 km away from Réunion that implies a long-distance plume-ridge interaction. Already Morgan (1978) suggested that a sublithospheric flow channel connecting the upwelling plume and the ridge is responsible for the creation of the Rodrigues Ridge, a rather eye-catching feature not at all parallel to the hotspot track or recent plate motions.

And there is one more noteworthy hypothesis associated with Réunion, based on extremely old zircons found at Mauritius; it postulates that the hotspot track has (coincidentally) been created on top of a Precambrian microcontinent (Ashwal et al., 2017).

The RHUM-RUM experiment (completely alcohol-free…)

Concerning the (present-day) state of the Réunion plume at greater depths, seismic tomography is the most promising tool to answer the question if it is indeed fed by a deep plume or not. But given that the island is rather remotely located and a classical plume tail is expected to be quite narrow, there are plenty of technical obstacles, and it was not until 2006 that Montelli published the first seismic image of a continuous plume conduit reaching into the deep mantle. More recent global tomography models also image the Réunion plume as a clearly resolved, vertically continuous conduit at depths between 1,000 and 2,800 km (French and Romanowicz, 2015).

In 2012-2013, the French-German RHUM-RUM project (Réunion Hotspot and Upper Mantle – Réunions Unterer Mantel) aimed at an even higher resolved image of the plume. Therefore, 57 German and French ocean-bottom seismometers were deployed at the seafloor around Réunion for about a year (Stähler et al., 2016) – still the largest seismological experiment to image a deep oceanic mantle plume so far.

 

RHUM-RUM seismic stations

All seismic stations related to the RHUM-RUM project, with the 57 ocean-bottom seismometer stations shown in red. More information on the project can be found here.

With all that in mind, and as part of the RHUM-RUM project, I set up a regional numerical model with some colleagues from the GFZ Potsdam in order to assemble Réunion’s entire dynamic history. We used time-dependent plate reconstructions and large-scale mantle flow as velocity boundary conditions as well as a laterally varying lithosphere thickness in order to specifically simulate the Réunion plume (for details, see Bredow et al., 2017). In short: altogether, we were able to reproduce a crustal thickness pattern that at first order fits the observed hotspot track (although the method is not suited to reproduce a continental LIP such as the Deccan Traps). Moreover, the interaction between the plume and the Central Indian Ridge explained both the genesis of the Rodrigues Ridge and the gap in crustal thickness between the Maldives and Chagos – both features that have not been dynamically modelled before.

After our models were published, the active long-distance plume-ridge interaction beneath the Rodrigues Ridge was additionally confirmed by seismological studies in the RHUM-RUM project: first in a three-dimensional anisotropic S-wave velocity model comprising the uppermost 300 km (Mazzullo et al., 2017), and second by SKS splitting measurements (Scholz et al., 2018). Overall, these interdisciplinary studies confirmed Morgan’s long-standing hypothesis – more than 30 years after its original publication.

 

Cross section geodynamic plume model of Bredow et al. 2017.

Cross section of the geodynamic plume model, showing the long-distance plume-ridge interaction as predicted by Morgan (1978). Figure after Bredow et al. (2017).

Surface wave tomography showing the Reunion plume.

Cross section of the surface wave tomography model, showing the low velocity signature of the plume rising toward the base of the lithosphere underneath Réunion and the sublithospheric flow toward the Central Indian Ridge (CIR). Figure after Mazzullo et al. (2017).

The whole-mantle P- and S-wave tomography models from the RHUM-RUM project have yet to be published, but the (almost final) results presented at this year’s EGU (Tsekhmistrenko et al., 2019) were quite intriguing: while the plume conduit can continuously be followed down to the LLSVP in the deep mantle, the conduit is not as narrow and not nearly as vertical as classically expected!

Therefore I think it is quite safe to say that we have not yet heard the last of the Réunion hotspot and I hope that the next time you hear this name, maybe you will remember it as a rather remarkable spot on our planet…

 

Ashwal et al. (2017), Archaean zircons in Miocene oceanic hotspot rocks establish ancient continental crust beneath Mauritius, Nat. Commun., 8, 14,086, doi: 10.1038/ncomms14086.

Bredow, E. et al. (2017), How plume-ridge interaction shapes the crustal thickness pattern of the Réunion hotspot track, Geochem. Geophys. Geosyst., 18, doi:10.1002/2017GC006875.

Courtillot, V. E. and P. R. Renne (2003), On the ages of flood basalt events, C. R. Geosci., 335(1), 113–140, doi: 10.1016/S1631-0713(03)00006-3.

French, S. W. and B. Romanowicz (2015), Broad plumes rooted at the base of the Earth’s mantle beneath major hotspots, Nature, 525, 95–99, doi: 10.1038/nature14876.

Mazzullo, A. et al. (2017), Anisotropic tomography around Réunion Island from Rayleigh waves Journal of Geophysical Research: Solid Earth, 122, doi: 10.1002/2017JB014354.

Montelli, R. et al. (2006), A catalogue of deep mantle plumes: New results from finite-frequency tomography, Geochem. Geophys. Geosyst., 7, Q11007, doi: 10.1029/2006GC001248.

Morgan, W. J. (1971), Convection plumes in the lower mantle, Nature, 230, 42–43, doi: 10.1038/230042a0.

Morgan, W. J. (1972), Deep mantle convection plumes and plate motions, AAPG bulletin, 56(2), 203–213.

Morgan, W. J. (1978), Rodriguez, Darwin, Amsterdam, ..., A second type of Hotspot Island, J. Geophys. Res., 83(B11), 5355–5360, doi: 10.1029/JB083iB11p05355.

Richards, M. A. et al. (1989), Flood Basalts and Hot-Spot Tracks: Plume Heads and Tails, Science, 246, 103–107, doi: 10.1126/science.246.4926.103.

Scholz, J.-R. et al. (2018), SKS splitting in the Western Indian Ocean from land and seafloor seismometers: Plume, plate and ridge signatures, Earth Planet. Sci. Lett., Volume 498, 169-184, doi: 10.1016/j.epsl.2018.06.033.

Stähler, S. C. et al. (2016), Performance report of the RHUM-RUM ocean bottom seismometer network around La Réunion, western Indian Ocean, Adv. Geosci., 41, 43-63, doi: 10.5194/adgeo-41-43-2016.

Torsvik, T. H. et al. (2013), A Precambrian microcontinent in the Indian Ocean, Nat. Geosci., 6(3), 223–227, doi: 10.1038/ngeo1736.

Tsekhmistrenko, M. et al. (2019), Deep mantle upwelling under Réunion hotspot and the western Indian Ocean from P- and S-wave tomography, Geophysical Research Abstracts, Vol. 21, EGU2019-9447, EGU GA 2019.

Searching for future directions in tectonic modelling

Searching for future directions in tectonic modelling

Geoscientists frequently use forward geodynamic simulations to test hypotheses derived from geophysical and geologic observations. While numerical simulations of lithospheric deformation have lead to key advances in our understanding of tectonic processes, in many cases it remains difficult to ascertain whether numerical models reproduce observations for the correct underlying regions.  This week, John Naliboff and Jolante van Wijk discuss this issue, and talk about a White Paper being prepared by the Computational Infrastructure for Geodynamics (CIG) Long-Term Tectonics Working Group on this topic.

John Naliboff. Assistant Research Scientist in the Department of Earth and Planetary Sciences, UC Davis.

In recent years, advances in numerical methodology, high-performance computing and elucidation of complex geologic observations have enabled 3-D simulations of long-term lithospheric deformation at kilometer-scale resolution and with complex non-linear material behavior. The lithosphere models generally include rheological and compositional layering that delineate a brittle upper crust, ductile (viscous) to brittle mid- to lower crust and a brittle to ductile lithospheric mantle overlying a purely viscous asthenosphere. Inherently, the rheological behavior of distinct layers varies temporally through complex feedbacks between temperature, grain-size, strain-rate, phase transitions, flexural stresses and finite deformation. Across a wide range of tectonic settings, numerical investigations incorporating this type of behavior have qualitatively and in some cases quantitatively reproduced key first- and second-order geologic observations.

Jolante van Wijk. Associate Professor in the Department of Earth and Environmental Science, New Mexico Tech.

Despite these successes, numerous significant challenges remain as the computational tectonics community looks toward investigations that account for physical processes acting across a wide range of spatial and temporal scales (Figure 1). Geodynamic model development currently evolves around modifying existing models to include surface processes, thermodynamically consistent melt and volatile transport, metamorphic reactions, and brittle failure to reproduce characteristic features of the seismic cycle.

While many of these processes or features are active areas of research and have been addressed on an individual basis, it often remains unclear at best as to how one should numerically validate even the simplest models of lithospheric deformation. In other words, one can ask whether numerical models of lithospheric processes are reproducing key observations for the correct underlying reasons. Significantly, this question of validation equally applies to observational studies: given that many geologic processes contain significant feedbacks across vast spatial and temporal scales, to what degree can a set of specific observations be interpreted to meaningfully reflect first-order processes?

Continued close collaboration between observational, experimental, and computational Earth scientists is needed to overcome these challenges. At present, the CIG Long-Term Tectonics Working Group is preparing a white paper draft that outlines a 5-10 year vision for collaboration between computational Earth scientists and experimental and observational communities. Given the vast series of topics and disciplines associated with lithospheric dynamics, the White Paper will be organized around Transitional Domains within the lithosphere. The Transitional Domains include Earth’s surface, the brittle-ductile transition, the Moho, the mid-lithosphere discontinuity, and the lithosphere-asthenosphere boundary. For each of these domains we will address the following questions:

  1. How are these domains characterized in the Earth’s lithosphere?
  2. (How) have these domains been modeled previously?
  3. What steps can we take to improve the characterization of these transitional domains in numerical models?
  4. What new methods need to be developed to implement the domains?

 

Figure 1. Spatial and temporal scales associated with distinct lithospheric processes, which was published in Cooper et al., 2015, GSA Today, v. 25(6), pp. 42-43 (Figure 1a).

 

A new generation of geodynamic models will be developed to include the transitional domains. These models will need to be validated, using datasets from the observational and experimental communities, with newly developed techniques.

The White Paper will also include sections on increasing the value of lithospheric models for other scientific communities, and on a pathway toward increasing societal relevance of our modeling efforts.

Please contact John Naliboff (jbnaliboff@ucdavis.edu) and Jolante van Wijk (jolante.vanwijk@nmt.edu) for suggestions or questions. A draft White Paper will be presented at CIG’s annual meeting at AGU 2019 in San Francisco.

It doesn’t work! (Asking questions about scientific software)

It doesn’t work! (Asking questions about scientific software)

Numerical modelling is not always a walk in the park. In fact, many of us occasionally encounter problems that we cannot directly solve ourselves, and thus rely on help from others. In this month’s Wit & Wisdom post, Patrick Sanan, postdoctoral researcher at the Geophysical Fluid Dynamics group at ETH Zurich, will talk about asking the right questions about scientific software. As an experienced scientific software developer, Patrick has often been at the “receiving end” of questions regarding numerical modelling and hopes to guide you through some important points that could make life for you as well as your ‘helper’ a lot easier.

 

Patrick Sanan is a postdoctoral researcher at the Geophysical Fluid Dynamics group at ETH Zurich

Numerical modelling is essential for geodynamics; since we cannot directly measure relevant phenomena, we partake in the magic of making a set of reasonable assumptions, setting up a model, and letting a system evolve to produce insight. It’s beautiful. We gain an understanding of the subtle-yet-fundamental processes which shaped our apparently-so-special Earth. We turn our eye beyond, to other planets.

I’m not here to talk about that, though. I’m here to discuss an ugly part of the job, which can bring all the profundity to a screeching halt: what to do when the code doesn’t work.

You know the situation. You’re stuck. There’s no output. Segmentation fault. Error Code -123. You didn’t sign up for this…

You don’t know where to start looking. Is it your model parameters? Your physical assumptions? Are you using the code as intended? Is it your compiler? Can it be the cluster? Your .bashrc? Is your keyboard plugged in?

You’re frustrated. No one around you has a good suggestion. Why are these cruel computers doing this to you?

What to do? Ask for help, in the right way. In this blog, I’ll point out some facts about scientific software, try to use them to formulate an effective email in which you ask for help, and then try to extract some guiding principles.

Some Points on Scientific Software

First, I’ll list some observations I, as a developer, have made about scientific software:

  • Scientific software projects are usually short on maintainers and time.
  • Software designers are not psychic, but they are often experts at deductive problem-solving.
  • Reproducible, bisectable problems are surprisingly easy to solve. Other problems are surprisingly difficult.
  • Working reference cases are very valuable.
  • Sufficiently-complex computing tasks must be treated like lab experiments: one must document the setup and control as many factors as possible.
  • The solutions to most problems are obvious, once found.
  • Numerical software is harder to test than generic software, because floating point arithmetic and parallel computing lead to acceptable differences in numerical output. Higher-level understanding of physics or (parallel) numerical methods is often required to know if something is a “real problem”.

With these as a guide, let’s consider how one might try to resolve a confusing issue.

No pretty pictures here, just error messages (if you’re lucky).

Asking for Help: The Bad Way and the Good Way

Email is a common way to ask for help, as often the person who can best help you is across the world. Let’s say that I’m using a regional lithospheric dynamics code called Rifter3D. I’ve come across an error I have no explanation for. I can’t figure it out, so I write an email to developers of the code. You might also write to a dedicated help address.

Hello - I'm using Rifter3D but on our local cluster I get this error:

/cluster/shadow/.lsbatch/1559505025.92314771: line 8: 951 Segmentation fault (core dumped) ./rifter3d -options_file options.opts

Do you know what I'm doing wrong?

The person on the other end wants to help, but has no information to work with and doesn’t know how much of their time you’re asking for. There is a better way:

Hello - I'm having some trouble diagnosing a problem using Rifter3D and was hoping
you could give me some pointers.

I've been working with Prof. XYZ and have modified a rifting scenario from XYZ et al. 2017 to study the effects of varying A on B. When I run a small case on my laptop, the simulation finishes as expected. I use the attached options_small.opts and run mpiexec -np 4 ./rifter3d -options_file options_small.opts. I'd like to run a bigger case (512 x 256 x 64) for 300 million simulated years, on 64 cores on our cluster.
However, my job fails, producing a segmentation fault almost immediately. I attached the job submission script (job.sbatch), and output from my run (lsf.o92314771). I am using the cluster's existing PETSc 3.11 module, and version 1.0.3 of Rifter3D. I can successfully run a simple isoviscous setup on this cluster - see attached job_small.lsf, option2.opts, and lsf.o92314681

Do you have any insight as to how I might be able to run this simulation?

Best,
E. Scholar

Attachments: options_small.opts options.opts job.lsf job_small.lsf options2.opts lsf.o92314771 lsf.o92314681

The recipient will likely respond with more questions as you work towards resolving the issue. Perhaps they’ll ask you for more information about how you built Rifter3D, or point out some unusual settings in your options file.

Why is this second email better? It is not simply longer, but

  • It clearly describes the problem. Just trying to precisely describe a problem has an almost-magical clarifying effect, and the solution will often appear. Software engineers call this rubber ducking.
  • It explains the true objective and how the problem is to be resolved. This is important to avoid the XY problem, describing a problem with a method to achieve a goal, without mentioning the goal itself.
  • It gives concrete information to reproduce the problem: the version of the code, input files, and launch commands/scripts.
  • It provides enough output to allow deductive reasoning, more than just a copy-and-pasted error message.
  • It’s polite
  • It shows some effort has already been put into investigating.
  • It notes similar, working cases.
  • It’s not too long, but it is detailed enough to allow for quick, intelligent follow-up questions. Supporting data are included as attachments or links.
  • It doesn’t make too many assumptions about the cause of the problem.

Boiling it Down: 3 Questions to Ask Yourself

When asking for help, consider these three questions. They will help with the central objective: clearly describing the problem.

  • Why do you need it to work?
    What is the context? What is the goal?
  • How do you show that it doesn’t work?
    What are the steps to reproduce your problem?
    How will you know that the problem is resolved?
  • What does work?
    How far are you from a working state? What similar cases work?

Here is a 1-page pdf which you can print out, with some of this advice:

Pdf and LaTeX source on GitHub

The “Real” Answer

To conclude, I will try to avoid an “XY problem” of my own. The most efficient way to resolve bewildering problems is to avoid them. To make an alpine analogy, the most important topic in avalanche safety is not how to dig someone out, it’s how to avoid risky terrain.

Real-life debugging. Avoid if at all possible (from Wikimedia commons)

First, borrow techniques from software engineering. Version control (e.g. git) will encourage you to save working states, amongst many other benefits. Next, leverage your intuition and experience as a geodynamicist. Always be able to quickly run and verify small, quick, simple cases, and test and visualize often. Look out for simple cases where you “know the answer ahead of time”: established benchmarks and analytical solutions.

The points in this article can help everyone save time (and not just running geodynamical models!): problems will be resolved more quickly, bugs will get fixed faster, and more time can be spent exploring more interesting questions than “Why doesn’t it work?”.