Advances in geodynamo modelling

ABSTRACT This paper reviews the remarkable developments in numerical geodynamo simulations over the last few years. Simulations with Ekman numbers as low as are now within reach and more and more details of the observed field are recovered by computer models. However, some newer experimental and ab initio results suggest a rather large thermal conductivity for the liquid iron alloy in Earth's core. More heat would then simply be conducted down the core adiabat and would not be available for driving the dynamo process. The current status of this topic is reported and alternative driving scenarios are discussed. The paper then addresses the question whether dynamo simulations obey the magnetostrophic force balance that characterises the geodynamo and proceeds with discussing related problems like scaling laws and torsional oscillations. Finally, recent developments in geomagnetic data assimilation are reviewed, where geomagnetic data and dynamo simulations are coupled to form a tool for interpreting observations and predicting the future evolution of Earth's magnetic field.


Introduction
Numerical geodynamo models have come a long way since the first attempts were published about 25 years ago (Glatzmaier andRoberts 1995, Kageyama andSato 1995). The main progress was driven by the always growing computer power, which allowed accessing increasingly realistic parameter regimes. The incorporation of new effects, like the heat flux pattern imposed by Earth's mantle, also contributed to the success story.
Detailed accounts of this development can be found in recent reviews by Roberts and King (2013) or Christensen and Wicht (2015), while the mathematical fundamentals are brilliantly covered by Braginsky and Roberts (1995). This paper provides a rather selective overview of some interesting recent developments. New simulations at very low Ekman numbers down to E = 10 −8 motivate a fresh look at the question whether the numerical dynamos operate in the "correct regime" or only yield seemingly "correct" results for the wrong reason. After all, the Ekman numbers and magnetic Prandtl numbers, even in these most recent and extremely costly simulations, are still many orders of magnitude away from realistic values.
The expression "correct regime" mostly refers to the correct force balance in the Navier-Stokes equation, where Coriolis force, pressure gradient, buoyancy, and Lorentz force contribute in first order, while viscous and inertial forces are much smaller. The too large Ekman numbers in the simulations translate into increased viscous forces. This in turn has to be counteracted by stronger convective driving, which can result in excessively large inertial effects. Intimately linked to the force balance is the question of the correct scaling for extrapolating the simulation results to realistic parameters. We provide updated discussions on these issues in section 4 and section 5.
Before doing so, however, we briefly introduce the mathematical formulation of the dynamo problem in section 2 and discuss the fundamental problem of driving the geodynamo process in section 3. Newer ab initio simulations of the thermal conductivity (de Koker andSteinle-Neumann 2012, Pozzo et al. 2012) in Earth's core suggest much larger values than previously assumed. This would mean that a lot of the heat emanating from Earth's core is simply conducted down the adiabat and therefore not available for driving the dynamo. At the moment, it remains unclear how the geodynamo, which operated during most of Earth's history, can be consolidated with our planet's thermal evolution. Section 3 provides an overview of the recent contributions to this interesting discussion.
A still relatively new development is the combination of dynamo simulations and geomagnetic observations in a data assimilation approach. This has several potential benefits and applications: First, there is the hope that the new framework would lead to dynamo simulations that better reflect the state of Earth's core. Second, the additional dynamical information brought in by the simulations could help to interpret and complement geomagnetic data. And third, similar to the wide-spread application in weather prediction, geomagnetic data assimilation allows predicting the future evolution of Earth's dynamo field. We provide an overview of the developments in this topic in section 6 before the paper closes with conclusions in section 7.

Mathematical formulation
The magneto-hydrodynamical equations for convection and magnetic field generation are formulated in a frame of reference rotating with rate about the z-axis. The rotation typically refers to the outer boundary and thus to Earth's mantle in the geodynamo context. Solved are equations for small disturbances around a hydrostatic, adiabatic, and well mixed background state (Braginsky and Roberts 1995). For planetary iron cores a homogeneous background is considered in the so-called Boussineq approximation. Background temperature and density, but also physical properties like thermal or electrical conductivities, are assumed to be homogeneous throughout the outer core. We will not discuss numerical methods here and refer to Glatzmaier (1984), Clune et al. (1999), Christensen and Wicht (2015), or Matsui et al. (2016).

Fundamental equations
The five defining equations are the Navier-Stokes equation

∂U ∂t
+ U·∇U + 2ẑ × U = −∇p + Ra (r/r o )Cr the induction equation the codensity equation the simplified continuity equation and the magnetic continuity equation Pressure p in equation (1) combines the non-hydrostatic pressure and centrifugal effects. C is the so-called codensity, which is further discussed below, andr andẑ are the unit vectors in radial direction and along the rotation axis, respectively. These equations have been non-dimensionalized by using the shell thickness D = r o − r i as a length scale, with r i the inner and r o the outer boundary radius, and Ω −1 as a time scale. Moreover, the magnetic field has been expressed in multiples of (ρμ) 1/2 ΩD, whereρ is the background density and μ the magnetic permeability. This non-dimensionalization avoids using diffusivities in the scaling, anticipating that quantities like the magnetic field strength B or the flow amplitude U should ultimately become independent of diffusive effects (Christensen and Aubert 2006).
Convection is driven by the gravity force acting on the density differences α C. Codensity C combines density differences due to deviations from the adiabatic background temperature and the homogeneous background composition, and the expansivity α = (1/ρ)∂ρ/∂C represents a combination of thermal and compositional effects. The latter comes into play because the light elements (sulfur, oxygen, carbon) mixed into the liquid core are not readily incorporated into the inner core matrix. In particular oxygen is expelled from the growing inner core front. Another example for compositional convection is the slow exsolution of MgO discussed in section 3.
Since the diffusivities for temperature and composition are vastly different, describing the evolution of both with only one evolution equation seems a gross simplification. However, since the simulations cannot resolve the small scale turbulence in Earth's core, the numerical diffusivities represent increased values that account for the unresolved turbulent mixing. And since the mixing acts similarly on all quantities, the turbulent diffusivities should also be comparable.
The variable ε on the right-hand-side of the codensity equation (3) models volumetric sinks or sources of different origins, for example radiogenic heating or the slow secular cooling.
The non-dimensional control parameters are the Ekman number the modified Rayleigh number the Prandtl number the magnetic Prandtl number and the aspect ratio Here g o is the outer boundary reference gravity. Parameter values estimated for Earth and used in dynamo simulations are listed in table 2 in section 4. Ekman number (6) measures the ratio of viscous to Coriolis forces. The modified Rayleigh number (7) can be related to the classical parameter defined for Rayleigh-Bénard convection in a non-rotating plane layer with an imposed codensity difference C instead of a temperature difference T: It measures the ratio of buoyancy effects to viscous and thermal diffusion effects, while the modified version (7) takes into account that rotation tends to suppress convection in a fast rotating system. Prandtl number and magnetic Prandtl number are ratios of the kinematic diffusivity ν, the thermal diffusivity κ, and the magnetic diffusivity λ = 1/(μσ ), where σ is the electrical conductivity. A more in depth discussion of the equations and involved approximations can be found in Braginsky and Roberts (1995).

Boundary conditions
In Earth's core, rigid flow boundary conditions with U = 0 at both boundaries seem appropriate. As we will further discuss below, however, the Ekman number is many orders of magnitude too high in the simulations. The viscous Ekman boundary layers, whose thickness scales like E 1/2 , may therefore remain too influential. Some authors therefore prefer to use stress-free conditions (for example Kuang and Bloxham 1997, Grote and Busse 2000, Simitev and Busse 2005, Aubert 2013).
Either constant temperature or constant heat flux boundary conditions can be assumed. Since the sluggish dynamics of Earth's mantle controls the heat allowed to escape the core, heat flux boundary conditions are thus more appropriate at the outer boundary. Sometimes patterns deduced from lower mantle seismic tomography are imposed. Lower than average seismic velocities are typically interpreted as higher than average temperatures, which in turn imply a reduced temperature gradient and thus a smaller heat flux from the core to the mantle , Aubert et al. 2007). However, since the variations in seismic velocity may also have a compositional origin, different patterns seem possible (Amit and Choblet 2009).
Because the light elements released from the growing inner core are not effectively accommodated by the mantle, using a vanishing compositional flux is most realistic. A boundary condition for codensity should reflect a mixture of both convection forms, representing the respective contributions to the overall convective driving.
Both fixed flux or fixed codensity conditions are commonly used at the inner core boundary. However, since the inner core also evolves over time and its dynamics is coupled to that of the outer core, it seems more appropriate to think in terms of a matching rather than a boundary condition. Two different end-member approaches interpret the inner core as the driver or the slave of the outer core dynamics. Models of inner core convection show that under certain conditions a specific mode is preferred, where the inner core grows on one side but melts on the other (Alboussière et al. 2010, Monnereau et al. 2010. The lower boundary conditions for the outer core dynamics should reflect this in the form of a differential codensity flux often modelled with a spherical harmonic mode of degree and order one. However, the outer core dynamics could also determine where the inner core is cooled more effectively and thus grows more rapidly. This can be cast into a dynamical equation for the inner core boundary that ties the local cooling rate to the codensity flux (Braginsky andRoberts 1995, Glatzmaier andRoberts 1996).
When using flux boundary conditions rather than imposing a codensity difference across the core, the Rayleigh number needs to be modified. Imposing, for example, a total outer boundary heat flux Q o , yields the condensity scale Q o /(4π r 2 oρ c p DΩ) and the modified Rayleigh number where c p is the specific heat capacity at constant pressure. Because of the low electrical conductivity in the mantle, the outer core field becomes a potential field at r o and thus has to obey a boundary condition that can easily be formulated in terms of spherical harmonics (Christensen and Wicht 2015). At the inner boundary, the outer core field has to match an inner core field that obeys a simplified induction equation (2) with U = sω iφ , where ω i is the inner core rotation rate, s the distance to the rotation axis andφ the azimuthal unit vector. For simplicity, some authors assume an electrically insulating inner core, so that the field also becomes a potential field at r i . The differences between solutions with conducting or insulating inner core seem minor (Wicht 2005) as long as the aspect ratio is not too large, say a ≤ 0.5.

Diagnostic parameters
While the dimensionless input parameters discussed above define the model, the results are often quantified by non-dimensional diagnostic parameters. The parameters can be understood in different ways, for example as ratios of time scales, as ratios of forces, or more generally as ratios of terms appearing in the dynamical equations. The Ekman number can be interpreted as the ratio of the rotation period Ω −1 to the viscous diffusion time D 2 /ν. Prandtl and magnetic Prandtl numbers can be interpreted as ratios of respective (inverse) diffusion times. Important diagnostic parameters are the Rossby number the Reynolds number the magnetic Reynolds number the Elsasser number and the Lehnert number Here U and B stand for characteristic flow velocity and magnetic field strength, typically the respective rms values in the core. These parameters provide different dimensionless expressions for either magnetic field strength (Λ, Le) or flow amplitude (Ro, Re, Rm). Dynamically more interesting, however, is the fact that they can serve as proxies for the ratio of certain terms in the Navier-Stokes equation (1) or the induction equation (2).
The Rossby number estimates the ratio of the non-linear inertial term U·∇U to the Coriolis term 2Ω × U and quantifies the rotation influence on the dynamics. The Reynolds number, on the other hand, provides a measure for the ratio of the non-linear inertial term to the viscous force ν∇ 2 U. Ro and Re are obviously related via Ro = ReE.
The magnetic Reynolds number measures the ratio of induction and Ohmic dissipation in the induction equation. The Elsasser number is often used to estimate the relative importance of the Lorentz force J × B to the Coriolis force. However, Soderlund et al. (2015) show that the dynamic Elsasser number provides a more direct and better measure for the true force ratio when faster variations in B are considered. The Lehnert number, called Lorentz number by Christensen and Aubert (2006), is the ratio of the rotation time τ Ω = 1/Ω to the magnetic Alfvén time when assuming a length scale of order one. Here the Alfvén velocity is the typical propagation velocity of magnetic waves along field lines. The ratio of the Alfvén time scale to the flow (turnover) time scale τ U = D/U is the Alfvén Mach number The Nusselt number is the ratio of the total heat flux Q to the (theoretical) purely conductive heat flux Q c = 4π r i r oρ c p C/D. Christensen and Aubert (2006) define a flux Rayleigh number that is based on the convective heat flux Q U : Note the similarity to definition (12). The flux Rayleigh number is closely related to the power P generated by buoyancy forces, which can be derived by multiplying the Navier-Stokes equation with U and integrating over the shell: where U r and C represent dimensional values. In non-dimensional form, given in units of ρΩ 3 D 5 , this reads The integral can be expressed in terms of Q U and Christensen and Aubert (2006) show that the power is roughly given by for the Earth's aspect ratio. The exact relation between power and Rayleigh number actually depends on the thermal boundary conditions (Christensen and Aubert 2006, Aubert et al. 2009, Christensen 2010, Aubert 2018, and we only provide the expression for an imposed temperature (codensity) jump here.

Driving the geodynamo
The convective flows in the Earth's core are powered by the slow secular cooling of the planet and the related release of light elements from a growing inner core. Paleomagnetic data suggest that the geodynamo has already operated 3.5 Gyr ago (Tarduno et al. 2010) and newer findings seem to extend the age of the dynamo to 4.2 Gyr (Tarduno et al. 2015). The early dynamo must have been powered by super-adiabatic temperature gradients in a hotter planet. Compositional convection only became available once the core temperature profile crossed the solidus temperature of the core alloy, which is not a pure iron/nickel mixture but also contains light elements like silicon, sulfur or oxygen. In particular oxygen is not readily incorporated into the solid iron matrix and is thus released at the growing inner core front (Alfè et al. 2007. Consolidating such a scenario with Earth's thermal history became questionable with the increased thermal conductivity estimates published by de Koker and Steinle-Neumann (2012) and Pozzo et al. (2012). These larger values meant that more heat is simply conducted along the adiabat and is thus not available for driving the dynamo. We will discuss the related re-evaluation of dynamo driving scenarios in section 3.2 after reviewing new results on the thermal and electrical conductivities in Earth's core in section 3.1.  Xu et al. (2018) All listed experiments are Diamond Anvil Cell (DAC) setups while older experiments used shock waves. Gomi et al. (2013) and the older experiments measured the electrical conductivity, highlighted in bold, and use the Wiedemann-Franz Law to calculate the thermal conductivity. Konôpková et al. (2016), however, infres the thermal conductivity from laser flash heating (McWilliams et al. 2015). Saha and Mukherjee (2016) use a similar method and confirm the low thermal conductivity. All experimental data are extrapolated to account for the effects of impurities by the light core elements and, if required, are also extrapolated in temperature and pressure.

Thermal and electrical conductivities in Earth's core
The new larger conductivity values by de Koker and Steinle-Neumann (2012) and Pozzo et al. (2012) are the result of quantum mechanical ab initio calculations that simulate the interactions between a sample of representative molecules and atoms. The interaction between the numerous electrons is typically reduced to a one particle problem by using an effective potential that depends on the electron density. Different approximations for the effective potential are in use in this so-called density-functional theory approach. With a Core-Mantle Boundary (CMB) value of k = 100 W m −1 K −1 , these new thermal conductivity values are about three times higher than classical ones. Table 1 compares the older low estimate of Stacey and Loper (2007) with the newer ab initio calculations and experimental data. Listed are not only thermal but also electrical conductivities. Since both electric currents and heat flux are mainly carried by electrons, the two conductivities are related. The Wiedemann-Franz law (WFL) states that the ratio of thermal to electrical conductivity grows linearly with temperature with a proportionality factor L called the Lorentz number: Recent diamond anvil cell experiments yield ambiguous results. Gomi et al. (2013) and Ohta et al. (2016) measure a high electrical conductivity and, using the WFL, predict a high thermal conductivity similar to the values suggested by the ab initio calculations. Konôpková et al. (2016), on the other hand, designed a diamond-anvil cell experiment where they deduce the thermal conductivity by flash heating one side of the probe and monitoring the temperature changes on the other (McWilliams et al. 2015). They find a particularly low thermal conductivity of only 27 ± 7 W m −1 K −1 that would pose no energetic problem for the geodynamo. Experiments by Saha and Mukherjee (2016), who measured the thermal conductivity with a similar method, confirm these low values.
The experimental values could in theory be brought in line when allowing for a strong pressure and temperature dependence of the Lorentz number. However, the ab initio simulations suggest that the dependence remains rather small (de Koker and Steinle-Neumann 2012, Secco 2017). Table 2. Input and output parameters for standard simulations, high-end simulations and Earth's core conditions. The Pr = 1 value cited for Earth refers to thermal diffusion. Composition diffusion would yield much larger values of up to Pr = 1000 (Braginsky and Roberts 1995).

Energy considerations and dynamo history
Because of the huge viscosity differences, the convective heat transport through the mantle is much less efficient than through the core. The lower mantle therefore forms a bottleneck that controls the heat flow, Q o , through the core-mantle boundary. In addition to the internal heat Q S lost by secular cooling, Q o has two contributions associated to the growth of Earth's solid inner core: the latent heat Q L and a compositional component Q C which accounts for the gravitational energy released by core differentiation.
Adding up the different sources yields where we have written the secular cooling as a volumetric term with source density expressed in terms of the changes in outer boundary background temperatureT o . Variables with a tilde generally express the adiabatic and well mixed background state. We have neglected the following minor contributions: (1) the change in gravitational energy due to core shrinking, (2) the heat of reaction that goes along with the change in outer core composition, and (3) radiogenic heating. Typically, a maximum of 300 ppm of potassium is considered realistic for the core which has only a limited impact on the energetics (Nimmo 2015).
Like the secular cooling, latent heat and light element release at the inner core boundary can be written in terms of a CMB temperature drop ∂T o /∂t with proportionality factors K S , K L , and K C that depend only on material properties (e.g. see Davies 2015): The inner core growth rate ∂r i /∂t is also directly related to ∂T o /∂t with, in first order, a proportionality factor K that depends on the gradients in core adiabat and melting curve (Braginsky and Roberts 1995): Given a time dependent CMB heat flux Q o (t) we could integrate equation (30) back in time to determine when the inner core would have started to nucleate . Neither kinetic or magnetic energy nor viscous or Ohmic heating appear in the global energy balance (30) because we consider a global slow evolution of the background state. On these time scales, convection and the dynamo action merely "serve" to transport heat to the outer boundary and to homogenise the composition.
Flow and magnetic field are brought into play when considering the entropy balance, where it matters at which temperature energy is infused or extracted from the system: The two last terms on the right hand side are dissipative contributions. The source q a = k(∇T) 2 /T describes diffusion down the adiabat, while ϕ is the total dissipative entropy production. Because of the small magnetic Prandtl number, Joule dissipation dominates so that This becomes somewhat more accessible when using an effective dissipation temperaturẽ T J so that where Q J is the total Ohmic or Joule heating for which at least some estimates are available. Since the adiabatic temperature varies by only about 35% throughout the core, the error introduced by choosing a meanT J for the entropy equation remains limited (Labrosse 2015). Energy equation (30) and entropy equation (32) formulate the energetic dynamo constraints with basically two unknowns: the heat flux Q o through the outer boundary and the dissipative entropy production ϕ. Two approaches are thus possible: For a given outer boundary flux Q o one could deduce the cooling rate and inner core growth rate via equation (30) and then calculate ϕ or Q J via equation (32) and equation (34). Alternatively, assuming a minimum ϕ allows deducing the cooling rate, inner core growth rate, and CMB heat flux. Typically, ϕ = 0 is assumed to determine the smallest required CMB flux Q min o for an operating dynamo.
Unfortunately, neither Q o nor the dissipation are well constrained for Earth. Christensen and Tilgner (2004) estimate that today's Joule heating Q J is rather small, somewhere between 0.2 and 0.5 TW. Estimates for today's CMB heat flux, on the other hand, range between 7 TW and 17 TW (Nimmo 2015). For the new higher thermal conductivity values by de Koker and Steinle-Neumann (2012) or Pozzo et al. (2012), about Q a = 15.5 TW are conducted down the adiabat, while the minimum heat flux is about Q min o ≈ 5-8 TW , Labrosse 2015, Nimmo 2015. From purely energetic considerations, Q o could thus be sub-adiabatic because the additional power from the compositional convection would suffice to drive the dynamo.
A sub-adiabatic CMB heat flux, Q o < Q a , translates into a temperature gradient that is stably stratifying and thus tends to suppress convection in a region underneath the CMB. Gubbins et al. (2015) estimate the thickness of the sub-adiabatic layer to about 750 km. When taking the radial dependence of the thermal conductivity into account, a "sandwich" stably stratified layer could develop at some distance from the CMB, even when Q o is slightly larger than Q a (Gomi et al. 2013, Labrosse 2015. Moreover, the layer developing for Q o < Q a would be much thicker than predicted by Gubbins et al. (2015).
Seismic observations indeed suggest that a stably stratified region may exist just underneath the CMB and estimate a thickness of up to 450 km (Kaneshima 2018). Magnetic field observations  and dynamo simulations (Christensen 2018) may be compatible with a thinner layer of say 200 km but not with the thicker layer of several hundred km suggested by, for example, some of the models discussed by Labrosse (2015).
Before the inner core started growing, the CMB heat flux had to be super-adiabatic. The minimum dissipation condition ϕ = 0 increases the required heat flux Q min o to about 1.1 Q a . Davies (2015) estimates that 3.5 Gyr ago the adiabatic CMB heat flux was Q a ≈ 18 TW while Q o amounted to at least Q min o ≈ 19.3 TW. Based on these considerations, it seems reasonable to assume a CMB heat flux that was at least 10% super-adiabatic before inner core nucleation started and possibly slightly subadiabatic thereafter. Integration in time provides the inner core growth history and yields a rather young inner core with ages ranging between 450 Myr  and 614 Myr (Labrosse 2015). Figure 1, modified from Davies et al. (2015), compares the evolution history that unfolds for different thermal conductivities, different amounts of radiogenic heating, and different density jumps δρ at the inner core boundary. The seismically inferred δρ value constrains the compositional convective power that becomes available at inner core freezing. Evolution models using a classical low value of thermal conductivity (green and blue) yield older inner cores with ages up to 1.8 Gyr and relatively low CMB temperatures and heat flux values 3.5 Gyr ago. The temperatures are in a range where at least partial melting of the lower mantle seems possible.
The revised large high heat flux values (red) lead to younger inner cores, higher CMB temperatures, and huge Q min o values. Adding 300 ppm of potassium reduces the CMB temperature, since less heat flux has to come from secular cooling. These values are compatible with recent coupled evolution simulations for mantle and core by Nakagawa and Tackley (2013), which suggest that CMB temperatures higher than 5000 K and Q o values larger than 40 TW may indeed have been possible during the first Gyr of Earth's history. With estimated lower mantle solidus temperatures around 4200 K (Andrault et al. 2011), this Figure 1. Core-mantle boundary temperature (y-axis) and CMB heat flux (numbers inside symbols in TW) 3.5 Ga ago against the inner core age (x-axis) in a study by Davies et al. (2015). Classical low thermal conductivity values (green and blue) predict an older inner core, lower CMB heat flux and lower CMB temperatures. High thermal conductivity models (red) predict an inner core younger than 600 Myr, high CMB heat flux and larger CMB temperatures. The latter comes down when radiogenic heat by 300 ppm of 40 K is assumed (open red symbols). Symbols connected by lines assume Q o = Q min o before inner core nucleation and Q o = Q a during the nucleation. Results from other studies than Davies et al. (2015) are shown in yellow (Nakagawa and Tackley 2013), pink (Nimmo 2015), orange (Labrosse 2015) and maroon (Driscoll and Bercovici 2014). Figure taken from Davies et al. (2015). (Colour online) means that a substantial part of the mantle would remain molten much longer than classical magma ocean models indicate (Abe 1997). Labrosse et al. (2007) and Ziegler and Stegman (2013) report that a magma ocean at the bottom of the mantle may have delayed the onset of dynamo action by up to 2 Gyr. This negative effect is owed to substantial radiogenic heating and latent heat from magma crystallization that would significantly contribute to the overall heat flow and isolate the core from more efficient cooling. The excess radiogenic heating is caused by the large amount of radioactive elements assumed to have accumulated in the magma ocean during mantle fractionation.
Estimates for the current CMB temperature range around 4000 K and are thus still close to the solidus temperature for lower mantle materials. Partial melting thus seems possible even today and may be responsible for the ultra low seismic velocity zones identified just above the CMB (McNamara et al. 2010). Figure 2 illustrates the power available in a model by Labrosse (2015) that assumes Q o = 1.15 Q a during the whole core evolution. The power-based magnetic field scaling further discussed in section 5 suggests the rms field strength indicated by the solid line in the lower panel. Dots show the diverse paleomagnetic estimates for the geomagnetic field strength. The figure illustrates that the dynamo can be supported before inner core growth started, albeit with a somewhat smaller rms field strength B. The increase in B due to the onset of inner core nucleation is relatively mild because B depends only on the third root of the power, as we will discuss in section 5.

Paleomagnetic evidence for the inner core age
An interesting question is whether paleomagnetic findings could constrain the onset of inner core growth. Biggin et al. (2015) analyze various paleointensity data and report a statistically significant increase in Earth's dipole moment about 1.3 Gyr ago. If indeed related to the start of inner core nucleation, the relatively old age would support a classically small thermal conductivity. However, Smirnov et al. (2016) argue that the dipole moment has been overestimated in many paleomagnetic studies covering the area after the inferred onset of inner core nucleation. The true variation in dipole moment over Earth's history may actually be surprisingly small.
Assuming evolution models similar to the ones discussed above, Landeau et al. history. The results, subsumed in figure 3, suggest that, while the internal field is indeed weaker before than after the onset of inner core nucleation, the surface field, which could potentially be sampled by paleomagnetic methods, shows the opposite behaviour. The reason is the predominant location of dynamo action. When the dynamo is exclusively driven by secular cooling, it operates more evenly throughout the core. In the presence of an inner core, however, dynamo action is preferentially located close to the inner core boundary where buoyancy driving due to latent heat and light element release is strongest. The retreat to greater depths upon the onset of inner core nucleation causes the surface field to decrease. Landeau et al. (2017) and Heimpel and Evans (2013) report another interesting effect that may help to constrain the onset of inner core growth. Dynamo simulations with an inner core show a particularly weak field at high northern and southern latitudes. This can be explained by the fact that the convective columns, which are important for creating the large scale magnetic field, preferentially array around the inner core equator. Flow into and down the (cyclonic) columns tends to advectively concentrate magnetic field where they touch the outer boundary . The resulting pronounced patches are also observed in the historic and present geomagnetic field. Above and below the inner core, convective flows are more plume-like than columnar. When approaching the outer boundary, the rising fluid is mostly diverted towards the tangent cylinder and the convective columns. In terms of a spherical harmonic representation, the high-latitude low in field intensity translates into an axial octupole field contribution that has the opposite sign to the axial dipole. When there is no inner core, however, somewhat stronger field production at depth leads to intensified field at high latitudes and thus to an axial octupole that has the same sign as the axial dipole.
The change in octupole contribution could potentially be observed in paleomagnetic directional data. Kent and Smethurst (1998) report statistically shallower magnetic field inclinations for paleomagnetic data prior to 250 Myr ago, which could indeed hint on a different octupole contribution. Dynamo simulations by Heimpel and Evans (2013) confirm the possible connection to inner core nucleation. This would argue for a young inner core and therefore support the higher heat flux values, but the timing based on the mild directional changes is certainly challenging.
Another interesting feature discussed by Landeau et al. (2017) is the possibility of two alternative solution branches when there is no inner core. Simulations started with a weak magnetic field can remain in a weak multipolar field configuration where convective flows strongly break the equatorial symmetry (see also Driscoll 2016). When initialised with a strong dipolar magnetic field, however, the field remains in a dipolar configuration with stronger fields. The latter is the more likely scenario for the geodynamo. The coexistence of two stable dynamo branches, often called bistability, has been reported by several other authors, for example by Simitev and Busse (2009) or Schrinner et al. (2012), and most recently by Petitdemange (2018).

Alternative scenarios
In order to explain the existence of the geomagnetic field during the time when a magma ocean at the bottom of the mantle may have delayed core dynamo action, Ziegler and Stegman (2013) suggest that the ocean itself may have harboured a dynamo during the first 2 Gyr of Earth's history. Estimates of the magnetic Reynolds number indicate that such a scenario is indeed conceivable.
An alternative driving mechanism for core convection is considered by O'Rourke and Stevenson (2016) and Badro et al. (2016). The solubility of MgO in iron increases with temperature. Today's equilibrium MgO concentration in Earth's core is about 1.1 wt% but it rises by about 2 wt% when the temperature is 1000 K higher. Large impacts during the late heavy bombardment have not only delivered considerable amounts of MgO but also raised the temperature. The slow exsolution of MgO during the subsequent cooling could have driven core convection and thus the dynamo. Badro et al. (2016) estimate that powering the geodynamo over the last 4 Gyr would require large impactors with combined masses comparable to the mass of the moon-forming impactor (2.5% to 20% of Earth's mass). Table 2 compares the non-dimensional parameters inferred for the Earth with those in typical and in the most advanced numerical simulations. Many parameters are several orders of magnitude smaller or larger than one, which translates into a huge span of time and/or length scales. While the magnetic field is dominated by the dipole component, the large Reynolds number indicates the presence of small scale turbulent motions in Earth's core. Compromises obviously have to be made, preferably in a smart way that preserves the key characteristics important for the magnetic field generation.

Fundamental considerations
The small Ekman and Rossby numbers imply that the dynamics in Earth's core is ruled by a balance between pressure gradient, Lorentz or Magnetic force (M), buoyancy or Archimedean force (A), and Coriolis force (C). This so-called MAC balance is often contrasted by other possible candidates that may arise because of a too strong Viscous contribution (VAC) or a too strong Inertial contribution (CIA) in the dynamo simulations. The pressure gradient is always a necessary contribution to guarantee the first order equilibrium. We refer to Christensen (2010), Jones et al. (2011), or Roberts and for recent reviews of the different possibilities and will mostly concentrate on the MAC balance in the following.
A primary balance between pressure gradient and Coriolis force means that the dynamics tries to minimise variations in the direction of the rotation axis so that ∂U/∂z ≈ 0. Flows obeying this Taylor-Proudman theorem are called geostrophic. The small Ekman and Rossby numbers therefore imply that the convection organises itself in the form of Taylor columns illustrated in figure 4. For smaller E values, the columns actually assume a more sheet-like structure, stretching away from the rotation axis. This suggests that there are at least two fundamental length scales: a large length scale in the direction parallel to the rotation axis and a small, predominantly azimuthal, length scale ⊥ . The sheet-like structure may warrant the introduction of a third scale in the direction perpendicular to the rotation axis that we ignore here.
While remains close to the system length scale D, ⊥ strongly depends on the parameters. Linear stability analysis for the onset of (non-magnetic) convection shows that ⊥ ∼ E 1/3 . Increasing the Rayleigh number leads to faster flows, larger Rossby numbers, less geostrophic geometries, and smaller ⊥ , as is also demonstrated in figure 4. When Lorentz forces enter the leading order force balance, however, ⊥ increases significantly, as we will discuss below.

Towards earth-like solutions
While it is desirable to decrease the Ekman number as much as possible, the affordable value is often dictated by the available numerical resources but also depends on the specific problem. When, for example, exploring the rarely happening geomagnetic field reversals, particularly long simulations are required, which are only possible at not too small Ekman numbers. Once the Ekman number has been decided on, the other parameters should be chosen such that the simulations yield dynamics as close as possible to the geodynamo process. To guarantee dynamo action, the magnetic Reynolds number has to be large enough so that the magnetic field production overcomes Ohmic dissipation. Christensen and Aubert (2006) report a minimum required value of Rm = 50 for their suite of models. For a given numerical model, dynamo action is guaranteed by increasing the Rayleigh number beyond a critical value Ra D where Rm becomes large enough. Close to onset, the dynamo is often sub-critical in the sense that, once the dynamo is established, Ra could be lowered to values somewhat below Ra D while still maintaining dynamo action (Morin andDormy 2009, Petitdemange 2018). When Ra exceeds a second critical value Ra m , a more geodynamo like dipole-dominated configuration is replaced by multipolar solutions where the axial dipole looses its dominant role.
Like Ra c and Ra D , Ra m depends on the other system parameters. The non-linear inertial effects are often instrumental in breaking symmetries and may also play a role in disrupting the ability of the flow to sustain a stable dipole-dominated field. Christensen and Aubert (2006) confirm this idea by showing that the transition to multipolar dynamos can be characterised by a modified Rossby number Including the typical flow length scale ⊥ in the definition is an attempt to account for the gradient in the advection term of the Navier-Stokes equation (1). Christensen and Aubert (2006) estimated ⊥ based on the kinetic energy spectrum: wherel denotes the energy-weighted mean spherical harmonic degreē E l being the energy of degree l flow contributions. The transition to multipolar fields happens at Ro ≈ 0.1 (Christensen and Aubert 2006). Note, however, that the original definition introduced by Christensen and Aubert (2006) uses a typical length scale that may slightly deviate from ⊥ . In order to establish how the transitional Ro depends on the control parameters, Olson and Christensen (2006) fit an empirical law to numerical dynamo results and find This is little intuitive but at least shows that decreasing E allows for larger Ra values and thus larger flow amplitudes while still maintaining a dipole-dominated field at Ro < 0.1. The true meaning of this threshold and the role of Ro is still debated. Several authors argue that viscous friction still plays a too large role in most numerical simulations (Soderlund et al. 2012, King and Buffett 2013, Roberts and King 2013, Oruba and Dormy 2014, Dormy 2016. Viscous effects could, for example, determine the typical flow length scale and could also influence the transition from dipole-dominated to multipolar dynamos (Soderlund et al. 2012, Oruba andDormy 2014). This may certainly be the case for the more typical runs at Ekman numbers larger than say E = 10 −5 but seems less likely for the most recent extreme simulations that we will further discuss below. Christensen et al. (2010) use four different criteria to more rigorously quantify the similarity between the geomagnetic field and simulation results: dipole dominance, equatorial symmetry, axial symmetry, and skewness, i.e. the tendency of the field to concentrate in "patches". Figure 5 subsumes their results in the parameter space spanned by the magnetic Ekman number E λ = E/Pm = τ Ω /τ λ (x-axis) and the magnetic Reynolds number Rm = τ λ /τ U (y-axis). The new low Ekman number cases of Aubert et al. (2017) and Schaeffer et al. (2017) have been added. The importance of the two parameters E λ and Rm suggests that the ratios of flow turnover time τ U = D/U, magnetic diffusion time τ λ = D 2 /λ, and rotation time τ Ω = 1/Ω play an important role. More Earth-like simulations fall into the wedge outlined by dashed lines in figure 5. Its location confirms that lower Ekman numbers and Rayleigh numbers, significantly, larger than Ra D but smaller than Ra m , are required to yield an Earth-like magnetic field. Figure 6 illustrates how closely the simulation indicated by a red triangle in figure 5 resembles the geomagnetic field. The similarity includes the weaker and even inverse field at high latitudes, the pronounced patches close to the tangent cylinder, and the strong patches around the equator.
The Pr dependence has, for example, been explored by Simitev and Busse (2005), who report that for small values Pr < 1 magnetic fields are generally not dipole-dominated.
The role of the magnetic Prandtl number is not well understood. There is a minimum magnetic Prandtl number, the critical Pm c , beyond which dynamo action becomes possible. Based on their suite of numerical simulations, Christensen and Aubert (2006) derive the empirical rule Pm c = 450 E 0.75 . Pm can thus only be decreased in conjunction with E. Close to Ra D , increasing either Ra or Pm leads to larger Rm values and can therefore promote dynamo action. However, there is a limit to this effect. When increasing Ra, the flow and thus also the magnetic field become increasingly small scale, which in turn increases Ohmic diffusion.
Since increasing Pm not only increases Rm but also decreases E λ , figure 5 suggests that this offers an alternative (and certainly numerically cheaper) path to more Earth-like simulations. Dormy (2016) and Dormy et al. (2018) report that a third dynamo branch can indeed be found when increasing the magnetic Prandtl number at low Rayleigh numbers. This branch is characterised by a particularly strong dipole-dominated magnetic field. The respective dynamos may therefore more closely obey the MAC force balance planetary dynamo are thought to operate in than most of the other dipole-dominated numerical dynamos at lower Pm.
However, increasing Pm may never lead to reversing dynamos ) and generally has a much smaller effect on the dynamics than increasing Ra. The complex dynamics and field structure of the geomagnetic field seems to require a sufficiently non-linear flow and thus larger super-critical Rayleigh numbers. More research is clearly required to further elucidate the role of Pm and to explore the properties of the interesting third strong-field dynamo branch. Dormy (2016) speculates that these strong-field solutions can be connected to the geodynamo by a uni-dimensional distinguished limit where both Ekman number and magnetic Prandtl number approach small values. Aubert et al. (2017) present such a path in parameter space, inspired by the MAC balance, that leads from one of their already Earthlike low-E simulations to even more realistic solutions. Uni-dimensional means that the path only depends on one parameter called here: The sub-index zero indicates the parameters of a reference simulation. A value of ≈ 10 −7 would correspond to Earth-like conditions. Keeping the Prandtl number fixed to Pr = Pr 0 = 1, Aubert et al. (2017) decrease step-wise from 1 to = 3.33 × 10 −4 . The path starts at the reference configuration with E 0 = 3 × 10 −5 and Pm 0 = 2.5, typical for many dynamo simulations, and currently reaches E = 10 −8 and Pm = 0.045, the lowest Ekman and Prandtl numbers achieved in a full dynamo simulation so far. Values below E = 3 × 10 −6 become feasible with a smart form of hyperdiffusion which progressively damps smaller flow and codensity scales beyond 3 ⊥ . The use of hyperdiffusion has been criticised for strongly influencing numerical results. However, Aubert et al. (2017) verify that their smart form has no drastic impact by comparing respective simulations with fully resolved results at Ekman numbers E ≥ 3 × 10 −6 . The path approach has the advantage that it preserves the already Earth-like properties in the reference solution, for example Rm, , and the similarity to the geomagnetic field according to the criteria by Christensen et al. (2010). Other properties, like M A or the presence of fast waves (Aubert 2018), become increasingly realistic, and the MAC balance is more and more closely obeyed.

Reproducing secular variation
Geomagnetic field variations cover a broad range of time scales. The fastest signal from the core are geomagnetic jerks, which may take some months. Slow variations in reversal frequency over some tens of million years form the other extreme of the spectrum. The time scale of more typical secular variations roughly obey the relation where the index denotes the spherical harmonic degree l and τ l is the typical time scale for magnetic field contributions B l of degree l: Lhuillier et al. (2011) and Christensen et al. (2012) report a typical secular variation time scale of τ SV ≈ 450 yr. This excludes the axial dipole contribution that is exceptionally stable and has a time scale of roughly a millennium. The inverse dependence of the time scale spectrum on degree l is consistent with an advective transport of magnetic features by a velocity field with τ U ≈ τ SV . Geomagnetic accelerations are characterised by a shorter time scale of τ SA ≈ 10 yr, nearly independent of the spherical harmonic degree, at least for l ≤ 10 (Holme et al. 2011, Christensen et al. 2012). Dynamo simulations not only reproduce the τ l and τ SA spectra but also the ratio τ SV /τ SA (Christensen et al. 2012, Aubert 2018.
Because of the unrealistic parameters, there are several alternatives for rescaling the dimensionless time in the numerical simulations to real time. For analysing secular variations, the flow overturn time τ U is the natural choice, while τ λ may be more appropriate for longer time scales. Both alternatives yield the same result when the magnetic Reynolds number Rm = τ λ /τ U assumes Earth-like values of about 1000.
The unrealistic Alfvén Mach number in the numerical models means that torsional oscillations and other magnetic waves are too slow compared to the overturn time. The simulations by Aubert (2018) demonstrate that magnetic waves indeed become faster and more pronounced at more realistic parameters. Aubert (2018) also reports that quasigeostrophic Alfvén waves are a promising candidate for explaining geomagnetic jerks. However, since the waves only account for a small fraction of the total secular variation, their unrealistic time scale seems of secondary importance. Similar inferences may hold for the 60 yr time scale that could be explained by magneto-gravitational waves (Buffett 2014) in a stably stratified layer underneath the CMB.
A typical secular variation feature is the westward drift of up to 14 km/year at low latitudes in the Atlantic hemisphere. Low latitude secular variations in the Pacific hemisphere, on the other hand, are particularly weak, as is illustrated in figure 7. Aubert et al. (2013) report that, while a standard dynamo simulation cannot reproduce this dichotomy, they succeeded with a specific setup. Their Coupled Earth model (CE) combines a CMB heat flux pattern inspired by seismic lower-mantle observations with an inner core that grows faster in one hemisphere than in the other. Moreover, the inner core and mantle are coupled via gravitational torques and the outer boundary condition is stress-free. This results in a dynamics where inner core and mantle rotate eastward while the flow has a mean westward direction in order to preserve angular momentum. Figure 7 compares the secular variation in a standard dynamo, the CE model, and from geomagnetic data. The flow pattern in the CE simulation is characterised by a large-scale, wave-number m = 1 gyre, which only comes close to the surface in the Atlantic hemisphere and thus explains the secular variation dichotomy. Recent simulations at E = 10 −7 (Schaeffer et al. 2017) demonstrate that such a feature may also form spontaneously for a standard setup with homogeneous boundaries, at least when the Ekman number is small enough. A very similar gyre has been inferred from geomagnetic observations (Pais and Jault 2008, Aubert 2013, Gillet et al. 2013, Bärenzung et al. 2018).

Force balances and scaling laws
Since the simulations cannot run at realistic parameters, they have to be scaled to planetary conditions for a more rigorous comparison. Several scaling laws have been suggested over the years, and recent reviews can be found in Christensen (2010) and Jones et al. (2011). Most concern simple key properties, like the rms magnetic field strength B, the rms flow velocity U, the typical length scale ⊥ , or the Nusselt number Nu, and describe how these properties depend on the system parameters, typically in the form of power laws. These laws can be purely empirical, using fits to simulation data, but gain more credibility when based on theoretical considerations, for example one of the different candidate force balances.
Ideally, the numerical data clearly support one of these candidates and thus point towards the force balance that rules the dynamical regime covered by the simulations. Another important test is whether the scaling laws also explain observations, for example the field strength of planetary dynamos. We discuss these different aspects for a scaling law based on the MAC balance in section 5.1.
In section 5.2 we then review examples for a direct analysis of the force balance and also illustrate the impact of the magnetic field on the flow structure. Indirect evidence, on the other hand, is provided by torsional waves and the Taylorization, discussed in section 5.3.

The MAC balance
In the following, we concentrate on the scaling laws in the MAC balance, following largely the proposal by Davidson (2013). When taking the curl of the Navier-Stokes equation (1) the pressure gradient drops out and the amplitude of the magnetic, Archimedean, and Coriolis terms can be estimated to Here we have used ∇ × (J × B) ∼ B 2 / 2 ⊥ , assuming that the dominant flow length scale ⊥ is also responsible for the dominant magnetic field gradients. The expression for the total power (24) suggests that we can approximate the energy input per mass by where V is the volume of the shell. The buoyancy term in equation (42) can thus alternatively be written as: The balance between Lorentz and Coriolis forces provides an estimate for the magnetic energy per mass, while the balance between buoyancy and Coriolis force yields an estimate for the flow length scale, Since roughly represents the system dimension (see above), only one additional expression is required to link the four unknowns B, U, ⊥ , and P . One possibility is the balance between the energy input and diffusion. In planetary dynamo regions, where diffusion is clearly dominated by Ohmic effects, the energy input is roughly balanced by Ohmic dissipation: The energy loss per mass can then be approximated by where Ohm the typical length scale characteristic for Ohmic dissipation. This suggests a second relation for the magnetic energy density: Here, f Ohm is the ratio of Ohmic dissipation to the total dissipated power. While f Ohm is likely very close to one in planetary dynamos, this is not necessarily true for the numerical simulations because of the much larger magnetic Prandtl number. The magnetic scale Ohm is likely the scale where magnetic induction and Ohmic diffusion balance, i.e. where the respective magnetic Reynolds number is around one: U Ohm Ohm /λ ≈ 1. This suggests where ω Ohm is the typical value of the flow vorticity at length scale Ohm . Turbulence theory seems to show that the smaller scale vorticity is also determined by the driving power P (Davidson 2013). Dimensional arguments for a pure power dependence then yield or Le = B/(ρ 1/2 μ 1/2 Ω ) ∼ f 1/2 in dimensionless form (Davidson 2013). Figure 8(b) demonstrates that dynamo simulations and also the geodynamo seem to roughly obey this scaling. A least-square fit of log(Le) to log(P) yields a best-fit exponent of 0.31, slightly smaller than 1/3. However, when considering the large spread of the data, the agreement seems convincing enough. Further support comes from the fact that the scaling also successfully predicts the magnetic field strengths of Jupiter and of fully convective stars (Christensen et al. 2009). Figure 8 includes geomagnetic values: Estimates for Earth's rms core field range from about 0.25 mT, inferred from downward continuing the observed field to the core-mantle boundary, to the 4 mT predicted from torsional oscillation frequencies . The velocity required for estimating Earth's Rossby, Reynolds, or magnetic Reynolds number is based on the typical secular variation time scale τ SV (see above). The power dissipated by Earth's dynamo has been assumed to range between the minimum estimate of 0.2 TW from Christensen and Tilgner (2004) and a 10 times higher value. The scalings shown in the top row of figure 8 successfully connect the simulations with the geomagnetic estimates.
Combining the two estimates for the magnetic energy, equation (45) and equation (51), with equation (46) for length scale ⊥ yields a scaling law for the flow amplitude: which reads in dimensionless form: This allows to eliminate the U-dependence in the ⊥ scaling: Ohm P 1/9 .
The three relations for Le, Ro, and ⊥ form the main result of the MAC scaling introduced by Davidson (2013). Panel (a) of figure 8 shows that the best-fit exponent for the P dependence of Ro in the numerical data is 0.41. This is close to the theoretical prediction of 4/9 = 0.4 and nicely connects the flow amplitude in the simulations with that inferred for Earth. Because of the limited P range in combination with the small exponent, we cannot expect the numerical simulations to support the scaling for the length scale (46). No related plot is therefore shown here. Typical simulations at E = 3 × 10 −5 reach a dimensionless power that is about 10 7 times larger than for Earth. This means that the typical length scale in these simulations would only be about a factor six too large, which seems good news for dynamo simulations. We have seen above that the length scale progressively decreases with growing Ra and decreasing E in non-magnetic convection. Once the Lorentz force contributes to the leading order force balance, however, this tendency is efficiently mitigated. While this is reminiscent of so-called magnetostrophic systems, where a strong enough imposed field enforces ≈ ⊥ , the situation is obviously somewhat more complex in full simulations where the field is dynamically created. Additional scaling relations can be derived when combining the above expressions. Equating equation (45) and equation (49) yield: When ignoring the weak dependence on the power, this agrees with the empirical fit to dynamo simulation results by Christensen and Tilgner (2004). The fact that dynamo simulations can reach realistic Rm values would mean that not only ⊥ but also Ohm are already reasonably well represented in current dynamo simulations . As already mentioned above, viscous effects may determine ⊥ in at least the larger Ekman number dynamo simulations. Several authors report that ⊥ indeed scales with E (King and Buffett 2013, Roberts and King 2013, Oruba and Dormy 2014, possibly with the E 1/3 dependence expected for non-magnetic convection. However, the interpretation of the ⊥ -scaling is particularly difficult. One reason is that ⊥ spans little more than one order of magnitude in typical current dynamo simulations. Another reason is that the results can depend on the specific method used for estimating ⊥ (Dormy et al. 2018). When using, for example, equation (36), the additional small scales that are excited at smaller Ekman numbers could already decrease the ⊥ estimate, even when the dominant flow scale remains unchanged. Aubert et al. (2017) therefore identify ⊥ with the scale at which the kinetic energy spectrum peaks. This definition, however, can also be somewhat ambiguous since the spectrum is relatively flat around the peak. More research on the parameter dependence of the flow length scales certainly seems in order.
A combination of the Le and Ro scalings, namely (52) and (54), provides an estimate for the Alfvén Mach number: Ohm P 1/9 . (57) Using once more the Le-scaling allows to establish a relationship between magnetic Reynolds number and Elsasser number: The weak dependencies on the power rationalises why the typical dynamo simulations are already quite realistic in terms of Rm, M A or , which scales like Le 2 . When plotting log-log scaling laws, there is a certain danger that the chosen nondimensionalization introduces a spurious trend when x and y axis are multiplied with different powers of a parameter or physical property. Christensen (2010) demonstrates that the simulations continue to support the scaling laws discussed above, independent of the chosen scales. Panels (c) and (d) of figure 8 show the scaling for flow and field amplitudes when using the viscous diffusion time τ ν = D 2 /λ as a time scale instead of the rotation period. For the U scaling, the best-fitting exponent 0.48 is now somewhat larger than the predicted 4/9, while the best-fitting exponent 0.33 for the B scaling nearly perfectly agrees with the prediction. The alternative non-dimensionalization also highlights the special role of the low Ekman number simulations by Aubert et al. (2017). Note, however, that the authors discuss slight deviations from the MAC scaling in an attempt to account for the possible effects of hyperdiffusion in some of their simulations.
The mild differences in the exponent for the different suggested scalings illustrates the still existing uncertainties. A discussion of alternative scalings can, for example, be found in Christensen and Aubert (2006), Christensen (2010), or Stelzer andJackson (2013). The numerical simulations may support an additional weak Pm dependence, and Christensen and Aubert (2006) show that the results more convincingly fall on the scaling lines when an additional factor of Pm 0.11 is included. Christensen and Aubert (2006) also assume a CIA balance rather than a MAC balance and use an empirical scaling for the Ohmic dissipation power (Christensen and Tilgner 2004) instead of equation (52). This leads to slightly smaller scaling exponents of 3/10 for the power dependence of Le and 2/5 for the power dependence of Ro. Figure 8 illustrates that the data are not sufficient to distinguish these alternatives (red lines) from the MAC scalings (blue lines). More rigorous statistical analysis may indicate a preference for one or the other scaling (Stelzer and Jackson 2013) but likely remains inconclusive unless additional extremely costly low Ekman number simulations will close the gap to the geodynamo parameters even further. Oruba and Dormy (2014) point out that power based scaling laws necessarily reflect the statistical balance between power input and dissipative output and are likely biased by the fact that viscous dissipation plays a too strong role in the numerical models. Scaling laws that only depend on the control parameters would certainly be preferable.

Accessing the true force balance
Several authors directly quantify the force balance in their simulations , Soderlund et al. 2015, Dormy 2016, Yadav et al. 2016, Schaeffer et al. 2017, Dormy et al. 2018. Calculating, storing, and analysing all forces can be demanding in terms of numerical resources and book keeping. Like in the MAC balance discussed above, the curl of the Navier-Stokes equation (1) and thus the vorticity equation is therefore often considered, which has the advantage that the pressure gradient and the conservative contributions of Coriolis and Lorentz forces drop out. This procedure will, however, emphasize the small scale contributions and may thus not be representative for the balance that rules at larger scales (Yadav et al. 2016). Yadav et al. (2016) therefore compare the rms of the sum of Coriolis force and pressure gradient with the rms of the individual other forces in the Navier-Stokes equation (1), leaving out the boundary layers where viscous effects necessarily become dominant in their rigid flow simulations. They illustrate that the MAC balance is gradually approached when decreasing the Ekman number from E = 10 −4 to E = 10 −6 . At E = 10 −6 , viscous and inertial forces are at least one order of magnitude smaller than the MAC contributors (including the pressure gradient). The effect of the Lorentz forces on the flow in their E = 3 × 10 −6 runs is illustrated in figure 9. The convective structures are clearly larger in the dynamo simulation (right) than in the non-magnetic sister simulation (left). Another effect of the Lorentz force is that the Nusselt number, a measure for the heat transport by convection, is amplified by up to a factor two. Aubert et al. (2017) adopt the force balance analysis proposed by Yadav et al. (2016) and in addition calculate the rms forces at different spherical harmonic degrees l. Figure 10 compares the respective spectra for two of their models, a standard case at E = 3 × 10 −5 and a lower Ekman number case with E = 3 × 10 −7 that uses smart hyperdiffusion. In both simulations, pressure gradient and Coriolis force form the primary balance, which guarantees a high degree of geostrophy in the flow. The residue is very effectively balanced by buoyancy and Lorentz forces for length scales smaller than about ⊥ . Inertial forces amount to about 10% of the Lorentz force at E = 3 × 10 −5 and to about 2% at E = 3 × 10 −7 . Viscous forces increase with growing l (decreasing length scale) but are at least another decade smaller than inertial forces for l < 30 and never enter the leading order balance. Buoyancy decreases significantly for scales smaller than ⊥ (l > 10), which suggest that these scales are predominantly driven by non-linear cascades that transport energy from larger to smaller scales. The smart hyperdiffusion progressively damps scales beyond l = 30 and clearly affects the rms forces beyond l = 50.
Other recent low Ekman number simulations at E ≥ 3 × 10 −7 (Soderlund et al. 2015) and E ≥ 10 −7 (Schaeffer et al. 2017) have not been analysed to the same level, but seem to obey the MAC balance to a similar degree.

The Taylor state
Dynamos that strongly obey the MAC balance operate in the Taylor state. This refers to the balance of azimuthal forces on so-called geostrophic cylinders that have the rotation axis as a common geometry axis. The only purely geostrophic motions that perfectly obey the Taylor-Proudman theorem can be described by the relative motions of these cylinders: U G (s) = sω(s)φ, where s is the cylindrical radius. When integrating the azimuthal component of the Navier-Stokes equation (1) over these cylinders, pressure gradient, buoyancy, and Coriolis force drop out. Since this leaves the Lorentz force as the only potential remaining candidate from the MAC balance, the integral must involve some cancelation that brings the axisymmetric azimuthal Lorentz force down to the level of the next integrated force, either inertia or viscosity. This cancellation, also called Taylorization, is quantified by the normalised integral over geostrophic cylinders where z − and z + refer to the bottom and top of the cylinders andF L is the axisymmetric azimuthal Lorentz force. The smaller the global Taylorization T = {Tay}, the clearer a dynamo obeys the MAC balance; curly brackets denote the mean over all geostrophic cylinders here. Disturbances of the Taylor state lead to Torsional Oscillations (TOs), one-dimensional Alfvén waves that rely on the magnetic coupling between geostrophic cylinders as a restoring force and travel away or towards the rotation axis. Their typical time scale is the Alfvén time based on the rms magnetic field B s perpendicular to the rotation axis: TOs have been observed in the geomagnetic field , Gillet et al. 2015 and in dynamo simulations , Teed et al. 2014, Schaeffer et al. 2017, Aubert 2018) and figure 11 shows prominent examples for both. They tend to have large wave lengths ≥ 0.1 D and a time scale of around 6 yr in Earth. For a more in-depth explanation of the Taylor state and TOs we refer to Roberts and King (2013) and Jault and Légaut (2005). Estimating the expected Taylorization for Earth is difficult. Because of the large TO wave numbers, the most effective viscous contribution to the force balance is the friction at the boundaries. The respective force can be estimated tō where d = E 1/2 D is the thickness of the Ekman shear boundary layer and the overbar indicates the azimuthal average. The integrated viscous force on the geostrophic cylinders is then roughly where we have used that the integral only contributes over the Ekman layer thickness. Lorentz force and the inertial advection term can be estimated to and Balancing the remaining integrated Lorentz force by viscous friction (61) thus suggest a global Taylorisation: Here we have assumed that the geostrophic flow amplitude U G is not much different from U and thatF L ≈ F L . Using ⊥ ≈ 10 −2 D, suggested by equation (55), yields a value of T ≈ 10 −8 for Earth. Balancing Lorentz force by advection, on the other hand, suggests where we have assumed that the typical advective force (64) also provides a reasonable approximation for the respective azimuthal force on geostrophic cylinders. Numerical simulations indeed show that the strong rotational constraint organises the convective columns in a rather efficient way that minimises the cancellation over geostrophic cylinders (Gastine and Wicht 2012). Formula (66) suggests T ∼ 10 −4 for Earth, which means that inertial effects clearly dominate over viscous effects. When the Taylor state is disturbed by torsional oscillations, the respective inertial term ∂U G /∂t provides the balance, suggesting Here we have used the Alfvén time scale τ A = D(ρμ) 1/2 /B and once more assumed that U G is not much different from the typical velocity U . For Earth, this formula yields T ∼ 10 −4 and thus a value very similar to that indicated by the advective balance (66). This suggests that torsional oscillations do not significantly reduce the Taylorization in Earth's core. Figure 12 compares T values from simulations with E ≤ 3 × 10 −5 from various authors. The data are rather scattered but seem to be more compatible with the linear dependence on M A suggested by scaling (67) than with the quadratic dependence implied by (66). This suggests that the Taylor state is significantly disturbed most of the time in the simulations, perhaps more so than in Earth. The fact that even the low Ekman number simulations by Aubert et al. (2017) and Schaeffer et al. (2017) deviate from the proposed scaling indicates Figure 12. Comparison of the Taylorization T in numerical simulations with E ≤ 3 × 10 −5 from various authors (filled symbols). The simulations from Aubert et al. (2017) are mostly compatible with the linear dependence on M A suggested by equation (67), but only when the additional dependence on the length scale ⊥ is ignored (solid vs. dotted line). The quadratic dependence expected in a Taylor state according to equation (66) seems little compatible with the simulation data. Open symbols show the prediction according to the viscous balance equation (65), which clearly underestimates T , in particular for the low Ekman numbers. Note that the cases by Aubert et al. (2017) use stress-free boundary conditions. The dashed-dotted line shows the predicted vicously limited T along the path defined by Aubert et al. (2017). (Colour online) that the role of the Taylor state and of torsional oscillations are not completely understood. Figure 12 also shows that the viscous limit (open symbols) clearly underestimates the Taylorization for most of the simulations.
The large interest in the Taylor state and torsional oscillations is not only motivated by the fact that they relate to the MAC balance. The torsional wave time scale (60) provides an estimate of the typical magnetic field strength in the core not accessible by other means (Gillet et al. 2015). Moreover, the transfer of angular momentum between torsional oscillations and Earth's mantle can potentially explain length-of-day variations not accounted for by other mechanisms (Gillet et al. 2015).
The Taylor state is also the basis for an alternative approach of solving the dynamo problem. When inertial and viscous forces are neglected altogether, the Navier-Stokes equation (1) becomes diagnostic. The flow can then be calculated from a given codensity distribution and magnetic field, provided these fulfil a solvability condition. In his original work, Taylor (1963) showed that this condition is what we nowadays call "Taylor's constraint", T = 0. He also devised a semi-analytical concept for calculating U. The induction equation (2) and the codensity equation (3) still have to be integrated in time, but potential numerical difficulties related to the Navier-Stokes equation (1) are avoided. This comes at the cost, however, that one has to find a flow field which induces a magnetic field that (sufficiently) obeys Taylor's constraint.
Following this scheme has proven difficult, and so far only highly parameterised meanfield approaches have been considered. More recently, Wu and Roberts (2015) succeeded in finding an axisymmetric mean field solution in a spherical shell. We refer to Roberts and Wu (2014) for a recent overview of the subject. One may, however, question whether it is desirable to completely get rid of inertia. Roberts and Wu (2014) show how Taylor's theory could be modified to at least incorporate torsional oscillations, but the revised approach would still lack Reynolds stress driven winds, turbulent energy cascades, and possibly also magnetic field reversals. Should the value of T ≈ 10 −4 be confirmed for Earth, the inertial effects may indeed prove too influential to be neglected.

Data assimilation
Data assimilation, a mathematical framework originally developed for weather forecasting, combines observations and physical models to provide a more complete view of a system (see Kalnay 2003, for a review). While being widely used in meteorology and oceanography, its applications to deep Earth geophysics like mantle convection or core dynamics has only started. Geomagnetic data assimilation confronts dynamo simulations with geomagnetic observations. This may lead to more Earth-like numerical dynamo models and improved reconstructions of the core dynamics in the period spanned by geomagnetic observations. In analogy to its most wide-spread application in weather forecasting, geomagnetic data assimilation also yields refined forecasts of the geomagnetic field evolution.

General formalism
We start by introducing the general data assimilation formalism. The state of the system of interest is defined by a generic model state vector x of size N x , which subsumes the unknown variables but can also contain model parameters. In the dynamo context, x can comprise the spectral representation coefficients of magnetic field B, flow U and codensity perturbations C: where subindex T corresponds to the transpose operator, l and m denote spherical harmonic degree and order, and k numbers the radial level r k in the dynamo region. The true state x t we want to solve for is constrained by the observation vector y o of size N o , where H is an observation operator that screens out the generally sparse observations from the full state vector. The observation error o subsumes different contributions like measurement errors, noise, or errors due to the observation operator H. The observation error statistics is described by the error covariance R = o oT and a mean typically assumed to cancel out: o = 0. Angular brackets indicate mean or expectation values. Because the number of observations N o is typically much smaller than the number of unknowns N x , and because o is often not well known, we cannot invert equation (69) directly. Data assimilation provides additional information in the form of a background state x b , which is related to the true state by where the statistics of the error b obeys the covariance matrix P b = b b T . Data assimilation is a general concept for combining observational data and background information in order to find a model x that best estimates the true state x t . The probability p(x | y o ) for a state x, given the observations y o , can be written according to Bayes theorem as This is usually referred to as the posterior distribution, while p(x) represents the background knowledge or prior information, p(y o | x) the likelihood of the observations conditioned to the model, and p(y o ) = p(y o | x)p(x) dx is a normalising factor known as the evidence (Carrassi et al. 2017). Assuming all distributions are Gaussian, one can write the likelihood as the prior distribution as and the posterior distribution as Here J is a cost function defined as Different criteria are used to find what is deemed the best estimate for the state x, also called the analysis x a . The maximum likelihood estimate is given by the mode of the posterior distribution in (74). An alternative is the minimum variance estimate, which is given by the minimum of equation (75) and coincides with the mean of the posterior distribution (Talagrand 1997). Both estimates coincide for truly Gaussian distributions.

Incorporating the dynamics
The extension of the formalism to the time domain is accomplished by incorporating the model dynamics, for example in the form of numerical dynamo simulations. The dynamical model M i,i−1 propagates the state vector from time t i−1 to t i : where the error M i represents the inevitable imperfect character of the model dynamics. Its statistics is given by the model error covariance Data assimilation techniques follow two main approaches of implementing the temporal domain: variational and sequential. In the variational approach, the cost function J is evaluated for the entire time window spanned by observations. The goal is to minimise J by finding the right initial condition x 0 that brings the whole model trajectory as closely to the observations as possible. Finding the minimum of J requires calculating the gradient of the cost function, ∇ x J(x), which proves to be prohibitively expensive at N x ∼ O(10 6 ). The adjoint method (Le Dimet and Talagrand 1986) offers an alternative but requires a backwards integration in time in order to establish the sensitivity of the observations to the initial conditions. Since variational methods have not been explored much in geomagnetic data assimilation so far, we focus on the sequential approach.

Sequential data assimilation
Sequential assimilation methods proceed by a succession of forecasts and analysis steps. A commonly used method is the Kalman Filter algorithm (KF, Kalman 1960), which assumes a linear model so that M i,i−1 in equation (76) is simply a matrix. The forecast step, starting with the initial state x a i−1 , reads Note that the forecast state x f represents a time dependent background state, and that in addition to the model error M i we also have to take into account the propagated initial error The model covariance is also forwarded in time, which proves to be computationally challenging since it involves N x model realizations.
Whenever observations are available, an analysis step is performed where the forecast is corrected in order to better explain the observations: Here is the Kalman Gain matrix, which formulates a compromise that depends on model and observation uncertainties (via the covariances) and propagates information from observables to the whole system state. More details on sequential filtering can be found in Cohn (1997). The analysis provides a new initial condition for the next forecast. Ideally the state is drawn closer to its true value with each cycle.
For a non-linear system like the geodynamo, we can use the full dynamics to forecast the state vector but would have to rely on the tangent linear model approximation M = ∂M/∂x for forecasting the covariances, as in the Extended Kalman Filter (EKF) formalism. In order to avoid the related numerical costs, however, many studies assume a stationary model covariance in an approach commonly known as Optimal Interpolation (OI).
The Ensemble Kalman Filter (EnKF, Evensen 1994) offers a more compelling method where an ensemble of models is used to approximate the model error covariance. Although the ensemble approach is still computationally demanding, the forecast phase is easily parallelised, neither equation (78) nor equation (80) have to be computed directly, and there is Figure 13. Example of an EnKF data assimilation for the surface magnetic field of degree l = 12 and order m = 8 in a dynamo model. The ensemble mean is given by circles and the error bars show the σ confidence interval based on the ensemble spread. After the initialisation, the ensemble assimilates noisy observations (stars) over a sequence of forecast (blue) and analysis (red) cycles and quickly approaches the reference true state from which the observations were created. Though the data assimilation is stopped at t = 0, the model continues to evolve close to the reference trajectory. For longer time scales, the predictive power naturally degenerates due to the chaotic nature of the system. (Colour online) no need to store the covariance matrix for the next forecast. An example of an EnKF based on an ensemble of N e = 512 dynamo models is shown in figure 13. The assimilation shown in the figure is from a synthetic twin experiment where different realizations of the same model are used to generate observations and to compose the background ensemble. The dynamo model is marginally Earth-like and with Rm = 300 and E/Pm = 10 −5 would be located in the middle lower part of the wedge in figure 5. The moderate E = 10 −4 makes integrating the large ensemble of models possible. In the figure, mean (dots) and spread (colored envelope) of the ensemble represent the best estimate of the model and its uncertainty respectively. The assimilation of a sequence of observations (stars) draws the dynamo models progressively closer to the true state. After the assimilation period, the ensemble is free to evolve unconstrained by observations, providing a forecast for the rest of the time window.

Peculiarities of geomagnetic data assimilation
Geomagnetic observations from magnetic observatories (INTERMAGNET 1 ) and satellite missions (e.g. Swarm, CHAMP, Ørsted) measure the field at different altitudes. The signal is a convolution of several sources. In addition to the Earth's core dynamo, induced and remnant crustal magnetisation, ionospheric and magnetospheric currents, and currents induced in the ocean also contribute (Hulot et al. 2015). Internal and external fields can be separated in the spherical harmonic domain based on the different functional form of the respective potential fields. However, since the internal field remains an entanglement of crustal and core signal, core field models are restricted to the range l ≤ 14 where the crustal contribution remains subdominant. Notable core field models are the historical gufm1 (Jackson et al. 2000), the more observatory-based COV-OBS (Gillet et al. 2013), and the more satellite-based CHAOS (Finlay et al. 2016).
To avoid the problematic separation of the different contributions, geomagnetic data assimilation has so far treated these spherical harmonic internal field models as "data". For l ≤ 14, this means that only N o = 224 spherical harmonic coefficients are available to constrain a state vector that even for a moderate dynamo simulation already comprises N x ≈ 10 6 unknowns . What comes to the rescue here are the symmetries and correlations favoured by the core dynamics, which statistically tie the unknowns to the observations. This crucial information is carried by the background model covariance matrix. Figure 14 illustrates the correlations between magnetic and flow field in a standard dynamo simulation. In spectral space, the covariance is dominated by correlations between alternating degrees but within the same order m. This reflects the preferred equatorial symmetry and the fact that the magnetic field is predominantly produced by non-axisymmetric flows acting on axisymmetric field. The correlations in grid space illustrate the preferred equatorial antisymmetry in the magnetic field but equatorial symmetry in the flow field. Also shown is the relatively strong correlation, with depth which is a consequence of the overall columnar flow structure.
For a successful data assimilation approach, it seems crucial to select a dynamo model that reproduces the fundamental aspects of the geomagnetic field as best as possible (see section 4). Coupled Earth has therefore been the model of choice in several studies so far, for example Aubert (2014). In order to combine model and observations, the magnetic field strength and time scale of the nondimensional simulations have to be rescaled to the geodynamo context, which is typically done by assuming an rms field strength of 10 3 µT and a secular variation time scale of τ ⊕ SV = 415 yr (see section 4.3).

Geomagnetic data assimilation attempts
Geomagnetic data assimilation studies can differ considerably in the background dynamo model, the data assimilation method, or their main objectives. Possible objectives are (i) testing the approach with synthetic experiments, (ii) assessing the impact of the different models in order to identify the most Earth-like, (iii) describing the dynamics over the time span where observations are available, and (iv) predicting the evolution of the geomagnetic field beyond that time span. Examples for each objective are discussed in the following.

Synthetic experiments
Synthetic experiments use noisy outputs from dynamo simulations as observations, as in figure 13. Since the true state is known, the assimilation errors can actually be quantified. Due to the complexity of self-consistent dynamo simulations, there are only very few geomagnetic data assimilations studies that use the variational approach. A recent application for a simplified, inertia-free, dynamo model can be found in Li et al. (2014). The authors show that the initial condition of the system, and hence its full trajectory, can be retrieved when not only the magnetic field but also the flow are observed at the CMB. Another noteworthy example, using the EnKF and a full self-consistent dynamo model, is the experiment by Fournier et al. (2013). They explore different ensemble sizes and report that 500 ensemble members suffice to assure convergence. The spin-up time until their model more closely follows the true state takes roughly one millennium.

Model testing
Model testing has been used by Aubert (2013) to establish which ingredient makes Coupled Earth a more Earth-like model, as already discussed above. While Aubert (2013) applies the method to individual snapshots (epochs), a sequential assimilation approach seems more appropriate to represent the dynamical evolution. Such an approach has been used for instance by Tangborn and Kuang (2015) who assumed a simple model covariance matrix where only correlations of the magnetic field with observations of the same degree l and order m are taken into account. Moreover, they use a simple prescribed depth dependence. Tangborn and Kuang (2015) report that the combination of strong convective forcing and a deeper reaching correlation improves the forecast. In an updated version, Kuang and Tangborn (2015) show that the additional use of the secular variation as observation considerably helps to constrain the internal state, although the spin-up time of the assimilation remains long.  Barrois et al. (2017) with corrections proposed by Barrois et al. (2018). On the left, the panels show the core flow streamlines, evidencing the gyre underneath the Atlantic Hemisphere, superimposed with the COV-OBS.x1 radial geomagnetic field. The middle panels display the combined contribution of diffusion and subgrid scale action to the secular variation signal from observations. On the right, the panels show the streamlines and colour-coded horizontal divergence of the flow that translates to up-and downwellings. The corrected figure has been produced by O. B. (Colour online)

Recovering geodynamo dynamics
Modeling the dynamics over the time span where data are available is usually referred to as reanalysis. Aubert (2014) use data from geomagnetic field models and secular variation in combination with the CE dynamo model to infer the state vector for the whole core for a succession of snapshots. The more complex sequential approach has so far only been attempted for simpler dynamical models. Instead of arising as a solution to the Navier-Stokes equation, the flow evolution obeys a stochastic process at the CMB, and the interaction with the magnetic field is computed by using the so called frozen-flux induction equation which neglects magnetic diffusion. The reduced size of the stochastic model allows using a large ensemble in an EnKF scheme. Barrois et al. (2017) and Bärenzung et al. (2018) apply this approach to the assimilation of COV-OBS.x1 (Gillet et al. 2015) magnetic field and secular variation models with slightly different setups. Barrois et al. (2017) implement a parameterisation of magnetic field diffusion and subgrid scale interactions in the frozen flux equation to infer the unmodelled interactions of field with rapid flow underneath the CMB. Figure 15 shows their core flow reanalyses for three different epochs. The figure differs from the original published in Barrois et al. (2017) because of corrections proposed by Barrois et al. (2018). A prominent feature is the large scale gyre already discussed above.
The study by Bärenzung et al. (2018) shows that the gyre has a strong non-geostrophic contribution, with higher amplitudes in the southern hemisphere but faster accelerations in the northern hemisphere. The study also tracks another secular variation feature that recently got more attention: the high latitude jet underneath Siberia (Livermore et al. 2017), which seems to have accelerated during the last four decades. The mechanism behind the gyre dynamics and its possible connection to the high latitude jet underneath Siberia remains unknown.

Predictions
Predictions of the geomagnetic field have several practical applications. They are particularly important because of the role the geomagnetic field plays for space weather, which can heavily impact Earth-orbiting satellites and space missions (Mandea and Purucker 2018). Every five years, different groups from the geomagnetism community make a joint effort to update the International Geomagnetic Reference Field (IGRF, Thébault et al. 2015) 2 . They provide model coefficients for the large scale internal magnetic field and its secular variation based on different methodologies and data selections. Based on an intricate weighted averaging, the models are combined to generate the reference field and its secular variation (see Thébault et al. 2015a for details on the calculation of the last IGRF). Assuming a linear time dependence of the coefficients, the IGRF also provides a magnetic field prediction for the next five years. Despite the simple method, the level of accuracy is more than adequate for the short term horizon (e.g. Bärenzung et al. 2018). On time scales longer than a decade, however, the linear extrapolation becomes increasingly less reliable and more involved models accounting for the flow dynamics begin to pay off (Aubert 2015).
The predictability of the nonlinear dynamo process is generally limited. A common measure is the e-folding time τ e , the time it takes for an initial error to be exponentially amplified. While the e-folding time of Earth's atmosphere is about 2 to 3 days, it amounts to τ e ≈ 30 yr for Earth's core processes ). This is a crude global measure, however, that may not apply to specific geomagnetic features.
One of the most important geomagnetic field features for space weather is the South Atlantic Anomaly (SAA). The SAA is a region of anomalously low magnetic field intensity where high energetic particles can more easily penetrate Earth's atmosphere. This, for example, increases the chances for electronic disruptions in low orbit satellites electronics (Heirtzler 2002). The SAA has intensified and expanded during the last century  and is also moving westward. A prediction exercise for the SAA evolution based on dynamo simulations by Aubert (2015) is shown in figure 16. The study focuses on a "single epoch" EnKF assimilation approach, in which the main field and secular variation from COV-OBS and CHAOS-5 are used to estimate initial conditions for an ensemble of Coupled Earth dynamo models, which is then integrated forward in time. Hindcast tests, in which past observations are used to forecast the present field, show that predictions of the SAA are trustworthy for about 50 years. The results from Aubert (2015) indicate that the SAA will intensify and split into two patches over the next few decades.
The decrease of the dipole component over the last roughly 400 years has raised the question whether this could indicate the beginning of a geomagnetic field reversal (Olson et al. 2009). Data assimilation predictions suggest that the decrease may continue at least for the next century at a nearly constant rate (Aubert 2015). Reversals are difficult to explore with full data assimilation approaches, however, since they happen very rarely and last about ten thousand years. Morzfeld et al. (2017) therefore use simplified low-dimensional models for the evolution of only the dipole component. They sequentially assimilate data from paleomagnetic intensity records, using a Particle Filter method instead of an EnKF, and report that this allows predicting the reversal likelihood for a 4000 yr period. Meduri and Wicht (2016) perform a statistical analysis of long dynamo simulations and of paleomagnetic records and show that the likelihood for a reversal only increases significantly once the dipole moment has decreased below 50% of its mean value. This means that today's dipole low has no statistical significance.

Conclusion
Our overview concentrates on selected issues that received particular attention during the last few years. The revised high thermal conductivity values indicated by ab initio simulations and some laboratory experiments continue to challenge the classical scenario for driving Earth's dynamo. Future experiments will hopefully resolve this issue or confirm the need for alternative ideas.
We also discussed the question whether recent numerical simulations operate in the MAC force balance expected to rule the geodynamo. The most advanced simulations run at E = 10 −7 , or even E = 10 −8 when using smart hyperdiffusion. While this is still seven orders of magnitude too large, the scaling behaviour and detailed analysis suggest that the primary force balance is indeed Earth-like. Secondary forces are still not as small as they should be, but may nevertheless have only a negligible impact on the dynamo process. Naturally, the very small scales and/or very short term dynamics are not represented correctly. The models are thus not adequate for studying the turbulence present in the Earth's core but seem to do a very good job in modelling the dynamo.
Since the numerical simulations do not exactly follow the predicted scaling laws, with slight differences in the second digit of the scaling exponent, there may also be room for improvement in theory and/or for even more extreme dynamo models. However, going the extra mile in parameter space will be extremely costly, and the results may not be very different from today's low Ekman number simulations.
Having said this, we should also stress that the recent low-E runs show some interesting phenomena that have not been reported before. They may not prove relevant for primary dynamo mechanism but are certainly interesting when it comes to capturing the detail of geomagnetic field dynamics. Examples are the fast quasi-geostrophic Alfvén waves discussed by Aubert (2018) or the large gyre that seems to form spontaneously in the simulation by Schaeffer et al. (2017). Both phenomena certainly deserve more in-depth studies.
Another interesting development is that numerical dynamo simulations are actually starting to have practical applications. As integral part of data assimilation frameworks, they can provide valuable missing information on the geomagnetic field dynamics and thus make advanced predictions of the future magnetic field evolution possible.
Some topics were only mentioned in passing or neglected altogether. We discussed torsional oscillation in some depth, but there are several other types of waves that may also play a role in Earth's core. First work goes back to Hide (1966) and we refer to  for an overview and to Hori et al. (2015) or Buffett et al. (2016) for recent studies.
While the codensity approach is still the method of choice for most dynamo simulations, there are also some attempts to take the different properties of thermal and compositional convection into account (Manglik et al. 2010, Nakagawa 2011. Numerical limitations usually restrict the respective studies to cases where the compositional diffusivity is only ten times larger than the thermal diffusivity, while differences in the order of 10 3 would be appropriate. Bouffard et al. (2017) present an interesting new approach where the evolution of the compositional field is described with a practically diffusion free particlein-cell method. First results promise interesting new solutions where, for example, the light material emanated from the growing inner core front could accumulate underneath the core-mantle boundary to form a stably stratified layer. It remains to be seen whether this could be an explanation for the layer suggested by some seismic studies (Kaneshima 2018).

Disclosure statement
No potential conflict of interest was reported by the authors.