101 Geodynamic modelling: How to design, carry out, and interpret numerical studies

Geodynamic modelling provides a powerful tool to investigate processes in the Earth’s crust, mantle, and core that are not directly observable. However, numerical models are inherently subject to the assumptions and simplifications on which they are based. In order to use and review numerical modelling studies appropriately, one needs to be aware of the limitations of geodynamic modelling as well as its advantages. Here, we present a comprehensive, yet concise overview of the geodynamic modelling process applied to the solid Earth, from the choice of governing equations to numerical methods, 5 model setup, model interpretation, and the eventual communication of the model results. We highlight best practices and discuss their implementations including code verification, model validation, internal consistency checks, and software and data management. Thus, with this perspective, we encourage high-quality modelling studies, fair external interpretation, and sensible use of published work. We provide ample examples from lithosphere and mantle dynamics and point out synergies with related fields such as seismology, tectonophysics, geology, mineral physics, and geodesy. We clarify and consolidate 10 terminology across geodynamics and numerical modelling to set a standard for clear communication of modelling studies. All in all, this paper presents the basics of geodynamic modelling for first-time and experienced modellers, collaborators, and reviewers from diverse backgrounds to (re)gain a solid understanding of geodynamic modelling as a whole.


Introduction
The term geodynamics 1 combines the Greek word 'geo' meaning 'Earth' and the term 'dynamics' -a discipline of physics that 15 concerns itself with forces acting on a body and the subsequent motions they produce (Forte, 2011). Hence, geodynamics is the study of forces acting on the Earth, and the subsequent motion and deformation occurring in the Earth. Studying geodynamics contributes to our understanding of the evolution and current state of the Earth and other planets, and the natural hazards and resources on Earth.
1 See the glossary in the Supplementary Material for definitions of the bold words that occur throughout the text inside the Earth where different spatial and temporal scales play a role ( Figure 1). So, geodynamics considers brittle and ductile deformation and plate tectonic processes on the crustal-and lithospheric scale as well as convection processes in the Earth's mantle and core. Lithospheric and crustal deformation operates on scales of tens to hundreds of kilometres and thousands to millions of years to study both the folding of lithospheric plates and the faulting patterns accommodating them (Watts et al., 2013). Mantle convection encompasses individual flow features like mantle plumes and subducted slabs on scales of hundreds 25 of kilometres and hundred-thousands of years as well as whole mantle overturns on thousand kilometre and billion year scales (Androvandi et al., 2011;King, 2015). In the outer core, there are both small convection patterns on hundred-metre and hour scales and large-scale convection processes taking place over tens of kilometres and ten thousands of years (Bouffard et al., 2019). Inner core dynamics is believed to act on similar scales as mantle convection, although it is has not yet been established whether or not its Rayleigh number allows for convection (Deguen and Cardin, 2011). In addition to the large-scale dynamics 30 of the different spherical shells of the Earth, there are other processes that play a role and operate on yet other different spatial and temporal scales. For example, earthquakes and -on larger time scales earthquake cycles -facilitate the large-scale deformation of the lithosphere, but their nucleation process can occur on scales as little as milliseconds and millimetres (Stein and Wysession, 2009;Beroza and Kanamori, 2015). Similarly, surface processes, such as erosion can play an important role in lithosphere dynamics, although it takes place on a smaller spatial and temporal scale (Pelletier, 2008). For both lithosphere 35 and mantle dynamics, processes including magma dynamics and grain dynamics have been shown to be important, which both have smaller spatial and temporal scales (Solomatov and Reese, 2008). Hence, in order to fully understand geodynamics, many processes on different temporal and spatial scales need to be taken into account simultaneously. However, this is difficult to achieve using only observational or experimental studies.
Indeed, since the study of geodynamics is predominantly occupied with processes occurring below the surface of the Earth, 40 one of the challenges is the limited amount of direct observations in space and time (Gerya, 2019). Engineering limitations are responsible for a lack of direct observations of processes at depth, with, for example, the deepest borehole on Earth being a mere 12 km (e.g., Ganchin et al., 1998), less than 0.2% of the radius of the Earth (6371 km). In addition, the available rock record grows increasingly scarce when going back in time. Other disciplines, such as seismology, geology, and geodesy, relate to geodynamics as they study the (indirectly) observable expressions of geodynamic processes at the surface and at depth. 45 Disciplines like mineral physics and rock mechanics relate to geodynamics by studying the physical properties and flow laws of rocks involved in geodynamic processes.
To compensate for the limited amount of data in geodynamics, many studies have turned towards modelling. Roughly speaking, there are two main branches of modelling based on the tools that they use: analogue modelling, where real-world materials are used in a lab, and physical or mathematical modelling, where physical principles and computational methods 50 are used. Analogue models have many advantages of which the most straightforward ones are that they are automatically subject to all the physical laws on Earth and have unlimited resolution. However, it is non-trivial to scale the models to realworld applications that occur on a different temporal and spatial scale. In addition, it is hard to determine variables such as stress and strain across the entire domain at any given time. Physical models aim to describe processes through equations. For relatively simple processes, analytical solutions or semi-analytical models can be used (e.g., Lamb, 1879; Samuel and 55 Bercovici, 2006;Montési and Behn, 2007). However, when the equations and their dependencies become increasingly complex in their description of geodynamic processes, numerical methods become necessary to solve the set of equations. Hence, numerical modelling is the use of numerical methods to solve a set of physical equations describing a physical process.
Numerical models are powerful tools to study any desired geometry and time evolution with full control over the physics and access to all variables at any point in time and space in the model, but they come with their own set of limitations and 60 caveats. Despite recent developments in numerical methods and modelling studies, it remains difficult to combine diverse spatial and temporal scales in one numerical model. This is due to the fact that numerical models are increasingly more difficult to solve computationally when large contrasts (e.g., in space or time) between the physical properties of the variables are present. Therefore, different modelling approaches have been developed to tackle the different natures of the processes outlined in Figure 1. Ideally, these modelling approaches would be combined in the future. In this article, we focus on the complexity beyond analytical solutions can be compared against observations and results from other codes, but again cannot be proven to be correct or, indeed, the only possible model producing such results.
The models we are concerned with here are all forward models, where we obtain a model result by solving equations that describe the physics of a system. These results are a prediction of how the system behaves given its physical state which 105 afterwards can be compared to observations. On the other hand, inverse models start from an existing data set of observations and aim to determine the conditions responsible for producing the observed data set. A well-known example of these kinds of models are the tomographic models of the Earth's mantle, which determine the 3D seismic velocity structure of the mantle based on seismic data sets consisting of e.g., P-wave arrival times or full waveforms (Bijwaard and Spakman, 2000;Fichtner et al., 2013). The possibility of incorporating data and estimating the model uncertainties in inverse models has recently led to 110 an increasing amount of efforts combining both inverse and forward methods in geodynamic modelling. For example, adjoint methods (Reuber et al., 2018b(Reuber et al., , 2020, pattern recognition Atkins et al. (2016), and sequential data assimilation (Bocher et al., 2016(Bocher et al., , 2018 have been incorporated in lithosphere and mantle dynamics models. One of the most common ways of solving the inverse problem in geodynamics is to run many forward models and subsequently test which one fits the observations best (Baumann and Kaus, 2015). This is difficult, because individual forward models are computationally expensive which results 115 in a limited amount of forward models that can be run. Additionally, there are too many parameters in the forward models to invert for. Hence, even in inverse models of geodynamics, a thorough understanding of the forward model is required, which we focus on in this work. However, the prospect of estimating uncertainties in geodynamic modelling and the increasing inclusion of data make the combination of inverse and forward modelling an exciting avenue in geodynamics.

120
A scientific modelling study encompasses more than simply running a model as illustrated in Figure 2. First and foremost, the modeller needs to decide which natural process they want to study based on observations to increase the understanding of the Earth's system. The modeller then needs to choose the equations that describe the physics of the processes they are interested in, i.e., the physical model (Section 2). Here, we confine ourselves to the conservation equations of mass, momentum and energy (Section 2.1). In order to solve the physical model, a numerical model is needed, which solves the discretised versions of the 125 equations using a specific numerical scheme (Section 3). Before applying the code to a particular problem, the code needs to be verified to ensure that it solves the equations correctly (Section 4). Once it has been established that the code works as intended, the to-be-tackled geological problem needs to be simplified into a model setup with appropriate geometry, material properties, initial and boundary conditions and the correct modelling philosophy (Section 5). This model setup should then be validated in terms of robustness and against observations (Section 6). After the model has successfully run, the results of the simulation 130 need to be interpreted, visualised, and diagnosed while taking into account the assumptions and limitations that entered the modelling study in previous steps (Section 7). The final step of a modelling study is to communicate the results to peers as clearly and reproducible as possible (Sections 8,9). Each of the above-mentioned steps of a modelling study correspond to sections in this paper, making this a comprehensive guide to geodynamic modelling studies ( Figure 2). Figure 2. The geodynamic modelling procedure. A modelling study encompasses everything from the assemblage of both a physical (Section 2) and a numerical model (Section 3) based on a verified numerical code (Section 4), to the design of a simplified model setup based on a certain modelling philosophy (Section 5), the validation of the model through careful testing (Section 6), the unbiased analysis of the produced model output (Section 7), the oral, written, and graphical communication of the modelling approach and results (Section 8), and the management of both software and data (Section 9).
From seismology, we know that on short time scales the Earth deforms like an elastic medium. Our own experience tells us that when strong forces are applied to rocks, they will break (brittle failure). But from the observation that continents have moved over the history of the Earth, and from the postglacial rebound of regions like Scandinavia, we know that on long time scales the mantle strongly deforms internally and behaves like a highly viscous fluid Gerya, 2019). The time that controls which deformation behaviour is dominant is the Maxwell relaxation time (Ricard, 2007), which is on the 140 order of 450 years for the Earth's mantle . The rocks in the Earth's interior are not the only material that deforms in this way. In fact, many other materials show the same behaviour: solid ice, for example, also flows, causing the constant motion of glaciers (Ricard, 2007).
In the following, we will focus on problems that occur on large time scales on the order of thousands or millions of years (i.e., 450 years, the Maxwell time, see Figure 1). Accordingly, we will treat the Earth's interior as a highly viscous fluid. 145 We will also treat Earth materials as a continuum. This implies that, on a macroscopic scale, there are no mass-free voids or gaps (Gerya, 2019). Under these assumptions and keeping the large uncertainties of Earth's materials in mind, we can apply the concepts of fluid dynamics to model the Earth's interior.

Basic equations
We have discussed above that setting up a model includes choosing the equations that describe the physical process one wants 150 to study. In a fluid dynamics model, these governing equations usually consist of a mass balance equation, a force balance or momentum conservation equation, and an energy conservation equation. The solution to these equations states how the values of the unknown variables such as the material velocity, pressure, and temperature (i.e. the dependent variables) change in space, and how they evolve over time (i.e. when one or more of the known/independent variables change). Even though these governing equations are conservation equations, geodynamic models often consider them in a local rather than a global 155 context, i.e., material and energy flows in or out of a system, or an external force is applied. Additionally, the equations only consider specific types of energy, i.e., thermal energy, and not, for example the potential energy related to nuclear forces. This means that for any system considered -or, in other words, within a given volume -a property can change due to the transport of that property into or out of the system, and thermal energy may be generated or consumed through a conversion into other types of energy, e.g., through radioactive decay. This can be expressed using the basic principle: 160 Influx + Generation = Outflux + Accumulation + Consumption.
Depending on which physical processes are relevant, different terms may be included in the equations, which is a source of differences between geodynamic studies. Since each term describes a different phenomenon and may imply characteristic time and length scales, being aware of these relations is not only important for setting up the model, but also for interpreting the model results and comparing them to observations. Mathematical tools like scaling analysis or linear stability analysis can be 165 helpful to determine the relative importance of each of the terms (see also Section 7.2), and for determining which terms should Figure 3. The governing equations: conservation of mass (Eq. 2, Section 2.1.1), conservation of momentum (Eq. 5, Section 2.1.2), and conservation of energy (Eq. 6, Section 2.1.3) with different types of rheology (Section 2.2.1. ⇢ is the density, t is time, v the velocity vector, the stress tensor, g the gravitational acceleration vector, Cp the heat capacity, T the temperature, k the thermal conductivity, H a volumetric heat production term (e.g., due to radioactive decay) and the term S = S1 + S2 + S3 accounts for friction heating, adiabatic heating, and the release or consumption of latent heat (e.g., associated with phase changes), respectively. Note that the plastic rheology depicted here is the geodynamic approximation of brittle failure.
be included in any given model. In solid Earth geodynamic modelling, we are usually interested in problems occurring on large time scales on the order of thousands or millions of years ( Figure 1). Accordingly, we usually model the Earth's interior as a viscous fluid. The corresponding equations are the Stokes Equations, describing the motion of a highly viscousFluid driven in a gravitational field, and the energy balance written as a function of temperature ( Figure 3). These equations are partial 170 differential equations, describing the relation between the unknowns (velocity, pressure, and temperature) and their partial derivatives with respect to space and time. We look at each equation in more detail in the following sections.

Mass conservation
The first equation (2) describes the conservation of mass: 175 where ⇢ is the density, t is time, and v is the velocity vector. The conservation of mass (2) states that in any given volume of material, local changes in mass @⇢/@t can only occur when they are compensated by an equal net influx or outflux of mass r · (⇢v). These local changes in mass are caused by the expansion or contraction of material, which implies a change in density ⇢. This usually happens as a response to a change in external conditions, such as a change in temperature (thermal expansion/contraction) or pressure. Other specific examples in the Earth are phase transitions or chemical reactions. 180 Because the first term is the only one in the equation that explicitly includes a time-dependence, it also introduces a characteristic time scale into the model (for details, see Section 4 of Curbelo et al., 2019). It describes how fast regions where the pressure is bigger (or smaller) than the hydrostatic pressure can relax by causing the material to expand (or compress) against viscous forces. This is also called the viscous isentropic relaxation time scale (Curbelo et al., 2019). In the Earth, this time scale is of the order of a few hundred years for the upper mantle to a few tens of thousands of years for the lower mantle.

185
The time scale of viscous relaxation is usually shorter than that of convective processes, and is often shorter than the time scales a model is focused on. In addition, these local changes in mass are often quite small compared to the overall mass flux in the Earth. Accordingly, many geodynamic models do not include this term (see also Jarvis and McKenzie, 1980;Glatzmaier, 1988;Bercovici et al., 1992), and instead Eq.
This means that the net influx and outflux of mass in any given volume of material is zero. The density can still change, e.g., if material is advected into a region of higher or lower pressure (i.e., downwards or upwards within the Earth), but these changes in density are always associated with the motion of material to a new location. Using this approximation still takes into account the largest density changes. For example, for the Earth's mantle, density increases by approximately 65% from the Earth's surface to the core-mantle boundary.

195
In some geodynamic models, particularly the ones that span only a small depth range, even this kind of density change is small. For example, within the Earth's upper mantle, the average density changes by less than 20%. This is the basis for another approximation: assuming that material is incompressible (i.e., its density ⇢ is constant). In this case, (2) becomes Assuming incompressibility also has implications for another property of the materials in the model, the Poisson's ratio ⌫.

200
Poisson's ratio is a measure frequently used in seismology to describe the volume change of a material in response to loading, or in other words, how much a material under compression expands in the direction perpendicular to the direction of compression.
Incompressibility implies a Poisson's ratio of ⌫ = 0.5. When converting material properties from geodynamic models to wave speeds for seismological applications, incompressible models result in unrealistic infinite P-wave speeds unless a different value for the Poisson's ratio is assumed during the conversion (van Zelst et al., 2019).

Momentum conservation
Eq. (5) describes the conservation of momentum, or, in the way we use it here, a force balance: where is the stress tensor, ⇢ is the density, and g is gravitational acceleration vector. The conservation of momentum states that the sum of all forces acting on any parcel of material is zero. Specifically, the first term r · represents the net surface 210 forces and the second term the net body forces. In mantle convection and long-term tectonics models, the gravity force ⇢g is generally the only body force considered, and the gravitational acceleration g is often taken as a constant of 10 m/s 2 for global studies, because it varies by less than 10% over the whole mantle depth (Dziewonski and Anderson, 1981), whereas for lithospheric-scale studies g is taken to be 9.8 m/s 2 . Other body forces like electromagnetic or Coriolis forces are negligibly small compared to gravity. Conversely, these forces are very important for modelling the magnetic field in the outer core.

215
Hence, depending on the geodynamic problem, additional terms become relevant in the force balance.
The surface forces are expressed as the spatial derivatives of the stress r · . Stresses can be perpendicular to the surface of a given parcel of material like the pressure; they can act parallel to the surface; or they can point in any other arbitrary direction.
In addition, the forces acting on a surface may be different for each surface of a parcel of material. Accordingly, the stress can be expressed as a 3 ⇥ 3 tensor in 3 dimensions (3D), where each entry corresponds to the component of the force acting in one 220 of the three coordinate directions (x, y, z) on a surface oriented in any one of the three coordinate directions, giving a total of 3 ⇥ 3 = 9 entries. One of the most complex choices in the design of a geodynamic model is the relation between these stresses and the deformation of rocks, i.e., the rheology (Section 2.2.1).
Together, Eqs. (2) and (5) are often referred to as the Stokes Equations. In their more general form, called the Navier-Stokes equations, Eq. (5) contains an additional term ⇢Dv/Dt that governs inertia effects, describing the acceleration of 225 the material. In the Earth's mantle, material moves so slowly that inertia (and, consequently, the momentum of the material) does not have a big influence on the motion of the material. Consequently, this term is very small and usually neglected in geodynamic models. In other words, we look at the behaviour of the material on such a long time scale that from our point of view, material accelerates/decelerates almost immediately to its steady-state velocity. Conversely, for models of seismic wave propagation (which cover a much shorter time scale), it is crucial to include the inertia term because it introduces the physical 230 behaviour that allows the formation of seismic waves.
Dropping the inertia term means that the Stokes Equations (2) and (5) describe an instantaneous problem if the mass conservation is solved using one of the common approximations (3) or (4)). The time does not appear explicitly in these equations, so the solution for velocity and pressure at any point in time is independent of the history unless this is explicitly incorporated in the material properties (Section 2.2). The model evolution in time is incorporated through Eq. (6), which 235 describes the conservation of energy. For more details on the Stokes equations and a complete derivation, we refer the reader to Schubert et al. (2001).

Energy conservation
Eq. (6) describes the conservation of thermal energy, expressed as ⇢C p T :

240
where ⇢ is density, C p is the heat capacity, T is the temperature, t is time, v is the velocity vector, k is the thermal conductivity, H is a volumetric heat production term, and S = S 1 + S 2 + S 3 are other sources of heat. Changes in thermal energy over time (⇢C p @T /@t) can be caused by several different processes. The term ⇢C p v · rT corresponds to advective heat transport, i.e., the transport of thermal energy with the motion of material. Consequently, it depends on the velocity the material moves with v, and on how much the temperature, and with that the thermal energy, changes in space, which can be expressed through 245 the temperature gradient rT . Heat is also transported by conduction, which is represented by the term r · (krT ). When the thermal conductivity k is larger, the curvature of any temperature variations as expressed by the second derivative r · rT becomes steeper and more heat is diffused. If the model includes no additional heating processes, Eq. (6) can be simplified by dividing it by ⇢C p : 250 which introduces the thermal diffusivity  = k/(⇢C p ).
The terms on the right-hand side of Eq. (6) describe other heating processes that can be significant in a geodynamic model.
Internal heating can be expressed as the product of the density ⇢ and a volumetric heat production rate H. Its most common source is the decay of radiogenic elements. Accordingly, this process is important in regions where the concentration of heatproducing elements is high, such as the continental crust. Frictional heating, also called shear heating and viscous dissipation, 255 describes the amount of energy that is released as heat when material is deformed. The more stress is required to deform the material, and the larger the amount of deformation, the more shear heating is produced. This can be expressed as S 1 = ⌧ :" (where ⌧ is the deviatoric stress, i.e. the part of the stress not related to the hydrostatic pressure, and : is the matrix inner product, so ⌧ :" = P i P j ⌧ ij"ij ). Adiabatic heating describes how the temperature of the material changes when it is compressed or when it expands due to changes in pressure, assuming that no energy is exchanged with the surrounding material during 260 this process. Hence, the work being done to compress the material is released as heat. This temperature change depends on the thermal expansivity ↵, which expresses how much the material expands for a given increase in temperature, and on the changes in pressure over time: S 2 = ↵T Dp/Dt. Since the dominant pressure variation in the Earth's interior is the effect of the lithostatic pressure increase with depth, the usual assumption is that pressure only changes due to the motion of material with velocity v along a pressure gradient rp, resulting in S 2 = ↵T (v · rp).

265
When material undergoes phase changes, thermal energy is consumed (endothermic reaction) or released (exothermic reaction) as latent heat. This happens both for solid-state phase transitions, such as from olivine to wadsleyite, and for melting of mantle rocks. For thin phase transitions, this can lead to abrupt changes in mantle temperature where phase transitions are located, such as around 410 km depth. The amount of energy released or consumed is proportional to the density ⇢, the temperature T , and the change in the thermodynamic quantity entropy across the phase transitions S, i.e., S 3 = ⇢T SDX/Dt 270 (Schubert et al., 2001, Chapter 4.1). Here, X is the fraction of the material that has already undergone the transition, sometimes also called the phase function. For mantle melting, X would be the melt fraction. The terms discussed above are generally the most important heating processes in the Earth's interior. While there may be other sources of heat or heat transport, they are assumed to be so small that they can be neglected in most cases, such as heat transport by radiation, or the energy required for reducing the mineral grain size and creating new boundaries between the grains.

Approximations of the governing equations
In the previous sections on the mass, momentum, and energy equations, we have already seen that there are different ways these equations can be simplified and that there is a choice in which physical processes to include. Based on an analysis of the equations, there are a number of different approximations that are commonly used in geodynamics (see also Gassmöller et al., 2020;Schubert et al., 2001) 280 The Anelastic Liquid Approximation (ALA, Jarvis and McKenzie, 1980) assumes that lateral density variations are small relative to a reference density profile and can be neglected in the mass (3) and energy conservation equations (6), which instead use the density on the depth-dependent reference profile. Only the buoyancy term in the momentum equation (5) uses a temperature-and pressure-dependent density, which is approximated using a Taylor expansion (Appendix A1). This approximation is commonly used in whole-mantle convection models, where the density changes substantially from the surface 285 to the core-mantle boundary, and compressibility is an important effect. Its use is appropriate as long as deviations from the reference profile remain small or are not important for the interpretation of the model. Processes that may invalidate this assumption are, e.g., substantial expansion or contraction caused by temperature changes that are not on the reference profile, such as cooling or heating in the thermal boundary layers caused by conduction, which can not be modelled with this approximation Phase transitions of different types of materials that cause large density or temperature changes also invalidate 290 this assumption, since only one type of material can be represented by the reference profile (see also Gassmöller et al., 2020).
The Truncated Anelastic Liquid Approximation (TALA, Jarvis and McKenzie, 1980;Ita and King, 1994) is based on the same assumptions as the ALA, but further assumes that the variation of the density due to pressure variations can be neglected in the buoyancy term as well. It has similar applications as ALA, but introduces an imbalance between energy dissipation calculated from the Stokes equation and heat dissipation in the energy equation. Consequently, it should be avoided when 295 energy dissipation is important in the model (Leng and Zhong, 2008;Alboussière and Ricard, 2013).
The Boussinesq approximation (BA, Oberbeck, 1879;Boussinesq, 1903;Rayleigh, 1916) assumes that density variations are so small that they can be neglected everywhere except in the buoyancy term in the momentum equation (5), which is equivalent to using a constant reference density profile. This implies incompressibility, i.e., the use of Eq. (4) for mass conservation.
In addition, adiabatic and shear heating are not considered in the energy equation (6). This approximation is valid as long as 300 density variations are small (see also Section 2.1.1) and the modelled processes would cause no substantial shear or adiabatic heating. The Boussinesq approximation is often used in lithosphere-scale models. Due to its simplicity, the approximation of incompressibility is sometimes also adopted for whole-mantle convection models, where it is only approximately valid since the density changes by 65% across the mantle  and it has been shown that compressibility can have a large effect of the pattern of convective flow (e.g. Tackley, 1996).

305
The extended Boussinesq approximation (EBA, Christensen and Yuen, 1985;Oxburgh and Turcotte, 1978) is based on the same assumptions as the BA, but does consider adiabatic and shear heating. Since it includes adiabatic heating, but not the associated volume and density changes, adiabatic heating leads to artificial generation of energy in the model. Consequently, the extended Boussinesq approximation should only be used in models without substantial adiabatic temperature changes.
Before deciding to use any of these approximations in a geodynamic model, it is important to make sure that the underlying 310 assumption are justified for that particular application. For a comparison between some of these approximations using benchmark models, see e.g., Steinbach et al. (1989); Leng and Zhong (2008); King et al. (2010); Gassmöller et al. (2020). Also note that technically, these approximations are all internally inconsistent to varying degrees (Section 6.2), since they do not fulfil the definitions of thermodynamic variables, but use linearised versions instead, and use different density formulations in the different equations. Nevertheless, many of them are generally accepted and widely used in geodynamic modelling studies, as 315 they allow for simplifying the equations to make it easier to obtain solutions.

The constitutive equations: Material properties
Using the equations discussed in the previous section 2.1 to model how Earth materials move under the influence of gravity requires some assumptions about how Earth materials respond to external forces or conditions. These relations are often called constitutive equations, and they relate the material properties to the solution variables of the conservation equations, like 320 temperature and pressure. Which of these relations is most important for the model evolution depends on the application. In regional models that impose the plate driving forces externally, the relation between stress and deformation, i.e., the rheology (Section 2.2.1) can be the most important parameter choice in the model, and can control most of the model evolution. Since buoyancy is one of the main forces that drive the flow of Earth materials, it is often most important to take into account how the density depends on other variables, and to include this dependency in the buoyancy term, i.e., the right-hand side of Eq. (5).

325
This is particularly important for mantle convection models, and is described by the equation of state (Section 2.2.5) of the material. In models of the lithosphere and crust, stress and state of deformation become more important.

Rheology of Earth materials
The rheology describes how materials deform when they are exposed to stresses. This relation between stress and deformation enters the momentum equation (5) in the term r· , which represents the surface forces, and is also used to compute the amount 330 of shear heating ⌧ :" in the energy balance. In many cases, the material reacts differently to shear stresses (acting parallel to the surface) and normal stresses (acting normal to the surface and leading to volumetric deformation). Correspondingly, the resistance to these different types of deformation is expressed by different material parameters: the shear viscosity or shear modulus, and the bulk viscosity or bulk modulus. 335 We will start by discussing a very simple type of rheology with the following features: (i) only viscous (rather than elastic or brittle) deformation, (ii) deformation behaviour that does not depend on the direction of deformation (corresponding to an isotropic material), (iii) a linear relation between the stress and the rate of deformation (corresponding to a Newtonian fluid), and (iv) a low amount of energy dissipation under compaction or expansion (corresponding to a negligible bulk viscosity). In this case, the (symmetric) stress tensor writes:

Viscous deformation
where p is the pressure, I is the unit tensor, ⌧ is the deviatoric stress, ⌘ is the (dynamic) shear viscosity and"(v) is the rate of deformation tensor (often called strain rate tensor), which is directly related to the velocity gradients in the fluid, and is defined aṡ

345
This strain rate tensor describes both volume changes of the material and changes in shape, and can be written as the sum of these two deformation components. The prime in Eq. (8) denotes the deviatoric part of the tensor, i.e., the part that refers to the change in shape. It follows that the case of incompressible flows, where the density is assumed to be constant and r · v = 0 (Eq. 4), also implies"(v) 0 ="(v). Hence, incompressibility implies that there is no volume change of the material, so this part of the strain rate tensor is not taken into account.

350
In the Earth's interior, viscous deformation is the dominant rock deformation process on long time scales, if temperatures are not too far from the rocks' melting temperature. Under these conditions, imperfections in the crystal lattice move through mineral grains and contribute to large-scale deformation over time. The physical process that is thought to most closely resemble this idealised case of a Newtonian fluid is diffusion creep, where single atoms or vacancies, i.e., the places in the crystal lattice where an atom is missing, migrate through the crystal. Diffusion creep is assumed to be the dominant deformation 355 mechanism in the lower mantle. Other creep processes exhibit different behaviour. Dislocation creep, which is considered to be important in the upper mantle, is enabled by defects in the crystal lattice that are not just a single atom, but an entire line of atoms. Dislocation creep is highly sensitive to the stress applied to a rock, which means that the relation between stress and strain rate is no longer linear, but a power law. Peierls creep, which has an exponential relation between stress and rate of deformation, becomes important at even higher stresses (Gerya, 2019). Grain boundary sliding describes deformation due to 360 grains sliding against each other along grain boundaries, which has to be accommodated by another mechanism that allows the grains themselves to deform (e.g., diffusion or dislocation creep). It becomes important at high temperatures and high stresses (Hansen et al., 2011;Ohuchi et al., 2015).
Consequently, the viscosity of rocks strongly depends on e.g., temperature, pressure, stress, size of the mineral grains, deformation history, and the presence of melt or water and it varies by several orders of magnitude, often over small length 365 scales. These flow laws can be expressed in a generalised form using the relation (Hirth and Kohlstedt, 2003): where d is the grain size, fH 2 O is the water fugacity and relates to the presence of water, is the melt fraction, and A, n, q, r, ↵ ? , E, V and R are constants. Each of these parameters varies depending on the mineral and creep mechanism, and which of these dependencies is most important depends on the type of model. For example, diffusion creep exhibits a linear 370 relation between stress and strain rate, so n = 1, whereas dislocation creep strongly depends on stress (typically n ⇠ 3 4), but there is no dependence on grain size (i.e., q = 0). Accordingly, temperature is probably the most important control in the lower mantle, while in the upper mantle, the stress-dependence plays a dominant role. Grain size d is another potentially strong influence (e.g. Hirth and Kohlstedt, 2003;King, 2016;Jain et al., 2018), but this relation is not well constrained, which makes it hard to include this effect in geodynamic models . Generally, several deformation mechanisms with 375 different flow laws (Eq. 10) can be active at the same time, and their individual strain rates add up to the total 'effective' rate of deformation.

Elastic deformation
On shorter time scales, the elastic behaviour of rocks becomes important in addition to viscous flow. In the case of elasticity, the stress tensor is related to the strain tensor through the generalised Hooke's law = C : "(u) where C is the fourth-order 380 elastic tensor and " is the strain tensor which depends on the displacement vector u as follows: Hence, for elastic deformation, the stress is proportional to the amount of deformation rather than the rate of deformation, as is the case for viscous processes. Due to the inherent symmetries of , ", and C, the tensor C is reduced to two numbers for isotropic media, and the stress-strain relation becomes: where is the first Lamé parameter and µ is the second Lamé parameter (also referred to as shear modulus), which describes how much the material deforms elastically under a given shear stress. These two moduli together define the Poisson's ratio of a material (also see Section 2.1.1).

390
For strong deformation and large stresses on short time scales, such as occurring in the crust and lithosphere, brittle deformation becomes the controlling mechanism. In this case, the crystal lattice breaks, and a fault or a network of faults (which can range from the atomic to the kilometre-scale) accommodates deformation. The relative motion of the two discrete sides of the fault is limited by the friction on the interface. However, in our continuum formulation, we cannot represent discrete faults and hence the deformation in geodynamics models localises in so-called shear bands of finite width, which can be interpreted as 395 faults on crustal and lithospheric scales (van Zelst et al., 2019). We approximate the brittle behaviour by so-called plasticity, characterised by a yield criterion that determines the yield stress that the maximum stress must satisfy. When locally built-up stresses reach this yield stress (also called yield strength), the rock is permanently deformed and the plastic strain rate is no longer zero. Note that in other Earth science communities the terms 'plastic deformation' and 'plasticity' are used to describe all non-recoverable, time-dependent deformation (e.g. solid-state creep, like dislocation creep) that is commonly also referred 400 to as ductile or viscous deformation (e.g., Karato, 2008;Fossen, 2016). Here we use 'plasticity' to refer to non-recoverable, nearly instantaneous, pressure-dependent yielding at high stresses and fracturing.
The expression for the yield stress yield of a rock is linked to the choice of plasticity yield criterion. One of the most well-known criteria is the Mohr-Coulomb criterion based on laboratory experiments (Handin, 1969;Kachanov, 2004): Here, C is the cohesion and f is the angle of friction. The strength of rocks is primarily characterised by these two parameters which are obtained from laboratory measurements. It is worth noting that the effective strength of the rock in Eq. (13) increases with the total pressure p. Another common way of denoting this yield criterion is where µ f represents the friction coefficient of the rock. If µ f = 0.6 is assumed, Eq. 14 reduces to Byerlee's law (Byerlee, . The numerical implementation of the Mohr-Coulomb yield criterion is convoluted. As a consequence, the related Drucker-Prager yield criterion (Drucker and Prager, 1952) is often preferred. In the case of a two dimensional (2D) model with a plane strain assumption and a Poisson's ratio of 0.5 (i.e., an incompressible model), the Drucker-Prager and Mohr-Coulomb yield criteria are equivalent for consistently chosen parameters. However, in general the Drucker-Prager and Mohr-Coulomb criteria 415 are not equivalent and the resulting strength of the rocks can differ by as much as 300% between them (Wojciechowski, 2018).
A common simplification of the yield criteria in Eqs. (13) and (14) is to use the lithostatic pressure p lith = ⇢gz instead of the total pressure, thereby ignoring contributions to the pressure by, e.g., the dynamic pressure and pore fluid pressure. This assumption effectively makes the yield criterion purely depth-dependent and is sometimes referred to as the depth-dependent von Mises criterion (Spiegelman et al., 2016). The constant von Mises yield criterion is another simplification of the  Prager yield criterion where it is assumed that f = 0 (von Mises, 1913). This in itself is also an approximation of the Tresca yield criterion, which is based on laboratory experiments like the Mohr-Coulomb yield criterion (Turcotte and Schubert, 2012).
In nature, the strength of rocks can change over time by, e.g., forming a fault gauge which lowers the strength of the rock.
To account for this variation in strength over time on tectonic time scales, the cohesion and friction coefficient of the rock can be made dependent on the strain or strain-rate in what is called strain-or strain-rate-weakening (also called softening):

425
when the strain rate increases, the strength of the rock is lowered. Similarly, strain-or strain-rate-hardening can be applied (e.g., Massmeyer et al., 2013). When geodynamic models are used to study fault and earthquake processes, they focus on a shorter time scale of rock strength variations ( Figure 1) and a more sophisticated way of strength variations needs to be used.
Experiments have shown that the strength of the rock is better represented as depending on slip velocity (i.e., the velocity with which the two sides of a fault slide past each other). In hybrid seismo-thermo-mechanical models, it is often assumed that the 430 slip velocity is equivalent to the strain-rate apart from a multiplication with an assumed constant fault width (van Dinther et al., 2013a). Complex friction laws based on the parameterisation of lab experiments have been developed to describe the change in friction coefficient µ f (Eq. 14) of rocks during an earthquake varying from simple linear slip weakening approximations (Ida, 1973) to complex rate- (van Dinther et al., 2013a;Dal Zilio et al., 2018;van Zelst et al., 2019) and rate-and-state dependent friction laws (Lapusta et al., 2000;Sobolev and Muldashev, 2017;Herrendörfer et al., 2018).

435
In the Earth, the relation between stress and deformation can be very complex, and deformation is a combination of elastic and viscous behaviour and brittle failure. In general, both the viscosity and the elastic moduli can depend on temperature, pressure, chemical composition, the presence of melt and fluids, the size and orientation of mineral grains, the rate of deformation or the deformation history of the material. Consequently, Earth materials are usually not isotropic, and the strength of the 440 material depends on the direction of deformation. To incorporate this behaviour into a geodynamic models, the viscosity and elastic moduli have to be expressed as tensors, and cannot be reduced to one or two parameters (Mühlhaus et al., 2002;Lev and Hager, 2008). This complexity is not taken into account in most geodynamic models. This is on the one hand because it is computationally expensive and requires substantial effort to implement in a model, but also because there is substantial uncertainty in how, specifically, rocks deform anisotropically. For models that include anisotropy, see for example Mühlhaus

The equation of state
The relation between density and temperature and pressure and sometimes other variables, like chemical composition is often called the equation of state. It describes the state of the material (density) under a set of conditions (temperature, pressure, 450 etc.). It incorporates material properties such as the thermal expansivity ↵, which describes how much the rock expands when the temperature is increased, and the compressibility, which describes how much the volume of a rock decreases when it is exposed to higher pressures. These relations should be chosen in a way that is thermodynamically consistent (see Section 6.2).
Depending on the model application, there are many different equations of state that can be used. For models that aim to capture the first-order effects of a given process, analyse the influence of a material property on the model evolution, or 455 develop a scaling law (see section 5.2.2), it is often appropriate to simplify the relationships between solution variables and material properties. In mantle convection models, the temperature usually has the largest influence on density variations that affect buoyancy. So in the simplest case, the equation of state may be a linear relationship between density and temperature.
However, the most important variable depends on the application. For example, chemical composition plays an important role in both the outer core, where temperature variations are very small, and in lithospheric deformation.

460
On the other end of the spectrum, there are models where it is less important to understand the influence of different material properties on the model evolution as the goal is to fit existing observations (e.g., seismic wave speeds, see also Section 5.2.1). In such a scenario, it is often desirable to include the material properties in the most realistic way that our current knowledge allows. This requires knowledge about mineral physics and thermodynamics, and is therefore often handled by external programs. Software packages such as Perple_X (Connolly, 1990), HeFESTo (Stixrude and Lithgow-Bertelloni, 2005) 465 or BurnMan (Cottaar et al., 2014), in connection with a mineral physics database, include the complex thermodynamic relations that determine how mineral assemblages evolve under Earth conditions. They compute material properties in dependence of the solution variables of the geodynamic model (usually temperature and pressure, but alternatively, specific entropy and volume, or a combination of these variables, can be used; see Connolly, 2009).

470
The mass, momentum, and energy conservation equations are visibly coupled since velocity (derivatives) and pressure enter the stress tensor, and thermal energy transport due to advection depends on the velocity. More importantly, the previous sections have highlighted how material properties such as ↵, ⇢, C p depend on pressure, temperature and other quantities such as composition and that deformation is partitioned between various mechanisms. Their pressure-, temperature-and strain ratedependence makes the exact partitioning complex and renders the conservation equations (2), (5) and (6) very nonlinear with 475 regards to the primary variables v, p and T . Hence, the equations usually contain nonlinear terms that imply a relationship between the solution variables that is no longer linear.
In addition to temperature, pressure, and velocity there may be other conditions in the model that change over time and are important for the model evolution, but not directly related to changes in temperature, pressure and velocity. A common example is the chemical composition of the material (which can refer to the major element composition, but may also relate to 480 e.g., the water content). In this case, a transport equation is required for every additional quantity that should be tracked in the model and moves with the material flow: where c i is the to-be-tracked quantity, t is time, (v) is velocity, and n c is the number of compositions, fields, or materials present. This advection equation assumes that at any given location, changes in composition over time @c i /@t are caused by  Other physical processes may require additional terms or additional equations. Examples are the generation of the magnetic field in the outer core (Jones, 2011), two-phase (McKenzie, 1984Bercovici et al., 2001) or multi-phase flow (Oliveira et al., 2018;Keller and Suckale, 2019), disequilibrium melting (Rudge et al., 2011), complex magma dynamics in the crust (Keller et al., 2013), reactive melt transport (Aharonov et al., 1995;Keller and Katz, 2016), (de)hydration reactions and water transport (Faccenda et al., 2012;Magni et al., 2014;Quinquis and Buiter, 2014;Wilson et al., 2014a;Nakagawa et al., 2015), the 495 evolution of the mineral grain size (Hall and Parmentier, 2003;Solomatov and Reese, 2008;Bercovici and Ricard, 2012;Cerpa et al., 2017;Mulyukova and Bercovici, 2019), the fluid dynamics and thermodynamics of a magma ocean (Labrosse et al., 2007;Solomatov, 2000;Ulvrová et al., 2012), the interaction of tectonic processes with erosion and other surface processes (Burov and Cloetingh, 1997;Roe et al., 2006;Thieulot et al., 2014;Ueda et al., 2015;Sternai, 2020), anisotropic fabric (Mühlhaus et al., 2002;Lev and Hager, 2008;Heister et al., 2017;Faccenda, 2014;Perry-Houts and Karlstrom, 2018;Király et al., 2020a), phase transformation kinetics (Bina et al., 2001;Tetzlaff and Schmeling, 2009;Quinteros and Sobolev, 2012;Agrusta et al., 2014) and inertial processes and seismic cycles (van Dinther et al., 2013b, a;Sobolev and Muldashev, 2017;Tong et al., 2019). Since the aim of this contribution is to introduce the general concepts of geodynamic modelling, we will not discuss these more complex effects here, and will instead focus on the main conservation equations presented above.
In the end, the partial differential equations (2), (5), (6) and potentially equation (15) must be supplemented by the geometry 505 of the domain on which they are to be solved (e.g., 2D/3D, Cartesian, cylindrical, spherical; Section 5), a set of initial conditions for the time-dependent temperature and compositions and a set of boundary conditions (see Section 5.3).

The numerical model
The conservation equations of Section 2.1 generally cannot be solved analytically. However, computers are capable of finding approximate solutions to equations (2), (5) and (6). To this end, we must define a numerical model, i.e. the mathematical 510 description of the physical model in a computer language, to solve numerically. As computers cannot by construction represent a continuum, the above equations must be discretised and then solved on a finite domain divided into a finite number of grid points, which together form the mesh of the computation. This process requires specific numerical methods (Section 3.1) and choices (Section 3.2-3.7) and verification of the resulting numerical code (Section 4). For a recent in-depth analysis and outlook of numerical methods in geodynamics, we refer the reader to Morra et al. (2020).

Numerical methods
The solutions of the continuum equations described in Section 2.1 need to be computed on a finite number of moments in time and locations in space. Rewriting the equations of Section 2.1 for a discrete number of points is called discretisation.
The three main methods of discretisation in geodynamics are the finite element method (FEM), the finite difference method (FDM) and the finite volume method (FVM). Often combinations of these methods are used to deal with time and space 520 discretisation separately. All three methods convert Ordinary Differential Equations (ODEs) or Partial Differential Equations (PDEs), which may be nonlinear, into a system of linear equations that can be solved by matrix algebra techniques. None is intrinsically better than the other, although there are differences that make a certain method more suitable for certain types of science questions. For instance the finite element and finite volume methods are capable of dealing with non-Cartesian geometries, such as the spherical shape of the Earth, or topography at the surface of a model, while the finite difference method 525 is not. On the other hand, the finite difference method is more intuitive in its design than the other two. An overview of the basic principles of these numerical methods is available in Ismail-Zadeh and  and in Zhong et al. (2015) and an in-depth exposé of the finite difference method in geodynamics is available in Gerya (2019). It is worth noting that spectral methods are also encountered in mantle dynamics modelling (Ismail-Zadeh and Petrunin et al., 2013), as is the boundary element method (BEM) (Morra et al., 2009(Morra et al., , 2010 and the Lattice Boltzmann Method (Mora and Yuen, 2017). Golabek et al., 2018) and radial basis function method (RBF, Arrial et al., 2014) have also been used in geodynamic modelling. In what follows we focus on the most popular methods, i.e. finite differences, finite elements, and finite volumes.

535
The discretisation concept for all three main methods (FEM, FDM, FVM) is identical. The domain is subdivided in cells or elements, as shown in Figure 4. In essence, the methods look for the solution of the equations of Section 2.1 in the form of combinations of polynomial expressions (also called shape functions in FEM terminology) defined on each element or cell.
To illustrate this concept, we provide a small example here for the conservation of energy using the finite difference method, which is based on a Taylor expansion keeping only first and second order terms (see Appendix A for the complete example). In 540 one dimension, and under the assumptions that there is no advection or heat sources, and that the coefficients are all constant in space, Eq. (6) becomes what is commonly called the heat equation: This is to be solved for T (x, t), for x between 0 and L x , and t between 0 and t final . The time discretisation describes how to use a temperature distribution that is known at time t n to compute a new temperature at time t n+1 = t n + t, with t being 545 the time step size. The discretisation in space means that this temperature distribution is computed at a finite number np of points. The simplest way of choosing the position of these points is such that x i = ih, with h = L x /(np 1) being the distance between points (which is often referred to as the resolution or grid size of the model) and the point indices running from 0  i  np 1. As shown in Appendix A, there are different ways of discretising Eq. (16), but the so-called explicit version where the subscripts refer to space indices and the superscripts refer to time indices, i.e., T n i = T (x i , t n ). This also means that T 0 i is the initial temperature condition and it is needed to start the calculation. In this example, we know all temperatures at time step n and we can easily compute the temperatures at time step n + 1 at all locations i inside the domain, while the temperature on the boundary T (x i=0 ) and T (x i=np 1 ) is prescribed by the boundary conditions. It also illustrates what the 555 name of the Finite Difference Method means: when going from the continuum to a finite number of grid points, derivatives are approximated by differences in temperature between these points.
While finite difference codes only need to specify the order of the approximation they rely on (and the associated stencil), finite element codes often specify the type of element they use (i.e., the specific polynomial expressions for the shape functions) using specific terminology. Earlier codes such as Citcom(S/CU) (Moresi and Solomatov, 1995), ConMan (King et al., 1990) 560 and SOPALE (Fullsack, 1995) relied on the computationally cheap (i.e. yielding the smallest possible linear system size) Q 1 ⇥ P 0 element. In this case Q 1 indicates continuous bi/trilinear velocity and P 0 (or Q 0 ) discontinuous constant pressure (Donea and Huerta, 2003). Modern codes such as ASPECT (Kronbichler et al., 2012) and pTatin3D (May et al., 2015) rely on the more accurate Q 2 ⇥ Q 1 and Q 2 ⇥ P 1 (sometimes denoted Q 2 ⇥ P disc 1 ) elements, respectively. In this case, Q 2 stands for continuous bi/triquadratic velocity, Q 1 stands for continuous bi/trilinear pressure and P 1 stands for discontinuous linear 565 pressure.
For all methods, the discretisation process results in a linear system of equations with its size being the number of unknowns, i.e., a multiple of the number of nodes/elements. This system of equations is written as where A is a large and very sparse matrix, X is the vector of unknowns, typically consisting of velocity components, pressure, 570 and temperature on the grid, and b is the known right-hand side vector. In each time step, the system is solved, new solution fields are obtained, post-processed, and analysed. If the model evolves in time, the domain and therefore the mesh may evolve, compositions are transported and a new system is formed.

Kinematical description
After the discretisation step, a kinematical description must be chosen to define how the material is going to move in the model.
Eulerian codes have a fixed mesh through which material flows. Since the evolution of the top boundary of the model is often of prime importance in geodynamical studies as it accounts for the generated topography, a feature that is directly and easily observable on Earth, the air above the crust must be modelled as well to allow for the formation of topography. This 580 air layer has been coined 'sticky air' (Crameri et al., 2012) because of its much higher viscosity (⇠ 10 18 Pa s) than real air (⇠ 10 5 Pa.s). This higher viscosity is necessary to avoid numerically tricky viscosity contrasts in the domain of more than 20 orders of magnitude. Although the sticky air viscosity seems high, it remains very low compared to the effective viscosity of the crust (⇠ 10 2 4 Pa s) bordering the sticky air, so that no shear forces are exerted on the crust surface by the air layer. Hence, topography can be developed. Finite difference codes typically use the Eulerian kinematical description and use the sticky air 585 approach. Problems associated with the use of sticky air and a possible solution are presented in Duretz et al. (2016).
In contrast to a Eulerian kinematical description, the mesh of Lagrangian codes deforms with the computed flow and therefore does not require the use of sticky air to model topography. This limits Lagrangian codes to small deformation. For example, subduction processes would quickly deform the mesh to such a point that it would not be suitable for accurate calculations.
PyLith (Aagaard et al., 2013), a finite element code for dynamic and quasistatic simulations of crustal deformation, uses a 590 Lagrangian kinematical description.
Finally, as its name implies, the Arbitrary Lagrangian-Eulerian (ALE) kinematical description combines features of both the Lagrangian and Eulerian formulations. In geodynamical codes, it often amounts to the mesh conforming vertically or radially to the free surface, while retaining its structure in the horizontal or tangential direction ( Figure 5). This approach forms the basis of the SOPALE (Fullsack, 1995), FANTOM (Thieulot, 2011), pTatin3D (May et al., 2015), and ASPECT (Rose et al.,

Solving the linear system
The discretisation process outlined in section 3.2 leads to a linear system as given in Eq. (18). There are many ways to solve these large linear systems of equations, which can be split into two families: direct methods and iterative methods. Direct methods exploit the sparsity of the matrix (i.e., the ratio of non-zero terms to the number of terms in the matrix is very small) 600 and arrive at the desired solution within a predictable number of operations. Because of the memory requirements of the direct solvers to store and process the matrix terms, these tend to be used for 2D applications (or low resolution 3D problems) only.
Iterative methods start with an initial solution (a guess) and incrementally improve on it until a convergence criterion is 605 met, i.e. the remaining error (residual) is small. Common iterative methods in geodynamic codes are the Conjugate Gradient method and the GMRES (Generalized Minimal RESidual) method (Saad, 2003;Langer and Neumüller, 2018). Note that the choice of iterative solving method is intrinsically tied to the properties of the matrix A. The most popular iterative solvers packages are PETSc (Balay et al., 1997) (used in pTatin3D, May et al., 2015) and Trilinos (Heroux et al., 2005) (used in ASPECT, Kronbichler et al., 2012).

Computer architectures & parallelisation
Early computing architectures of the 1970s were quite limited by today's standards and predominantly relied on sequential programming where one task is performed after the other ( Figure 6). Hence, early models in the 1970s were confined to a few thousands of unknowns at most (Stephansson and Berner, 1971;Beaumont and Lambert, 1972). The computer architectures and processor speeds on which calculations are performed have since vastly increased in the past decades following Moore's 615 law (Eijkhout, 2013) resulting in models with several hundreds of millions of unknowns now possible (Alisic et al., 2012).
Accommodating such enormous computational loads was only made possible by the use of parallel architectures ( Figure 6) and their growing availability. Such architectures can be multi-core processor-based commodity-grade desktop computers, home-made (so-called Beowulf clusters) or super calculators counting up to a few millions of processors as listed in the top 500. The majority of the state-of-the-art geodynamics codes run in parallel on distributed memory systems and rely on the 620 MPI (Message Passing Interface) standard ( Figure 6). In MPI, domain decomposition is used, i.e. each computer thread of tasks only 'knows' about a subset of the mesh for which the solution is obtained. This allows for a virtually unlimited number of threads but it requires that the whole code is built with parallelism in mind to achieve maximum efficiency (Saad, 2003).
Other codes, such as the surface processes code FastScape (Braun and Willett, 2013) or the geodynamics code I2/I3ELVIS (Gerya and Yuen, 2007), rely on a different approach based on OpenMP which is limited to running on shared memory 625 systems. OpenMP can be added later on to an already existing sequential code and targets areas of the code which take the most time. Although appealing at first, this level of parallelism is limited by the number of CPUs attached to the single memory (typically a few dozen at the maximum). It is worth noticing that many codes or the linear algebra libraries that they link to use a combination of both MPI and OpenMP. When documenting the parallel performance of a code, one often talks of strong and weak scaling (Kronbichler et al., 2012;630 May et al., Kaus et al., 2016). Strong scaling is defined as how the solution time varies with the number of processors for a fixed problem size. Optimal strong scaling is achieved if the solution time is inversely proportional to the number of processors. Conversely, when looking at weak scaling, both the number of processors and the problem size are increased by the same factor. This also results in a constant workload per processor and one therefore expects a constant solution time for optimal weak scaling. In practice, all geodynamic modelling software have different limits for their strong and weak scaling, 635 limiting the size of the model they can effectively compute.

Dealing with nonlinearities
As mentioned in Section 2.3, the set of partial differential equations to be solved is likely to be nonlinear, i.e. the coefficients of the equations depend on the solution. For example, the viscosity depends on the velocity via the strain rate tensor, which makes the term of the momentum equation that is the product of viscosity and strain rate a nonlinear term. Special techniques 640 must then be employed to solve this system and so-called nonlinear iterations are carried out on the linearised equations until a When an Arbitrary Lagrangian-Eulerian approach is used, the domain width is often kept constant in geodynamic applications, such that the mesh only deforms vertically to accommodate the topography. In the Lagrangian formulation, the mesh deforms with the velocity computed on its nodes. converged solution is obtained. Note that these nonlinear iterations are distinct from the iterations taking place in the iterative method employed to solve the system. One type of nonlinear iterations are the fixed-point ('Picard') iterations, where a guess of the solution is made and it is used to compute the solution-dependent coefficients. The linear system is solved, and the coefficients are updated with the new solution. This process is repeated until the left-hand side and the right-hand side of the 645 linear system match up to a given tolerance. Then, the iterative process has reached convergence. This technique is sub-optimal in its convergence rate as it often requires dozens or hundreds of iterations for strongly nonlinear problems, i.e., the linear system must be solved as many times (Spiegelman et al., 2016;Glerum et al., 2018).
State-of-the-art codes now all rely on some form of Newton iterations based on Newton-Raphson's method, which guarantee quadratic convergence towards the solution in contrast to the linear convergence of the Picard iteration, provided the 650 initial guess is close enough to the solution of the nonlinear problem (May et al., 2015;Kaus et al., 2016;Rudi et al., 2015;Spiegelman et al., 2016). It is worth noting that the implementation of Newton's method is substantially more complex and therefore difficult than a Picard iteration (Fraters et al., 2019a).

Tracking materials
Aside from the conservation of mass, momentum, and energy, special care must also be taken when solving the advection 655 equation (15)  throughout the domain at the beginning of the simulation with an assigned material identity, e.g., crust, mantle lithosphere, or sediments. These properties are then projected onto the mesh on which they are needed for solving the partial differential 660 equations. Once a new velocity field has been computed, these markers are advected with the flow velocity field, and the process is repeated at each timestep. This method is often found in finite element (Tackley and King, 2003;Deubelbeiss and Kaus, 2008;Thielmann et al., 2014;Gassmöller et al., 2018Gassmöller et al., , 2019 and finite difference codes alike (Gerya and Yuen, 2007;Pusok et al., 2016). Other employed techniques are the level-set method (Samuel and Evonuk, 2010;Hillebrand et al., 2014;Braun et al., 2008), the grid-based advection method (sometimes referred to as compositional fields method, Figure  Each method comes with its own strengths and weaknesses, and needs to be assessed in the context of the problem at hand. For example, particle-based methods are relatively easy to implement and not subject to numerical diffusion (a form of undesirable smearing of sharp gradients) but take up a lot of memory and are difficult to parallelise. Field-based methods are 670 notoriously difficult to advect and require numerical stabilisation, but are straightforward to parallelise and much lighter in terms of memory usage.

Multiphysics
On Earth, the lithosphere interacts with the cryosphere, biosphere, atmosphere, hydrosphere, magnetosphere, and other systems, and the deformation of the lithosphere is related to many natural hazards such as earthquakes, volcanic eruptions, land-675 slides and tsunamis. These systems are often inherently multi-scale with different processes occurring on vastly different time and/or length scales ( Figure 1) and/or multi-physics where the coupling of different physical processes is important. Such multiphysics or multidisciplinary research often sees researchers coupling a geodynamics code with another existing one, be it for surface processes (Braun and Yamato, 2010;Thieulot et al., 2014;Ueda et al., 2015;Beucher and Huismans, 2020), atmosphere evolution (Gillmann and Tackley, 2014), planetary impact (Golabek et al., 2018), plate reconstruction modelling 680 (Bull et al., 2014;Gassmöller et al., 2016), elastic flexure (Naliboff et al., 2009), dynamic rupture models and tsunami models (van Zelst et al., 2019;Madden et al., 2020). This coupling often takes place in the boundary conditions of the geodynamics code at the surface, or at the bottom of the lithosphere or mantle. This coupling can be one-sided when, e.g., plate velocities are prescribed. The coupling can also be two-sided when, e.g., the uplift rate and a current topography of a geodynamic model are used in a surface process model which returns an eroded surface which forms the new top boundary of a geodynamic model.

685
Coupling codes and the clever use of boundary conditions (Section 5.3.1) are promising avenues to incorporate geodynamic models in a host of multidisciplinary research.

Code verification
At this point, we have described the model in terms of governing and constitutive equations, and discretised and solved the system using appropriate numerical methods. Hence, an application code has been obtained at this point. However, before any 690 scientific study can be performed with confidence, the code must be tested to do what it is intended to do. This process has two components: verification and validation (Section 6). Verification means 'solving the equations right', while validation means 'solving the right equations' (Roache, 1997). In this section, we will briefly focus on code verification, which involves using tests and benchmarks.
In the software engineering community, the importance of tests is well acknowledged: "Without tests every change is a 695 possible bug" (Martin, 2008). While there are only a few equations to solve in geodynamics (Section 2; Lynch (2005)), the coefficients and the constitutive equations (i.e., rock rheology, multi-component processes) to describe the physical system often turn them into complex nonlinear problems, which can be difficult to test. Tests should verify various parts or whole functionality of the code. It is also desirable to test as many features of the code as possible. Of course, one can think of an infinite number of tests for even a simple geodynamics application, but because code verification can be time-consuming, a 700 balance has to be found between test coverage and production. Without doubt, verification is an absolutely necessary stage of code development as its purpose is to test the robustness of the code in a broad range of situations and applications.
The recommended approach is to implement only what is needed, test it, and then refactor and expand the system to implement new features later on, and test again. Implementing an automatic testing framework can speed up these steps. Good tests also follow the FIRST acronym: Fast (run quickly), Independent (do not depend on each other), Repeatable (in any envi-705 ronment), and Self-validating (have a boolean output, pass or fail) Tests (Martin, 2008). Without proper testing, model results may look beautiful, but are probably wrong; and having only simple tests may not catch bugs in complex geodynamic applications (Gerya, 2019). Moreover, just inspecting the output and saying 'it looks right' is not sufficient. Instead, checking error norms and expected convergence rates or other numerical diagnostics is more robust. Tests should verify the code on a broad range of challenging cases relevant to the scientific question. For example, in lithosphere dynamics modelling, the code should 710 correctly handle time-dependent problems with a variable viscosity that depends on temperature and strain rate, as well as complex geometries with a free surface, and the transport of chemical heterogeneities.
The types of tests one can do to benchmark codes relate to numerical implementation and functionality. In general, this means comparing the numerical solution obtained from solving the system of equations to analytical solutions, results of analogue experiments, numerical results from other codes and general physical considerations (Gerya, 2019). However, one 715 can also make smaller unit tests that verify correct functionality of individual routines, i.e., not necessarily routines that solve for physics, and make sure they do not change when modifications are done somewhere else in the code. For example, parallelisation is not a physical process but poses new problems that can be tested with unit testing. In a parallel code, one can verify that the parallel communication routines work correctly with a simple hello world test. Recently, unit tests can be automatically incorporated in the development process through platforms like GitHub, Gitlab, BitBucket, and Travis CI. Analytical solutions for code verification can also be used in the form of method of manufactured solutions (MMS) (Roache, 1997(Roache, , 2002. In principle, the user manufactures an exact solution or coefficients, without being concerned about their 730 physical meaning. The method of manufactured solutions primarily verifies correct implementation and the convergence order of a particular numerical method, i.e., the accuracy of the solution with increasing spatial and temporal resolution. Even if the method does not focus on the physical meaning of the solution, some physical considerations still need to be taken into account when using manufactured solutions, e.g., coefficients corresponding to physical parameters need to be positive. The method of manufactured solutions can be a powerful technique as it can be applied to arbitrarily complex partial differential equations 735 (e.g., linear, nonlinear, time dependent) and to arbitrarily complex domains and boundary conditions. Recent computational advances make this method easy to use by creating manufactured solutions using symbolic algebra packages such as Maple, MATLAB, Mathematica or SymPy. The method of manufactured solutions is not as commonly used in geodynamics as in engineering or applied mathematics, but it provides potential for more robust testing frameworks (Thieulot, 2014;May et al., 2013;Donea and Huerta, 2003;Burstedde et al., 2013). When analytical solutions are not possible, numerical experiments of the same model setup (i.e., equations, boundary conditions, geometry, and parameters) can be tested with a number of different codes within the community. These are called community benchmark and allow for testing complex problems by comparing similar diagnostics (e.g., non-dimensionless parameters such as the Nusselt number; Section 7.2) across codes. The results are then compiled, and the best average behaviour among codes is taken as the benchmark for the test. Examples of community benchmarks include thermo-mechanical 745 convection (Blankenbach et al., 1989;Zhong et al., 2008;King et al., 2010;Arrial et al., 2014;Tosi et al., 2015), subduction (Schmeling et al., 2008;van Keken et al., 2008), Rayleigh-Taylor instabilities (van Keken, 1997;Louis-Napoléon et al., 2020), and motion of the free surface (Crameri et al., 2012). Important to note is that community benchmarks are not bulletproof against bugs and do not necessarily provide insights into the numerical or physical behaviour of the system.
Comparison with analogue experiments is important for calibrating numerical models for complex processes, required for 750 numerical modelling of large-scale tectonic processes (e.g., plastic failure). They can be used for both verification and validation of numerical models. For example, modelling sandbox experiments (Buck and Sokoutis, 1994;Buiter et al., 2006Buiter et al., , 2016 (Davaille, 1999;Davaille et al., 2011), indenter block experiment (Tapponnier et al., 1982;Peltzer and Tapponnier, 1988), and subduction dynamics (Schellart, 2005;Duarte et al., 2013;Király et al., 2020b). Since there are fundamental differences between numerical and laboratory experiments (e.g., the numerical model is discrete, whereas the laboratory model has 'infinite resolution'), the model results are often not identical, and instead certain characteristic features of the solutions need to be compared.

760
As the importance of testing is revealed, more software engineering practices are required to keep codes clean, testable, and robust (see Section 9). Designing a suite of complete tests is just as important as building efficient, fast, realistic, highresolution numerical 3D codes (Gerya, 2019). There are many excellent books on analytical solutions (Turcotte and Schubert, 2012), keeping codes clean (Martin, 2008;Hunt and Thomas, 1999) and building robust testing frameworks for geodynamics (Ismail-Zadeh and Tackley, 2010;Simpson, 2017;Morra, 2018). As a final note on code verification, complex systems of 765 equations (i.e., multi-phase, multi-component, nonlinear rheology) can still be badly posed and are difficult to test. For this reason, verification of simpler systems is important, and complex solutions should be validated using scaling analysis and against natural observations (see Section 6).

From the real world to a model setup
Designing a model is not straightforward. Before starting to design a model, it is important to understand the difference between 770 a code and a model. What the code is to a model, is an electric engine to an electric car. Since a working electric engine (/code) does not guarantee a working electric car (/model), the car has to be designed and tested as carefully as the engine itself.
The car should neither be too simple (the addition of seats can indeed prove necessary), nor should it be too complicated Figure 7. Different model complexities for the heart (top) and the Earth (bottom). A simple model can be more useful: the basic shape of the heart has likely become the most successful model, indeed a true icon, only because it was neither too complex (it can be reproduced easily), nor too simple (its characteristic shape is still recognisable). To find the right level of complexity is challenging and must repeatedly be considered carefully by the modeller for each new modelling task at hand.
(two steering wheels might be an overkill). For subtle design choices, simplifying is by no means simple and has to be done carefully. Furthermore, the car should be customised for the purpose of driving: the intention to drive a racing circuit as fast 775 as possible is different from the intention to drive for shopping as reliably as possible. The underlying philosophy therefore has to be clear from the onset of the model design as it builds the basis for, and gives purpose to, any setup. Lastly, driving should start from a suitable point in time (a warm summer day might be better suited than a cold winter night) and in a suitable surrounding (flat asphalt might be more appropriate than a bumpy road). Hence, both initial and boundary conditions also need special consideration.

780
How to design a simple but not oversimplified geodynamic model that is based on a certain modelling philosophy and applies suitable initial and boundary conditions, and how to tackle this challenge in detail, is therefore outlined in the following sections.

Simplifying without oversimplifying
A model is, by definition, a simplified representation of a more complex natural phenomenon. This is a simple and obvious truth 785 that is easily forgotten when geodynamic models are interpreted, presented, and reviewed. It is the modeller's responsibility to not only constantly remind themselves, but also others about this key underlying fact.
The complexity of the planet Earth as a whole is vast. It is therefore challenging to reconcile such true complexity with a desired simplicity. A model can easily become too complex and, just as easily, oversimplified ( Figure 7). A geodynamic modeller should always strive to meet the appropriate middle ground in between a too complex and too simple model. This is 790 one of the first major tasks of a modeller when setting up a new modelling approach to answer a certain hypothesis: everything should be made as simple as possible, but not simpler.
So, how simple should a model optimally be? The answer to this question is not an easy one, as it strongly depends on the purpose of the model, the capabilities to diagnose and understand it, and the hypothesis that it will test. It is clear though that a more complex model does not necessarily mean a better model. In fact, a simpler model is often better than a more complex 795 model. A simpler model is clearer to understand, clearer to communicate, and, by making fewer assumptions, more likely right (Occam's razor; Thorburn, 1915), although not necessarily directly applicable to the real world.
There are various ways to reduce model complexity (Figure 8). Simplifying the physical model is one of them. Any given physical complexity of the natural phenomenon in question has to be evaluated and a decision has to be made to either reproduce it, parameterise it (i.e., mimic the phenomenon with a simplified approach), or neglect it in the model. This decision is based 800 upon the model's purpose and the spatial and temporal extent of the process under investigation.
Further model simplification is achieved through numerical adjustments. For example, all the following studies model plate tectonics, but the geometry of the model can be complex (e.g., a 3D spherical domain like Arnould et al. (2020)) or simple (e.g., a 1D or 2D domain like Bercovici and Ricard (2014)). One can choose the temporal extent of the model to be more complex (e.g., time-dependent with a long evolution as in Crameri et al. (2012)) or simply instantaneous as done by Stadler et al. (2010).

805
For the same model, such simplifying choices make for a generally simpler model. However, a model with simpler geometry (or time evolution) can, of course, also feature more complicated processes than one with a more complex geometry. Indeed, a simpler model geometry (or time evolution) often enables the modelling of more complex physical processes.
The numerical model complexity can also be adjusted by changing the initial and boundary conditions (heterogeneous like Conrad and Gurnis (2003) or homogeneous like Foley and Becker (2009)) and the imposed forcing (space and time dependent 810 forcing like Quere and Forte (2006) or self-consistent like Rolf and Tackley (2011)). Overall, complexity should be decided based on the scientific question addressed and the focus and scope of the study.

Modelling philosophies
There are two overarching geodynamic modelling philosophies: specific modelling and generic modelling (Figure 9). The first philosophy attempts to reproduce the specific state of a certain geodynamic system (e.g., based on a specific observation) with 815 a specific model to better understand the system's specific state. In contrast, generic modelling attempts to produce different regimes of behaviour of a certain geodynamic system (e.g., based on a general observation) to better understand the system's general behaviour.
Both overarching modelling philosophies can result in instantaneous and time-dependent studies (Tackley et al., 2005).
Instantaneous models are focused on resolving a certain state and usually rely heavily on a comparison with measured geo-820 physical, geological, and geochemical data. Time-dependent models tend to focus more on the evolution or the natural state of a system. For general modelling, time-dependent studies can be performed either focusing on the transient evolution from a known initial state or the statistical (quasi-)steady state.

Specific modelling procedure
Modelling aimed to compare and understand a specific state of a geodynamic system necessitates the following procedure.

825
First, a specific observation in a certain region has to be defined. Second, a hypothesis about the control mechanism(s) has to be outlined. Third, a model setup needs to be designed considering three key aspects. The model needs to be able to produce the observed feature, include the hypothetical control mechanism(s), and physically link the control mechanism to the observed feature. Fourth, the model has to be simplified to be easily understandable without being oversimplified.
For specific modelling in particular, the modeller needs to keep in mind that there is no guarantee that the suspected control 830 mechanism is the actual, nor the only, controlling mechanism (see Section 1.2). A specific modelling philosophy is often used to understand the circumstances that facilitated natural hazard events, like earthquakes, in order to improve hazard analyses. For example, Smith and Sandwell (2006)  topographic history of a certain subduction zone. These investigations first build the best informed model of a subduction zone (e.g., the Cocos subduction zone in Southwestern Mexico, or the Southern Alaskan subduction zone in North America) and then use it to reproduce and study how this subduction zone impacts the topographic evolution of the adjacent topographic heights (e.g., the Trans-Mexican Volcanic Belt) as done in Gérault et al. (2015) or intra-continental shear zones (e.g., the Denali fault) as done in Haynie and Jadamec (2017). A global specific modelling example is the global lithosphere-asthenosphere modelling 840 study of Osei Tutu et al. (2018), which obtains plate boundary friction coefficients and the minimum asthenospheric viscosity cutoff value by optimising the fit with observed global plate motions.

Generic modelling procedure
Modelling aimed to understand a general behaviour of a geodynamic system necessitates the following procedure. First, a general, first-order observation has to be defined. Second, a hypothesis about the controlling parameters and their possible 845 range has to be outlined. Third, a model setup needs to be designed considering two key aspects. The model needs to include the control mechanism(s), and it needs to be built on a set of assumptions for simplification.
For generic modelling in particular, the assumptions that go into designing the geodynamic model are key and need to be specified and described clearly (Section 8.1).
When a generic modelling philosophy is applied, a general geodynamic feature is investigated via a parameter study, where 850 a certain parameter space is mapped out and can be represented by a so-called regime diagram (e.g., Lourenço et al., 2020).
For example, surface elevation changes due to plate subduction can be investigated without specifically relating to one particular subduction zone, but to subduction zones in general. One way of doing that is to model subduction and its surface response, and vary all key subduction parameters over their individual Earth-like ranges as done in e.g., Crameri et al. (2017). Other generic modelling examples are quantifying crustal thickness at mid-ocean ridges for various spreading velocities (e.g., Katz, 2008), 855 reproducing general magma dynamics (e.g., Spiegelman, 1993) or magma transport behaviour (e.g., Yamato et al., 2012), the onset of convection in a planetary mantle (Turcotte and Schubert, 2012, section 6.19), testing for what Rayleigh numbers and Clapeyron slopes a phase transition induces layered convection (e.g., Christensen and Yuen, 1985), quantifying the amount of entrainment of a dense layer into mantle plumes (e.g., Lin and van Keken, 2006a, b;Jones et al., 2016), investigating the plate tectonic regimes of a planet, which might range from a stagnant lid with only one plate to a mobile lid, similar to modern plate 860 tectonics (e.g., Petersen et al., 2017;Lourenço et al., 2016;Lourenço et al., 2020), investigating under which conditions the flow in the Earth's outer core would cause an Earth-like dynamo (Christensen, 2011;Christensen and Wicht, 2015;Wicht and Sanchez, 2019), and mapping out the dominance of inner-core convection, rotation, or translation depending on its viscosity and conductivity (e.g., Deguen et al., 2013).

Boundary and initial conditions 865
After choosing the equations that will be solved and the model geometry, both initial and boundary conditions are needed to solve the numerical model. The solution of the numerical model will depend on the used initial and boundary conditions, so it is important to choose them carefully (Quinquis et al., 2011;Chertova et al., 2012).

Boundary conditions
The at the boundary. This is typically used to mimic the 660 discontinuity at the bottom boundary of asthenospheric-scale models.
(5) Prescribed stresses are Neumann boundary conditions as the stresses relate to the velocities through the derivatives. They can be used to mimic e.g. plate push. (6) An open boundary (Neumann) is a special case of prescribed stresses, where material can freely leave the domain (Chertova et al., 2012). In addition, infinity-like or 'external' boundary conditions can be applied where a boundary is modelled as being far away (Gerya, 2019). This is typically applied in lithospheric-scale models, where an 900 external free or no slip boundary is applied to the bottom boundary, mimicking a free or no slip boundary at a distance L from the actual boundary. Similar to this, a Winkler boundary condition assumes an isostatic equilibrium at the bottom boundary, analogous to applying a free surface boundary condition at the bottom of a lithospheric-or crustal-scale model (Burov et al., 2001;Yamato et al., 2008). In combination with a free slip boundary condition, a layer of less dense (⇢ ⇡ 1) 'sticky air' (or 'sticky water') material is often used to model topography in methods that use an undeformable mesh and can therefore not 905 employ a free surface boundary condition (Section 3.3, Figure 4) (Schmeling et al., 2008;Crameri et al., 2012).
Another type of commonly used boundary conditions in geodynamics are the periodic boundary conditions. They are different in nature from the purely mathematical boundary conditions listed above, as they do not explicitly prescribe any part of the solution. Instead, they 'link' boundaries together to approximate a larger (or infinite) system of which the model setup is Boundary conditions can be used to drive the system by, e.g., prescribing the velocities of plates resulting in lithospheric extension or convergence. Hence, the modeller could assimilate data on plate motions from the geologic record into the model through the boundary conditions to improve the predictive power of the model (Colli et al., 2015). Similarly, stress boundary conditions can be implemented to simulate a load on the model, such as an icesheet. The data assimilation into the models via 920 boundary conditions is typically used for specific modelling studies (Section 5.2).
Considering boundary conditions for the energy equation for models of the Earth's mantle and crust, it is common practice to prescribe fixed temperatures (Dirichlet) at the Earth's surface and at the core-mantle boundary, and prescribed heat fluxes On the inflow part of the boundary, it is useful to prescribe the temperature (e.g., slab age), whereas on the outflow part of the boundary, insulating Neumann boundary conditions can be used.
It is also possible to constrain degrees of freedom inside the domain. For example, in lithospheric-scale models, a velocity 935 prescribed within the slab (i.e., not at the boundary) is an example of an internal driving force to prescribe subduction in the absence of initial slab pull (e.g., van Dinther et al., 2013b).
The choice of boundary conditions can alter the modelling results with, e.g., different mechanical boundary conditions on the sidewalls of the mantle affecting the resulting subduction behaviour (Chertova et al., 2012;Quinquis et al., 2011). It is also important to choose boundary conditions consistent with the rest of the model (Section 6.2). Hence, the modeller should be 940 careful when selecting the boundary conditions and keep in mind how they affect the produced model results.

Initial conditions
Initial conditions are required for time-dependent equations. Together with boundary conditions, they control how the model evolves over time and are required to set up and drive the model. Since the conservation of energy (6) contains a time derivative @T /@t, we need to define an initial temperature field. This is difficult, since we have even less knowledge about what the Earth 945 looked like at depth in the past than the present-day. Luckily, there are some useful strategies to come up with reasonable initial conditions. In the case of oceanic lithosphere, the modeller can choose the age of the lithosphere and calculate the corresponding temperature profile according to well-established cooling models, such as the half-space cooling model or the plate cooling model (Turcotte and Schubert, 2012). For the continental lithosphere, an initially linear temperature profile or slightly more complex geotherm can be prescribed (Chapman, 1986). In the regions of the mantle that are not part of the thermal 950 boundary layers, it is often a reasonable assumption to start from an adiabatic temperature profile, which is the temperature path along which material in a convecting mantle moves if it does not exchange energy with its surroundings (see also adiabatic heating). As the conservation equations of mass and momentum do not contain time derivatives of the velocity or pressure, we do not need to provide an initial velocity or pressure field. However, from a numerical point of view, a reasonable initial guess for the pressure or velocity can reduce the amount of iterations needed in iterative solvers and hence speed up computations.

955
This becomes particularly relevant when using pressure-and strain-rate dependent rheologies.
The initial conditions also include the initial compositional geometry and material layout and/or history in the model, since the transport equation (15) contains a time derivative as well. For general parametric studies of the dynamics of a system, simple geometric blocks could be used to set up the initial geometry, representing for example the lithosphere. For more complex models, the initial conditions could be inferred from (regional) tomographic models. For complex models of specific regions, it 960 is often difficult to manually create detailed geometries that correspond with geologic or tomographic observations. Therefore, several tools have recently been developed, such as geomio (Bauville and Baumann, 2019), the SlabGenerator (Jadamec and Billen, 2012), and the Geodynamic World Builder (Fraters et al., 2019b), to automate and simplify setting up models with complex, 3D geometries. Another choice the modeller has to make concerns the initial chemical heterogeneity present in the Earth's mantle. The simplest choice is to assume that the mantle is homogeneous or has been mixed so well that heterogeneities 965 are on such a small length scale that they do not influence mantle dynamics. However, we know from geochemical data that the mantle is chemically heterogeneous. In addition, subduction zones continuously introduce new chemical heterogeneities into the mantle. Hence, another option is to initialise the model with chemical heterogeneities on a given length scale, for example representing a model like the marble-cake mantle (Allègre and Turcotte, 1986), or to include distinct chemical reservoirs in the model setup. The initial state of deformation in terms of accumulated plastic or viscous strain or the state variable in rate-and-970 state friction laws can also be prescribed as initial conditions. Such initial strain can represent preexisting heterogeneity in the crust or lithosphere formed through deformation prior to the model start time. The top of the model can be used to assimilate topographic data into the model for a certain region or from another model (Section 3.3).
To initially drive the model in the absence of driving boundary conditions, specific initial conditions inside the model domain, so-called 'initial perturbations', can be applied. A frequently used example in mantle modelling is a density or temperature 975 perturbation in the middle of, or distributed throughout, the mantle. Such a perturbation ensures that subsequent deformation starts immediately, rather than after the accumulation of numerical rounding errors over time, and it localises the deformation in the region the modeller is interested in, e.g., a plume rises in the middle of the model.
Another common example of initial conditions is the so-called weak seed in numerical models: a small zone with artificially lowered strength used to localise the deformation in the model in the desired region. Without a weak seed, numerical inaccu-980 racies will determine the location of instabilities and subsequent deformation, typically at the boundaries of a model. This is undesirable, because this results in irreproducible, random models and deformation influenced directly by boundary conditions.
Hence, weak seeds are necessary in models and can represent naturally occurring heterogeneity in the Earth such as a previous fault zone or a region of melting (Buck et al., 1999). Weak seeds are commonly used in lithospheric-or crustal-scale models of rifting (Huismans et al., 2005;Allken et al., 2011), strike-slip (Schierjott et al., 2020), subduction (Erdos et al., 2015), and 985 continental collision to localise the deformation and force the model to behave in such a way that the relevant process can be studied (Gray and Pysklywec, 2012;Yang et al., 2020). They can take numerous shapes and sizes and the lower strength of the weak seed can be achieved through many different methods, including a mechanically weaker seed (i.e., lower friction coefficient (Allken et al., 2012;Erdos et al., 2019)), initial strain or a temperature anomaly such as a locally raised thermal lithosphere-asthenosphere boundary (Brune et al., 2017). The choice of weak seed affects the numerical results, as it could for 990 instance lead to different modes of rifting (Dyksterhuis et al., 2007). However, most studies argue that the weak seed does not significantly alter their conclusions. The effect of weak seeds and robustness of the models can vary between codes.
The modeller should always keep in mind that the initial conditions determine the model outcome. That is, after all, their purpose, since otherwise there would be no localised deformation or initial drivers. Efficiently starting the model is solely at the discretion of the modeller, who aims to artificially mimic a process they are interested in. Since these initial conditions, 995 in combination with the boundary conditions, are critical for the model development, the choices the modeller makes are sometimes referred to as the hand of god that helps the models along in the beginning. It is therefore important that initial conditions and their effects on the model are acknowledged and described alongside other important details of the modelling approach (Section 8.1).
After the code has been successfully tested and benchmarked (Section 4), every individual model setup with its particular modelling strategy (Section 5.2) and initial and boundary conditions (Section 5.3) should be carefully validated to make sure that it contains no detectable flaws, is internally consistent, and represents the geological problem to the best of our knowledge.

Common numerical problems
Despite successful code verification (Section 4), nonphysical artefacts can still arise from the specific model setup, for example from incorrect initial or boundary conditions (Section 5.3). These issues can usually be spotted through monitoring solver convergence behaviour and visual inspection of the solution throughout the model evolution, with model breakdown (i.e. a crash of the program) and unexpected behaviour being the most obvious red flags. In this section, we describe a number of common problems and their potential solutions.
A resolution test should be standard to check whether a certain model, or model aspect, is mesh-dependent or not. So, the 1010 modeller should check the change in model results with higher mesh resolution. In the ideal case, from a certain resolution onward the solution no longer changes significantly or the spatial discretisation error becomes smaller than other errors like that resulting from time discretisation, i.e., the numerical solution has converged. It is desirable to use a grid spacing where, for example, the thermal boundary layers or crustal compositions are well resolved and the resolution does therefore not affect the model evolution anymore. However, it is not always possible to completely avoid grid-dependency. For example, 1015 most implementations of brittle (plastic) deformation (e.g., Duretz et al., 2020) do not include an internal length scale and are therefore grid-dependent (Figure 10c). This grid dependency causes shear zone width to keep decreasing with increased resolution. There is an active research effort to include internal length scales that can limit grid-dependency (Lavier et al., 2000;Choi and Petersen, 2015;Duretz et al., 2020). The practical solution is to have one fixed resolution for the affected domain throughout the modelling study after assessment of changes in the overall model behaviour with resolution. Another problem that can occur when mesh deformation is allowed in finite elements (i.e., in Lagrangian and ALE methods; Figure 5) is distortion of the mesh elements. The quality of the mesh decreases when element aspect ratios change too much from 1 for triangular as well as quadrilateral elements (e.g., the ratio of the width and height of a 2D quadrilateral element (Figure 4), and similarly in 3D), and the accuracy of the computations therefore decreases. When the element distortion is too  Lemiale et al., 2008;Kaus, 2010;Spiegelman et al., 2016;Glerum et al., 2018). The angle of the shear bands is lower for a coarser resolution, while thinner shear bands at a higher angle form for a resolution that is 4 times higher. large, the computation will crash or at the very least become extremely inaccurate. To avoid such a distortion of the elements, local or global remeshing (i.e., mesh regeneration) or mesh smoothing can be applied (e.g., Zienkiewicz et al. (1995) Burov and Cloetingh, 1997;Sternai, 2020), but remeshing of the free surface also introduces some numerical diffusion (Jammes and Huismans, 2012) that could help remedy instabilities.
Visual inspection of the modelling results can uncover other issues. For one, smaller features (e.g. a subduction interface) can be seen to spread out over time and disappear; this can be due to diffusion or smearing of the advected field in the grid-based 1040 advection method. On the other hand, steep gradients of advected fields can lead to oscillations of these fields normal to the gradients. Mitigating such under-and overshooting requires more diffusion or different stabilisation algorithms of advection (e.g. Lenardic and Kaula, 1993;Brooks and Hughes, 1982;Guermond and Pasquetti, 2011).
Unstable yet very popular finite element pairs such as the Q 1 ⇥ P 0 element (Section 3.2) are prone to spurious oscillations and an element-wise chequerboard pressure pattern (Figure 10b as shown in Figure 18 of Thieulot et al. (2008) or in Donea 1045 and Huerta (2003)). Before this pressure can be used in rheological expressions it must be postprocessed and smoothed so as to remove the chequerboard mode (Lee et al., 1979).

Internal consistency of the model
Internal inconsistencies can arise from disagreements in the modeller's choices in terms of boundary conditions, density formulations in the different governing equations and the equations' coefficients. Not all inconsistencies are easily detectable or 1050 manifest themselves as numerical problems. For example, when the net prescribed in-and outflow through the model boundaries is not (close to) zero, while a model is assumed incompressible, volume is no longer conserved and the solver for the Stokes equations might crash or return a nonsensical solution. When a free surface is used, this problem might be overseen, as the surface can rise or fall in response to a net increase or decrease of volume, respectively. This physical inconsistency is also harder to detect in compressible models. Another example are prescribed surface velocities based on, for example, plate 1055 reconstructions models, which can add unrealistic amounts of energy into the modelled system.
Care should also be taken that the assumptions made to simplify the treatment of density in the governing equations (see Section 2.1 and specifically 2.1.4) agree with one another. The simplest accepted combination of simplifications is the Boussinesq Approximation.
The thermodynamics of Earth materials are very complex, especially in multi-phase, multi-material systems, hence they are 1060 often simplified in the numerical model. For example, in nature the thermal expansivity varies both smoothly and abruptly with depth (e.g. Figure 7 of Stixrude and Lithgow-Bertelloni (2007)), but in models it is often taken as constant or merely smoothly increasing with depth. At first, the material properties described in Section 2.2.5 may appear to be independent. However, the definition of properties like density, thermal expansivity, specific heat, and compressibility need to satisfy thermodynamic relations in order for them to be consistent. These thermodynamic relations can be derived through thermodynamic potentials.
For example, the thermodynamic potential of the thermal expansivity is defined as : where ↵ is thermal expansivity. It defines the percentage increase in volume of a material per degree of temperature increase at constant pressure.
The thermodynamic potential of the isothermal compressibility is defined as : where T is the isothermal compressibility and the potential describes the percentage increase in density per unit change in pressure at constant temperature. For the isobaric heat capacity, the thermodynamic potential is defined as :

1075
where C p is the isobaric heat capacity and the potential is defined as the ratio of the increment of heat q added to the material to the corresponding change in temperature T . The thermodynamic potentials also imply that: All these relations have to be fulfilled at all temperatures and pressures. Consequently, it is often not immediately apparent if a given material description is thermodynamically consistent or not, and equations of state used in the geodynamic literature do 1080 not always take this into account , Section 6.10). For example, it would be thermodynamically inconsistent to take thermal expansivity, compressibility at constant entropy (along an adiabat) and isobaric specific heat all as constants, if the temperature varies in the model, because, from Eq. 22, any temperature dependence in ↵/⇢ leads to a pressure dependence of the specific heat. fact, a model of a natural system like the geodynamic models described here can never truly be validated (Oreskes et al., 1994). This is because natural systems are never truly closed. For example, neither the rocky planetary surface, nor the core-mantle boundary are true closed boundaries devoid of any temperature, compositional, or mechanical exchange with the outside world.

Testing the model against observations
As such, all model results are always non-unique and the model cannot be validated even if compared to a natural observation.

Model analysis 1100
After ensuring the accuracy, consistency and applicability of the model results, these can now be used to address the hypothesis the modeller set out to test according to a particular modelling philosophy (Section 5.2). The raw model output therefore requires analysis. While analysing the geodynamic model results, the modeller has to keep in mind all simplifications made during model setup (e.g., what forces and processes were included), the initial and boundary conditions (e.g., whether subducting plates are free to move laterally or are attached to the boundaries), the resolution (e.g., whether the resolution is high enough 1105 to resolve a certain process), and all other numerical and physical model assumptions and uncertainties. Most importantly, the model cannot be mistaken for the real Earth (or any other real planet).
Model analysis includes visual (qualitative) diagnostics and quantitative diagnostics. These two important, partly overlapping, aspects are discussed in detail below.

1115
To cover the wide range of potential visualisation products, a multitude of visualisation programs is available. Some of the commonly used software packages are gnuplot for creating 2D and 3D graphs, GMT for making maps and plotting data on maps (Wessel and Luis, 2017), MatLab for scripted analysis and plotting of model data, matplotlib for static and interactive plotting in Python, ParaView for looking at 2D and 3D data sets both interactively and in batch jobs (Ahrens et al., 2005;Ayachit, 2015), and similarly VisIt (Childs et al., 2012). There are even fully immersive visualisation options like CAVE (Billen et al., 1120(Billen et al., 2008 and virtual reality display technologies (Kreylos and Kellogg, 2017) that overcome the necessity to project 3D data to 2D surfaces for visual inspection. However, no matter the software, a modeller has to pay careful attention to circumvent the most common visualisation pitfalls for effective graphics as is outlined in Section 8.2.

Quantitative model diagnostics
Mere visual inspection of the model results is not sufficient to analyse and interpret the outcome of the simulations; a quan-1125 titative analysis of the results is also required. Deciding what specific post-processing is to be done can be a time-consuming process. There are a range of non-dimensional numbers that can be calculated to characterise the flow of fluids in a range of geodynamic environments. If multiple physical processes influence the behaviour of a system, the non-dimensional numbers derived from the governing equations (Section 2.1) can be used to analyse the relative importance of each of these processes.
Examples are the Rayleigh number, Ra, the Knudsen number, Kn, the Mach number, M, Argand number, Ar, the Ekman 1130 number, Ek, the Reynolds number, Re, the Peclet number, Pe, the Prandtl number, Pr, and the Nusselt number, Nu. The Rayleigh number, in particular, is a widely used non-dimensional number characterising the vigour of buoyancy-driven flow.
Many of these non-dimensional numbers can be defined differently for different environments, but they can all be insightful diagnostic metrics. With increasing model complexity, newer diagnostic numbers have been defined, such as plateness (Tackley, 2000;Alisic et al., 2012;Stein et al., 2014).

1135
Further analysing and diagnosing a model then varies with the modelling approach that has been taken (see Section 5.2).
For specific modelling (i.e., model approaches directly comparing to an observation to understand the origin of a specific state into what causes a specific state of the system. For generic modelling (i.e., models used to reproduce basic fluid dynamics to understand how a specific system works), a modeller should diagnose whether there are individual models that exhibit similar behaviour within the model parameter space and whether there are controlling parameters that can, when changed, cause a switch from one regime to another. In general, often-used diagnostic quantities to define regimes are the root-mean-square (RMS) and maximum values of certain 1150 model parameters like velocity or temperature. For example, the root-mean-square velocity over time can be used to check whether a model has reached steady-state or shows signs of periodicity (e.g. Tosi et al., 2015). For subduction models, the slab dip, slab tip depth, slab sinking velocity, plate age at the trench, and the interaction of the slab with the mantle transition zone are often monitored (e.g. Garel et al., 2014;Crameri and Lithgow-Bertelloni, 2018). Specific isotherms, fault length and offset, rift symmetry, and the migration of fault activity are often tracked in extensional models (e.g. Huismans and Beaumont, stagnant lid versus mobile lid), and the degree of layering are often extracted (Heyn et al., 2018;Foley, 2018;O'Neill et al., 2016).

Automated model diagnostics
While some model analysis can be done by hand, more elaborate post-processing that becomes increasingly popular nowadays 1160 needs to be automated using open-source, testable, and extendable algorithms and shared as user-friendly software (Greene and Thirumalai, 2019). To achieve such next-generation post-processing, like plate boundary tracking, extraction of lithosphere thickness, or computing the dynamic or isostatic topography (Crameri, 2017b), the output of geodynamic codes should ideally use a standard, widely-accepted format and include metadata that can also be accessed by machines. Accessing individual subsets of the data, like individual time steps or parameter fields, should be straightforward. Model output should therefore be 1165 independent from computational details, like the number of computational cores the results have been produced on.
A few software packages that allow for automated post-processing and diagnosis of geodynamic models are available to support geodynamicists with analysing their increasingly complex models, and the large datasets originating from them. However, such tools are rare, because while most individual researchers spend a large amount of time in coding post-processing scripts, they often do not share those scripts with the geodynamics community. Moreover, scripts that are shared in the context Scientific colour maps (Section 8.2). Moreover, StagLab is versioned and extendable, and it can process output of more than one mantle convection code, while also offering low-effort compatibility for any other machine readable output.
Other recent developments include the automated comparison of observations to model predictions to find the smallest misfit between the two. Such statistical and probalistic inversion methods help determine the model parameters, e.g., mantle viscosity or crustal density, that result in the best fit of the model solution with the observed quantity through forward geodynamic 1185 modelling (e.g. Baumann et al., 2014;Baumann and Kaus, 2015;Li et al., 2017;Reuber et al., 2018a;Morishige and Kuwatani, 2020;Ortega-Gelabert et al., 2020).
Scientific results are only of value if they are communicated to the wider scientific community. No matter whether they are spoken or written, the first aspects to get right when communicating science concerns letters, words, and phrases. Since geodynamic modellers, like most other life forms, tend to learn most effectively by observing and copying other fellows, it is no surprise that we tend to speak and write in a similar way to our mentors, peers, and friends. While there is generally nothing wrong with that process, it does, however, make for an excellent breeding ground for problems related to semantics that can lead to serious miscommunication (e.g., Crameri, 2018a). Written and spoken communication therefore needs to be held dear and handled carefully by a geodynamic modeller.

1195
The semantics behind a modelling publication or presentation need to be in tune with the approaches taken in the modelling itself. If a modelling study is suitably communicated, there will be less misunderstanding about what the presented model stands for, what it does not stand for, and what the drawn conclusions mean.
Because geodynamic models are per definition simplifications of a natural system (Section 1.1), their individual features should not be mistaken for an exact replica of their natural counterparts. When communicating modelling aspects, semantic 1200 differentiation between the feature in the model and in nature helps to avoid confusion. For example, refer to 'the modelled slab' instead of 'the slab' or 'the subducted Nazca plate'. On a similar note, distinguish between 'thermo-chemical piles', which are collections of material with different thermal and chemical properties as the surrounding mantle, and 'LLSVP's', which are observed regions of low shear wave velocity along the core-mantle boundary.
In addition, care has to be taken with absolute statements, like 'X on the Earth is due to Y', when drawing conclusions from the model results. As discussed in Section 5.2, models can only demonstrate a certain likelihood of an hypothesis, and in particular specific modelling studies should acknowledge that there is no guarantee that the suspected control mechanism is the actual, nor the only, controlling mechanism. Statements like 'X on the Earth is likely due to Y', 'Y is a potential explanation for X', and 'if our assumptions A, B, and C are fulfilled, then XYZ will probably happen in the Earth' are more correct and prevent misconceptions.

1210
Communicating a geodynamic modelling study, however, goes beyond semantics. The suitable words and phrases are most effective when combined with an appropriate manuscript structure and effective still and even motion graphics. Combined, these forms of communication make a new scientific insight accessible.

Structure of a geodynamic modelling manuscript
Peer-reviewed scientific papers are essential to disseminate relevant information and research findings. In particular, it is im-1215 portant to make those results understandable and reproducible in the methods and results sections. Reviewers will criticise incomplete or incorrect method descriptions and may recommend rejection because these sections are critical in the process of making the results reproducible and replicable. In this section, we briefly review the structure of a manuscript and highlight the parts required in a geodynamic numerical modelling study (Pusok, 2019).
the structure of scientific papers evolved from mainly single-author letters and descriptive experimental reports to a modernday comprehensive organisation of the manuscript known as 'theory-experiment-discussion' (Audisio et al., 2009;Liumbruno et al., 2013). The formal IMRAD structure (i.e., Introduction, Methods, Results, And Discussion/conclusions) was adopted in the 1980s, and at present, is the format most widely used and encouraged by scientific journals. The IMRAD structure facilitates modular reading and locating of specific information in pre-established sections of an article (Kronick, 1976). In 1225 geodynamics, the general structure of a manuscript follows the IMRAD structure ( Figure 11).
A good introduction should answer the following questions: what is the problem to be solved? What previous work has been done? What is its main limitation? What do you hope to achieve? How do you set up your investigation? One major mistake is to attempt to do an extensive literature review in the introduction, which often goes off-topic. The introduction serves as the stage to lay out the motivation for the study, and any background reading should focus on the question being addressed.

1230
The methods section is considered one of the most important parts of any scientific manuscript (Kallet, 2004). A good methods section allows other scientists to verify results and conclusions, understand whether the design of the experiment is relevant for the scientific question (validity), and to build on the work presented (reproducibility and replicability) by assessing alternative methods that might produce differing results. Thus, the major goals of the methods section are to verify the experiment layout and allow others to reproduce the results. Here, we outline standards for reporting the methods in numerical geodynamic 1235 modelling.
First, the methods should be plain and simple, objective, and logically described and should be thought of as a report of what was done in the study. Unstructured and incomplete methods can make the manuscript cumbersome to read or even lead the reader to question the validity of the research. Generally, journals have guidelines on how the methods should be formatted, how many words can be used, but not necessarily what it should contain because it varies from field to field. The 'who, what, when, 1240 where, how, and why' order proposed by Annesley (2010) breaks the methods section down in the following questions: who performed the experiment? What was done to answer the research question? When and where was the experiment undertaken?
How was the experiment done, and how were the results analysed? Why were specific procedures chosen? The answers to these questions should be adapted to every field (i.e., in geodynamics, 'results were obtained using code X on cluster Y'). Here, we focus on methods that have primarily theoretical (mathematical and physical) and numerical (computational) components, but 1245 geodynamics studies may have other aspects such as a data component (e.g., collection and post-processing of data) or an analogue component (e.g., laboratory experiments). Figure 11 shows a breakdown of the most important elements of a manuscript and of the methods section in particular.
The methods should start with a brief outline (one paragraph) describing the study design and the main steps taken to answer the scientific question posed in the introduction. The outline should be logical and go from theoretical elements, to numerical aspects, to analysis and post-processing. First to be described is the theoretical framework. This includes the mathematical and physical concepts used in the study including the governing equations, constitutive equations, initial and boundary conditions, and any other relevant theory describing the model setup. Figure 11. Manuscript structure for a geodynamics numerical modelling study following the IMRAD structure. In particular, the methods section has to include a description of the physical and numerical model, the design of the study, and of any techniques used to visualise and analyse the numerical data. This is followed by a section on the computational approach explaining how the theory and the model are translated into computer language (Section 3). This includes details on numerical methods (discretisation and other numerical techniques, After the model setup has been explained, the methods should contain a section describing the design or layout of the study in detail. What is being tested/varied? How many simulations were performed in terms of model and parameter space? For example, one can use different model setups (i.e., lithosphere-scale and mantle-scale subduction models) with varying parameters in the same study. Why perform those simulations and vary those parameters? A summary table is handy to indicate all simulations that were run and which parameter was varied in which run.

1265
Analysis, visualisation and post-processing techniques of numerical data must also be described in the methods section. This is a step generally ignored, but it is important to be open about it, e.g., 'visualisation was performed in ParaView/Matlab, and post-processing scripts were developed in python/matlab/unicorn language by the author'. If the post-processing methods are more complex, the author can provide more details (i.e., statistical methods used for data analysis). It is also good exercise for peer-reviewing and community engagement to provide these post-processing scripts (Section 9).

1270
Information should also be given on code and data availability. This has originally been part of the methods section, but recently journals have introduced data management requirements (Section 9) and this information may have a designated location in the manuscript template. However, it is good practice to write this information in the methods section. The authors should indicate code availability, code verification, and input files or other data necessary to reproduce the simulation results (e.g., installation guides). Additional questions to be answered in the methods section are: where were the simulations per-1275 formed, and on how many cores? What was the average runtime of a simulation? Can the model results be reproduced on a laptop/desktop computer or is access to a cluster required?
Before moving to other sections, model assumptions need to be stated clearly, in either the description of the theory or the numerical approach. Geodynamics is a field in which we take a complex system like the Earth or another planetary body and simplify it to a level that we can extract some understanding from (Section 5). In doing so, we rely on a physically consistent 1280 set of assumptions. It is important to bear in mind that this set of assumptions may not always be obvious to the reader. As long as assumptions are explicit and consistent (i.e., through clear and honest communication), the reviewers and readers will find fewer flaws about the study. It is always good practice to write a complete methods section for every manuscript, such as the one described here. However, some journals will ask for a short version to be included in the main manuscript and have the complete methods section in a separate resource (i.e., in the appendices, supplementary data, supporting information, online 1285 repository).
Complementary to the methods section, the results section should be a report of the results obtained. The main goal of the results section is to answer the question of what you found. Any interpretation of the results or reference to other studies should be reserved for the discussion. For example, results in a mantle convection model might show that dense material accumulates at the bottom of the domain (i.e., core-mantle boundary). The interpretation of these results is that they provide a mechanism 1290 to explain how LLSVPs (large low shear velocity provinces) have formed. Illustrations, including figures and tables, are the most efficient way to present results (Section 7). However, authors must only include material and information that is relevant to demonstrate the scientific arguments discussed in the next section. Therefore, to avoid distraction, writers should present additional data as supplementary materials, e.g., movies of the whole simulation, whereas only a few snapshots are provided in the main body of the manuscript.

1295
The discussion section relates all the questions in the manuscript together: how do these results relate to the original ques- At this point in preparing the manuscript, the authors have all the necessary elements to write the abstract, conclusions, and come up with a descriptive title. Both the abstract and conclusion summarise the entire publication, but in a different way; one as a preview and one as an epilogue, respectively. It is crucial to focus a paper on a single key message, which is communicated in the title.
In the end, every scientific manuscript has additional components such as the references, acknowledgements, supplementary 1310 material, software and data availability, and author contributions ( Figure 11) that contain further information about how the study was funded, conducted, and shared with the community.
In this section, we have primarily referred to scientific articles, but scientific manuscripts can also be reviews, editorials, and commentaries. The structure and contents of these manuscripts differ for each type. Each publisher and journal has its own style guidelines and preferences, so it is good practice to always consult the publisher's guide for authors. Finally, even though 1315 scientific manuscripts may have a rigidly defined structure due to journal guidelines, there is still plenty of flexibility. In fact, the best manuscripts are creative, tell a story that communicates the science clearly, and encourage future work.

Effective visualisation
There are many different ways to visualise geodynamic models and it is challenging to figure out how to do so most effectively.
However, avoiding the most-common visualisation pitfalls is the best start for any modeller looking into visually communicat-1320 ing important results across the research community and possibly beyond. The key aspects to remember when creating figures, thereby preventing misleading visual impressions, are the following: (1) Scales, like graph axes and colour bars, must always Figure 12. Effective visualisation through a scientific use of colours. Non-scientific colour maps (a,b) like rainbow always misrepresent data, are often not intuitive, and inaccessible to a large fraction of readers, while scientific colour maps (c,d) like lajolla or vik (Crameri et al., 2020) ensure unbiased and intuitive data representation and are inclusive to readers with colour-vision deficiencies and even colour blindness. be included to allow quantification of data values.
(2) Bar plots must always have a zero baseline, to not mislead the reader with altered relative bar heights.
(3) Pie diagrams should be avoided as angles and areas are not easily quantifiable by the human brain and therefore are not directly comparable to each other. These problems are exaggerated when pie charts are displayed as 1325 3D structures, which causes the values underlying the pieces closest to the viewer to appear artificially larger than the others.
(4) Heatmaps (i.e., plots with differently coloured tiles) should have numbered tiles that include the data value, as surrounding colours heavily distort the impression of a given colour, which can mislead the viewer's perception of the underlying data values significantly (Crameri et al., 2020). (5) Colours must be applied in such a way that data is reflected correctly and is inclusive to all readers (Crameri et al., 2020). Scientifically-derived colour maps exist, like Colorbrewer, MPL, Cividis, CMO-1330 cean, CET, and Scientific colour maps (Kovesi, 2015;Thyng et al., 2016;Nuñez et al., 2018;Crameri, 2018b), and must be chosen over unscientific default colour maps like rainbow, jet, seis, and others ( Figure 12). (6) Visualisation should be subject to the same scientific scrutiny as other research methods to communicate data and concepts truthfully, accessibly, and clearly.
All aspects of a figure need to be explained and understandable. While filtering, resampling or otherwise postprocessing model results instead of plotting raw data can improve the message purveyed by the figure, such actions should be mentioned in 1335 the figure caption. Some numerical codes work, for example, in dimensionless numbers (Section 7.2) and require scaling of the model output before they can be related to observations. However, too much information jammed in a figure can easily render the figure unusable to the reader. Again, the modeller has to simplify enough to reach the sweet spot, without oversimplifying the figure (compare with Section 5.1). Everything that can be removed without losing key information should be removed.
Unnecessary and/or duplicated axes labels, e.g., those repeated across multiple panels, can be removed. The same applies 1340 to other figure aspects like colour bars. Colour bars in multiple figure panels applied to the same parameter should always maintain the same range (i.e., map the same colour to the same data value) and, if they do, displaying just one colour bar is sufficient.
Displaying 3D models effectively is challenging, and somewhat arbitrary, as the third dimension is often difficult to convey in a still image. Given the current dominant occurrence of non-interactive, two-dimensional canvases (e.g., the pdf format), 2D 1345 slices of parameter fields often represent the model more effectively than 3D volumes. The combination of various data sets, like flow characteristics on top of a temperature field, can be effective, but is also challenging. Velocity arrows, for example, should not overlap, nor distract from the remaining important content of the figure. If the velocity in a 3D visualisation is displayed using arrows, they should be coloured according to their magnitude because their lengths are distorted by the 3D perspective. Stream lines and a coloured contour plot of the velocity field often provide a more suitable solution to display the 1350 flow direction and patterns, and its velocity magnitudes, respectively.
An uninformed, unscientific use of colours can not only exclude a large fraction of the readership, for example through hardly distinguishable colour combinations for readers with a colour-vision deficiency (like the most common red-green colour blindness), but also significantly distort the underlying model data visually (e.g., Crameri, 2017a, 2018c, and Figure 12). In fact, the distortion can be more than 7% of the displayed data range (e.g., a temperature variation of 20 degrees could look 1355 like 24 degrees in one part along a colour axis that ranges between zero and a hundred degrees, while it looks like 17 degrees in another part), which in most cases easily amounts to the biggest error in a modelling study. Therefore, modellers should avoid using default colour maps blindly, but instead look for scientific colour maps (Crameri et al., 2020). If a software does not offer them, or at least allows to import them, users should contact the software developers. Hence, scientific colour maps should be the default in common geodynamic codes and visualisation programs.

1360
Scientific colour maps are perceptually uniform to prevent data distortion, and they are perceptually ordered to make data intuitively readable, colour-vision deficiency friendly and optimally readable in black and white to include all readers. Suitable scientific colour maps (e.g., Crameri, 2018b) that offer fair data representation and are fully inclusive have become readily usable (either as built-in options or after importing them) with all major modelling, visualisation, and graphics software, like MatLab, Python, R, GMT, QGIS, ParaView, VisIt, Mathematica, gnuplot, GIMP, and Inkscape. Indeed, the scientific colour 1365 maps suite (e.g., Crameri, 2018b) offers custom colour map classes for any data type. The sequential class (e.g., batlow) is best suited for monotonically increasing data like a temperature field. The diverging class (e.g., vik) is best suited for zero-centred data like a horizontal velocity field. A multi-sequential class (e.g., oleron) should be used with zero-centred, non-diverging data like a combination of ocean bathymetry and land topography and the cyclic class (e.g., romaO) is tuned for circular data like azimuthal direction. Moreover, these scientific colour maps are offered as continuous type (e.g., batlow) to visualise small-1370 scale data variation, discrete type (e.g., batlow10) to highlight similar data values more clearly, and as categorical type (e.g., batlowS) to clearly differentiate not ordered data for line or scatter plots. More in-depth information, a clear user guide to using colour and colour maps, and a concise list of the major currently available scientific colour maps are provided in Crameri et al. (2020).
9 Software, data, and resource management 1375 Just like any other study, numerical modelling studies should be reproducible (i.e., the same researcher doing the same computation on the same machine produces the same measurement at stated precision), and most importantly replicable (i.e., different researchers can obtain similar enough results from independent development) (Bollen et al., 2015;Goodman et al., 2016). Note that the actual terms used for these concepts vary across the sciences (Association for Computing Machinery, 2016; Karlsen, 2018;Plesser, 2018), with the terms reproducible and replicable for example used interchangeably. Repro-1380 ducibility and replicability imply that the software used to conduct the study as well as the specific model setups, installation environment specifications, and post-processing workflow should be available (Perkel, 2020) to interested peers, and preferably, everybody. More and more, scientific journals are requesting or even requiring the publication of data and software along with the manuscript. Although the requirements vary per journal, it is good practice to adhere to these principles for every publication.

1385
Before development starts, software developers and modellers involved in the development of the software they use should consider setting up a software management plan (SMP) (see this SMP template of the Software Sustainability Institute). This includes, but is not limited to, the following questions: what software will be produced? Who will use the software? How will the software be made available to the users? What support will be provided to the users? How will the software contribute to availability?
In the hero codes development strategy, one or only a few people are responsible for development and maintenance of a modelling software package (Thieulot, 2017). Community software efforts are often managed through version control software like git and svn and corresponding platforms like GitHub, GitLab, and Bitbucket (e.g., ASPECT, Kronbichler et al.

1395
These platforms provide open as well as private repositories for software development, issue tracking, automated testing, and project management, greatly simplifying the points addressed in the SMP. They also facilitate contributing the modeller's own development efforts to the main repository of the software such that they are available to other researchers. Moreover, extraction of statistics on the number of downloads, users, and contributors is made easy. However, these platforms themselves do not provide persistent identifiers and are not considered archiving facilities according to the FAIR Data Principles discussed below. input files, software, and machine specifications are properly described, identified, and made available, numerical data does not have to be archived, as the study can be reproduced from those elements. However, accessible model results can save the computational resources needed to recreate the model results, serve reviewer assessment and educational purposes, and post-processed results can be of interest to the general public (e.g. movies on YouTube or graphics on figshare; Pusok, 2020).

1420
To help modellers with the implementation of the FAIR Data Principles, publishers and data repositories formed the coalition COPDESS. The COPDESS website explains why and how to archive and publish software and accompanying data and includes FAIR author guidelines, links to the guidelines of associated journals, and a FAIR-aligned repository finder from the Enabling FAIR Data project (see re3data for a larger repository portal). An extensive explanation of the FAIR Data Principles research data management starts with addressing the following points (after Science Europe, 2018): (1) description of data, its collection/production and the reuse of existing data sets; (2) data metadata and management; (3) storage, backup and security during the project; (4) legal and ethical requirements, codes of conduct; (5) what, how and where will data be stored, accessed and identified?; (6) who will deposit and maintain the data?
Data repositories can be subdivided into institutional repositories, domain-specific repositories (e.g. EarthChem, IRIS, PAN-1430 GAEA and HydroShare), thematic data repositories (which differ from domain-specific repositories by having to transform the data to the repository's format yourself) and general repositories, like figshare, Dryad and Zenodo (Stall et al., 2020). The library of a modeller's institute can explain what repositories they support and what workflows already exist for archiving data with persistent identifiers. In general, institutional and domain-specific repositories provide more support and quality control in submitting the data, while general repositories do not set specific requirements for the data. Also, by depositing data in 1435 a domain-specific repository, it is more likely to be found by the target audience. Useful repositories also provide you with Not only data can have a persistent identifier, researchers can also create persistent digital identifiers, like an ORCID iD, for themselves. These identifiers can be connected to all types of research products, such as articles, grants, and peer reviews, providing acknowledgement of the work and preventing identity mix-ups. The article processing portals of many journals 1445 already allow researchers to link their ORCID iD to their profile, and these iDs are then linked to in the published manuscript.
Another development in article publication is the rise of reproducible articles, where not only the output, but the whole code used to produce the output is interactively included in the manuscript, using for example R Markdown or Jupyter Notebooks.
For massively parallel geodynamic computations, such a reproducible article is currently not feasible, but we would like to invite both authors and reviewers to critically assess the reproducibility of their geodynamic modelling studies (Pusok, 2020).

1450
One last thing modellers should consider is that numerical modelling does not come for free. As a community, we have to acknowledge the environmental impact of especially high-performance computing and data storage. In a busy year (e.g., 1 million CPU hours), computations of one researcher can emit up to 6 tons of CO 2 on a modern high-performance machine (e.g., a power draw of 400 W per processor; see Lannelongue et al. (2020) and their computing emissions calculator).
For comparison, a trans-Atlantic round-trip flight (e.g. Berlin-San Francisco) produces about 2 to 3 tons of CO 2 per person 1455 (myclimate CO 2 calculator, Ivanova et al. (2020)). A conscious effort should therefore be made to not waste computing power.
For example, short, low-resolution models should be favoured to test new implementations and setups and weak and strong scaling test results (Section 3.5) should be used to determine the optimum amount of processes for large production runs.
Furthermore, code developers should care about code optimisation from the start and could additionally provide both debug Computing and data centre cooling systems can be improved and waste heat should be used for example in district heating or direct heating of the surrounding buildings (Wahlroos et al., 2018;Oltmanns et al., 2020), although these solutions are typically beyond the modeller's control.
In addition to environmental costs, there are non-negligible financial costs to modelling. Access to high-performance machines can be expensive and a heavy entry in a modeller's budget. Moreover, the often big data that results from running 1465 numerical models needs to be stored, diagnosed, visualised, and shared. Large local or remote storage solutions, software licenses, and powerful personal computers are expensive. These financial modelling costs need to be acknowledged not only by modellers themselves but also by others, such as funding agencies. With conscious management of resources, software and data, we can ensure a fairer, more efficient and greener geodynamic modelling community.

Conclusions & outlook 1470
Geodynamic modelling studies provide a powerful tool for understanding processes in the Earth's crust, mantle, and core that are not directly observable. However, for geodynamic modelling studies to make significant contributions to advancing our understanding of the Earth, it is of paramount importance that the assumptions entering the modelling study and their effects on the results are accurately described and clearly communicated to the reader in order for them to be well understood. These assumptions are made at numerous stages of the numerical modelling process such as choosing the equations the code is 1475 solving, the numerical method used to solve them, and the boundary and initial conditions in the model setup.
Apart from acknowledging the assumptions made and their implications, it is important to view a modelling study in light of its intended philosophy. Generic modelling studies, usually characterised by extensive parameter variations, aim to understand the general physical behaviour of a system. Specific modelling studies on the other hand aim to reproduce a specific state of a specific geodynamic system and therefore rely more heavily on data comparisons.

1480
In order to make the geodynamic modelling process transparent and less prone to errors, good software management is necessary with easily repeatable code verification to ensure that the equations are solved correctly. Additionally, it is important that the model is internally consistent with regards to thermodynamics and boundary conditions. Then, for individual models, the results need to be validated against observations.
When communicating the results of a geodynamic modelling study to peers, it is important to provide both quantitative 1485 and qualitative analyses of the model. Fair presentation of the results requires clear, unbiased, and inclusive visualisation. The results should first be objectively described and presented, after which they can be interpreted in the discussion.
In addition to outlining these best practices in geodynamic numerical modelling, we have shown how to apply them in a modelling study. Taking these best practices into account will lead to clearly communicated, unambiguous, reproducible geodynamic modelling studies. This will encourage an open, fair, and inclusive research community involving modellers, col-1490 laborators, and reviewers from diverse backgrounds. We hope to set a standard for the current state of geodynamic modelling that scientists can build upon as future research develops new methods, theories, and our understanding of the Earth. Geody-namic modelling is bound to increasingly link with a growing number of disciplines, and we trust that the perspective presented here will further facilitate this exchange. In this appendix, we provide an example of how to translate a physical model into a numerical model through discretisation, such that it can be coded up into an algorithm. More in-depth details on numerical modelling can be found in Gerya (2019); Turcotte and Schubert (2012). For this example, we consider a simplified version of the energy conservation equation (6), i.e., the one-dimensional, transient (i.e., time-dependent) heat conduction equation without additional heat sources: where ⇢ is the density, C p the isobaric heat capacity, k the thermal conductivity, and T the temperature (Section 2.1.3). If we assume the thermal conductivity, density, and heat capacity to be constant over the model domain, the equation can be simplified to the heat equation: (16).
where  = k/⇢C p is the thermal diffusivity (Section 3.2). We want to solve this partial differential equation in time and

1505
1D space with appropriate boundary conditions using the finite difference method (Section 3.1) on the domain from x = 0 to x = L x . Note that this appendix is a simple introduction to and example of the finite difference method and the reader is referred to Gerya (2019) for a more thorough presentation of the method.

A1 Taylor series of functions
Before discretising the heat equation, we need to approximate all its terms, after which we can discretise these approximations.

1510
Both first-order and second-order derivatives are present in the heat equation (A2), which we can approximate by Taylor series in the finite difference method. Here, we briefly show how to approximate a general function f (x) that is continuous and differentiable over the range of interest. We assume that we know the value f (x 0 ) and all the derivatives @ n f /@x n at x = x 0 .
Then, in Section A2 and A3 we apply this to the specific terms of the heat equation (A2).

A1.1 First-order derivative 1515
The forward Taylor-series expansion of f (x), away from the point x 0 by a small amount h, is given by This can be rewritten as where the truncation error O(h) indicates that the full solution would require additional terms of order h, h 2 , and so on like in Eq. (A3). Hence, we have an approximation for a first-order derivative. Eq. (A4) is often called the forward derivative, as we can also expand the Taylor series backward (i.e. looking 'left' of x 0 , at location x 0 h) and the backward finite difference derivative then writes: To approximate second-order derivatives, we start with the Taylor expansions of function f (x) at locations x 0 + h and x 0 h: Adding these two equations together yields 1530 which is an approximation of a second-order derivative.

A2 Space discretisation
Now that we know how to approximate first-order and second-order derivatives, we can apply these approximations to Eq. A2.
First, we start with the spatial discretisation for which we need the approximation of the second-order derivative. We want to solve the heat equation on a 1D domain that is divided into separate parts (i.e., discretised). For simplicity, we will focus on 1535 three consecutive, discretely spaced points ( Figure A1). Using Eq. (A8), we can compute the second-order derivative of T at point x i assuming we know the values of T at x i 1 and x i+1 : where where the node spacing, or resolution, h is assumed to be constant.

A3 Time discretisation
The next step is to discretise the first-order time derivative in Eq. A2, using the approximation of the first-order derivative (Section A1.1). To discretise time we dividing it into discrete intervals of time, i.e., the time step t is the time between two consecutive measurements ( Figure A1). The time step is the equivalent of the grid size h in the spatial discretisation (Section A2). In order to distinguish the indices relative to space and time, in what follows we adopt the convention that the forward derivative of T at x = x + i and at time t n = n t as an approximation: This forward finite difference derivative is called first-order accurate, which means that a very small t is required for an accurate solution. The backward derivative from Eq. (A5) is then:

A3.1 Explicit formulation
Both n and i are integers; n varies from 0 to nstep 1, where nstep is the total number of time steps, and i varies from 0 to np 1, where np is the total number of grid points in x-direction. When the forward derivative of Eq. (A10) is used for the time derivative and coupled with the spatial derivative of Eq. (A9), the following approximation of Eq. (A2) is obtained: which can be rearranged to Hence, we have found an expression to compute the temperature T n+1 i at point x i for the next time step n + 1 from all known values at the current time step n. Such a scheme is called an explicit finite difference method which we arrived at through our 1560 choice of evaluating the temporal first-order derivative with forward differences (Section A1.1).

A3.2 Implicit formulation
An alternative approach to deal with the time discretisation is an implicit finite difference scheme, where we use the backward difference for the time derivative (Eq. A11). Together with the spatial derivative of Eq. A9, we then obtain This is often rewritten as follows in order to deal with the unknowns of time step n + 1 instead of the known time step n: The main advantage of implicit methods over their explicit counterpart is that there are no restrictions on the time step, since the fully implicit scheme is unconditionally stable. Therefore we will use the backward (implicit) scheme for the rest of this Eq. (A15) can be rearranged as follows: In contrast to the explicit formulation, we no longer have an explicit relationship which allows us to compute T n+1 i one by one by looping over index i. In other words, Eq. (A16) contains more than one unknown. Therefore, we need to combine these expressions for all unknown points in space and solve the resulting linear system of equations.
A4 Obtaining the linear system of equations 1580 We discretise the domain of length L x with 4 cells, i.e. i = 0, . . . 4 (np = 5). Since we have a second-order derivative in space, we need to prescribe two boundary conditions. We choose T (x = 0) = T 0 = 0 and T (x = L x ) = T 4 = 100. For simplicity, we assume that they do not change over time, so we omit the superscript. Finally, we assume that we know the initial temperature T 0 i for all locations i, i.e. initial conditions have been provided. We want to compute the solution at time t = t, or T 1 i . To simplify notations we define the dimensionless parameter s = ( t)/h 2 such that we get the following 5 equations, i.e., one 1585 for each point of the grid: where A is the coefficient matrix, T is the unknown solution vector, and b is the known right-hand-side vector. As opposed 1595 to the explicit approach, the linear system has a size given by the total number of nodes or points np. A is a sparse matrix their continued support. We thank the attendees of the short courses for their constructive feedback. Thanks to all our proofreaders for their  van Dinther, Y., Gerya, T., Dalguer, L., Corbi, F., Funiciello, F., and Mai, P.: The seismic cycle at subduction thrusts: 2. Dynamic implications of geodynamic simulations validated with laboratory models, J. Geophys. Res., 118, 1502Res., 118, -1525Res., 118, , 2013a van Dinther, Y., Gerya, T., Dalguer, L., Mai, P., Morra, G., and Giardini, D.: The seismic cycle at subduction thrusts: Insights from seismo-Glossary adiabatic heating A heat-producing mechanism that describes how the temperature of the material changes when it is compressed or when it expands (due to changes in pressure), assuming that no energy is exchanged with the surrounding material during this process (i.e., the process is adiabatic). 11, 37 analogue modelling Modelling performed in a laboratory with real-world materials such as sand or silicone acting as substitutes for Earth materials (Schellart and Strak, 2016). 2 analytical solution Exact solution obtained by mathematically solving a system of equations. 4, 28 Anelastic Liquid Approximation (ALA) An approximation to the equations of mass, momentum, and energy, commonly used in compressible mantle convection models. It assumes that lateral density variations are small relative to a reference density profile and can be neglected in the mass and energy conservation equations. 11, 12 Argand number Non-dimensional number (Ar) representing the ratio of the stress arising from crustal thickness contrasts (vertical stress) to the stress required to deform the material at ambient strain rates (horizontal stress) (Houseman and England, 1986;England and McKenzie, 1982). It is commonly used in mountain building dynamics as a measure of the tendency of an orogen to collapse under its own gravitational potential energy. 44 diffusion creep A mechanism of viscous deformation of rocks in the Earth's interior, facilitated by the motion of imperfections in the crystal lattice -single atoms or vacancies (i.e., places in the crystal lattice where an atom is missing) -through mineral grains, which accumulate to large-scale deformation over time. Diffusion creep is assumed to be the dominant deformation mechanism in the lower mantle. 7, 14 direct method Class of methods used in geodynamics to solve the linear system of equations in a finite number of operations obtained after discretising the partial differential equations. Typically robust, these methods give an exact solution barring numerical rounding errors, but are often memory intensive Eijkhout, 2013). 22

article-in-cell
Dirichlet boundary condition A type of boundary condition that specifies the value of the solution of an equation at the boundary, such as prescribing a zero velocity (no slip) boundary condition. 7, 9, 35 discrete element method (DEM) Family of numerical methods for computing the motion of a large number of small particles which are in contact with each other (Matuttis and Chen, 2014). 19 dislocation creep A mechanism of viscous deformation of rocks in the Earth's interior, facilitated by the motion of imperfections in the crystal lattice -lines of atoms or line defects -through mineral grains, which accumulate to large-scale deformation over time. This migration is highly sensitive to the stress applied to a rock. Dislocation creep is assumed to be an important deformation mechanism in the upper mantle. As the creep behaviour can be described by a power law dependence of the strain rate on the differential stress (e.g., Hirth and Kohlstedt, 2003;Karato, 2008) extended Boussinesq approximation (EBA) An approximation to the equations of mass, momentum, and energy, commonly used in incompressible models of convection or lithosphere dynamics. It is based on the same assumptions as the Boussinesq approximation, but additionally considers adiabatic and shear heating. As such, the reference temperature varies with depth. 12 FAIR Data principles that stand for and promote Findability, Accessibility, Interoperability, and Reuse of data. 54 finite difference method (FDM) Widely used numerical method that solves ordinary differential and partial differential equations by approximating the derivatives with finite differences in both space and time. 19 finite element method (FEM) Widely used numerical method for solving problems in geodynamics, engineering, and mathematics. The FEM subdivides a large system into smaller parts called elements and formulates the ordinary differential or partial differential equations on these, which ultimately results in a linear system of equations. 19 finite volume method Numerical method to solve partial differential equations in which volume integrals in a partial differential equation which contain a divergence term are converted to surface integrals using the divergence theorem. These terms are then evaluated as fluxes at the surfaces of each finite volume. 19 fixed-point ('Picard') iteration An iterative method for dealing with nonlinearities in the governing equations. The method starts by guessing an approximate solution. It then improves the solution in each nonlinear iteration by using the current approximate solution to compute the solution-dependent coefficients and solving the linear system, obtaining the approximate solution for the next iteration. The process is repeated until the left-hand side and the right-hand side of the half-space cooling model Analytical expression for the temperature as a function of time and depth in the case that a semiinfinite constant temperature lithosphere is suddenly subjected to a different temperature at the surface (Jaupart and Mareschal, 2011). 37 hand of god A modelling approach where a specific feature is imposed artificially by the modeller, rather than developing selfconsistently in the model (e.g., Glišović et al., 2012). Examples are prescribing plate velocities at the surface boundary, or weak zones at plate boundaries. 38 hello world A computer program that outputs or displays the message 'Hello, World!'. It is a simple sanity test to illustrate the basic syntax of a programming language and check that it is correctly installed and used. A parallel 'Hello world' test outputs a message from each processor. 28 hero code A code predominantly written and maintained by a single person. 54 IMRAD Acronym for Introduction, Methods, Results, And Discussion/conclusion, which stands for an organisational structure of scientific articles adopted by the majority of scientific journals. 47 initial condition Imposed physical starting conditions required to solve time-dependent equations like the temperature and other advection equations. They specify the initial values of the solution parameter, e.g., temperature, throughout the model domain. 37 initial perturbation Perturbation of a field imposed by the user to trigger a physical instability (typically for buoyancy-driven models). This is often applied to the initial temperature field in the form of a harmonic perturbation or noise with a user-defined amplitude. 38 inner core dynamics The motion within the Earth's solid iron-nickel inner core and the forces causing this motion. 2 internal heating A heat-producing mechanism that is intrinsic to a specific material, and not related to its motion. In geodynamics, this would be the decay of radiogenic elements, i.e., radiogenic heat production. 11 iterative method A class of mathematical methods that solves a problem by starting from an initial guess, and then improving the solution iteratively, i.e. by repeatedly applying a given solution method to the problem with each new approximate solution being based on the solution of the previous step. In numerical geodynamic modelling, iterative methods are often used to solve the large linear system that arises from discretising the equations. 4, 7, 22 Knudsen number Non-dimensional number (Kn) representing the lower limit in length scale where the Navier-Stokes equations are relevant and the transitions to other forms of transport in micro-to nanoscale flow channels. It is defined as the ratio of the mean-free-path of molecules (i.e., the average distance travelled by a moving particle between successive impacts) to the macroscopic length scale of the pore space. 44 level-set method Numerical technique that allows for tracking interfaces. The basic idea is that the location of the interface is defined as the zero level set of a field defined over the domain of interest. 26 linear stability analysis A method that can be used to analyse the behaviour of nonlinear systems, specifically, how they evolve from a given initial condition when small perturbations are added to it. The analysis involves linearising the equations (sometimes around an equilibrium), so that they can be solved analytically. 7 Mach number Non-dimensional number (M) representing the ratio of the velocity of material flow to the speed of sound in a  Neumann boundary condition A type of boundary condition that specifies the value of the derivative of the solution. Specifying the stress is a Neumann condition, as stress is related to the velocity through its derivatives. 9, 35 Newton iteration An iterative method for dealing with nonlinearities in the governing equations (e.g., Fraters et al., 2019a).
It introduces the residual r = A · X b of the linear system, which needs to be minimised. This is done by finding the roots of its derivative and requires computing the Jacobian J of the residual (a matrix containing all its first-order partial derivatives with respect to the solution variables), and solving J X = r for X. Note that this implies that the linear system is being solved for the update X to the solution of the original problem in each nonlinear iteration (and not the solution itself). 26 Newtonian fluid A viscous fluid with a linear relation between the stress and the rate of deformation (strain rate). In geodynamic modelling, the rheology of mantle rocks is sometimes approximated as that of a Newtonian fluid, because there is a linear relation between stress and strain rate if rocks deform by diffusion creep. 13 no slip A type of Dirichlet boundary condition for the Stokes equations in which the velocity is fixed to zero at the boundary. It is a special case of prescribed velocities. This boundary condition is often used to mimic the 660 km depth discontinuity at the bottom boundary of asthenospheric-scale models. 35 nonlinear In a nonlinear system, the reaction or output of the system is not proportional to the input, i.e., the relation between input and output is not linear. In a system of (partial differential) equations, as considered in geodynamic modelling, this means that any of the unknown variables (in our case, velocity, pressure, temperature, and in some cases other advected quantities) appear in the equations in a form that is not of polynomial degree one (e.g., a higher polynomial degree, or a different relationship, like an exponent). This often occurs when material properties depend on the unknown variables.

18, 19, 23
numerical code Software that in the context of geodynamics solves the governing equations (usually conservation of momentum, mass, and energy) based on a numerical model. 4, 6 numerical implementation The process of putting the physical/mathematical model into computer language as a series of numerical routines. 16, 28 numerical model The representation of a mathematical or physical model on a computer. 4, 6, 7, 19 numerical modelling Simplified representation of a natural system based on physical principles that is performed on a computer. This involves using numerical methods to find solutions to, or analyse, the equations that make up a mathematical or physical model. 3 Nusselt number Non-dimensional number (Nu) representing the relative efficiency of convective heat transport. It describes the ratio of the total heat flux in the presence of convection (i.e., heat advection and conduction combined) to the heat flux in the absence of convection (i.e., heat conduction only). Earth's mantle is characterised by Nu ⇡ 10 100, meaning that the convection of heat is much more efficient than the conduction of heat. 44 open boundary In the geodynamics context, a type of Neumann boundary condition for the Stokes equations in which the lithostatic pressure used to prescribe the normal component of the traction on the boundary, allowing for material entering and leaving the domain based on the over-pressure. 36 OpenMP (Open Multi-Processing) An application programming interface that consists of a set of compiler directives, library routines, and environment variables that support the parallelisation of codes running on shared memory systems. 22 parameter space The parameters and their range of values that are tested in a -in this case geodynamic modelling -study. 8,

34
parameter study A geodynamic modelling study focused on determining the dependency of the model solution and evolution on the parameters under investigation. Typically, the value of one or more parameters is varied systematically (see parameter space) and the resulting model behaviour is analysed for corresponding trends. 9, 34 Peclet number Non-dimensional number (Pe) representing the ratio of heat diffusion to heat advection time scales in a fluid transporting mass and heat. 44 Peierls creep A mechanism of viscous deformation of rocks in the Earth's interior with a very strong, exponential dependence on the stress applied to a rock (also called 'low-temperature plasticity'). It is inferred to apply when creep is controlled by glide of dislocations (line defects) through mineral grains. 14 periodic boundary condition Type of boundary condition that does not explicitly prescribe any part of the solution. Instead, it 'links' boundaries together to approximate a larger (or infinite) system of which the model setup is merely a part: any materials or flows passing through one boundary interface re-enter the model domain through the opposite boundary interface. 36 physical model A conceptual model of a natural process or system, based on the laws of physics. These models generally use equations to approximate the natural behaviour. 4, 6-8, 19 physical or mathematical modelling The process of approximating a natural system or process by developing a mathematical or physical model. In geodynamics, the mathematical descriptions follow the laws of physics, and physical and mathematical modelling can be considered to mean the same. 2 plane strain assumption Assumption that the strain in the direction perpendicular to the applied stresses equals zero. It is often used in 2D Cartesian models of the lithosphere. 16 plasticity Non-elastic and non-recoverable deformation of materials characterised by a yield criterion. Once local stresses reach the yield strength of a material, it starts to deform plastically. In geodynamic modelling, plasticity is commonly used to describe brittle deformation, specifically, tectonic fracture/failure in shear. There are a number of failure criteria used to model brittle behaviour, in particular, Byerlee's law (which is empirically-determined, Burov, 2011) and the pressure-dependent Mohr-Coulomb failure criterion. 2, 10, 15 plate cooling model Analytical expression for the temperature as a function of time and depth in the case that the lithosphere is of finite thickness and subjected to a constant temperature at the bottom and the top (McKenzie, 1967). 37 Prandtl number Non-dimensional number (Pr) representing the ratio of momentum diffusivity (⌘/⇢) to thermal diffusivity in a fluid. It depends on the properties of the fluid, namely its viscosity, specific heat, and thermal conductivity. The Prandtl number is approximately 10 25 for Earth's mantle. Consequently, a common approximation for mantle and lithosphere models is the 'infinite Prandtl number' approximation, which eliminates inertia effects from the governing equations.

11, 44
prescribed stresses A type of Neumann boundary condition for the Stokes equations in which the stress at the boundary is prescribed, for example to mimic ridge push on a plate moving into the domain. 35 prescribed velocities A type of Dirichlet boundary condition for the Stokes equations in which the velocity is prescribed, for example to define the plate motions along the top the model. Also called kinematic or in-/outflow boundary conditions.

7, 35
processor Electronic circuitry in a computer that executes instructions given by a computer program. Also called central processing unit (CPU). 5, 10-12, 22, 25 quadtree A data structure which partitions a 2D model domain by recursively subdividing cells into four quadrants or regions ('children'). How often each cell is subdivided can vary cell by cell, so that relevant areas are subdivided more often (leading to smaller cells). An octree is the 3D analogue where each cell is subdivided into 8 children cells. 23 radial basis function method (RBF method) Numerical method based on real-valued functions whose values depend only on the distance between the input and some fixed point. 20 Rayleigh number Non-dimensional number (Ra) representing the ratio of heat conduction to heat convection time scales. For a system with given boundary conditions and geometry, a value for a critical Rayleigh number, Ra critical , exists that marks the transition from pure conduction to fluid convection as the dominant mode of heat transfer. 2, 44 regime diagram Commonly the synthesis of analysing the results of a parameter study that shows how the individual models' behaviour can be grouped into several regimes (e.g. slab retreat, slab advance and slab stagnation in a subduction model), with the studied parameters on the axes of the diagram. 34 remeshing Manipulating the mesh in numerical methods to improve its quality by minimising skewness and/or element aspect ratios. 41 replicable Results of a prior study can be duplicated with new data and tools (e.g., codes) (Bollen et al., 2015;Goodman et al., 2016). Not to be confused with reproducible. 10, 46, 53 reproducible Results of a prior study can be duplicated with the same data and tools (e.g., codes) (Bollen et al., 2015;Goodman et al., 2016). Not to be confused with replicable. 9, 46, 53 resolution test Investigation of the changes in the model outcome with varying resolution. Ideally, the model outcome (i.e., solution) converges with the differences in solution decreasing with increasing resolution (i.e., smaller grid size). 39 Reynolds number Non-dimensional number (Re) representing the ratio of inertial forces to viscous forces in a flow. (Very) Low Reynolds numbers are typical of fluid flow in the mantle. 44 rheology The study of the deformation of matter in response to an applied force (see Ranalli, 1995;Jaeger et al., 2007). 13 Robin boundary condition A type of boundary condition that consists of linear combinations of the values and of the derivatives of the solution. 35 scaling analysis A technique to simplify equations by using non-dimensional numbers that combine multiple dependent variables and characterise the physical system. The goal of scaling analysis is to obtain dynamic similarity in which nondimensional numbers can predict the behaviour of the system (i.e., systems characterised by similar non-dimensional numbers behave similarly). Also known as dimensional analysis (Barenblatt, 1996). 7, 29 scientific colour maps Scientifically derived colour gradients that are perceptually uniform (i.e., do not distort the data visually), perceptually ordered (i.e., are intuitively readable), and colour-vision deficiency friendly (i.e., universally readable) to map data values to colours (Crameri et al., 2020). 53 shared memory systems A computer architecture with a (large) block of memory that can be accessed by several processors.
This allows for fast communication between the different threads of a parallel program. 8, 22 shear band Narrow zones of extreme shear deformation, resulting from strain localisation. In numerical geodynamic modelling, where discrete faults cannot be represented, and plasticity is used to approximate brittle failure, shear bands are often interpreted as faults on crustal and lithospheric scales, resulting from coalescence of preexisting micro-cracks and fractures into a single frictional shear band (fault) at sufficiently high strains (Burov, 2011). 15, 40 smoothed particle hydrodynamics method (SPH method) A mesh-free Lagrangian computational method used in mechanics of continuum media, such as solid mechanics and fluid flows, astrophysics, and ballistics. The method is based on dividing the material/fluid into a set of particles with mass that interact with each other. 19 software management plan A well-defined plan that describes the goals in software development and management: what software will be produced and for whom, and how will this software be managed and maintained. 53 sparse matrix A matrix in which most of the entries are zero, which is common for matrices obtained from discretising partial differential equations. Which numerical methods are most efficient for solving a linear system depends on the properties of a matrix, including its sparsity. 21, 61 specific modelling Modelling approach that aims to reproduce and understand a specific state of a certain geodynamic system (e.g., based on a direct set of observations, or with the aim of modelling a specific geographic area). 32, 44 spectral methods Numerical methods to solve ordinary or partial differential equations, generally involving the use of the fast Fourier transform. The solution is expressed as a sum of certain 'basis functions' (for example, as a Fourier series), and the coefficients in the sum are computed such that they satisfy the differential equation as well as possible. 19 stencil a geometric arrangement of neighbouring nodes that relate to the node of interest on which finite fifference derivatives are formulated. 20 test Steps or programs intended to establish the quality, performance, or reliability of a code, especially before it is taken into widespread use. 27 thermal diffusivity Thermal property of a material representing the rate of transfer of heat of a material from hot to cold areas, defined by the thermal conductivity divided by density and specific heat capacity at constant pressure. 11, 57 thermodynamic potentials Scalar quantities used to describe the thermodynamic state of a system, e.g., the internal energy, Helmholtz free energy, Gibbs free energy, enthalpy, and the Landau potential. 41 Truncated Anelastic Liquid Approximation (TALA) An approximation to the equations of mass, momentum, and energy, used in compressible mantle convection models. It is based on the same assumptions as the Anelastic Liquid Approximation, but further assumes that changes in density due to pressure variations can be neglected in the buoyancy term in the momentum equation as well. 12 unit test Generally automated tests to ensure that a section of an application (known as the 'unit') meets its design and behaves as intended. The unit can be an individual function or procedure, but can also extend to the whole module or object. 28 validation Testing of the equations against observations to confirm that the right equations are solved. Not to be confused with verification. 12, 27 verification Testing of whether the equations are implemented and solved correctly. Not to be confused with validation. 6, 12, 27 version control software System responsible for managing and tracking changes to codes, enabling easier and less errorprone collaboration on codes. 4, 11, 54 viscous fluid A material that deforms as a continuum, without the formation of gaps (e.g., bubbles), with a predominantly viscous rheology. For a viscous rheology, the stresses are proportional to rate of deformation (strain rate) in the material, and deformation is irreversible. In geodynamic models of the mantle and lithosphere, it is often assumed that rocks deform like viscous fluids (although other deformation mechanisms also play an important role), and Stokes Equations can be used to describe their mass and momentum balance. 6-8, 11 viscous isentropic relaxation time scale The time necessary to expand (or contract) regions of excess (or deficit) pressure against viscous constraints (e.g., the characteristic time scale of postglacial rebound). In the Earth, this timescale is on the order of a few hundred years for the upper mantle to a few tens of thousands of years for the lower mantle. 9 visualisation pitfalls Inaccurate representation of data, thus leading to miscommunication of data, by altering the impression of actual local data magnitudes and gradients, or making the data graphically inaccessible to parts of the readership. 43 volume-of-fluid method A Eulerian numerical technique for tracking and locating the free surface or fluid-fluid interface. 26 weak scaling A method for quantifying the parallel performance of a code. It describes how the solution time varies when both the number of processors and the problem size (the number of unknowns) are increased by the same factor. The ideal result for this scenario is a constant solution time. 23 weak seed Area/volume in the domain prescribed by the user to have weak mechanical properties. It ensures that the deformation initially localises on the weak seed rather than on numerical noise and/or near boundaries. 38