Inferring internal properties of Earth ’ s core dynamics and their evolution from surface observations and a numerical geodynamo model

Over the past decades, direct three-dimensional numerical modelling has been successfully used to reproduce the main features of the geodynamo. Here we report on efforts to solve the associated inverse problem, aiming at inferring the underlying properties of the system from the sole knowledge of surface observations and the first principle dynamical equations describing the convective dynamo. To this end we rely on twin experiments. A reference model time sequence is first produced and used to generate synthetic data, restricted here to the large-scale component of the magnetic field and its rate of change at the outer boundary. Starting from a different initial condition, a second sequence is next run and attempts are made to recover the internal magnetic, velocity and buoyancy anomaly fields from the sparse surficial data. In order to reduce the vast underdetermination of this problem, we use stochastic inversion, a linear estimation method determining the most likely internal state compatible with the observations and some prior knowledge, and we also implement a sequential evolution algorithm in order to invert time-dependent surface observations. The prior is the multivariate statistics of the numerical model, which are directly computed from a large number of snapshots stored during a preliminary direct run. The statistics display strong correlation between different harmonic degrees of the surface observations and internal fields, provided they share the same harmonic order, a natural consequence of the linear coupling of the governing dynamical equations and of the leading influence of the Coriolis force. Synthetic experiments performed with a weakly nonlinear model yield an excellent quantitative retrieval of the internal structure. In contrast, the use of a strongly nonlinear (and more realistic) model results in less accurate static estimations, which in turn fail to constrain the unobserved small scales in the time integration of the Correspondence to: J. Aubert (aubert@ipgp.fr) evolution scheme. Evaluating the quality of forecasts of the system evolution against the reference solution, we show that our scheme can improve predictions based on linear extrapolations on forecast horizons shorter than the system e-folding time. Still, in the perspective of forthcoming data assimilation activities, our study underlines the need of advanced estimation techniques able to cope with the moderate to strong nonlinearities present in the geodynamo.


Introduction
The Earth's fluid core is a dynamic, yet sparsely observed system.Direct or indirect measurements of the planet magnetic field are the main source of data used to probe the dynamical state of the core, giving access only to the largescale image of the magnetic field poloidal component over the external surface of the system.Integral constraints based on geodetic data also provide a useful secondary source of data.Under the assumption that the large-scale field temporal variations are dominated by a diffusionless induction process (Roberts and Scott, 1965), the radial magnetic field and its rate of change (called the secular variation) have been used over the past decades to estimate the fluid flow below the core surface, at the origin of the temporal variations (see e.g.Finlay et al., 2010a, for a review).In conjunction with length-of-day data and dynamical models of torsional waves, the knowledge of these core flows can be used to probe the magnetic field strength deep in the core (most recently Buffett et al., 2009;Gillet et al., 2010).Other dissipative constraints on the magnetic field strength can also be derived from short timescale (daily) measurements of the Earth's nutations, as recently done by Buffett (2010).Strategies aiming at inferring the properties of Earth's core dynamics from surface observations commonly encounter problems of non-uniqueness, and spatial resolution problems.The errors induced by these problems now tend to exceed by far the measurement errors, which have become very small at the age of satellite magnetic observation.As an illustration of a non-uniqueness problem, closure constraints are needed in order to perform core flow inversions (most recently the quasi-geostrophic assumption, Pais and Jault, 2008).The nonlinearity involved in the core flow problem is an instance of a spatial resolution problem, in the sense that the unresolved small scales of the magnetic field are responsible for part of the observed secular variation at large scale, in a way which is difficult to predict.This, in turn, complicates the evaluation of the flow large scales (Eymin and Hulot, 2005;Gillet et al., 2009).One elegant, yet not much explored way to handle both problems is through enforcing dynamical consistency of the solutions, that is, solving for a velocity field with a time evolution consistent with first principle evolution equations.The problem of determining core flow then becomes an inverse problem where an initial condition is sought, the evolution of which will subsequently explain and predict the observations at various points in time, hereby all dynamically connected together.Along these lines, the now expanding field of geomagnetic data assimilation aims at optimally combining physical laws and observations of Earth's core dynamics (see Fournier et al., 2010, for a recent review).This expansion capitalizes on the progress made over the last twenty years by data assimilation techniques in other fields of research, most importantly (from the viewpoint of core dynamics), atmospheric dynamics and physical oceanography.Regarding the specific three-dimensional Kalman filter which is central to the following, the interested reader is referred to the authoritative monograph by Evensen (2009) for a detailed description of its theoretical foundations, and its extension to nonlinear problems in the form of the so-called ensemble Kalman filter.For further reading, Kalnay (2010), Brasseur (2006), Elbern et al. (2010) and Houser et al. (2010) review recent applications of the Kalman filter to the analysis of the atmosphere, of the ocean, of air quality and of land surfaces, respectively.
First principle equations suitable for the direct modelling of core dynamics (Braginsky and Roberts, 1995) are now routinely solved numerically, and have had considerable success in reproducing the first-order features of the geomagnetic field: morphology and dipole dominance of the field, secular variation and reversals (main advances recently reviewed by Christensen, 2011).These equations include the induction equation, Navier-Stokes equation with convection described in the Boussinesq approximation, and an equation for the transport of a buoyancy field (which, in the Earth's core, is of both thermal and chemical origin).The main difficulty faced by these three-dimensional, self-consistent simulations is the current impossibility to reach numerically the physical parameters of natural dynamos.This is related to the great disparity between the diffusion coefficients of the thermal, chemical, magnetic and velocity fields.As the situation is not likely to improve in the foreseeable future, progress in the field has been achieved over the recent years by identifying and scaling phenomena where some or all of these diffusivities play only a secondary role.For instance, a large set of numerical models has revealed that the magnetic field strength does not depend on any diffusion coefficient, only on the available power to drive the dynamo (Christensen and Aubert, 2006;Aubert et al., 2009).A connected study (Christensen et al., 2010) showed that a morphological similarity can be obtained between the geomagnetic field and the output of numerical dynamos if only three time scales are in reasonable proportion when compared to their Earth counterparts: the rotation period of the planet, the characteristic time scale of advection of the magnetic field by the fluid flow, and the characteristic time for the diffusion of the magnetic field.
Here we wish to use the information and dynamically consistent solutions provided by numerical geodynamo models in order to carry out inverse modelling.Our long-term aims are (i) to estimate the dynamical state of Earth's core from surface observations, (ii) to assess the extent to which such estimations are affected by (or immune to) non-uniqueness and spatial resolution problems, and (iii) to determine the magnitude of the associated errors.Linear estimation of the system state invisible parts (also called Kalman filtering or stochastic inversion) is a method of choice for this type of problem.Its efficiency is classically tested through the procedure of synthetic (twin) experiments: a reference solution is first computed, and used to generate a catalog of surface data, which are in turn used to recover the solution, starting from a wrong initial guess.Previous attempts (Liu et al., 2007) used parameterised, ad-hoc covariance properties to perform such estimations, and focused mostly on the evolution of the observed, surficial part of the system.The novelty of our approach stands in a preliminary numerical computation of the system multivariate covariance properties.This approach has already been used in a companion paper (Fournier et al., 2011), addressing the two-dimensional core flow problem described above.Here we proceed to the determination of the three-dimensional internal structure, which can then be used as an initial condition for time evolution, opening the way to data assimilation practice.In the following, Sect. 2 presents the numerical model and inversion technique.Section 3 presents the results of the numerical experiments, which are then discussed in Sect. 4.

Numerical geodynamo model
Our numerical model is formulated as in Aubert et al. (2009).We solve for the velocity field u, magnetic field B, and codensity (or density anomaly) field C in a spherical fluid shell between radii r i and r o , of aspect ratio r i /r o = 0.35 rotating about an axis e z , using the equations Here r is the radius vector.The fundamental scales underlying the dimensionless scheme are the inverse shell rotation rate −1 for time, the shell gap D for length, (ρµ) 1/2 D for magnetic induction, where ρ is the fluid density and µ the magnetic permeability of the fluid.The kinematic, magnetic and thermal Ekman numbers are defined as Here ν, λ, κ are respectively the viscous, magnetic, and thermochemical diffusivities of the fluid.As detailed in Aubert et al. (2009), the distribution of boundary mass anomaly fluxes can be determined from a parameterised thermodynamical model of Earth's core evolution.Here we choose an idealised situation which is thought to be representative of Earth at present (Lister, 2003), where the mass anomaly flux F originates entirely from the inner boundary, and the outer boundary has null mass anomaly flux.The Rayleigh number Ra Q is thus where g o is gravity at radius r = r o .As shown in Aubert et al. (2009) the dimensionless volume sink term for mass anomaly corresponding to this situation is S T /ξ = −3/(r 3 o − r 3 i ).The other boundary conditions at both boundaries are no-slip for velocity, and insulating for the magnetic field.The numerical implementation PARODY-JA is used in this study (Dormy et al., 1998;Aubert et al., 2008).The fields u, B are expanded into toroidal and poloidal scalars, which, together with the scalar C, are described using a finite-difference scheme in the radial direction with up to 160 grid points, and a spherical harmonic decomposition in the lateral directions up to degree and order 133.
Table 1 summarizes the properties of the two models which have been integrated for this study.The models have been chosen so as to provide end-members in physical complexity and semblance to the geomagnetic field.Model  Christensen and Aubert (2006).The mean harmonic degrees l u,B in the velocity and magnetic field are as defined in Christensen and Aubert (2006).The median harmonic orders m med u,B are the orders which, on average, separate the velocity and magnetic power spectra in two domains of equal energy.The e-folding time τ e and secular variation time τ sec are as defined in Lhuillier et al. (2011a,b).The morphological semblance to the geomagnetic field χ 2 is defined according to Christensen et al. (2010).Bottom two rows: parameters relevant to the determination of the model covariance matrix P. Unless otherwise noted, the numerical experiments use a covariance matrix which is determined from the number n of free run samples reported in the table.Also reported are the number of radial nodes used for the determination of the matrix, the degree and order l P max up to which this matrix is determined, the size of the state vector (complex coefficients) and the number of coefficients involved in the determination of P (or reduced state vector size).

Model 2
Fig. 1.Dynamical magnetic fieldline imaging (DMFI) representation of the numerical models magnetic field line structure.The magnetic field lines are rendered as grey tubes, with thicknesses proportional to the local magnetic energy density.The inner core surface is color-coded according to the radial magnetic field strength.The outer surface is color coded similarly, and made transparent with an opacity proportional to the radial magnetic field strength (see Aubert et al., 2008, for further details).

25
Fig. 1.Dynamical magnetic fieldline imaging (DMFI) representation of the numerical models magnetic field line structure.The magnetic field lines are rendered as grey tubes, with thicknesses proportional to the local magnetic energy density.The inner core surface is color-coded according to the radial magnetic field strength.The outer surface is color coded similarly, and made transparent with an opacity proportional to the radial magnetic field strength (see Aubert et al., 2008, for further details).Christensen et al. (2010).The main morphological difference arises from the high concentration of magnetic flux into a small number of patches, and also in the lack of equatorial field features (which was not taken into account in the ratings of Christensen et al., 2010).Model 2 fares much better with respect to semblance to the geomagnetic field.It can also be considered a physically more suitable model because the three time scales mentioned in the introduction (magnetic advection, magnetic diffusion, rotation period) are in better proportion when compared to the Earth than for model 1.
The ratio of the magnetic diffusion and magnetic advection time scale is the magnetic Reynolds number Rm, which is very close to the value of about 800 which can be expected in the Earth's core if surface velocity estimations are representative of the deep flow (Christensen and Tilgner, 2004).
The ratio of the length of day and the magnetic diffusion time scale is the magnetic Ekman number E λ , which in model 2 is one order of magnitude closer to the value E λ ≈ 10 −9 expected in the core (e.g.Christensen and Aubert, 2006).As a consequence of the higher Reynolds numbers, model 2 also has stronger nonlinearities than model 1, and the effect of a stronger magnetic advection results in an apparent decorrelation between the surface and internal magnetic structures (see Aubert et al., 2008, and also Fig. 1).The increased temporal complexity of model 2 can be quantified using the ratio of the e-folding time of the system τ e , or time constant for the exponential divergence of two infinitesimally close solutions (Hulot et al., 2010b;Lhuillier et al., 2011a), to the characteristic time scale for secular variation τ sec (Christensen and Tilgner, 2004;Lhuillier et al., 2011b).The increased spatial complexity can be evaluated through the mean harmonic degrees l u,B , as defined in Christensen and Aubert (2006), and median harmonic orders m med u,B in the power spectrum of the velocity and magnetic field.It is important to mention that although they are quite different, both models have more than half of their energy within the harmonic order range m = 0 − 13, which means that there is reasonable hope that surface observations supplied within the same spectral range could constrain well the internal structure.

Rescaling the model output
Although this is not fundamental when only synthetic data are used, any attempt to integrate time-dependent geomagnetic data into a numerical model will require the dimensionless model output to be rescaled to the geophysical world.If the model operated at the same parameter values as the Earth's core, it would be enough to use the canonical scales presented in the last section.However, as the model operates far from Earth's core conditions, we have to resort to units underlain by scaling principles known (or thought) to hold both in the model and in the Earth's core, so that the various quantities, once presented in these new units, should have similar values in the model and in the core.
Following previous work on the secular variation time scale (Christensen and Tilgner, 2004;Lhuillier et al., 2011b;Fournier et al., 2011), time will be presented in units of τ sec , which is roughly 500 yr in the Earth's core.Velocity will be presented in units of D/τ sec , which is roughly 4.4 km yr −1 in the core.Following the Christensen and Aubert (2006) scaling, magnetic field will be presented in units of f 1/2 ohm (ρµ 3 p 2 D 2 ) 1/6 , where p is the convective power density in the shell and f ohm the fraction of this power which is dissipated through Ohmic effects.Using the high-power scenario presented in Aubert et al. (2009), one magnetic field unit then amounts to 1.7 mT for the present Earth.Finally, the relationship between the convective power, codensity and velocity (Eq.A1 in Christensen and Aubert, 2006) prompts to present the co-density field in units of pτ sec /g o D. Using again the high-power scenario, one codensity unit amounts to 10 −5 kg m −3 .

Best linear unbiased estimate of the internal structure from surface data
For a given discrete instant t i in a numerical dynamo simulation, we define a (column) state vector (superscript T denotes the transpose) which contains the complex values of the poloidal (superscript p), toroidal (superscript t) harmonic scalars of the fields u,B and the co-density C, for each harmonic degree and order l,m on the nodes j of the radial grid.The full size of x is on the order of one to ten million elements (see Table 1).Covariance matrix calculations presented in this study however use a version of x which is decimated by taking a subset of radial nodes and harmonic coefficients (see Sect. 3.1 and Table 1), yielding a typical size on the order of ten to a hundred thousand elements.The state vector is centered and normalised for the time series to have zero mean and unit variance.We assume that the various elements in x have a Gaussian distribution, with probability density function (pdf) P(x) ∝ exp(−x P −1 x/2), where P is called the covariance matrix of the model and the prime symbol denotes the transpose complex conjugate.The validity of the Gaussian assumption is explored in Fournier et al. (2011).Some deviations from a Gaussian behavior can be expected, in which case the best linear estimate about to be derived is not optimal anymore in the sense of maximum a-posteriori pdf, but still remains the estimate of minimum variance.
The state vector has an observable part and a hidden part.Our goal is to provide an estimate of the hidden part from the observable part.We define an observation operator H which extracts its observable part from x.The observation operator is thus a rectangular matrix with a number of rows equal to the size of the observations y and a number of columns equal to the size of the state vector x.Because we are dealing with the equivalent of global field models at the core-mantle boundary, and given their current resolution limits (e.g.Olsen et al., 2010), we define the observable part of the state vector as the poloidal magnetic field B p lm (r o ) at the outer boundary up to degree and order 13.The corresponding observation operator contains ones in the entries corresponding to an observed quantity and zeros otherwise.The operator contains an additional sub-block when the rate of change of B p lm (r o ) is observed as well (up to degree 13).Considering the radial part of the induction Eq. (2) on the fluid side of the outer boundary, where the velocity field vanishes, we obtain Prescribing the time derivative of the outer boundary poloidal magnetic field is thus equivalent to prescribing the radial component of ∇ 2 B p lm .The sub-block dedicated to ∂B p lm (r o ,t i )/∂t thus contains a discrete Laplacian operator written on the fluid side of the outer boundary.The vector Hx comprises up to 210 elements.
Our inverse problem seeks a state vector x such that where y is a set of observations, statistically centered and normalised using the same means and variances as those used for x.In a general context, the observations bear some error o , with a covariance matrix R = E( o o ), where E stands for the expected value.In other words, the likelihood of y if x is realised is Given the above mentioned sizes of x and y, the problem posed by Eq. ( 12) is vastly underdetermined.Our preferred estimate of x is the best linear unbiased estimate, which minimizes the functional This estimate is the most likely given the data and model covariance properties (see Fournier et al., 2011, for details).
Looking for the extrema of J (x) one finds the best linear inverse solution, which takes a simple form when x is centered: with the Kalman gain matrix Equation ( 15) is ubiquitous in geomagnetic field modelling (Gubbins, 1983) and core flow modelling (review in Finlay et al., 2010a).In both cases it is usually referred to as stochastic inversion, but we also note that it is formally identical to one of the Kalman filter equations (e.g.Fournier et al., 2010).The stochastic inversion will be more efficient when P contains strong correlations between the observed and unobserved parts of the state vector.Our primary goal being to test the accuracy and prediction power of the inversion with synthetic data, we consider the data error-free i.e. o = 0 and R = 0, in an attempt to isolate errors resulting from the inversion from all other error sources which can arise in a realistic context.
www.nonlin-processes-geophys.net/18/657/2011/ Nonlin.Processes Geophys., 18, 657-674, 2011 In a time-dependent context, the stochastic inverse can be used to initialise the numerical model and perform forecasts of the system evolution.This forms the backbone of sequential data assimilation.At a given later time where the numerical model has updated the state vector to a value x f (the forecast), an analysis of the system can be performed by comparing the observed part of the forecast Hx f to the observations y available at the analysis time.The state vector is then corrected to the new value x a such that The analysed state vector is then used as a new starting condition, completing an assimilation cycle.Each time the system is corrected, P should be updated according to a corresponding evolution equation (the set of equations updating P and the state vector and describing their time evolution is called the Kalman filter, Kalman, 1960).As this is a computationally very intensive task, P is often assumed to be time independent, which amounts to assuming that although the analysis tends to reduce the error on the system state knowledge, the nonlinear dynamics operating between two analyses restores this error back to its natural (free run) value.When using such a frozen covariance matrix, the assimilation technique is referred to as optimal interpolation (or OI, see e.g.Kalnay, 2003, for a review).Atmospheric and oceanic data assimilation usually resort to matrices operating in physical space (Kalnay, 2003, §5.4).There, the choice of an a-priori defined correlation length results in sparse banded structures which are easy to process, even with state vectors with sizes comparable to what is reported in Table 1.The accuracy of such an approach depends on whether in-situ measurements are available with sufficient quantity.The geomagnetic assimilation case is completely different because of the lack of in-situ measurements.In that context, the information cannot be efficiently and accurately propagated radially downward past the correlation length if the above approach is employed.We thus need to process a covariance matrix with a full structure, obtained in spectral space, in order to perform an accurate and efficient propagation of this information.It can then be understood that practical computational considerations set a limit on the size of such a matrix, which can thus update only a subset of the state vector.This limitation was not present in our previous study (Fournier et al., 2011), where the two-dimensional character of the problem permitted high-resolution inversions.When used with this limitation, scheme ( 16) can be unstable due to the deleterious influence of the uncorrected variables.Here we control these instabilities using a slightly modified version of the OI scheme, where only a fraction β (0 ≤ β ≤ 1) of the forecast is re-injected at analysis stage If β = 0 then at each analysis time, the system is set to x = 0 (time average of the dynamo simulation) before the analysis Table 2. Summary of the time-dependent assimilation scheme used in this study (noise-free data).

Preliminary computations
Covariance matrix P determined from a free run (Sect.3.1), frozen for the entire duration of the assimilation run, Observation operator H (Sect. 2.3), Kalman gain matrix K = PH HPH −1 .

Initialisation step
x f (t 0 ) = 0 (time average of the simulation).
Analysis step is performed; each analysis is thus an inverse of the corresponding data with no memory from previously inversed data.If β = 1 the analysis corrects the full forecast resulting from the previous time integration; the whole forecast is thus re-injected for the next analysis cycle.A value of β between 0 and 1 will help mitigate the two possibilities.There is no theoretical justification for the introduction of β in the analysis.The justification is practical, as we observed that the use of the regular Kalman filter analysis (β = 1) resulted in overenergetic estimates of those variables not directly impacted by the observations and the truncated covariance matrix P.

Computation and structure of the model covariance matrix
A correct determination of P is central to the quality of the inversion ( 14)-( 15).Here we approximate P using the multivariate statistics of the numerical model.This matrix is thus computed during a preliminary "free run" of the model, where a numerical integration is performed and a large number n (see Table 1) of state vector snapshots x(t i ) are extracted, with a typical time lag between the snapshots on the order of the e-folding time of the system to ensure decorrelation between snapshots.In terms of classical dynamo time scales, the duration of the free run is 17.5 magnetic diffusion times D 2 /λ.In terms of the advective time rescaling used in this study, the duration is about 580τ sec , amounting to 290 000 yr if τ sec = 500 yr.Once each time series of the state vector components is centered and normalised to unit variance, if the vectors x(t i ) are stored as columns into a matrix X, then P can be estimated through For reasons of storage and cpu time limitations, about half of the numerical grid radial nodes are used to determine P, the remaining nodes (see Table 1) being computed by linear interpolation.We have checked that this has no impact on the quality of the inversion.Likewise, only harmonic coefficients up to degree and order l P max are retained.In order to capture the correlations involving the dominant scales of the system, we set l P max so that it exceeds both m med u,B and l u,B (see Table 1).The convergence of the coefficients defining P is checked by monitoring the effect of doubling and decimating the number of samples used to build the matrix (see also Fig. 8a).
According to Eq. ( 15), not all coefficients of P are actually needed to compute the matrix K: only the correlations involving one observed quantity have an effect on the result.We thus represent on Fig. 2 sub-blocks of P displaying the correlations between B p lm (r o ) and one of the fields (here u p lm at various radii; the structure is similar for other fields, with some differences detailed below).Harmonic coefficients lm have been ordered along a one-dimensional lexicographic scheme gathering all possible degrees (up to l P max = 15 in this case) for each given order.
In the case of model 1, Fig. 2 shows that the internal structure is linearly coupled with the surface observations, as could be expected from the visualisation presented in Fig. 1.Correlation coefficients are indeed up to 0.9 when correlations between B p lm (r o ) (or its time derivative) with another field at depth are considered.Correlations between harmonic coefficients of same order m result from the linear coupling terms present in Eqs.(1)-(3).In particular, the Coriolis force (which is a dominant force acting on the fluid) produces correlations between harmonic coefficients which have different degrees l, provided again that they share the same harmonic order.Correlations between harmonic coefficients of different orders are almost non-existent, an indication of weak J. Aubert and A. Fournier: Inferring internal properties of Earth's core dynamics nonlinear couplings.The dominant azimuthal wavenumber of the surface magnetic field thus reflects that of the internal convection flow, as expected from a dynamo mechanism (e.g Olson et al., 1999;Aubert et al., 2008) where the magnetic energy is sustained through stretching, twisting and folding of the magnetic field lines by a columnar convection flow.The correlations tend to peak at or around the dominant azimuthal wavenumber of the dynamo, as measured by m med u,B in Table 1.The influence of the leading equatorial symmetry properties of the flow and magnetic field is seen in the checkerboard pattern of the correlations: if l obs and l deep are respectively the harmonic degree of B p lm (r o ) (or its time derivative) and the internal field, then l obs + l deep needs to have odd parity if the internal field is u p lm , C lm or B t lm , and even parity if the internal field is B p lm and u t lm .As the radius decreases towards the inner-core boundary, the correlation matrix displayed in Fig. 2 exhibits an upper triangular structure: coefficients with l obs > l deep tend to preserve their correlation with depth (track 1 of Fig. 2) while coefficients with l obs < l deep tend to lose their correlation as depth increases (track 2).We ascribe this effect to the dominance of the Coriolis force, leading to a strong correlation between magnetic field patches at the surface (large l obs ) and convection columns at depth (lower l deep ).
The covariance matrix of model 2 (not shown) has a similar visual structure, with less marked correlations peaking at about 0.7 for the magnetic field and 0.3 for its secular variation.The increased nonlinear dynamics indeed tends to blur the linear relationships between the surface and deep fields.One could expect to see increased correlations between harmonic coefficients of different order m as a result of the same nonlinear dynamics.This is indeed the case but the signal remains small, with cross correlations peaking at less than 0.1.In general, nonlinear dynamics is thus not beneficial to the correlation between the surface and deep fields, but a reasonable predictive power of the deep structure from surface observations can still be expected.

Synthetic inversion tests, model 1
Once the matrices P and K are computed (see Table 2), we then proceed to the computation of a reference time series of the model which will be used to benchmark the efficiency of the inversion.About one hundred to a few hundreds of snapshots spaced by 0.01, 0.05, 0.1 and 0.2 time units (respectively equivalent to 5, 25, 50 and 100 Earth years if τ sec ≈ 500 yr) are extracted (in selected cases, a spacing of 0.02 time units or 10 yr is also used).The surface poloidal magnetic field and secular variation harmonic coefficients in this reference time series will be subsequently referred to as the "data".A twin run is next initialized from a wrong initial guess (the time average of the free run, see Table 2).The data are then injected in the assimilation algorithm and the quality of the reference trajectory recovery is evaluated.In addition to being started from different initial conditions, the reference run and its twin may have different physical parameters, in order to simulate the effects of modelling errors arising from an imperfect physical description of the system (Liu et al., 2007).Here, however, we wish to isolate the errors associated with the inversion for the deep structure from all other sources of errors, and we thus use the same set of physical parameters for the free run determining P, the reference and the assimilation runs.
A first qualitative evaluation of the static and timedependent inversions is presented in Fig. 3.The first column represents the reference state at an arbitrary time t of the reference time series.The quality of the internal structure retrieval depends on the stage that the assimilation has reached at time t (see Table 2).If the assimilation is initialised exactly at that time, its state (second column) is the time average of the model, the best guess one can make in the absence of data.If the assimilation is initialised and analysed at that time, the resulting static inversion (third column) considerably improves the estimation.The knowledge of P allows for instantaneous propagation of the information to all fields throughout the whole shell, giving very reasonable (but visibly underpowered) estimates of the internal fields.The low strength of the estimated fields is a general property of the linear estimation based on correlations (see for instance Fournier et al., 2011).The detailed agreement between the reference and the recovery is further improved if the assimilation has already performed several cycles when time t is attained (fourth column): this shows that the dynamics has a beneficial effect in determining the amplitude of the scales which are neither observed, nor correlated to the surface observations.A quantitative assessment of the recovery quality is presented in Fig. 4, where we compute the energetic misfit M u and correlation coefficient C u between the estimated (est) and reference (ref) velocity fields, Here V is the shell volume.The corresponding quantities M B and C B are also presented for the magnetic field.For the sequence shown in Fig. 3, misfits are moderately low, with the energy of the difference between estimated and reference amounting to about 20 % to 50 % of the reference energy.
Correlation coefficients are very high, on the order of 0.7 to 0.9.The recovery is better for the velocity field than for the magnetic field, because of the richer small-scale content of the latter (due to a magnetic Prandl number larger than one).
It is worth noticing that the recovery quality undergoes sizeable fluctuations and does not monotonically increase with the number of analyses.In the present case, strong fluctuations are experienced as the reference time series happens to pass through a state of abnormally low kinetic energy (lowest energy at time 1.97),only experienced a few times in the free  2).Column 3 is the assimilation state if only the initialisation step and the first analysis step as been performed at time t.Column 4 is the assimilation state if it has been initialised at time t = 0, and if 14 analyses and forecast steps have been subsequently performed with a time lag 0.1τ sec (about 50 yr in Earth time) until time t = 1.3τ sec (note that this column then represents the forecast, not the analysis).For this run we used β = 0.75.All fields are presented in the units proposed in Sect.2.2, with their tentative rescaling to Earth's core values in parentheses.
run computation used to build P. The internal structure estimates are thus worse in the vicinity of this statistical outlier.
The fact that the recovery quality fluctuates stands in contrast with the ideal Kalman filter which, when used in conjunction with a linear model and when updating a properly initialized error covariance matrix of the model, statistically reduces the misfit between the recovery and the reference (see e.g.Fournier et al., 2010).Our assimilation scheme loses this property because of its imperfections: a timeindependent covariance matrix is used and a part of the solution length scales spectrum (for l and m above l p max ) is left untouched at analysis time and remains only determined by the dynamics.These uncorrected, dynamically determined scales tend to backreact on the other scales through the nonlinear couplings present in the dynamical equations, thus diverting the system trajectory away from the reference.The fluctuating recovery quality is then the result of a balance between the error decrease due to the introduction of new information and the error increase due to the scheme imperfections.Using a factor β (Eq.17) smaller than 1, we can increase the correction brought by the data, while decreasing the part of dynamically determined scales which is reinjected at analysis time.Comparing assimilation sequences performed with β = 0.75 and β = 1 in Fig. 4, we see that the recovery quality remains similar at the first few analyses and subsequently becomes better, with smaller fluctuations if β = 0.75.Note that this does not result from the removal of the energy contained in the dynamically determined scales (this represents a small part of the total energy), but reflects a better determination of the whole state vector.As we shall see later, there exists an optimum value for β (which is close to 0.75 for model 1), meaning that it is important to maintain dynamically determined scales, although not at the level they tend to freely reach.It should be noted that additional Fig. 4. Correlation and misfit coefficients after analysis, for the same series of twin experiments as in figure 3, for which β = 0.75 (black), and for a similar series where β = 1 (grey).Time is in units of τsec, as detailed in section 2.2, with a tentative rescaling in real time in parentheses.The assimilation is initialised and first analysed at time t = 0, the values of M and C reported at negative times are those obtained between the first reference sample and the model time average (which is the initialisation step estimate, second column of figure 3).The vertical dashed line is the time at which the fourth column of figure 3 has been obtained.28 Fig. 4. Correlation and misfit coefficients after analysis, for the same series of twin experiments as in Fig. 3, for which = 0.75 (black), and for a similar series where β = 1 (grey).Time is in units of τ sec , as detailed in Sect.2.2, with a tentative rescaling in real time in parentheses.The assimilation is initialised and first analysed at time t = 0, the values of M and C reported at negative times are those obtained between the first reference sample and the model time average (which is the initialisation step estimate, second column of Fig. 3).The vertical dashed line is the time at which the fourth column of Fig. 3 has been obtained.experiments with variable noise levels added to the observations (not shown here) show that the benefits of setting β = 0.75 are preserved when observations are imperfect.
To illustrate the issue of observed, estimated and dynamically determined scales, we next focus on energy spectra of the surface magnetic field (Fig. 5).Since the data are considered perfect, performing an analysis always results in a perfect match between the observed part of the system and the observations, as illustrated by the vanishing residual between the spectra of the analyses and the reference up to degree 13.The first analysis (red curve in Fig. 5) uses the correlations contained in P to additionally estimate the unobserved surface field coefficients between degree 14 and degree l p max = 15.Coefficients with degrees larger than l p max are not estimated by the first analysis.In contrast, the 15th analysis (green curve, performed immediately after the forecast presented in the fourth column of Fig. 3) benefits from a dynamical determination of these coefficients, and reduces the misfit to the reference by a factor of 2 for harmonic degrees between 15 and 30.
An important aspect of data assimilation is the evaluation of the forecast quality.When dealing with a real system, it is indeed impossible to evaluate the recovery quality of the internal structure as it was done in Fig. 4. It is however possible to use surface data in order to evaluate a-posteriori how well the system has predicted a given time evolution.This .Top: energy spectra of the magnetic field at the model surface for the reference solution (black), the first analysis (red) and the 15th analysis (green).The procedure through which these analyses are obtained is described in figure 3. Bottom: energy spectra of the differences between the reference, the first and the 15th analyses.Units as described in section 2.2, with their tentative rescaling in parentheses.
29 Fig. 5. Top: energy spectra of the magnetic field at the model surface for the reference solution (black), the first analysis (red) and the 15th analysis (green).The procedure through which these analyses are obtained is described in Fig. 3. Bottom: energy spectra of the differences between the reference, the first and the 15th analyses.Units as described in Sect.2.2, with their tentative rescaling in parentheses.
provides a combined assessment of the performance of the assimilation scheme, of its possible biases and of the suitability of the model for describing the real system (dealing with synthetic data, our study does not cover this last aspect).The quality of the i-th forecast x f i can only be assessed from the standpoint of the observer and its limited access to the system.It is thus evaluated on the observed part of the system only, using the instantaneous innovation (or instantaneous forecast error) vector d i = y i − Hx f i (which indeed represents the difference between the observable part of the forecast and the data), or more precisely through its norm Here our norm definition is adapted such that each harmonic coefficient is multiplied with l(l + 1)/r o prior to the evaluation of the norm, such that the result is a rms value of the radial magnetic field at the outer boundary (or its time derivative).One important property of the innovation vector d i is that its statistical expected value should be zero for an unbiased assimilation scheme (e.g.Talagrand, 2003).Computing the cumulative mean innovation provides a quantitative way to test this prediction.Figure 6 presents both quantities for various assimilation sequences.We recall that in our case, the observed part of the system  3 and other sequences with same access to observable quantities and assimilation parameters but varying forecast horizon (which is the same as the time lag between analyses), from 0.01τsec to 0.2τsec (solid black lines, each circle represents an analysis, forecast horizon increases vertically).The innovation is computed only with the magnetic field (not the secular variation).Also plotted are an assimilation sequence with β = 1 (solid grey line, forecast horizon 0.1τsec) and the root-mean squared radial magnetic field amplitude at the outer boundary for the reference sequence (dashed line).Units as in previous figures.30 Fig. 6.Instantaneous innovation (or instantaneous forecast error) d i and cumulative mean innovation d k as a function of assimilation time, for the sequence presented in Fig. 3 and other sequences with same access to observable quantities and assimilation parameters but varying forecast horizon (which is the same as the time lag between analyses), from 0.01τ sec to 0.2τ sec (solid black lines, each circle represents an analysis, forecast horizon increases vertically).The innovation is computed only with the magnetic field (not the secular variation).Also plotted are an assimilation sequence with β = 1 (solid grey line, forecast horizon 0.1τ sec ) and the root-mean squared radial magnetic field amplitude at the outer boundary for the reference sequence (dashed line).Units as in previous figures.
comprises the poloidal magnetic field up to degree and order 13 and its secular variation up to a variable degree and order (13 in the case of Figs. 3 to 6).However, to facilitate comparison between sequences where secular variation is assimilated to a variable degree (see Fig. 7 below), we will evaluate the forecast quality only on the observed magnetic field poloidal coefficients up to degree 13.As expected from the discussion of Fig. 4, d i does not decrease with the number of assimilation cycles, but oscillates about a slowly increasing baseline.This long-term trend does not mean that the assimilation gets worse over time, but is simply the consequence of an increase of the reference sequence magnetic energy through time (dashed line in Fig. 6).In line with Fig. 4, d i is significantly reduced when β is decreased from 1 to 0.75.Regardless of the forecast horizon (which is the same as the time lag between analyses), d k decreases sharply after one system overturn time (which is about 0.3 in units of τ sec , a value which is quite independent of the chosen model, Christensen and Tilgner, 2004;Lhuillier et al., 2011a).At very long assimilation times, d k ceases to decrease, revealing the existence of a forecast bias in our scheme.The bias significantly decreases when β is decreased from 1 to 0.75.This shows that it is connected to the influence of the uncorrected variables of the state vector on the corrected variables.Although it would certainly be desirable to implement a bias removal strategy in the assimilation scheme, the present bias is not likely to be a limiting factor in practical applications, its level being typically one order of magnitude lower than the intrinsic forecast error introduced by the assimilation.Furthermore, long assimilation times (here 8500 yr) are needed to reveal the presence of this bias.Our synthetic experiments, where high-quality data are assimilated in such long sequences, obviously do not represent a practical situation, given the presently available record of Earth's magnetism on centennial to millenial timescales (e.g.Hulot et al., 2010a).
To really become meaningful, the absolute forecast quality should be compared to that obtained with other prediction strategies, the simplest of which is the "no-cast", or use of the present magnetic field as a forecast for the future.Another still simple strategy is the linear extrapolation of the existing surface magnetic field, making use of the secular variation data.The results of series of no-casts, linear extrapolations and assimilations are reported in Fig. 7, for various forecasts horizons never exceeding the system e-folding time 0.6τ sec .The average forecast error is defined as where n a is the total number of assimilation cycles.The average forecast error d should not be confused with the cumulative mean innovation, d k .It can be first noted that d approximately follows a power-law of the forecast horizon, for all the prediction strategies which we have studied.The left panel of Fig. 7 shows that regardless of the forecast horizon, an assimilation of surface magnetic field coefficients always makes a better forecast of these coefficients (by about a factor of 2) than a no-cast performed with the same amount of data.This is also true (middle panel) if the linear forecast is refined using the secular variation data up to degree 8 (the corresponding assimilation then also uses this additional data).Finally, if secular variation coefficients up to degree 13 are used, linear extrapolations will perform better than the assimilation for forecast horizons shorter than 0.03τ sec (to confirm the robustness of this result, an additional assimilation sequence has been carried out with spacing 0.02τ sec ).Here, the errors introduced in the retrieval of the internal structure are too large to provide a forecast which can match the high accuracy of the linear extrapolation.In this latter case, the typical forecast bias of our scheme (5 10 −4 mT for an 0.01 τ sec horizon) is less than a third of the gap separating the assimilation and the linear forecast (about 1.7 10 −3 mT).Removing the forecast bias would thus not entirely bring the assimilation in line with the linear forecast.Fig. 7. Average forecast error d of surface poloidal magnetic field coefficients up to degree and order 13, for various prediction strategies and variable forecast horizon.From left to right: The no-cast (poloidal magnetic field coefficients up to degree 13 at a given time are used to forecast at a later time), the linear forecasts (poloidal magnetic field coefficients (B) up to degree 13 and secular variation (SV) up to degree 8 or 13 are used for a linear extrapolation) are compared with the assimilation using the exact same amount of data.For each panel, the average forecast error is computed on an equal number of cycles for the assimilation and the linear prediction.All sequences have β = 0.75 as in figure 3. Units as in previous figures.The averaging time interval for all sequences presented in this figure is from 0 to 6τ sec .
31 Fig. 7. Average forecast error d of surface poloidal magnetic field coefficients up to degree and order 13, for various prediction strategies and variable forecast horizon.From left to right: the no-cast (poloidal magnetic field coefficients up to degree 13 at a given time are used to forecast at a later time), the linear forecasts (poloidal magnetic field coefficients (B) up to degree 13 and secular variation (SV) up to degree 8 or 13 are used for a linear extrapolation) are compared with the assimilation using the exact same amount of data.For each panel, the average forecast error is computed on an equal number of cycles for the assimilation and the linear prediction.All sequences have β = 0.75 as in Fig. 3. Units as in previous figures.The averaging time interval for all sequences presented in this figure is from 0 to 6τ sec .
In a final series of experiments (Fig. 8), we use the forecast quality as an indicator to perform a number of checks on our data assimilation framework.The first of these is an evaluation of the adequate number n of free run snapshots needed for building a robust estimate of P, as well as their spacing.Figure 8a shows that more than 1000 samples with a spacing at least equal to the e-folding time of the system are needed for a reasonable determination of the covariance matrix.Conversely, using a too small set of samples leads to forecast errors which may exceed the error made with a no-cast strategy.The next test evaluates the impact of the factor β which controls the amount of correction brought by the data at analysis time. Figure 8b shows that an optimum in the forecast quality is reached for β = 0.75, which underlines again the potentially deleterious effect of re-injecting all the dynamically determined scales into the system after analysis.Re-injecting some of these dynamically determined scales is however beneficial to the forecast quality, as seen from the regular decrease of the average error from β = 0 to β = 0.75.Finally, we have performed experiments (Fig. 8c) where the correction at analysis time has been turned off for one or several fields.Obviously, correcting only the magnetic field harmonic potentials results in a better forecast than not assimilating anything, but the forecast quality will outperform that of a no-cast strategy only if the flow and buoyancy harmonic potentials are also corrected.Note that the best improvement comes from updating the buoyancy potential, presumably because the flow driven by the buoyancy anomaly is then also well estimated.This last test emphasizes the benefit of using multivariate statistics for the estimation, and also illustrates the difficulties encountered by monovariate, modelled covariances (Kuang et al., 2009) in forecasting the field with better accuracy than that of a linear extrapolation.

Synthetic inversion tests, model 2
The satisfying results obtained with model 1 can be understood through the strong linear couplings existing between the observed and unobserved part of the system.To evaluate the impact of stronger nonlinearities, we now turn to synthetic experiments performed with model 2. Following the prescriptions obtained through the analysis of model 1, the co-variance matrix for model 2 was built using n = 978 samples in the free run, with a spacing between samples of 0.125τ sec , which is about three times longer than the system e-folding time τ e = 0.04τ sec .The duration of the free run is 0.65 magnetic diffusion times or 122τ sec , amounting to 61 000 yr if τ sec = 500 yr.It should be noted from Fig. 8 that if we were to use n = 978 also in model 1, this would not substantially degrade the quality of the inversions.The differences present in the following results should thus not be ascribed to the size of the error space (which is the rank of the covariance matrix and is equal to n).The matrix coefficients were computed up to degree and order l P max = 30.All assimilation experiments have been performed using β = 0.75.
Figure 9 presents an equivalent to Fig. 3 for model 2, with an additional line presenting the observable part of the surface radial magnetic field, which is now clearly different from the complete magnetic field (for model 1 most of the 8. Average forecast error d for series of experiments where the surface magnetic field is assimilated up to ee 13 (the secular variation is not assimilated here), and the forecast horizon is 0.1τ sec .a.: the factor β is 0.75 and the number n of samples used to build P is varied by decimating an ensemble of 8155 snapshots ned during the model free run, the spacing between each snapshot being 0.13τ sec .Two vertical dashed indicate the values of n for which the spacing between snapshots is equal to the e-folding and overturn s of the system.b.The number n of samples is set to 4098 and the factor β is varied from 0 (the system et to its time average state prior to each analysis) to 1 (the analysis is performed on the forecast resulting the previous time integration).c.The correction at analysis time has been turned off for some, or all fields 4098, β = 0.75).Same time averaging interval as in previous figure .32 Fig. 8. Average forecast error d for series of experiments where the surface magnetic field is assimilated up to degree 13 (the secular variation is not assimilated here), and the forecast horizon is 0.1τ sec .(a) the factor β is set to 0.75 and the number n of samples used to build P is varied by decimating an ensemble of 8155 snapshots obtained during the model free run, the spacing between each snapshot being 0.13τ sec .Two vertical dashed lines indicate the values of n for which the spacing between snapshots is equal to the e-folding and overturn times of the system.(b) The number n of samples is set to 4098 and the factor β is varied from 0 (the system is reset to its time average state prior to each analysis) to 1 (the analysis is performed on the forecast resulting from the previous time integration).(c) The correction at analysis time has been turned off for some, or all fields (n = 4098, β = 0.75).Same time averaging interval as in previous figure .surface magnetic field was observable).With the exception of this additional line, all other lines have the same color bars as the corresponding lines of Fig. 3, in order to highlight the relevance (and possible shortcomings) of the scaling procedure which we have adopted (see Sect. 2.2).From an order-of-magnitude standpoint, the rescaling is indeed satisfactory, but there are variations within the order of magnitude which the scaling theory fails to describe.For instance, model 2 has a slightly larger internal magnetic energy than model 1 (third line of Fig. 9) but less magnetic flux escapes at the surface (second line).The large-scale velocity and buoyancy anomaly fields have roughly the same amplitude as in model 1, but more powerful small-scales, located mostly within the strong plumes emerging from the inner boundary.
When compared to model 1, the static stochastic inversion of one data sample (second column of Fig. 9) captures less of the internal structure of the dynamo.The system is indeed less observable, in the sense that a reduced state vector which is now ten times larger (see Table 1) needs to be constrained by the same amount of observation.Furthermore, the stronger nonlinearities present in model 2 tend to even out correlation peaks resulting from the linear couplings in the covariance matrix of model 2, resulting in less accurate estimations.Still, at the exception of the deep magnetic field (third line), which is too small-scaled to be well constrained by the surface observations (first line), a number of gross details of the internal solution are retrieved.The first estimate is underpowered as it was for model 1.If more data samples have been previously assimilated through integration of the time evolution scheme (third column), the estimate reaches the same power as the reference.After a time of 0.7τ sec (corresponding to about two overturns), the deep fields are only qualitatively recovered in a morphological sense.The gross details of the convection flow and thermal plumes are roughly into place.However, small scale details created by the nonlinear dynamics, which are especially prominent in the deep magnetic field map, are clearly not constrained by the surface observations and the scheme.The frozen covariance matrix which we have obtained for model 2 lacks the ability to reliably estimate small-scale details from the large-scale observations.Indeed we have seen in Sect.3.1 that couplings between different harmonic orders are nonexistent, and that couplings between harmonic degrees are prominent under the condition that the solution is rather well controlled by the linear part of the equations (including the Coriolis force).The strong nonlinear dynamics is thus left free to populate these small-scales in a rather unconstrained way.This deleterious effect is also enhanced by the fact that the size of the uncorrected part of the state vector is ten times larger in model 2 when compared to model 1 (Table 1).
Although the deep structure of model 2 is not well retrieved, the assimilation seems to make (at least visually, see the first line of Fig. 9) decent forecasts of the system evolution.On Fig. 10 we again compare the quality of these forecasts with that of other prediction strategies.In contrast to Fig. 9. Results of a twin experiment with model 2, where the surface radial magnetic field (first line) and its secular variation (not shown), both given up to degree and order 13, are used to retrieve the internal structure (second to fifth line).The first column is the reference.Similarly to Fig. 3, column 2 is the estimate if only the initialisation and the first analysis steps have been performed, and column 3 is the estimate obtained after performing a full data assimilation sequence where 7 data samples spaced 0.1τ sec time units each are injected in the model.Note that for this last column, the state of the system is represented after the forecast (not after the analysis).For this run we used β = 0.75.For comparison purposes, units and color bars for lines 2 to 5 are the same as in Fig. 3. Fig. 7, the error for the various strategies no longer follow a power-law of the forecast horizon as this horizon approaches the system e-folding time τ e = 0.04τ sec .Indeed it is expected (Hulot et al., 2010b;Lhuillier et al., 2011a) that the error starts to become of macroscopic (same order of magnitude as the reference) amplitude beyond this point, as confirmed by our experiments.The left panel of Fig. 10 shows that the assimilation never performs significantly better than the nocast.Moreover, when secular variation data up to degree 13 are also used, the linear extrapolations are shown to be much better than the assimilation forecasts, for all horizons below the e-folding time where both strategies again yield the same error.Given that Re = 858 for model 2, the secular variation at the outer boundary is now mostly controlled by flow advection underneath the outer boundary.The flow responsible for the secular variation is nonlinearly coupled to the observations, and thus not well grasped by our linear estimation technique.This explains the worse results obtained by the assimilation.
Here, our assimilation scheme faces a problem of relevance: it is less efficient than the linear forecast for horizons below the system e-folding time, and for longer horizons, prediction strategies can only be envisioned if the analysis done before the time integration has very low errors in the determination of the internal structure.As, for the Earth, we may not be able to observe enough of the system to perform such low error analyses, this reduces the prospect of implementing a reasonably complex and realistic geomagnetic data assimilation framework in an operational context.However, assuming that the errors of the assimilation and linear extrapolation strategies are decorrelated, and if their statistics are known, the predictions from both strategies can be combined in order to provide a third prediction which is better than the two original ones.Such an additivity of accuracies is a concept at the root of the linear estimation theory.Here again, a Kalman filter is the adequate tool to perform this task.If we assume that we can estimate the error covariance matrix R of linear extrapolations, and if, for a given forecast horizon, we have both the linear extrapolation forecast x f l and the assimilation forecast x f , then the combined forecast x f c is such that In a more readable presentation which involves the observable parts y f l = Hx f l , y f = Hx f and y f c = Hx f c , and the model covariance matrix reduced to the observable part O = HPH , the combined forecast writes: The combined forecast is thus simply an optimal interpolation of the two other forecasts.Figure 11 illustrates this concept on a forecast horizon 0.05τ sec , roughly equal to the system e-folding time.Diagonal coefficients of the matrix R (Fig. 11a) are estimated following the same procedure as in Sect.3.1, using 100 linear forecasts performed in a preliminary run of the model.Cross-correlations in the error statistics are neglected (non diagonal coefficients are set to 0).Errors are naturally statistically centered, and according to the formalism outlined in Sect.2.3, they are normalised with the corresponding diagonal variances used for P.An error equal  11. a.: empirically-estimated diagonal coefficients Rlm,lm of the error covariance matrix R for linear forecasts of model 2 at the horizon 0.05τsec.Coefficients are normalised with the corresponding variances of the free run, and are ordered according to the same one-parameter lm ordering scheme as in figure 2. Non-diagonal coefficients are set to 0. b.: instantaneous innovation (or instantaneous forecast error) di (circles) for the assimilation sequence shown in figure 9 with forecast horizon 0.05τsec.Also represented are the results of a linear extrapolation using the exact same amount of data, i.e. the surface magnetic field and secular variation coefficients up to degree 13 (squares), and the results of a combined forecast (diamonds) as defined in equation ( 25), using the empirically-estimated matrix R.   9 with forecast horizon 0.05τ sec .Also represented are the results of a linear extrapolation using the exact same amount of data, i.e. the surface magnetic field and secular variation coefficients up to degree 13 (squares), and the results of a combined forecast (diamonds) as defined in Eq. ( 25), using the empirically-estimated matrix R.
to one thus means that the linearly forecast harmonic coefficient varies as much as what was observed for the coefficient itself during the model free run.From Fig. 11.a it appears that the linear forecast performs better on large scales than on small scales.The optimal interpolation outlined above has the effect to mitigate this error by injecting more of the assimilation forecast when the linear forecast error is large.The overall root-mean-squared forecast error thus decreases (Fig. 11b).The quality of the combined forecast outperforms that of the linear extrapolation by about 10 %, and that of the assimilation by about 40 %.More importantly, the combined forecast is shown to be always better than the best of the two strategies.This clearly underlines the interest of having several independent prediction strategies at hand when attempting to perform forecasts, especially if each strategy is far from perfect.

Discussion
In this study, we have explored how a linear estimation technique, the stochastic inverse (or Kalman filter), together with www.nonlin-processes-geophys.net/18/657/2011/ Nonlin.Processes Geophys., 18, 657-674, 2011 a sequential assimilation scheme, can provide inferences on the internal structure of numerical dynamos from surface observations only.The way the stochastic inverse handles the non-uniqueness problem is through the use of prior information, obtained by directly computing the multivariate statistics of the numerical models from a suitably large set of snapshots with enough spacing between each other.Strong correlations were found between the surface observations and the internal fields, which arise from the linear part of the dynamical equations and the leading influence of the Coriolis in the dynamics.This led to excellent synthetic recovery results with a weakly nonlinear model (model 1).Stronger nonlinearities (model 2) however defeat the linear estimation technique to some extent, even though nonlinearity itself can be used with some success in estimating the part of the underlying field spectrum which is neither observed nor correlated to surface observations (Fig. 5).This could provide a possible path towards improving spatial resolution problems of the existing inversion strategies mentioned in the introduction.
With kinematic Reynolds numbers of order 10 9 and magnetic Reynolds number of order 10 3 (e.g.Christensen and Tilgner, 2004;Aubert et al., 2009) it is certain that nonlinearities are important in Earth's core dynamics.Our results are encouraging but highlight the need for estimation techniques specifically designed to handle a large amount of nonlinearity.When dealing with such magnetic Reynolds number (as is the case in model 2), advection by the flow underneath the outer boundary is responsible of most of the secular variation in the observed range.The linear estimation of core flow performed in this study is clearly insufficient and could be replaced by a classical inversion of the radial induction equation (see for instance Finlay et al., 2010a), which handles the nonlinear nature of core flow advection correctly.Another important point is to implement a time evolution algorithm for the model covariance matrix, in order to have an instantaneous matrix which, for each analysis time, should be more adapted to the evaluation of nonlinear correlations than the generic, frozen covariance matrix in which, as we have seen, only the correlations subsequent to the linear part of the system arise.A promising method is for instance the ensemble Kalman filter (Evensen, 1994), where the model error statistics are evolved through the use of an ensemble of models states created around the actual model trajectory with the help of a Monte-Carlo method.Such a method is however much more costly in terms of computer power than the method which we have presented here, which had requirements on disk space (400 GB for model 2) and random access memory (30 GB) only at the time of the computation of the frozen model covariance matrix.
Geomagnetic data assimilation is however still in its infancy, and more acute problems need to be solved before geomagnetism can integrate the advanced data assimilation techniques routinely used in atmospheric and ocean sciences.One of these issues is the choice of a numerical model for performing an operational inversion for the internal structure of the geodynamo.As discussed in the introduction, the Earth's core has an extraordinary disparity in the diffusivities of the involved fields, as measured for instance by the magnetic Prandtl number P m ≈ 10 −6 of liquid iron.Numerical models are not able to handle such diffusivity contrasts and operate at Prandtl numbers close to unity.Even if we can reach an acceptable level of magnetic turbulence (Rm = O(1000)), current computer power limitations make it impossible to reach a level of hydrodynamic turbulence which approaches that of the Earth's core.As a result, even if the magnetic induction phenomenon is correctly simulated, we should be cautious about the relevance of the large-scale velocity output of our simulations and inversions.In the absence of a clear path towards improving this situation, our understanding of the physical grounds underlying the morphological semblance between numerical dynamos and the geomagnetic field (Christensen et al., 2010) suggests that we should select a model which has the lowest possible ratio between the length of the day and the magnetic dissipation time scale (the magnetic Ekman number) and a ratio between the magnetic advection and diffusion time scale approaching 1000 (the magnetic Reynolds number).In that sense model 2 should be much better suited than model 1 for the task of inverting real geomagnetic data, but here a trade-off should be made between the physical relevance of the model and the ability to linearly estimate hidden state variables, which, as we have seen, is precisely hampered by the fact that the magnetic Reynolds number is large.
It is interesting to briefly review how geomagnetic data assimilation strategies currently in development deal with the issue of this improper rendering of the real physics.Kuang et al. (2010) propose that these modelling errors evolve on a time scale longer than that at which observation data can be made.In that case they can be mitigated by performing two closely spaced assimilation sequences.The method is promising and already resulted in a secular variation model for the latest generation of the IGRF (Finlay et al., 2010b).It is however limited to surface field forecasting activities.Progress in determining the internal dynamical structure could be achieved through combining data assimilation with asymptotic assumptions on core dynamics, for instance the quasi-geostrophic assumption (Canet et al., 2009), or building Taylor states (Livermore et al., 2010) compatible with surface observations.In any case, using several internal structure modelling approaches within their spatial and temporal range of validity could help overcome the intrinsic limitations of each strategy.For instance, a quasi-geostrophic framework is more appropriate for rapid (decadal) flow variations (Jault, 2008) while a three-dimensional numerical dynamo is suitable for describing long-term, ageostrophic flows such as thermal winds (Aubert et al., 2010).
With the help of the scaling procedure presented in Sect.2.2, it is useful to recast the time axis of our models to the dimensional world, using the secular variation timescale τ sec ≈ 500 yr, in order to set some bounds on what could be expected from an operational implementation of our method for short-term geomagnetic data assimilation.Given the macroscopic errors obtained when inverting for the internal structure, our scheme cannot provide reasonable forecasts at horizons longer than the system e-folding time.For the Earth's core, this time was evaluated (Hulot et al., 2010b;Lhuillier et al., 2011a) at about τ e = 0.04τ sec ≈ 25 yr (similar to model 2).At very short horizons such as 0.01τ sec = 5 yr, we have also seen that linear extrapolations of the magnetic field evolution usually perform better than the assimilation.Still, as the two prediction strategies have complementary strengths and weaknesses, we have applied the principle of "additivity of accuracies" in order to show how the two prediction strategies can re-inforce each other to provide a better third forecast, provided the error statistics of each strategy is known (or can be estimated).Here again it should be stressed that the current generation of numerical dynamo models lacks the short timescale dynamics of Earth's core, such as the 6-yr oscillation studied by Gillet et al. (2010).The extent to which data assimilation based on these numerical models could improve our knowledge of this fast dynamics needs to be clarified.Still, it should not be forgotten that beyond the classical exercise of forecasting the surface magnetic field evolution, data assimilation could bring estimates and forecasts for other hidden and internal geodynamical properties for which we currently do not have any other estimation and forecasting strategy at hand.
Evaluating the typical forecast error of our assimilation scheme at the surface of the Earth, we find errors on the order of a microtesla for a 50 yr forecast, and of a few tenths of nanoteslas for a 5 yr forecast.These should be put in perspective with the typical observation errors present in geomagnetic field modelling, which amount for instance to tens to hundreds of nanoteslas in the historical period (Bloxham and Jackson, 1989) and a few nanoteslas in the satellite observation period (Olsen et al., 2010).An assimilation performed on real data will thus need to include a proper estimation of the observation error for the assimilated quantities.It might then be impractical to work with the harmonic coefficients supplied by geomagnetic field models such as gufm1 (Jackson et al., 2000), as the evaluation of the observation error covariance matrix is intricate.To tackle this issue, future data assimilation frameworks should be able to directly integrate pointwise direct or indirect measurements, along with their corresponding observation errors, or progress should be encouraged towards better evaluating the error associated with each harmonic coefficient in geomagnetic field models, as done most recently in the archeomagnetic context by Korte and Constable (2011).

Fig. 2 .Fig. 2 .
Fig.2.Representation of the coefficients of the covariance matrix P involved in the determination of the stochastic inverse matrix K.At five different radial levels (vertical axis), the colored squares map the modulus of the correlation coefficients between the surface poloidal magnetic field harmonic coefficients (first horizontal axis) and the poloidal velocity field coefficients (second horizontal axis).Harmonic coefficients are ordered according to a one-dimensional scheme where all admissible values of l are grouped together for each given value of m.These correlation coefficients are computed from a free model run (here from model 1) where n = 4098 instantaneous state vectors spaced by 0.3τ sec (half an e-folding time) each are extracted.The coefficients are computed up to degree and order l P max = 15, which corresponds to the one-dimensional parameter lm = 120.Two vertical tracks are drawn, representing the evolution of the correlation coefficients with depth for m = 4, between harmonics l obs = 5 of the observed field and l deep = 4 of the deep field (track 1), and between harmonics l obs = 5 and l deep = 12 (track 2).

Fig. 3 .
Fig.3.Results of a twin experiment with model 1, where the surface radial magnetic field (first line) and its secular variation (not shown), both given up to degree and order 13, are used to retrieve the internal structure (second to fourth line).Column 1 is the state of the reference at time t = 1.3τ sec .Column 2 is the assimilation state if only the initialisation step has been performed at time t (it thus represents the time average of the simulation, see Table2).Column 3 is the assimilation state if only the initialisation step and the first analysis step as been performed at time t.Column 4 is the assimilation state if it has been initialised at time t = 0, and if 14 analyses and forecast steps have been subsequently performed with a time lag 0.1τ sec (about 50 yr in Earth time) until time t = 1.3τ sec (note that this column then represents the forecast, not the analysis).For this run we used β = 0.75.All fields are presented in the units proposed in Sect.2.2, with their tentative rescaling to Earth's core values in parentheses.
Fig.5.Top: energy spectra of the magnetic field at the model surface for the reference solution (black), the first analysis (red) and the 15th analysis (green).The procedure through which these analyses are obtained is described in figure3.Bottom: energy spectra of the differences between the reference, the first and the 15th analyses.Units as described in section 2.2, with their tentative rescaling in parentheses.

Fig. 6 .
Fig.6.Instantaneous innovation (or instantaneous forecast error) di and cumulative mean innovation dk as a function of assimilation time, for the sequence presented in figure3and other sequences with same access to observable quantities and assimilation parameters but varying forecast horizon (which is the same as the time lag between analyses), from 0.01τsec to 0.2τsec (solid black lines, each circle represents an analysis, forecast horizon increases vertically).The innovation is computed only with the magnetic field (not the secular variation).Also plotted are an assimilation sequence with β = 1 (solid grey line, forecast horizon 0.1τsec) and the root-mean squared radial magnetic field amplitude at the outer boundary for the reference sequence (dashed line).Units as in previous figures.

Fig. 10 .Fig. 10 .
Fig. 10.Average forecast error d of surface poloidal magnetic field coefficients, for model 2, various prediction strategies and various forecasts horizons, computed as in figure 7.
Fig.11.a.: empirically-estimated diagonal coefficients Rlm,lm of the error covariance matrix R for linear forecasts of model 2 at the horizon 0.05τsec.Coefficients are normalised with the corresponding variances of the free run, and are ordered according to the same one-parameter lm ordering scheme as in figure2.Non-diagonal coefficients are set to 0. b.: instantaneous innovation (or instantaneous forecast error) di (circles) for the assimilation sequence shown in figure9with forecast horizon 0.05τsec.Also represented are the results of a linear extrapolation using the exact same amount of data, i.e. the surface magnetic field and secular variation coefficients up to degree 13 (squares), and the results of a combined forecast (diamonds) as defined in equation (25), using the empirically-estimated matrix R. 35

Fig. 11 .
Fig. 11.(a): empirically-estimated diagonal coefficients R lm,lm of the error covariance matrix R for linear forecasts of model 2 at the horizon 0.05τ sec .Coefficients are normalised with the corresponding variances of the free run, and are ordered according to the same one-parameter lm ordering scheme as in Fig. 2. Non-diagonal coefficients are set to 0. (b): instantaneous innovation (or instantaneous forecast error) d i (circles) for the assimilation sequence shown in Fig.9with forecast horizon 0.05τ sec .Also represented are the results of a linear extrapolation using the exact same amount of data, i.e. the surface magnetic field and secular variation coefficients up to degree 13 (squares), and the results of a combined forecast (diamonds) as defined in Eq. (25), using the empirically-estimated matrix R.

Table 1 .
Properties of the numerical models used for the study.First row: input parameters (see main text for definitions).Second row: output parameters.Earth's core values are estimated in