Geodynamic diagnostics, scientiﬁc visualisation and StagLab 3.0

. Today’s geodynamic models can, often do and sometimes have to become very complex. Their underlying, increasingly elaborate numerical codes produce a growing amount of raw data. Post-processing such data is therefore becoming more and more important, but also more challenging and time-consuming. In addition, visualising processed data and results has, in times of coloured ﬁgures and a wealth of half-scientiﬁc software, become one of the weakest pillars of science, widely mistreated and ignored. Efﬁcient and automated geodynamic diagnostics and sensible scientiﬁc visualisation preventing common pitfalls is thus more important than ever. Here, a collection of numerous diagnostics for plate tectonics and mantle dynamics is provided and a case for truly scientiﬁc visualisation is made. Amongst other diagnostics are a most accurate and robust plate-boundary identi-ﬁcation, slab-polarity recognition, plate-bending derivation, surface-topography component splitting and mantle-plume detection. Thanks to powerful image processing tools and other elaborate algorithms, these and many other insightful diagnostics are conveniently derived from only a subset of the most basic parameter ﬁelds. A brand new set of sci-entiﬁc quality, perceptually uniform colour maps including devon , davos , oslo and broc is introduced and made freely available (http://www.

Abstract.Today's geodynamic models can, often do and sometimes have to become very complex.Their underlying, increasingly elaborate numerical codes produce a growing amount of raw data.Post-processing such data is therefore becoming more and more important, but also more challenging and time-consuming.In addition, visualising processed data and results has, in times of coloured figures and a wealth of half-scientific software, become one of the weakest pillars of science, widely mistreated and ignored.Efficient and automated geodynamic diagnostics and sensible scientific visualisation preventing common pitfalls is thus more important than ever.Here, a collection of numerous diagnostics for plate tectonics and mantle dynamics is provided and a case for truly scientific visualisation is made.Amongst other diagnostics are a most accurate and robust plate-boundary identification, slab-polarity recognition, plate-bending derivation, surface-topography component splitting and mantle-plume detection.Thanks to powerful image processing tools and other elaborate algorithms, these and many other insightful diagnostics are conveniently derived from only a subset of the most basic parameter fields.A brand new set of scientific quality, perceptually uniform colour maps including devon, davos, oslo and broc is introduced and made freely available (http://www.fabiocrameri.ch/colourmaps, last access: 25 June 2018).These novel colour maps bring a significant advantage over misleading, non-scientific colour maps like rainbow, which is shown to introduce a visual error to the underlying data of up to 7.5 %.Finally, STAGLAB (http: //www.fabiocrameri.ch/StagLab, last access: 25 June 2018) is introduced, a software package that incorporates the whole suite of automated geodynamic diagnostics and, on top of that, applies state-of-the-art scientific visualisation to produce publication-ready figures and movies, all in the blink of an eye and all fully reproducible.STAGLAB, a simple, flexible, efficient and reliable tool made freely available to everyone, is written in MATLAB and adjustable for use with geodynamic mantle convection codes.

Overview
The first basic numerical geodynamic models were developed in the early 1970s (e.g.Minear and Toksöz, 1970;Torrance and Turcotte, 1971).Since then they have become more powerful and often more complex (see e.g.King, 2001;Gerya, 2011;Lowman, 2011;Coltice et al., 2017).Indeed, dynamically self-consistent geodynamic models used to reproduce the first-order characteristics of the complex platemantle system, like mobile surface plates, single-sided subduction and mantle plumes, that need a certain complexity (e.g.Crameri and Tackley, 2015).However, this complexity often inhibits a simple understanding of the full interplay between all individual physical aspects of these models: the models become too complicated to be easily explained or even fully understood.In addition, the more elaborate numerical codes powering these models (e.g.Zhong et al., 2000;Gerya and Yuen, 2007;Moresi et al., 2007;Tackley, 2008;Davies et al., 2011;Thieulot, 2014;Kaus et al., 2016;Heister et al., 2017) and the still increasing computational power available for their execution produce more and more raw data.The resulting amount of data to be processed easily exceeds the capability of a human scientist.Today, efficient, automated and intelligent geodynamic diagnostics are thus more important than ever to keep up with the advances made in numerical models.
Moreover, conveying new findings to the community critically depends on data visualisation as figures are pivotal to make raw data tangible, understandable and explainable (e.g.Gerya, 2010).However, it has also become increasingly difficult to visualise research given the increasing complexity Published by Copernicus Publications on behalf of the European Geosciences Union.
F. Crameri: StagLab 3.0 of models (e.g.high-resolution and 3-D geometry) and new and improved visualisation techniques (e.g.coloured figures and movies).Worrisome visualisation pitfalls arise that make figures confusing, unreadable or even misleading and hence unscientific.The rainbow colour map strongly deteriorates, for example, the underlying scientific data (Pizer and Zimmerman, 1983;Ware, 1988;Tufte, 1997) due to the inhomogeneous colour sensitivity of the human retina (Thomson and Wright, 1947).On top of that, it is confusing to people with some of the most common colour-vision deficiencies.Even though all this has been known for awhile (see Rogowitz andTreinish, 1996, 1998;Light and Bartlein, 2004;Borland and Ii, 2007, and #endrainbow), the rainbow colour map is still commonly used by scientists and commonly accepted by scientific journals.In addition to the current need for more advanced geodynamic diagnostics, it is thus clear that, today, efficient scientific quality visualisation has also become more important than ever.
Nevertheless, there is currently no tool available to perform key geodynamic diagnostics and scientific visualisation efficiently and reliably.STAGLAB (http://www.fabiocrameri.ch/StagLab, last access: 25 June 2018) is the post-processing and visualisation software that aims to fill this gap and fulfil all the necessities mentioned above.Previous versions of STAGLAB (Crameri, 2013(Crameri, , 2017a) ) were designed to specifically handle the raw data produced by the finite-difference, finite-volume multi-grid mantle convection code StagYY (Tackley, 2008).STAGLAB 3.0 (Crameri, 2017b), however, offers potential compatibility with other geodynamic codes as well and has, for example, already been tested successfully with output from the finite-element code Fluidity (Davies et al., 2011).
In the following, I will present key geodynamic diagnostics and how to automate them, discuss scientific and nonscientific colour schemes and visualisation approaches, and introduce STAGLAB, the all-in-one software package that makes powerful geodynamic diagnostics and state-of-the-art scientific visualisation easily accessible to everyone.

Geodynamic diagnostics
In science, it is desirable to rely on quantitative measures rather than on visual impression, also for reasons outlined in Sect.3. Here, I list a variety of key diagnostics for the dynamics of the plate-mantle system that all can be automatised and applied to digital data.While some of these diagnostics are already well known, others are significantly improved or newly introduced here.The geodynamic diagnostics covered here are listed in Table 2 and explained in detail below and can be predominantly applied to 2-D data or vertical slices of Cartesian 3-D data.All of the diagnostics are tested and implemented in the software STAGLAB 3.0 (see Sect. 4).

Generic flow diagnostics
A suite of generic flow diagnostics is often calculated and exported to file directly by fluid-dynamics codes (see e.g.Zhong et al., 1998;Tackley, 2000).These include, for example, root mean square (RMS) flow velocity, the Nusselt number, which is the ratio of the total heat flux across a fluid layer to its conductive value (Chandrasekhar, 1961), heat flow and mean temperature.These are often either given as a function of depth or as depth-average values.Other generic flow diagnostics like plateness and mobility, both characterising the surface boundary layer of mantle convection, are mentioned below in more detail in addition to a list of different approaches to derive residual mantle temperatures.

Plateness
One of the key characteristics of ocean-plate tectonics is wide, almost rigid plate interiors bounded by localised, weak plate boundaries (Crameri et al., 2018).Plateness is a measure of how localised the surface deformation is and thus for how well a surface represents ocean-plate tectonics (Weinstein and Olson, 1992;Tackley, 2000).Plateness is derived from the second invariant of the strain rate on an xy plane spanning the surface from which the total integrated εsurf and the surface-area fraction f 80 , in which the highest 80 % of that deformation occurs, can be calculated (Tackley, 2000).For perfect plates, f 80 would be zero, with deformation only taking place within infinitely narrow zones, while for uniformly distributed plates, f 80 would be 0.8.A non-dimensional plateness parameter, P , is finally derived by where f c is the characteristic surface-area fraction for a flow with given vigour and heating mode, and is, for example, f c = 0.6 for internally heated convection with a Rayleigh number of Ra H = 10 6 (Tackley, 2000).

Mobility
Another key characteristic of ocean-plate tectonics is the motion of the surface plate.Mobility is a measure of this motion across a specified surface.It is simply defined by where v rms,surf and v rms,global are the RMS velocities averaged over the surface and over the whole domain, respectively (Tackley, 2000).Values of M > 1 (e.g.M = 1.3) indicate plate-like motion, whereas a value of M = 1.0 would be measured for isoviscous mantle convection with no stiff top boundary layer.x 3 z 3 − z(x) dx a x 1 = 0, x 2 = 400 km, x 3 = x(z = z 3 ) and x 4 = x TR .b z 1 = 0 km is the sea level, z 2 = z min (x > x FB ), z 4 = z mean (x < x 2 ) and z 3 = z 4 − 2(z max − z min ) for x < x 2 .c See graphical representation in Fig. 2.

Residual mantle temperature
Extracting the residual temperature in a given domain is often very insightful, as geodynamic flows are often strongly temperature dependent and hence driven by the temperature anomalies.The residual mantle temperature can be defined in different ways.Most commonly, residual mantle temperature is defined as the temperature anomaly after normalising the temperature at each depth to the corresponding global horizontal mean (e.g.Labrosse, 2002).This definition is, however, less useful to distinguish local anomalies in global wide-aspect-ratio models.To distinguish a local anomaly, the field has to be normed to a regional instead of a global mean.Therefore, I outline here different ways of defining a residual mantle temperature.
1. Horizontal residual: the field is normalised to the horizontal mean at each depth.
2. Global residual: the field is normalised to the global mean.
3. Horizontal-band residual: the field is normalised to the mean across a finitely thick, horizontal band at each depth.4. Regional residual: the field is normalised to the regional mean surrounding each point in space.
The different definitions of mantle temperature anomalies can further be used to diagnose upwellings and downwellings and track mantle plumes (see Sect. 2.3.1 and 2.3.2).

Surface-plate and slab diagnostics
A first suite of more specific diagnostics focuses on the structure and dynamics of the top boundary layer that forms at and mostly operates within the uppermost part of a planet's convecting mantle.A variety of simple, general diagnostics, such as mean plate thickness, maximum plate stress and strain rate, upper-mantle temperature, density and viscosity, are extracted from only a few key parameter fields.Major and more complex diagnostics are outlined below and focus on long-wavelength surface topography as well as plate, plateboundary and slab dynamics.

Regional subduction topography
In the case of ongoing ocean-plate tectonics, the topography above a subduction zone typically displays a suite of characteristic regional features.These regional topographic characteristics can be individually measured using automated diagnostic algorithms as explained in Fig. 2 and Table 1.Specific algorithms are available for the following regional characteristics.
1.The viscous fore-bulge (or outer rise) is the upward deflection of the subducting plate outboard of the subduction trench.This transient, viscous uplift is mainly caused and controlled by the downward bending of the plate at the subduction zone into the low-viscosity mantle (de Bremaecker, 1977;Crameri et al., 2017).Here, the fore-bulge height, FB, is defined as the difference between its maximum elevation, FB 0 , and the minimum plate surface elevation, z 2 , at its side away from the subduction trench (see Fig. 2).
2. The subduction trench is the downward deflection that is located at the plate boundary precisely indicating the interface at Earth's surface between the upper and lower plate.Studies like Zhong and Gurnis (1994) and Crameri et al. (2017) suggest that it is likely of dynamic origin and continuously controlled by a multitude of factors.Here, the trench depth, TR, is defined by the maximum depth of the depression (TR 0 ) relative to the model's sea level (z 1 = 0 km).
3. The island arc (or volcanic arc) is the collisional high caused by horizontal plate compression.The plate strength has a major control on its resulting elevation (Crameri et al., 2017) and apart from its dynamic origin, it is often strongly affected by volcanism (Karig, 1971).Here, the island-arc height, IA, is defined by its maximum elevation (IA 0 ) relative to the characteristic upper-plate elevation, z 4 , outboard of the arc.4. The back-arc depression (or basin) is the upper-plate depression following the island arc further away from the trench.The depression's origin is often likely twofold and a combination of upper-plate extension or even spreading (Karig, 1971) and dynamic coupling (via the mantle wedge) with the sinking slab below (Crameri et al., 2017).Here, the back-arc depression depth, BAD, is defined as the difference of the maximum depression on the upper plate, BAD 0 , with the characteristic undeflected upper-plate elevation, z 4 , outboard of the back arc.This is necessary due to variable, isostatically induced elevation differences between the upper and lower plate.
Additional diagnostics of regional surface topography at subduction zones can be measured.The maximum horizontal extent of the back-arc depression, BAD extent , can, for example, be defined by the distance between the trench (x 5 ) and the far end of the back-arc deflection, x 3 , away from the trench.The latter point can be approximated by scanning the deflection for the one point that is still lower than twice the maximum vertical variation occurring in the undeflected reference upper-plate portion.Another diagnostic is given by the volume of the back-arc depression, which on a 2-D plane corresponds to a vertical area.It can be approximated by the area constructed by the maximum basin horizontal (x 3 -x 4 ) and vertical extent (z 4 -BAD 0 ).Yet another useful diagnostic is the tilt of the upper plate towards the subduction trench, as it has been shown to vary significantly during subduction evolution (Crameri and Lithgow-Bertelloni, 2017).The tilt can be tracked by a measure taken at a certain critical distance away from the trench.The tilt angle measurement is significantly improved (i.e. made more robust and less fluctuating) by taking a mean of multiple measurements taken just next to each other.
It is generally good practice to normalise the vertical amplitude of the characteristic topographic points mentioned above to a characteristic plate thickness defined, for example, by the thermal lithosphere thickness.This allows for scaling the obtained results to systems with different Rayleigh numbers and hence different plate thicknesses.

Topography components
In addition to the absolute surface topography, its isostatic and the remaining, non-isostatic residual components are useful to understand the various and diverse sources of longwavelength surface elevation.To derive these two topography components, the plate thickness, d p , has to be tracked along the horizontal extent of the model (see Sect. 2.2.3).As introduced in Crameri et al. (2017), the isostatic topography component for each vertical column along the model extent can be calculated using the base of the plate (e.g. as defined by a 1700 K isotherm) as compensation depth.Depending on whether the plate is denser or lighter than the mantle, it is then given by  with ρ p as the vertically averaged plate density at each horizontal point in space, ρ m the horizontal mean upper-mantle density just below the plate and away from any sinking slab, ρ air the air density, d LAB the variable thermal lithosphereasthenosphere boundary (LAB) depth as defined, for example, by a 1700 K isotherm along the model, and d p the variable plate thickness at each horizontal point in space that includes surface topography and so is the thickness between the rock-air interface and the base of the plate.Crucially, the resulting isostatic topography component has to be normalised throughout the model to produce a mean topography that corresponds to the sea level according to with z topo,iso as the mean of the model-wide isostatic topography component and x as the horizontal coordinate.This therefore ensures that both the conservation of volume in an incompressible mantle and the coherence of the surface plate are accounted for.
The residual topography component corresponds then simply to the non-isostatic part of the topography and is given by with z topo,total as the surface topography of the model.The residual topography component can be considered as the part of the topography that cannot be explained by the plate's iso-  static dynamics.It is important to point out that it can, however, be caused not only by dynamic sources from within the convecting mantle below, but also from sources (e.g.horizontal tectonic forces) within the plate itself.

Plate thickness and plate-core depth
An estimate of the plate thickness can be derived using temperature isotherms (generally the 1600 K isotherm) according to the thermal-lithosphere definition.A maximum depth limits the derived plate thickness at subduction zones, where the plate and hence the isotherms bend downward into the mantle.By definition, this approach simultaneously allows us to track the lithosphere-asthenosphere boundary and its topography.
Another more widely applicable option is to derive the plate thickness automatically from individual radial profiles of temperature (or viscosity).Since a typical plate cools via diffusion, the local temperature profile (nearly) linearly decreases throughout the boundary layer (from top to bottom) and becomes (nearly) isothermal in the convecting mantle below.The transition from diffusion-dominated to convectiondominated cooling is marked by a characteristic sharp bend in the profile.The bend in the profile can be tracked automatically by scanning for it from top to bottom.This gives an estimate of the local plate thickness, or a mean plate thickness, when performed on a horizontal mean (i.e.root mean square) temperature profile.The plate-core depth lies then simply about half-way down between the surface and the plate base.This depth is particularly useful to robustly track plate velocities and other similar characteristics.

Stagnant-lid diagnostics
Checking for a stagnant lid (i.e. the absence of subduction) can be done using a combination of tests.Key indicators for a stagnant lid are a low overall plate-thickness variation, a low maximum plate thickness throughout the model domain (i.e.no incipient subduction zone), the lack of cold lithosphere in the upper mantle (i.e.no mature subduction zone) and a low relative motion within the surface plate (a stagnant lid is mostly rigid).
If a stagnant lid is present, various diagnostics can be performed.Stagnant-lid diagnostics include, for example, the maximum plate stress, the maximum depth of plastic failure (i.e.brittle or ductile yielding) and the lithosphereasthenosphere boundary topography.The two latter physical complexities have been shown to be a reliable indicator of subduction initiation (Crameri and Tackley, 2016).

Plate-boundary tracking
Using advanced routines, converging and diverging plate boundaries can be tracked robustly and fully automatically.A basic, first-order plate-boundary tracking can already be achieved by diagnosing temperature and velocity fields only.A plate-boundary tracking routine can, however, be improved when additional parameter fields are considered (see Table 3).
To robustly find converging and diverging plate boundaries (i.e.subduction trenches and spreading ridges) within a 2-D model at a given point in time, a quite elaborated procedure is necessary (see Fig. 1).Several checks need to be performed initially to exclude the presence of a stagnant lid and thus the absence of subduction or spreading zones (see Sect. 2.2.4).If a stagnant lid is ruled out, potential plate boundaries can be located.Finding plate boundaries robustly necessitates finding all sharp and diffuse plate boundaries, finding the exact location at the very plate surface (e.g. the outcropping of the subduction fault) while accounting for plate topography and preventing multiple tracking of one and the same plate boundary.
The first derivative of the plate velocity indicates both convergent and divergent plate boundaries quite robustly by clearly detectable peaks.However, plate boundaries can sometimes become quite diffuse.Spreading ridges in coarsely discretised numerical models are, for example, often diffuse due to various secondary ridges forming in the young and weak plate around the main ridge.To find both Figure 2. Deriving the regional topographic characteristics at a subduction zone, which are from the subducting plate (right) towards the overriding plate (left) viscous fore-bulge (FB), subduction trench (TR), island arc (IA) and back-arc depression (BAD).Additional diagnostics explained here are the shallow-depth slab dip, θ , taken at the depth z 5 , and the bending radius R B derived by fitting a circle to the point of maximal plate bending.See text and Table 1 for more details (figure reproduced from Crameri et al., 2017).sharp boundaries accurately and not miss diffuse boundaries, critical values for the velocity change across plates and the horizontal distance between the measurement points need to be quite conservative (i.e. a low-velocity change measured far apart from the boundary).To not lose accuracy on the location, best practice is to start with a conservative check and gradually make critical values more restrictive to the point shortly before the boundary is not detected any longer.
The plate velocity at different depth levels inside the plate (e.g. at the plate surface and plate core; see Sect.2.2.3) can be used to distinguish between shallow and deep plate boundaries, and their horizontal offset additionally serves as an indicator of the polarity of the subduction system.

Plate velocities
By knowing the location of plate boundaries, the velocities of the converging plates can be diagnosed.Plate velocity is best measured in the cold, strong core of the plate(s) (see Sect. 2.2.3).The surface-plate RMS velocity is, for example, measured along the plate core.The individual plate velocities at convergent and divergent boundaries are then additionally most representative when measured close but not too close to the trench and ridge, respectively.These measurements then yield convergence and divergence rates.Once the subduction polarity is found (see Sect. 2.2.8), the upper and lower plate can be distinguished, and the trench velocity can be measured simply by the velocity of the upper plate just next to the trench.For comparison to the effective trench retreat, a theoretical trench velocity can be calculated during any given subduction phase.Assuming a non-deformable slab (which is not always the case), the current trench velocity can be approximated by where v Stokes is the vertical sinking velocity of the slab and θ is the shallow-depth slab dip angle (e.g.Capitanio et al., 2007).

Plate age
A theoretical subducting plate age, a Plate,theoretic , at the subduction trench can be derived by making use of the theory of the half-space cooling model as where h LP is the thickness of the subducting plate at the trench and κ is the thermal diffusivity (Turcotte and Schubert, 2014).

Slab-dynamics diagnostics
The current slab-tip depth can simply be tracked using the deepest point in a temperature contour outlining the cold material of the sinking plate.It can be refined by taking not the deepest point within the contour, but the point inside the contour that is the farthest away from the trench, as the slab tip usually is.
The subduction polarity of an asymmetric subduction zone can be found by checking for cold material in the uppermost mantle at two depth levels.A subduction polarity should only be given if the coldest spots at these two depth levels are shifted horizontally with respect to one another over a certain threshold.
If the subduction polarity is found, the shallow-depth slab dip, θ , can additionally be measured between these two depth levels, z 5A = 1.25d p and z 5B = 7/4z 2 , where d p is the mean lithosphere thickness away from the subduction zone.This results in a measuring depth, z 5 , for the slab dip angle located just below the lithosphere (see Fig. 2) according to The minimum, maximum and mean slab viscosity can be automatically derived using a specified area inside the slab.This area can be set in such a way that it spans the points outlining the sinking plate that lie within a square box around the point of maximal slab bending (see Sect. 2.2.9).If for some reason this method fails, the slab viscosity can instead be derived by spotting the coldest portion inside a sinking slab at a depth level below the surface plate.There, a small region surrounding the coldest spot is checked for the mean, the minimum and maximum viscosity.Depending on the task, either the mean value (e.g. for the indication of absolute slab viscosity) or the minimum value is more useful (e.g. for the calculation of bending dissipation; see Sect.2.2.10).
The slab-mantle viscosity difference can be derived from the minimum viscosity found in the slab (as outlined above) and the viscosity found in the surrounding upper mantle.The upper-mantle viscosity can be approximated by the horizontal median of viscosity found in the upper mantle below the plate.As such the anomalously high viscosity of the few cold regions (i.e. the sinking slabs) are weighted much less than the viscosity of the more common hot mantle surrounding them.The slab-mantle viscosity difference, η LA , is then the difference between slab and mantle viscosity.
The slab sinking rate, v Slab , is simply extracted using the vertical velocity measured at z 5B .An approximation of the total amount of water transported to the mantle can be calculated using the scaling law of Magni et al. (2014).The theoretical total water retention, W (kg m −2 ), of the slab is then where v Slab (cm a −1 ) is the sinking velocity of the slab, a Slab = a Plate,theoretic (Ma) is the age of the slab and T Mantle ( • C) is the potential mantle temperature.

Plate bending
The subduction bending radius, R B , can be fully automatically calculated using the current plate geometry.A spline method with temperature contours can be used to outline the subducting plate shape (see Petersen et al., 2017, for other methods).A certain temperature threshold thereby extracts cold plate portions including a sinking slab if present.Extracting the lowermost points in every horizontal column along the numerical grid then marks the down-going plate and simultaneously removes any non-subducting plate portions.The resulting band of points then needs to be fitted using a smoothing spline.This provides a line to represent the down-going plate geometry.The local curvature, R curv , along the line over its whole width is given by with dx and dy as the incremental spatial change in the line at the horizontal location, x, in the x and z direction, respectively.The minimum plate-bending radius at the kink of the subduction zone, R B (see Fig. 2), corresponds then simply to the maximum bending found along the considered plate portion and is thus

Viscous bending dissipation
The viscous bending dissipation within the lithosphere at a subduction zone can be calculated using the above definition of the subduction bending radius.Conrad and Hager (1999) provide an approximation for the viscous bending dissipation for the case of a purely viscous plate.Including the additions outlined in Buffett (2006) to consider a visco-plastic plate and neglecting dissipation in the subduction channel, the viscous dissipation is finally given by with C l = 1/6 as a constant, v p as the lower-plate velocity, d s the slab thickness and σ y,p as the maximum yield stress within the bending portion of the plate.The lithospheric bending dissipation, φ vd L , can be normalised using a certain value, φ vd L,char in W m −1 , to a normed value

Mantle-flow diagnostics
This second suite of more specific diagnostics focuses on the mantle dynamics operating in a planet's interior.The major diagnostics are outlined below and focus on the discrimination between active and passive upwellings and downwellings and the tracking of active mantle plumes.

Upwellings and downwellings
As highlighted in Fig. 3b, diagnosing thermally passive and active upwellings and downwellings is useful to distinguish whether material in a certain region is rising or sinking.Moreover, information can be extracted about whether this regional flow is thermally self-driven (i.e.active) or induced (i.e.passive).This task can be achieved using the actual flow field (i.e.velocity) and one kind of regional mantle residual temperature (see Sect. 2.1.3).The latter provides information on whether a patch of material is more or less buoyant than its direct surroundings.

Mantle-plume tracking
Mantle plumes can be tracked using the information outlined above about active upwellings and downwellings and a plume-tracking algorithm based on the one described in Labrosse (2002).Mantle plumes are here defined (and tracked) as hot or cold upwellings or downwellings that emerge and are connected to either of the two (hot or cold) boundary layers of the convecting flow (i.e.mantle).Hot and cold temperatures are marked as anomalies when the temperature at a given location exceeds a certain threshold (f hot or f cold ) in the range between the horizontally averaged temperature (T mean ) and the maximum (T max ) or minimum (T min ) temperature at a given depth level (z) according to where for hot anomalies, a threshold of f hot = 1 (or f hot = 0) defines anomalies that are 100 % (or 0 %) hotter than the horizontal average in the possible range between the mean and the maximum temperature.Once all the anomalies are located, they are checked for their connection to their respective boundary layer by a classical image processing procedure searching for connected pixels in a matrix (e.g.Kovesi, 2000).The hot and cold thresholds can thereby be chosen separately.

Field-variation diagnostics
Histogram plots of the variation in a specific field along a specified horizontal surface can prove very useful to get insight into the statistical physical behaviour of the state of a geodynamic system like the Earth's mantle (e.g.Fig. S1 in the Supplement).These statistics can provide mean and median values as well as the standard deviation of the parameter field under consideration.This kind of diagnostic can not only be applied to a perfectly horizontal surface, but also along a vertically slightly variable surface like along the surface-plate core.This enables, for example, the diagnostics of the strain-rate distribution within the surface plate(s).

Scientific visualisation
It is important to note that the key purpose of scientific visualisation is not about making data look pretty or entertaining; it is about creating comprehension and about delivering insight (e.g.Tufte et al., 1998).Unfortunately, there is not one single right way of visualising scientific data, but there are certain important ways to improve the presentation of scientific data and -crucially -there are several severe pitfalls that have to be prevented when doing so (e.g.Rougier et al., 2014).Here, I outline the unpleasant implications when using a non-scientific colour map (like rainbow), introduce a novel set of reliable, scientifically tested colour maps and outline additional ways to improve the positive impact of scientific figures.

Colours
Displaying and even printing scientific figures in colour has become standard for most journals.Random or even phys- ically based colour schemes that disregard the human eye's uneven colour perception (or its most-common mutations) most likely have multiple serious drawbacks.Therefore, they should be fully avoided by the science community (including authors and journals).Unfortunately, such colour maps are still widely and frequently used, even in high-impact scientific publications.I will therefore outline the most severe problems of unscientific colour schemes below and provide a novel set of ready-to-use, scientific quality alternatives.

Unscientific colour schemes
Using fancy colour schemes incorporating the whole colour range is appealing: they look peppy and have a lively appearance with their varying contrasts and multiple colours.Additionally, their main representative, the rainbow (a.k.a.jet) colour map, was and in some visualisation programs still is the default.These unscientific colour maps are thus widely blindly applied by authors while rarely criticised by reviewers and editors.
A colour scheme is unscientific as soon as it features one of the following aspects.
1.Both red and green colours: various forms of colour deficiencies can exist in human eyes, some of which make, for example, green and red undistinguishable.
2. Multiple different colours with similar lightness: colours like red, green and blue with similar lightness cannot be readily ordered against one another (e.g. from low to high values).
3. No gradual lightness gradient: a lack of a constant lightness gradient (from light to dark or vice versa) makes a colour map unreadable when printed in black and white.
4. No perceptual uniformity: perceptually non-uniform colour maps cause different parts of the data to be weighted differently (see Fig. 5).The green-cyan part of the colour spectrum has a lower contrast to the human eye than the yellow-red part.The greenish colours therefore hide low-amplitude data variation compared to reddish colours that amplify them.
Whether colours make a figure confusing or even unreadable to colour-blind people or whether colours introduce dramatic visual artefacts (see e.g.Light and Bartlein, 2004;Borland and Ii, 2007, and Fig. 5), the figure surely cannot be considered suitable for science any longer.In fact, applying the most moderate form of the commonly used rainbow colour map introduces an estimated error (calculated from the change in CIE76 lightness along the colour map) to the represented underlying data of up to the staggering amount of 7.5 % (see Fig. 5).
The misleading widespread view that visualisation is not an important part of science and thus also not worth spending time and money on (e.g. for external visualisation expertise) is therefore fundamentally wrong.(Crameri, 2018) used in STAGLAB is freely available online (http://www.fabiocrameri.ch/colourmaps,last access: 25 June 2018) in all major data formats and including the individual colour-map diagnostics that are performed using MATLAB software by Peter Kovesi (Kovesi, 2017).The colour maps, including devon, davos, oslo and broc, are all perceptually uniform and prevent distorting the data visually.The suite additionally includes lapaz and roma, scientific versions to replace the widely used unscientific rainbow and seis colour maps, and oleron, a scientific Earth-topography colour map to be used zero centred at sea level.

Scientific colour schemes
Useful and clear guidelines on what colour schemes to use and how to judge a given colour map have already been provided elsewhere in detail (e.g.Healey, 1996;Kelleher and Wagener, 2011;Silva et al., 2011).The most important points for choosing a suitable colour scheme can be summarised as follows.
1. Perceptual order: the different colours of a colour bar should be perceived as having the same order as the represented numerical values.A temperature scale should, for example, be represented by using the notions of cold and warm colours (and their proportional mixtures).
2. Uniformity and representative distance: two colours should convey the distance of numerical value between them, and colours representing equally differing values should also seem equally different.Clearly separated values must additionally be represented by clearly dis-  2015).CIE76 lightness variations along the standard rainbow colour bar (Kovesi, 2017) introduce a significant error of 7.5 % across the colour bar range to the underlying data, while the error of the perceptually uniform colour map davos (Crameri, 2018) is only 1.6 % locally.The impact of the much higher visual error is dramatic and can be seen by a slight shift of colour bar limits: while the same data look factually the same with the scientific davos colour map, the unscientific rainbow colour map introduces strong artificial boundaries and distortion to the underlying data.
tinguishable colours, and closer values must be represented by colours perceived to be closer.

3.
No artificial boundaries: if there are no boundaries in the represented numerical values, the colour scale should not create boundaries, but should rather look continuous.
4. Separation of bivariate information: two (or more) parameter fields represented in the same figure should be clearly separated by two clearly different colour schemes with no repeated colours.
A suitable scientific colour map makes a figure more intuitive and easier to understand and does not distort the underlying data.This can be done by adjusting a colour map to the parameter's nature and/or to the kind of parameter visualisation.Adjustment to a parameter's nature might be to plot temperature with a blue to red colour scheme as it is intuitively linked to a human's conception of hot (red) and cold (blue) as mentioned above.Adjustment to the kind of visualisation means, for example, that low to high values might be represented by varying the colours from white (low) to red (high), while a variation between two similar colours (e.g.blue and green) might clarify displaying positive and negative values around a given zero level (e.g.white).
Here, I introduce a novel set of fully scientific, perceptually uniform colour schemes (Crameri, 2018).Crucially, these intuitive colour schemes do not distort the data; they are readable even by people with a colour-vision deficiency; they are readable after being printed in black and white; they add no significant error to the underlying data; the scientific diagnosis of each individual colour map is provided alongside the colour map; they are available in all of the most common data formats; and they are freely available.Included in the perceptually uniform and visually appealing suite of novel colour maps shown in Fig. 4  that needs to be used zero centred at sea level.All of these colour maps perform significantly better in scientific tests for, for example, uniform perception and local variations in colour contrast (Kovesi, 2015(Kovesi, , 2017)).The full suite of these novel perceptually uniform colour maps including their individual scientific test results is freely available from http://www.fabiocrameri.ch/colourmaps(last access: 25 June 2018).

Figure design
Good scientific visualisation is accurate and clear.Plots need to show all, and only, the relevant aspects like the data itself, axis ticks and labelling, clear colour bars and instructive titles.Removing unnecessary clutter like duplicated labelling or unnecessary graphical forms and choosing a clear sansserif font and a lighter (e.g.grey) colouring of axis labels does shift the visual focus onto the most important part of the plot: the data.A visual focus is especially useful when plotting multiple subplots next to each other.Subtle visual guides (that for example group subplots or relate them to each other) can additionally help the reader to understand the figure by improving its clarity.Magnifying panels can be added to provide a view on the details while not losing the big picture.Adjusting the background colour is useful to adjust the figure to fit seamlessly into its surroundings whether this is a white conference poster or a dark presentation slide.Such visual refinements improve scientific figures and allow them to convey the precious scientific data accurately in an easily understandable and visually appealing manner.
4 StagLab: the software STAGLAB is the all-in-one, easy-to-use software that combines all geodynamic diagnostics (see Sect. 2) and the stateof-the-art scientific visualisation (see Sect. 3) outlined above and makes them accessible to everyone, including students and inexperienced modellers.Here, I provide an overview of various aspects of the software, including some of the powerful features, its thoughtful code design and the helpful external contributions to it, and outline how easily it can be used.

Supported model data
STAGLAB is optimised for the geodynamic finite-difference code StagYY (Tackley, 2008) but is easily made compatible with other codes.STAGLAB has, for example, already been used with 2-D output from the finite-element code Fluidity (Davies et al., 2011) that, in contrast, employs an unstructured numerical discretisation (see Supplement Fig. S2).
The data input for parameter fields simply has to be imported to MATLAB and subsequently needs to be adjusted to match the format outlined in STAGLAB's input routine f _readOther.

StagLab features
The STAGLAB user receives constant support from a builtin, friendly artificially intelligent operator, fAIo.It facilitates finding input data across unspecified folders and file numbers, facilitates saving figures and movies while preventing overwrites of existing files, prevents unnecessary errors during execution, ensures unbroken forward compatibility of previously used parameter files and keeps STAGLAB itself up to date.In the following, I list the key features that lie at the feet of such a streamlined user experience that make STAGLAB truly boost existing geodynamic models and the research behind them.

Dimensional scaling
For performance reasons, numerical models are often run using non-dimensional numbers.If a geodynamic code can be run in both dimensional and non-dimensional mode, STAGLAB will account for that by checking if the data files are dimensional or not and adjusts the dimensions fully automatically.Moreover, STAGLAB offers the possibility to convert non-dimensional values into sensible dimensional numbers using its Dimensional mode (SWITCH.DimensionalMode).A given set of dimensionalisation parameters can be defined in the function file f _Dimensions and then referred to by setting the corresponding flag to the variable IN.Parameter in the parameter files.

Automated geodynamic diagnostics
STAGLAB's incredible diagnostic capabilities (see Sect. 2) decipher geodynamic models within a fleetingly brief time span.In fact, performing the whole suite of plate-tectonics diagnostics listed in Table 2 with STAGLAB 3.0 within less than 2.2 s for a high-resolution (512 × 256) 2-D model on a power-efficient laptop (1.3 GHz processor; 8 GB RAM) is record breaking if not revolutionary.During this fleetingly short time period, real-time output is provided listing the most important model details and diagnostics and, in the case of problems, instructive warnings and error messages.
The real power of STAGLAB's geodynamic diagnostics lies, however, not only in its speed itself, but rather in the combination of both speed and robustness.Given the enormous variety occurring in a geodynamic system and related numerical models, STAGLAB's diagnostic routines have been trained excessively to become incredibly robust (see e.g.Fig. 6).

Plot and figure design
STAGLAB determines the model geometry automatically from the structure of the input data.It can handle 2-D Cartesian and 2-D cylindrical geometries with different aspect ratios as well as 3-D Cartesian and even 3-D spherical geome- tries.The 3-D data can be represented with both 2-D slices (in any direction normal to a side boundary) or with 3-D isosurfaces.For Cartesian models these two methods can even be combined into one figure.For 3-D spherical model data, STAGLAB additionally offers a large variety of map projections.
STAGLAB's visualisation routine is trimmed for accuracy, clarity and simplicity.Plots produced with STAGLAB show all the relevant data like axis ticks and labelling, clear colour bars and instructive titles.On top of that, STAGLAB offers two plotting modes, Analysis mode (SWITCH.AnalysisMode) to carefully examine the data and Publication mode (set as default) to clearly present the data (Fig. 9).Analysis mode offers detailed information and has refined axis ticks and more labels.Publication mode shifts the focus from the information surrounding the data to the data itself and labels subplots automatically to be easily referred to.This is achieved by removing unnecessary (e.g.duplicated) labelling, the choice of a well-readable and larger sans-serif font and a less distracting grey colouring (see Sect. 3.2).
STAGLAB makes plotting multiple subplots straightforward, while keeping the figure clear and focused.To compare different experiments, STAGLAB fully automatically adds subtle visual background guides (SWITCH.BackgroundGuides), areas that visually combine subplots of the same experiment, while separating them visually from the other experiments (e.g.Fig. 10).A similar example is the Time-evolution mode (SWITCH.TimeEvolutionMode) that adds time arrows to highlight the temporal evolution from one subplot to another (Fig. 6).To highlight smaller areas in a subplot and enlarge important details, magnifier panels (SWITCH.Magnifier) can be added to plots (e.g.Fig. 7c).STAGLAB has a convenient option to use a discrete colour map (SWITCH.DescreteColormap) that can help to outline regions of similar values more clearly (Fig. 8).Finally, a stunning Dark mode (SWITCH.ColorMode) can be switched on, which inverts relevant colours of the figure to be presented on a black background (Fig. 7).Dark mode is particularly useful to display STAGLAB figures on screens and via projectors.It is worth mentioning that all of these options and modes can be individually put into action with one single, simple switch.Moreover, options like these crucially enable STAGLAB figures to convey scientific data accurately and clearly, while still being visually appealing to the reader.

Plot types
Apart from plotting various parameter fields, STAGLAB can further produce various useful spatial or temporal graph plots (with data extracted from parameter fields) and a suite of special plots.The graph plots either show the horizontal plate velocity (PLOT.PlateVelocity; only for 2-D) or a fieldcontour topography (PLOT.PlateBaseTopography; e.g. for the lithosphere-asthenosphere boundary).Additionally, StagLab data graphs (PLOT.CustomGraph) can be used to visualise a large variety of previously processed and www.geosci-model-dev.net/11/2541/2018/Geosci.Model Dev., 11, 2541-2562, 2018 saved geodynamic diagnostics against each other or time (e.g.trench location, plate velocities, etc.; see Sect.4.2.5).Moreover, STAGLAB produces a variety of additional useful plot types that are listed below.
1. Grid (PLOT.Grid): the grid can be plotted separately or as an addition to a parameter field.This is particularly useful to check physical features against grid resolution or to highlight the grid's spatial variation.The grid lines can be plotted in actual resolution or coarser if the grid is too fine to be resolved with the given figure resolution.

Output files
STAGLAB produces post-processed data files in a variety of data formats.This might be useful to save a post-processed parameter field or simply to convert a field from one format to another.Available file formats are here .mat,.datand .txt.Particularly useful is the option to save a large variety of geodynamic diagnostics to data files.The diagnostics that can currently be saved are listed in Supplement Table S1.Having these data files is particularly useful to plot temporal (or other) graphs including some of the diagnostic data using the option to plot this STAGLAB data (see Sect. 4.2.4).

Software design
STAGLAB is written in MATLAB and compatible with all the latest MATLAB versions including MATLAB 8.4.0 (2014b) and newer (see Sect. 4.4.1 for more details).STAGLAB has its roots in STAGPLOT, a plotting routine introduced in Crameri (2013), and has been developed and extended further ever since with a few externally contributed routines.It now makes use of image processing and other advanced routines to provide highly accurate geodynamic diagnostics and scientific colour maps to maintain high accuracy throughout data visualisation.STAGLAB consists of an incredibly flexible parameter file that executes one of the three core applications STAGPLOT, STAGRPROF or STAGTIMEDAT.While STAGPLOT mainly handles two-or three-dimensional data of field variables, STAGRPROF and STAGTIMEDAT visualise data of radial profiles of horizontally averaged vari-ables (e.g.Supplement Fig. S3) and time-evolution graphs of globally averaged variables (e.g.Supplement Fig. S4), respectively.The parameter file allows users to change chosen parameters and switches to be different from the default setting.As such, a STAGLAB procedure is fully reproducible by saving a used parameter file, even after updates to the core routines as the parameter files in STAGLAB are forward compatible.Moreover, STAGLAB is extensively tested and heavily optimised for efficiency, which speeds up both computation and, crucially, also the research process as a whole (see Sect. 4.3.6).It also makes elaborated diagnostics and stateof-the-art visualisation easily accessible to everyone (including students and inexperienced modellers).And, last but not least, the geodynamic post-processing routines are fully open source.

Accuracy
STAGLAB is built for accuracy.Its automated diagnostics (see Sect. 2) offer more accuracy than interpretations by hand as they are based on quantitative measures rather than on visual impression.In order to maintain the high accuracy all the way into the data representation, perceptually uniform colour maps (see Sect. 3.1.2)are applied to additionally minimise the visually introduced error.

Flexibility
STAGLAB is built for flexibility.Its fully customisable parameter files contain a wealth of options (i.e.switches) to adjust post-processing and visualisation.Data from all model geometries, ranging from simple 2-D Cartesian to 3-D fully spherical, can be processed (see Sect. 4.2.4).Various plot additions like velocity arrows, streamlines or isolines can be added on top of other parameter fields (see Sect. 4.2.4).An automatic subplot arrangement enables direct comparisons between outputs from different time steps or even experiments.

Reproducibility
STAGLAB is built for reproducibility.The customised parameter files, from which any STAGLAB procedure is executed, are forward compatible and can be stored in and run from any possible directory.Given that a previously used parameter file is safely stored, it can therefore be reused at any time with any forthcoming version of STAGLAB and so reproduce previous geodynamic diagnostics and visualisations.Therefore, any work done with STAGLAB is and always will be fully reproducible.

Continuity
STAGLAB is built for continuity.Old parameter files are, on the one hand, always updated to be compatible with the latest version of STAGLAB, fully automatically and fully effortlessly.On the other hand, STAGLAB itself will be kept compatible with the latest versions of MATLAB.

Efficiency
STAGLAB is built for efficiency.A wealth of most common post-processing and visualisation tasks are just a click (and a few seconds or less) away.Crucially, STAGLAB diagnostic and visualisation tasks can be quickly reproduced when going through the common, improving iterations during the research progress.This keeps the time-consuming coding effort for producing, adjusting and updating the core software to a minimum.
In addition, STAGLAB itself is trimmed towards efficient computation.Time-consuming parts of the software are optimised by, for example, the vectorisation of loops as well as the efficient reading and transferring of large data structures.It can thus quickly handle the huge amount of data resulting from high-resolution 2-D and 3-D models.A useful performance switch is built in that allows users to reduce the file sizes of extremely large data files for quicker visualisation (SWITCH.ReduceSize).Another option, Quick mode (SWITCH.QuickMode), allows for an extra quick execution to ensure efficient post-processing throughout multiple time steps and files.

Reliability
STAGLAB is built for reliability.It has been coded carefully and tested extensively to prevent foreseeable problems.Although it is always optimised for the latest MATLAB version, STAGLAB's compatibility with previous MATLAB releases is also tested carefully.It is constantly being debugged and tested, and a growing number of users accelerates the exposure of hidden problems.Bug reports are very welcome  and should be sent to the author to ensure future reliable releases of STAGLAB.However, it has to be noted that since a fully bug-free software cannot be guaranteed, the scientific quality check always remains with the user.

Simplicity
The user interface (UI) design is a pivotal part of any software to simplify its application.STAGLAB therefore uses parameter files combining all the important switches in one place.The parameter file is clearly structured and can be reduced to only the routinely used switches thanks to having defaults to every individual switch.Using the parameter files ensures a streamlined usage even for complex figures.STAGLAB returns selected, clear and informative real-time output in MATLAB's terminal window during its operations.The display messaging allows users to monitor STAGLAB's progress, receive both scientific diagnostics and, in the case of problems, clear warnings thanks to its advanced error messaging system.The advanced error handling implemented in STAGLAB often prevents interruption during execution and displays warnings instead.An optional Verbose mode (SWITCH.Verbose) allows users to display more detailed information during the STAGLAB execution, which simplifies the debugging procedure.STAGLAB crucially simplifies the reading and saving of data files: an advanced file finder automatically checks for other possible file directories or the latest file number if the specified files are not found initially.Finally, STAGLAB itself is written carefully with a unified code structure and descriptive comments to simplify further code development.

Using StagLab
STAGLAB's software design aims towards an effortless user experience (see the Supplement user guide).Its application is simple and hence accessible to experienced numerical modellers as well as fresh beginners.Using STAGLAB is an incredibly efficient and motivating way, for students especially, to enter the world of numerical modelling, but also opens up new exciting doors for more experienced scientists through its incredible flexibility.

Prerequisites
STAGLAB is compatible with all three major operating systems running on Mac, Linux and Windows PC.It has been updated to the latest graphical improvements made to MAT-LAB and therefore performs best with the latest MATLAB version.Although it can function with older versions, MAT-LAB version 8.4.0 (i.e.MATLAB 2014b) or higher is necessary for many core functions and thus highly recommended.
Although a set of parameter-field data is necessary to use all the post-processing routines built into STAGLAB, it can already function with just a temperature field or any other single field that has to be processed.To complete more demanding tasks like plate-boundary tracking, more fields like temperature, velocity and topography are necessary.
STAGLAB performs a local directory search for the specified input file in case the file cannot be found directly.It also checks for the specific folder structure (based on StagYY's native folder structure) consisting of an image folder ".../<folder>/+im" containing the image files and an output folder ".../<folder>/+op" containing the (e.g.binary) output data files.Using default settings, STAGLAB will, in that case, look for binary data to import in "+op", while it will save the figures to "+im".This behaviour can, however, also be adjusted manually in the parameter file.

Download and installation
STAGLAB is available from http://www.fabiocrameri.ch/StagLab (last access: 25 June 2018).An included README file (http://www.fabiocrameri.ch/resources/README.pdf, last access: 25 June 2018) provides up-to-date, detailed instructions on installing and running STAGLAB.Once downloaded, it can then be installed by adding all of its files to the MATLAB search path.This can conveniently be done by running the included installation routine, f _I N ST ALL (Algorithm 1), which additionally checks for possible file duplications, or manually from within MATLAB (i.e.version 2014b or later).

Testing
STAGLAB has a built-in testing routine, f _T EST (Algorithm 2), to make sure it performs as expected on the current system and to highlight some of its capabilities.It performs a suite of automated tests for STAGLAB's core tasks and produces a suite of test figures from some included data files.
Algorithm 2 Testing command for StagLab.

Running StagLab
Once it has been added to the MATLAB search path, STAGLAB can be run via one of the provided parameter files (or a copy thereof).Included example parame-otherwise can mount up to a staggering 7.5 % across the displayed data range, as is shown to be the case for the commonly used rainbow colour map; this is an error that easily dominates in most data sets.The novel suite of scientific colour maps, including devon (http://www.fabiocrameri.ch/devon, last access: 25 June 2018), davos (http://www.fabiocrameri.ch/davos,last access: 25 June 2018), oslo (http: //www.fabiocrameri.ch/oslo, last access: 25 June 2018), broc (http://www.fabiocrameri.ch/broc, last access: 25 June 2018) and others, presents the precious scientific data undistorted and without excluding certain readers.In an unprecedented manner, the whole suite is made freely available including all scientific diagnostics and in all major data formats (see http://www.fabiocrameri.ch/visualisation, last access: 25 June 2018).
Thirdly, all geodynamic diagnostics and the state-of-theart scientific visualisation are packed into one single software package, STAGLAB, which is introduced here.STAGLAB is an easy-to-use MATLAB software that significantly facilitates post-processing and visualisation, the two crucially important aspects of research.It provides a powerful, fully automated and incredibly robust implementation of the geodynamic diagnostics outlined above.STAGLAB's efficiency turns laborious days of post-processing towards revealing hidden model secrets into an exciting and effortless 2.2 s of pure revelation (according to the measurement outlined in Sect.4.2.2).In the same breath, these revelations can be finely packed into a publication-ready figure or movie using fully reproducible, forward-compatible parameter files, while applying state-of-the-art visualisation techniques like its unique suite of scientifically tested, perceptually uniform colour maps.
STAGLAB is currently compatible with two of the widely used geodynamic codes, StagYY (Tackley, 2008) and Fluidity (Davies et al., 2011).In combination with StagYY output, it is capable of handling all different geometries (2-D and 3-D Cartesian, 2-D partial and full cylindrical, and 3-D spherical) and output (parameter fields, radial profiles and time-evolution data).With little effort, STAGLAB can also be adjusted to be compatible with additional mantle convection codes.The latest version of STAGLAB, freely available at http://www.fabiocrameri.ch/StagLab(last access: 25 June 2018), is a flexible, efficient, reliable and simple software that produces state-of-the-art, reproducible diagnostics and visualisation for upcoming and groundbreaking geodynamic models.

Figure 3 .
Figure 3. STAGLAB's plots showing a mantle convection model in 2-D cylindrical geometry for (a) temperature and (b) diagnosed active and passive upwellings and downwellings.

Figure 4 .
Figure 4.The novel set of scientific colour maps(Crameri, 2018) used in STAGLAB is freely available online (http://www.fabiocrameri.ch/colourmaps,last access: 25 June 2018) in all major data formats and including the individual colour-map diagnostics that are performed using MATLAB software by Peter Kovesi(Kovesi, 2017).The colour maps, including devon, davos, oslo and broc, are all perceptually uniform and prevent distorting the data visually.The suite additionally includes lapaz and roma, scientific versions to replace the widely used unscientific rainbow and seis colour maps, and oleron, a scientific Earth-topography colour map to be used zero centred at sea level.

Figure 5 .
Figure5.The visually introduced error to the underlying data.A common rainbow (a.k.a.j et) colour map hides low-contrast variations in the cyan-green part and amplifies it unnaturally elsewhere as highlighted by low-amplitude ripples in the two large colour bars (afterKovesi, 2015).CIE76 lightness variations along the standard rainbow colour bar(Kovesi, 2017) introduce a significant error of 7.5 % across the colour bar range to the underlying data, while the error of the perceptually uniform colour map davos(Crameri, 2018) is only 1.6 % locally.The impact of the much higher visual error is dramatic and can be seen by a slight shift of colour bar limits: while the same data look factually the same with the scientific davos colour map, the unscientific rainbow colour map introduces strong artificial boundaries and distortion to the underlying data.

Figure 6 .
Figure 6.Time evolution of 2-D Cartesian mantle convection using STAGLAB's Time-evolution mode highlighting the accuracy and robustness of the plate-boundary tracking algorithm even throughout a dramatic subduction-polarity reversal (figure adjusted from Crameri and Tackley, 2015).

Figure 7 .
Figure 7. STAGLAB's Dark mode and geodynamic diagnostics.While the default Light mode is intended for publications and posterpresentations, Dark mode is particularly useful for digital presentations on screens.STAGLAB's diagnostics highlighted here include slabtip depth (red cross), isostatic and residual topography components (dashed and dotted grey lines, respectively), plate-boundary tracking (white and grey rectangles), minimum plate-bending radius (red-white dashed circle) and resulting plate-bending dissipation.Geodynamic diagnostics like these foster a better understanding of complex models.

2.
Tracers (PLOT.Tracers): individual tracer information, like the tracer position and type, can be visualised in a separate plot.3. Streamfunction (PLOT.Streamfunction): there is an option to visualise the flow field in a separate plot or on top of any other field in the form of flow contours.This is a useful way to show the pattern and spatial extent of mantle flow cells.4. Streamlines (PLOT.Streamline): instantaneous streamlines can be plotted on top of another parameter field or as a separate plot to highlight the flow pattern of individual particles.5. Quiver (PLOT.Quiver): the option to plot distributed velocity arrows on top of any parameter field adds the possibility to highlight the flow direction and relative strength.The amount and scaling of velocity arrows can be adjusted manually if needed.

Figure 8 .
Figure 8. STAGLAB's continuous (a, b) versus its discrete (c, d) colour-map option that is useful to better convey the actual field value in a certain region shown for a model in partial-cylinder geometry.
4.3.1 External contributions to StagLab 3.0 STAGLAB calls the routine f _readStagY Y.m that was originally written by Boris Kaus to read StagYY's binary output directly into MATLAB.The routine f _Y Y toMap, which was originally written by Paul Tackley, is used to produce horizontal maps of fully spherical yin-yang data.The original routine f _readF luidity to read Fluidity data was provided by Fanny Garel.It further uses the figure-saving routine export_f ig, which was originally written by Oliver Woodford, the routine f lowf un originally written by Kirill K. Pankratov to derive the stream function, the routine MinV olEllipse by Nima Moshtagh to fit a minimumvolume ellipse around a point cloud, the routine plotboxpos by Kelly Kearney to derive the plot position more accurately, the routine hatchf ill2 originally developed by Neil Tandon to fill areas with a specific texture and a few routines including equalisecolourmap, sineramp2, normalise, show and strendswith developed by Peter Kovesi to provide the scientific colour-map diagnostics.

Figure 9 .
Figure 9. STAGLAB's Analysis mode (left) and Publication mode (right).While in Analysis mode as much detail as possible is provided to facilitate the examination of the data, the data (i.e. the key result) are put into focus and fully annotated in Publication mode.

Figure 10 .
Figure 10.STAGLAB's figure design routine turns overloaded data representation (left) into clear and focused visualisation (right) using various figure simplifications, while subtle visual guides make comparisons between a large number of models or model time steps easily understandable (figure adjusted from Crameri and Lithgow-Bertelloni, 2017).
source STAGLAB is open access and currently only requires a valid MATLAB license.Apart from the design routines, all STAGLAB files are additionally fully open source.Free usage and redistribution of STAGLAB and its individual routines and colour schemes, however, fall under the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/,last access: 25 June 2018).
-Necessity for diagnostics.-Improves the diagnostics without being a necessity.* -At the time of submission.