Dynamic responses of the Earth’s outer core to assimilation of observed geomagnetic secular variation

Assimilation of surface geomagnetic observations and geodynamo models has advanced very quickly in recent years. However, compared to advanced data assimilation systems in meteorology, geomagnetic data assimilation (GDAS) is still in an early stage. Among many challenges ranging from data to models is the disparity between the short observation records and the long time scales of the core dynamics. To better utilize available observational information, we have made an effort in this study to directly assimilate the Gauss coefficients of both the core field and its secular variation (SV) obtained via global geomagnetic field modeling, aiming at understanding the dynamical responses of the core fluid to these additional observational constraints. Our studies show that the SV assimilation helps significantly to shorten the dynamo model spin-up process. The flow beneath the core-mantle boundary (CMB) responds significantly to the observed field and its SV. The strongest responses occur in the relatively small scale flow (of the degrees L≈30 in spherical harmonic expansions). This part of the flow includes the axisymmetric toroidal flow (of order m=0) and non-axisymmetric poloidal flow with m≥5. These responses can be used to better understand the core flow and, in particular, to improve accuracies of predicting geomagnetic variability in the future.


Introduction
Geomagnetic field observed at the Earth's surface varies significantly in time: its temporal scales range from minutes to geological time scales.Though it was first noticed by mankind over 5000 years ago (Roberts, 1992), and its origin was sought as early as 800 years ago (Dibner Library 1980), the modern theory that the geomagnetic field is generated and maintained by convective flow in the Earth's outer core (geodynamo) was originated from the seminal work of Larmor (1919).Successful numerical simulation of the geodynamo was first carried out by Glatzmaier and Roberts (1995), and then followed by Kageyama and Sato (1997), and by Kuang and Bloxham (1997).Christensen et al (2010) provided a comprehensive summary of numerical geodynamo solutions and their relevances to geomagnetic observations.Assimilation of geomagnetic observations with numerical geodynamo models started less than a decade ago.Sun et al (2007) and Fournier et al (2007) used simplified magnetohydrodynamic (MHD) systems and synthetic data tested the applicability of assimilation of sparse magnetic data.Liu et al (2007) first used observation system simulation experiments (OSSEs) with a full dynamo and demonstrated clearly that one could use assimilation of magnetic field at the surface to estimate the dynamo state deep in the fluid core.Kuang et al (2008) published the first working geomagnetic data assimilation system MoSST − DAS in which the Gauss coefficients of various geomagnetic and paleomagnetic field models are assimilated with their MoSST geodynamo model (Kuang and Chao, 2003;Jiang and Kuang, 2008) for estimation of the core state and prediction of geomagnetic field.Kuang et al (2009) then used this assimilation system and 100 years of the Gauss coefficients from GUFM1 (Jackson et al 2000) and CM4 (Sabaka et al 2004) to understand the responses of the core state to surface geomagnetic observations, and their implications to core state estimation and SV prediction.We refer the reader to Fournier et al (2010) for a comprehensive review of the data assimilation algorithms for geomagnetic data assimilation (GDAS) and some of the early results.
Rapid advances have occurred in multiple facets of GDAS.Several independent assimilation systems have been developed to understand better the core dynamical state.For example, Aubert and Fournier (2011), and Fournier et al (2011, 2013) carried out OSSEs with synthetic observations and numerical dynamo models to examine possibilities of core state determination.Aubert (2013Aubert ( , 2014) ) investigated possibilities of inverting core state properties using the observed field and SV.In addition to the sequential data assimilation systems mentioned above, there are also efforts in developing GDAS systems based on variational data assimilation techniques.For example, Li et al (2011Li et al ( , 2014) ) have been continuing their effort on a new combined system of forward and adjoint systems.Encompassed application is the contributions of assimilation results to international geomagnetic reference field (IGRF) (Kuang et al, 2010), and efforts to determine field model error statistics (Gillet et al, 2013).
Despite these advances, GDAS is still in an early stage similar to that of early numerical weather prediction (NWP) (for a more comprehensive review, see, e.g.Kalnay 2003).Many important questions are still to be fully answered, such as comprehensive assessment of numerical dynamo system biases, observation and core state covariances and error statistics, and the dynamic responses of dynamo state to the observed geomagnetic field.The latter is of in particular importance to the spin-up processes of the numerical models which, in turn, determine how fast and how close the numerical solutions can be pulled to the true state of the core.
Concerns on the spin-up of the numerical models can be examined from the time scales of the observed field and of the numerical models.Global field model results from the past 400 years of geomagnetic data (e.g.Jackson et al, 2000;Sabaka et al, 2004Sabaka et al, , 2015;;Olsen et al, 2006Olsen et al, , 2014) ) show that the typical time scales τ l of the degree l components (Stacy 1992;Hulot and Le Mouël 1994;Olsen et al 2006) (1) varies from over 1000 years for the dipole (l = 1) to less than 100 years for higher degrees (see Figure 1).In ( 1), (g m l , h m l ) are the Gauss coefficients of the field, and ( ġm l , ḣm l ) are their first order time derivatives, i.e. the Gauss coefficients of the SV.Currently, the longest record for low degree (l ≤ 5) field coefficients is from the paleo/archeo magnetic data (e.g.Korte et al, 2011;Nilsson et al, 2014).The high quality coefficients for up to degree l ≤ 8 could be obtained from historical and observatory data (Jackson et al, 2000).Very high quality coefficients for degrees l ≤ 13 are obtained in the past 50 years with satellite magnetic data (Sabaka et al, 2004(Sabaka et al, , 2015;;Olsen et al, 2006Olsen et al, , 2014)).In summary, the data record is no more than 10 times of the typical time scales of the geomagnetic field.This brings the very concern on whether the observational record is sufficient to spin up numerical dynamo models.The model spin-up also has direct consequence on estimation of the core state.
How could we improve geomagnetic data assimilation systems within the observational limit?There are several areas for improvements.For example, improvements in global geomagnetic field modeling are needed since the Gauss coefficients from various field models have been used in most of the previous GDAS studies.Currently there are many field models covering different epochs (e.g.Jackson et al, 2000;Korte et al, 2011;Gillet et al 2013;Olsen et al, 2014;Sabaka et al, 2015).A unified field model covering the longest possible period could certainly reconcile differences in these models, and thus help greatly GDAS systems.There is an ongoing effort on constructing a unified global field model of the past millennium (private communication with Korte).The field model error statistical information of such unified field models, such as those in Gillet et al (2013), is also necessary for GDAS.
Improvement in the assimilation algorithms could also help data utilization.Some efforts were made by Kuang et al (2010) in which a subset of the Gauss coefficients (of lower degrees) with much longer records are assimilated first to speed up the model, followed by assimilating those of higher degrees for the past 100 years.Tangborn and Kuang (2015) showed, via a set of experiments, that such assimilation methodology can have positive impact on core state, and improve accuracies of predicting the subset of the Gauss coefficients not assimilated.Another example is employment of ensemble Kalman Filtering (EnKF) approach (Evensen, 1994).Fournier et al (2011Fournier et al ( , 2013) ) used OSSEs to show the potential to speed up the transfer of information from geomagnetic data to the core state.But such speedy transfer depends on model errors (that are in general very large due to limitations of numerical dynamo models) not considered in their studies.It should also point out that GDAS is computationally very expensive.Such expense needs to be considered in the algorithm improvement.
Another improvement is on exploiting and utilizing further geodynamic information embedded in surface geomagnetic measurements.An immediate candidate for such exploitation is the geomagnetic secular variation (SV), described by the first order time derivative ( ġm l , ḣm l ) of the Gauss coefficients since, as we will describe in the next section, they provide additional constraints on the core flow beneath the CMB, and on the radial variation of the magnetic field.The former is not new, as there is a long history of, started from Roberts and Scott (1965), core flow inversion from observed SV at the Earth's surface via the "frozen-flux" approximation (in which the Ohmic dissipation beneath the CMB is ignored).However, this approximation comes with the price: the core flow cannot be uniquely inverted (e.g.Roberts and Scott 1965;Backus 1968).Thus additional constraints on core flow properties are necessary in such core flow inversion studies (for more complete reviews, please read, e.g., Holme, 2007;Kuang and Tangborn, 2011).If the Ohmic dissipation is retained (no "frozen flux" approximation), then the observed SV imposes the constraints on the radial variation of the field in the core, as the latter is part of the magnetic induction.Since both field advection and Ohmic dissipation are included in geodynamo modeling, both kinds of constraints can be examined in MoSST − DAS or any other GDAS system without mathematical difficulties.
Therefore, a natural expansion of data utilization in GDAS is to assimilate both the field and its SV, so that the embedded geodynamic constraints can be used to make more optimal analysis, thus speeding up the transport of information from the surface geomagnetic observations to the dynamical state in the outer core.Since the SV is not included in the state vector of numerical geodynamo models, it will be connected through a non-linear observation operator, H, which transforms the model state space to the observations space.Obviously H will depend on, among others, fundamental physical properties of the magnetic field.
It should be pointed out here that assimilating the rate of change of geodynamic observables has been routinely used in numerical weather prediction (NWP).For example, precipitation rate, measured from a variety of satellite instruments is assimilated, despite not being a state variable in a GCM (Hou et al, 2000).It should also be pointed out that, in addition to core flow inversion (Roberts and Scott, 1965), there are also attempts to invert core dynamical state with both the surface observations and the dynamo models (Aubert, 2013(Aubert, , 2014)).The latter will benefit the SV assimilation.
In this paper, we describe in detail the results from our recent effort on assimilation of both the field and its SV.These results, from a series of experiments, will demonstrate the improvement in prediction, and knowledge on core flow responses to the SV assimilation.The results also provide valuable information for further development in this direction.
This paper is organized as follows: the numerical model details and the mathematical formulation for SV assimilation will be given in the next section.Followed are the experimental results we have with this assimilation approach.Discussions and plans for further improvements are presented in the last Section.

Mathematical Description
The mathematical formulation for SV assimilation depends on the numerical geodynamo models and the assimilation algorithms, in addition to the physics controlling the time variation of the magnetic field.In this section, we provide the mathematical methodologies used in MoSST − DAS employed in this study (Kuang et al, 2008;Sun and Kuang, 2015).But, with some modifications, they can be applied to other GDAS systems.

Dynamo State Vector and Geomagnetic Observation
MoSST − DAS utilizes the MoSST core dynamics model for time integration of the magnetic field (Kuang et al, 2008;Sun and Kuang, 2015).In this system, the state includes the velocity field v and the density anomaly δ in the outer core r i ≤ r ≤ r c (r i and r c are the mean radii of the ICB and CMB, respectively); and the magnetic field B in the outer core, the electrically conducting inner core r ≤ r i and the D"layer r c ≤ r ≤ r d (r d is the mean radius at the top of the layer).The superscript "T " in (2) implies the transpose.The solid mantle above the D"-layer r d ≤ r ≤ r s (r s is the mean radius of the Earth's surface) is electrically insulating.The whole system is defined in the reference frame fixed with the solid mantle.
The velocity field v and the magnetic field B are decomposed into the poloidal and toroidal components, with the scalars described via spherical expansions where 1 r is the unit radial vector, θ is the co-latitude, φ is the longitude, Y m l are the fully normalized spherical harmonic functions of degree l and order m, L M is the truncation order, and C.C. implies the complex conjugate part.P and T in (3) are called the poloidal and toroidal scalars.It is therefore convenient to write where the subsets are defined with the relevant spectral coefficients in (4), e.g., for the poloidal magnetic field.( 5) and ( 6) can be different for other dynamo models.
In geomagnetic field modeling, geomagnetic measurements are used to obtain the magnetic field B o originated from the core (simply called the geomagnetic field hereafter) that is described as where P m l is the Schmidt normalized associate Legendre polynomial of degree l and order m, (g m l , h m l ) are the Gauss coefficients (slightly different from the standard notation), and L o is the maximum degree (L o ≤ 13 in general).Since these Gauss coefficients (g m l , h m l ) are provided by different field models over the past 10000 years (e.g.Jackson et al, 2000;Korte et al, 2005Korte et al, , 2011;;Gillet et al, 2013;Olsen et al, 2014;Sabaka et al, 2015), they are used as the "observations" in our study.By (3), ( 4), ( 7) and ( 8), we can obtain the relationship between (g m l , h m l ) in ( 8) and b m l in (4) via the radial component B r of the magnetic field B With the definitions of Y m l and P m l , (9) requires that The spectral coefficients of the SV are the time derivatives of (10): where (˙) means the time derivative.

SV and Core State
Geomagnetic observations only provide the time series of (g m l , h m l ).The SV coefficients ( ġm l , ḣm l ) are actually derived.Assimilation of the SV thus raises two major concerns: could the SV be approximated as "instantaneously" measured, and whether it is redundant to the assimilation of the field?
Answers to the first concern depend on the significance of numerical errors in SV calculation.Consider, for example, a central difference scheme is used, Then the relative numerical error is of order where τ o is the typical time intervals of data series, and τ l , defined in (1), is the typical time scales of the observed geomagnetic field.In general, τ o ≤ 1 month in the field models using modern observatory and satellite data (e.g.Sabaka et al, 2004Sabaka et al, , 2015;;Olsen et al, 2006Olsen et al, , 2014)), while τ l ≤ 70 years (see Figure 1).Thus n ≈ 10 −6 , which leads to an order 10 −4 nT/year error in SV.On the other hand, the external field is several tens of nT at the Earth's surface (Sabaka et al, 2015), and changes on the solar cycle (∼ 11 years) and shorter time scales.Thus, n is negligible compared to those arising from, e.g.separation of the external and the internal magnetic signals.One could then argue that both the field and its SV are "concurrently" measured.
The redundancy is not an issue because the observed SV brings different knowledge of the core state x compared to the observed field.To see this, let us consider the magnetic induction of the poloidal magnetic field beneath the impenetrable and "free-slip" CMB (r = r − c ) and in the D"-layer In ( 12), the subscript "h" implies the horizontal components of the velocity field v, and η is the magnetic diffusivity of the outer core fluid; η d in ( 13) is the magnetic diffusivity of the D"-layer (η ≤ η d in general).These two equations show clearly that the observed ḃm(o) l will impose the constraint on v, and on the non-potential part of the poloidal field.
The latter, i.e. ( 13 Therefore, SV assimilation is not redundant to the field assimilation.Indeed, our earlier assimilation results in Figure 3

New Assimilation Approach
We have been using the sequential assimilation approach in MoSST − DAS (e.g.Kuang et al, 2008;Sun and Kuang, 2015).It can be summarized as follows: at the analysis time t a when the observation y is made, a new initial condition x a (called the "analysis") is made from the forecast x f and the observation y, future forecast for t > t a can then be made with the following initial value system: If there is a linear observation operator H that projects x to the observation space (where y is defined), then the analysis x a is of the form where K is called the gain matrix, P f and R are the error covariances of the forecast x f and of the observation y, respectively.( 15) is obtained to minimize the error |H • (x t − x a )| 2 between the analysis x a and the truth x t .In our previous studies, P f is calculated from an ensemble of x f (Sun et al, 2007;Sun and Kuang, 2015), or with some empirical formulations (Kuang et al, 2009;Tangborn and Kuang, 2015).The process is repeated again at the next analysis time t a + Δt (Δt is called the "analysis cycle").If only the observed field is assimilated, then By ( 2) and ( 5), H is linear and very simple where H b corresponds to the subset x b , and has only non-zero entries for b m l (r d ) with 0 ≤ m ≤ l ≤ L o .If the observed SV is also assimilated, then However, by ( 12) and ( 13), transformation between y ḃ and x f is a differentialfunctional projection and is denoted as H x f .One could of course construct an independent projection system which evaluates H x f directly (e.g.Kalnay, 2003).Alternatively, a linearization approximation H(x f ) ≈ H • x f could be made so that (15) can still be used., and a non-potential field that accounts for ḃm(o) Obviously, at the top of the D"-layer The relative errors of (21) are, via the Taylor expansion, of order [(r d − r c )/r c ] 3 .For a 20 km layer thickness, it is smaller than 10 −6 .( 21) allows us to extend the surface observations to the CMB.The observation vector y is now of the form The observation projection is again linear: with H defined in (18).However, H b now includes non-zero entries on all grid points in the D"-layer r c ≤ r i ≤ r d .We can use this approach to further construct an effective observed velocity field v o beneath the CMB (r = r − c ).Since ḃm l is continuous across r = r c the CMB, by ( 12), ( 13) and ( 21), we have Obviously, ( 24) is an under-determined system, since both b If the effective observed velocity field v o (r − c ) is included, then the observation vector y is where y ω includes, as shown in (3) and (4), the spectral coefficients ω Again, the linearized observation projection ( 23) is achieved.However, H includes additional subsets: where which can be calculated from those of the Gauss coefficients g m l and h m l .In this section, we only describe a formal procedure without going into the details.
In geomagnetic field modeling (Jackson et al, 2000;Sabaka et al, 2004;Korte and Constable, 2005;Olsen et al, 2006;Gillet et al, 2013), the Gauss coefficients, e.g.g m l , can be described in general as where S is the vector describing deterministic, model specific base functions in the time domain, e.g.B-spline functions, and α is the coefficient vector which describes the observation error statistics.
For illustrative purpose, we use the simplest error statistics for our derivation.Assume that geomagnetic observations (and thus α) are unbiased, and with known error covariances: where α t lm ≡ α lm is the truth (expectation) and C α is the observation error covariance matrix of α lm .Thus, by (28), Similar formulation applies to h m l as well.By ( 10) and ( 30), we have One can use this equation to evaluate the covariance at any location in the mantle, including r = r d the top of the D"-layer.If S in (30) is replaced by Ṡ, then we can obtain the covariance R lm ġ of the SV, and therefore the variance of ḃm(o) The full error covariance of b m(o) l (r) can then be determined from ( 21), ( 34) and (36).

Results
In this study, we focus only on ( 22 at the top of the D"-layer, mainly for two goals: to explore improvements of the assimilation system with the observed SV, such as the model spin-up process and rms of the observed minus forecast (O-F) of the magnetic field; and to understand responses of the core state x to the observed SV, in particular changes of the velocity field v beneath the CMB.Both are critical for determination of the effective velocity field v in (25), and thus for implementation of the more comprehensive observation (26).
We consider only the observations for the time period 1900 − 2000 simply because modern observatory and satellite data provide very high quality (g m l , h m l ) and ( ġm l , ḣm l ).These coefficients are from GUFM1 (Jackson the al, 2000) for 1900-1962 and CM4 (Sabaka et al, 2004)  , and the differences between the solutions of Case II and Case III are due to the assimilation of the observed SV ḃm(o) l .These allow us to understand clearly the responses of the core state to surface observations, and their dynamical consequences.
We use a modeled observation error covariance, since the actual error covariances of the field models are not yet available.The model error covariance R is assumed diagonal, with the diagonal elements defined as where 0 and 1 decreases linearly in time: 0 decreases from 0.01 in 1900 to 0.001 in 2000, and 1 decreases from 0.3 in 1900 to 0.1 in 2000.These imply that the relative errors in (38) decreases in time, but increases with the degree l.
We would like to point out here that Gillet et al (2013) provided a global field model which includes a full error covariance of the Gauss coefficients.This model and any future model with specified error statistic knowledge are more appropriate for GDAS.However, we conjecture that ( 38) is sufficient for our current objectives.

Responses of the Magnetic Field to SV Assimilation
The quantities used to understand the responses of the magnetic field are the (O-F) of the radial magnetic field B r and its SV Ḃr .Instead of using traditional (O-F), we prefer the following modified definition  From this figure, we can observe clearly that their magnitudes in Case III are approximately 30% smaller than those in Case II over the entire assimilation period, demonstrating a substantial improvement in forecast accuracies with the SV assimilation ( 21) and ( 22).
The SV assimilation also helps accelerate the dynamo model spin-up process.For example, we can observe from Figure 2 that the time variations of (O-F) B are nearly identical in both cases: they decay nearly monotonically over much of the assimilation period before leveling off in the last 20 years (from 1980 to 2000).But (O-F) Ḃ , as shown in Figure 3, are very different in the two cases: in Case II, it increases first from 1900 to 1940; and only starts to decay continuously in the last 20 years.In Case III, however, (O-F) Ḃ decays almost monotonically in time, except two small surges around 1940 and 1980.This implies that the dynamo core state x f responds stronger to the SV assimilation.In other words, the SV assimilation helps to accelerate the model spin-up process.
To better understand how do the forecasts b m(f ) l and ḃm(f) l respond to the observations y b in ( 17) and y b in ( 22), we examine first the (O-F) for individual degrees.In Figure 4 are (O-F) B for the degrees l ≤ 6. Improvements are clearly shown in all 6 degrees, as all values are smaller in Case III than in Case II.But we can also observe significant differences in individual degrees.For example, (O-F) B for the odd degrees (l = 1, 3, 5) increase in magnitude from around 1980.But those for the even degrees (l = 2, 4, 6) do not show either visible increases or increases far less significant than those for the odd degrees.
As shown in Figure 5, the difference between the odd and even degrees of (O-F) Ḃ is even more significant.There is still a strong surge in magnitude for l = 3 around 1980 in the both cases.But the reduction for l = 5 is minimal.In particular it does not decay monotonically in time in either case.These differences may indicate potential inconsistencies between the core dynamics of the model and the time variation of the Gauss coefficients.We will discuss this again later in this paper.

Responses of the Velocity Field to SV Assimilation
Why does the dynamo model respond faster and stronger in Case III than in Case II?We can find at least partial answers from the difference between the free-running model solutions x M (Case I), and the forecasts x f in Cases II and III, in particular the differences in the velocity field v beneath the CMB, because they are the direct consequences of the magnetic induction (12).The knowledge is also very important for obtaining the "effective" observed velocity field (25) for future studies.
Since in our geodynamo model, the CMB is impenetrable and is free-slip, the radial velocity v r = 0 and, by (3) and (4), the horizontal velocity v h depends on ∂v m l /∂r and ω m l at r = r c .Therefore, it is very convenient to examine the following two variables beneath the CMB: where v r is poloidal and describes the up-and-down welling, and ω r is toroidal and describes the differential rotation.The rms differences (M-F) of these two variables between the free-running model solutions (Case I) and the forecasts (Cases II and III) can be used to quantify the responses of the core flow to the assimilation of surface observations: In the above equations, • 2 is the L 2 −norm (or rms) over the CMB.
In Figure 6 are the non-dimensional (with the scaling factor 5 × 10 −6 year −1 for dimensional values) v r 2 (red) and ω r 2 (blue) of the free-running model (Case I).As shown in the figure, v r increases slightly in magnitude in the assimilation period, and ω r remains flat.But, the rms differences (M-F) v P (shown in Figure 7) and (M-F) v T (shown in Figure 8) increase in time, i.e. a growing divergence between the forecast state x f and the free-running model state x M .
From Figures 7 and 8, we can also observe that (M-F) of Case III (the solid lines) are slightly larger than those of Case II (dashed lines), implying that x f moves away from x M faster with the SV assimilation ( 22), another demonstration of improved model spin-up with the SV assimilation.However, the differences are much less significant than those of the magnetic field.This suggests the need for the effective observed velocity field v o to increase further (M-F) of the velocity field, and thus to expedite the model spin-up process.
To aid the future study of determining the effective core flow from the observed SV via (25), we need to understand better the details of (M-F), e.g.their distributions in the spectral space defined by the spherical harmonic degrees l and orders m.We shall pay special attention to their distributions in l, i.e. the summation of the terms in (42-43) with 0 ≤ m ≤ l for a given degree l; and their distributions in m, i.e. the summation of the terms in (42-43) with m ≤ l ≤ L M for a given order m.Since, as shown in Figures 7 and 8, the differences between the two cases are very small, we can focus only on Case III without loss of generality.
In Figure 9 is the distribution of (M-F) v P in the degree l, and in Figure 10 is its distribution in the order m.From the figures we can find that (M-F) v P varies substantially in the spectral spaces.As shown in Figure 9, the differences for the degrees 15 ≤ l ≤ 35 increase the fastest in time, and their magnitudes are the largest at the end of the assimilation period, with the peak at l = 20.The differences are much smaller and grows slower in time for the degrees l ≤ 5 and l ≥ 40.But, as shown in Figure 10, the distribution in m is more broad band: the differences for 5 ≤ m ≤ 35 increase rapidly in time and reach comparable values in magnitude at the end of the assimilation period.However, (M-F) for m ≤ 4 are very different: they remain small and nearly unchanged throughout the entire assimilation.These suggest that the responses of the poloidal velocity is dominantly non-axisymmetric (m > 0).
The distribution of (M-F) v T of the toroidal velocity, as shown in Figures 11  and 12, displays both similar and distinct characteristics.Its distribution in l is very similar to that of (M-F) v P , except that it peaks at a higher degree l = 30.But its distribution in m (Figure 12) is very different: the differences for m ≤ 20 remain comparable in both the magnitude and the time increasing rate.But they decay rapidly for larger m.It should be pointed out in particular that, opposite to (M-F) v P (in Figure 10), (M-F) v T of the axisymmetric toroidal velocity (m = 0) remains very large, implying that the axisymmetric toroidal flow is very sensitive to the surface observations.

Conclusions
In this study we have examined the consequences of assimilating the observed SV on geomagnetic forecasts and on the responses of the dynamo core state.We argued that, because geomagnetic data sampling frequencies are several orders of magnitude higher than those of the SV, the geomagnetic field and its SV are concurrently measured.We further demonstrated that the observed SV provides unique knowledge of the magnetic field and the velocity field in the core.Thus assimilations of the observed field and of the observed SV are necessary and are not redundant.
In this study, we incorporate the observed SV into the observation vector y via introducing the effective poloidal field b m(o) l (21) in the D"-layer, which is then used in the sequential assimilation algorithm (15).We designed three experiments (37) to identify the impact of SV assimilation: a free-running model dynamo simulation (Case I); an experiment with the assimilation of the observed field (Case II), and an experiment with the assimilation of both the observed field and its SV.The relative (O-F) of the field and SV, defined in (39), at the top of the D"-layer are used to measure the forecast accuracies; the (M-F) of the poloidal velocity field (42) and of the toroidal velocity field (43) beneath the CMB are used to characterize the responses of the core state to the SV assimilation.
The results of our experiments demonstrate clearly that the SV assimilation with (21) improves significantly the geomagnetic forecast accuracies since, as shown in Figures 2 and 3, both (O-F) B and (O-F) Ḃ in Case III are more than 20% smaller than those in Case II.In particular, the improvements occur to all degrees, as shown in Figures 4 and 5.The nearly monotonic decay in time of (O-F) Ḃ in Case III (Figure 3) shows clearly that the SV assimilation accelerates the spin-up of the dynamo model.
The improvement by the SV assimilation can be also seen from the differences (M-F) v P and (M-F) v T between the free-running model state and those of the assimilations.As shown in Figures 7 and 8, these differences grow rapidly in time, showing an accelerated departure of the core state with assimilation from the freerunning model state.The differences in Case III are slightly larger than those in Case II, further demonstrating the improvement brought by the SV assimilation, though such increment is less significant that those in the (O-F) of the magnetic field (Figures 2 and 3).
Our results have further implications.First, even with the help of ( 21), the dynamo model is still not fully spun up.For example, though (O-F) Ḃ decreases monodically in time, the SV forecast is still very far away from the observations, as (O-F) Ḃ ≈ O(1) for all degrees (see Figure 5).This can be shown further by the continuously growing differences (M-F) v P and (M-F) v T between the forecast velocity field and that of the free-running model (see Figures 7 and 8).In addition, the differences at the end of the assimilation period are still very small, approximately 10% in magnitude of the velocity field of the free-running model (Figure 6).
These suggest that much larger velocity differences (M-F) v P and (M-F) v T are needed over a shorter assimilation period for expediting the model spin-up.Assimilation of the effective observed velocity (25) could be an answer.However, as discussed earlier, ( 25) is an underdetermined system, since both v An alternative answer could be the core states inverted from the surface observations and dynamo solutions, such as those of Aubert (2013Aubert ( , 2014)).These can be used as the analysis of the assimilation system.But cautions should be taken with this approach.For example, the inverted velocity field beneath the CMB is actually derived with the observed field and SV and the magnetic diffusion of the dynamo state (Aubert 2014).This could potentially lead to dynamical inconsistencies, as well as uncertainties in error statistics.
Our results also show several new features that may have implications to field modeling and to core flow inversion.One new knowledge is from the time variation of (O-F) Ḃ .As shown in Figure 5, (O-F) Ḃ of the odd degrees (l = 1, 3, 5) are significantly different from those of the even degrees (l = 2, 4, 6): the values of the even degrees decay nearly monotonically in time; but those of the odd orders show either spikes (for l = 1, 3) during the assimilation, or even increase over time (for l = 5).These even-odd degree disparities suggest inconsistencies between the model and the observations.These inconsistencies could be entirely due to numerical dynamo model which may include a magnetic induction different from those in the Earth's outer core, or may include some mechanisms resulting in different symmetry properties of the core state.But the inconsistencies could also come from possible biases in the field models that are not included in the observation error covariances.For example, ionospheric ring current generated field (an external field component) contributes dominantly to the Gauss coefficients of degrees l = 1, 3, 5, and varies on time scales comparable to those of SV, e.g. the solar activity cycles (Sabaka et al, 2015).Model biases exist if this part of the signals is not well separated from those of the core field.This is potentially an area for application of geomagnetic data assimilation.
The core fluid flow responses, i.e. the differences (M-F) v P and (M-F) v T between the forecast velocity field v f h and the v M of the free-running model (see Figures 9-12), from our experiments could also help inversion of core flow from the observed SV.For example, the different characteristics in (M-F) v P and (M-F) v T suggest that the poloidal velocity field and the toroidal velocity field could be treated separately in the core flow inversion.It should be pointed out that purely toroidal core flow approximations were used in previous studies (e.g.Bloxham et al, 2002;Olsen and Mandea, 2008).The strong responses of the high degree core flow (l ≈ 20 for the poloidal flow and l ≈ 30 for the toroidal flow) to the observed SV (for l ≤ 8) indicate that higher degree velocity field should be included in the core flow inversion.For example, one would normally expect that, due to nonlinear effects, i.e. the quadratic terms in Navier-Stokes equation and the induction equation, the core flow up to the degrees twice as much as that of the SV should be sufficient for the core flow inversion (as in Aubert, 2013).But our results show that time evolution of the core flow leads to the strongest responses for the degrees more than triple of the maximum degree of the SV.Therefore, inversion of time-dependent core flow from the observed SV (up to degree 13) should include high degree (l > 40) spectral coefficients.
Again we should point out that our results could be improved with more sophisticated assimilation algorithms and field models with more accurate error statistics.For example, we anticipate more accurate estimation of (O-F) for both the field and the SV, and better assessment of the core state responses if a full ensemble approach is used for the covariance P f , and a more appropriate observation error covariance, e.g.those determined by Gillet et al (2013), than (38) used in this study.Regardless, our assimilation experiments have shown clearly the importance of SV assimilation, and the improvements that the SV assimilation brought to forecast accuracies and to model spin-up processes.

Figures legends
Figure 1 The time scales τ l (1) derived from CM4 for the period 1960-2000: the solid line is for the dipole (l = 1), with τ l ≈ 1500 years; the dashed line is for the non-dipolar components (l ≥ 2), with τ l ≈ 70 years.), showing a clear improvement in forecast accuracies.
), implies that, at the top of the D"-layer (r = r d ), one could use a purely potential field b m(p) l to match the observed field b m(o) l .However, b m(p) l can not recover the observed ḃm(o) demonstrate clearly that assimilation of b m(o) l could not reduce the differences between the forecast SV ḃm(f) l and the observed SV ḃm(o) l , called (O-F) of the SV, although that of the field is reduced very rapidly in the first few analysis cycles, a strong indication for the need of SV assimilation.
There are different means to linearize H(x f ).In our current study, we create an effective observed field b m(o) l defined in the D"-layer that matches both b m(o) l and ḃm(o) l .In this approach, b m(o) l comprises of a potential field that accounts for b m(o) l are unknown at r − c .But one can find the "best-fit" v o and b v and H ω include only non-zero entries for v Since the gain matrix K in (16) depends on the observation error covariance R, we need to determine the effective error covariance R for b m(o) l ), i.e. assimilation of the effective observed field b m(o) l which matches both the observed field b m(o) l and the observed SV ḃm(o) l 39) at r = r d .Replacing b m l by ḃm l in (39), we have (O-F) Ḃ of the SV.This modified (O-F) can tell us more accurately how close is the forecast to observation, because it eliminates the effect of changes in magnitude of the individual spectral coefficients.

Figure 2
Figure2are the (O-F) B and (O-F) Ḃ of Case II (dashed lines) and Case III (solid lines).From this figure, we can observe clearly that their magnitudes in Case III are approximately 30% smaller than those in Case II over the entire assimilation period, demonstrating a substantial improvement in forecast accuracies with the SV assimilation (21) and (22).The SV assimilation also helps accelerate the dynamo model spin-up process.For example, we can observe from Figure2that the time variations of (O-F) B are nearly identical in both cases: they decay nearly monotonically over much of the assimilation period before leveling off in the last 20 years(from 1980 to 2000).But (O-F) Ḃ , as shown in Figure3, are very different in the two cases: in Case II, it increases first from 1900 to 1940; and only starts to decay continuously in the last 20 years.In Case III, however, (O-F) Ḃ decays almost monotonically in time, except two small surges around 1940 and 1980.This implies that the dynamo core state x f responds stronger to the SV assimilation.In other words, the SV assimilation helps to accelerate the model spin-up process.To better understand how do the forecasts b Figures11 and 12) are needed to determine v (o) h .For example, as shown in the two figures, the non-axisymmetric (m > 0) poloidal velocity v m l around the degree l = 20, and the toroidal velocity ω m l around the degree l = 30 and order m ≤ 20 should be given more attention, as they are most sensitive to the surface observations.An alternative answer could be the core states inverted from the surface observations and dynamo solutions, such as those ofAubert (2013Aubert ( , 2014)).These can be used as the analysis of the assimilation system.But cautions should be taken with this approach.For example, the inverted velocity field beneath the CMB is actually derived with the observed field and SV and the magnetic diffusion of the dynamo state(Aubert 2014).This could potentially lead to dynamical inconsistencies, as well as uncertainties in error statistics.Our results also show several new features that may have implications to field modeling and to core flow inversion.One new knowledge is from the time variation of (O-F) Ḃ .As shown in Figure5, (O-F) Ḃ of the odd degrees (l = 1, 3, 5) are significantly different from those of the even degrees (l = 2, 4, 6): the values of the even degrees decay nearly monotonically in time; but those of the odd orders show either spikes (for l = 1, 3) during the assimilation, or even increase over time (for l = 5).These even-odd degree disparities suggest inconsistencies between the model and the observations.These inconsistencies could be entirely due to numerical dynamo model which may include a magnetic induction different from those in the Earth's outer core, or may include some mechanisms resulting in different symmetry properties of the core state.But the inconsistencies could also come from possible biases in the field models that are not included in the observation error covariances.For example, ionospheric ring current generated field (an external field component) contributes dominantly to the Gauss coefficients of degrees l = 1, 3, 5, and varies on time scales comparable to those of SV, e.g. the solar activity cycles(Sabaka et  al, 2015).Model biases exist if this part of the signals is not well separated from those of the core field.This is potentially an area for application of geomagnetic data assimilation.The core fluid flow responses, i.e. the differences (M-F) v P and (M-F) v T between the forecast velocity field v f h and the v M of the free-running model (see Figures 9-12), from our experiments could also help inversion of core flow from the observed SV.For example, the different characteristics in (M-F) v P and (M-F) v T suggest that the poloidal velocity field and the toroidal velocity field could be treated separately in the core flow inversion.It should be pointed out that purely toroidal core

Figure 2
Figure 2 The rms (O-F ) B of the magnetic field in Case II (dashed line) and Case III (solid line).In both cases, (O-F ) B 1 and decays monotonically after the first 3 analysis cycles, and then levels off in the last 20 years.This shows the continuing improvement in the forecast accuracies.In addition, the (O-F ) results in Case III (with the assimilation of b m(o) l and ḃm(o) l ) are in general

Figure 3
Figure 3 Similar to Figure 2, but for (O-F ) Ḃ of the SV.In Case II (dashed line), (O-F ) Ḃ = O(1) for much of the assimilation period before decays gradually in the last 20 years, implying that there is no similarity between the forecasted SV ḃm(f) l

Figure 4
Figure4The (O-F ) B of the first 6 spherical harmonic degrees in Case II (dashed lines) and Case III (solid lines).

Figure 6
Figure 6The non-dimensional v r 2 (red) and ωr 2 /10 (blue) beneath the CMB r = r − c from the free-running model solutions (Case I).The dimensional values can be obtained with the scaling factor 5 × 10 −6 year −1 .

Figure 7
Figure7The (M-F ) of the poloidal velocity field v r as defined in (42).The dashed lines are the results without SV assimilation (Case II) and the solid lines are those with the SV assimilation (Case III).

Figure 8
Figure 8 Similar to Figure 7, but for The (M-F ) of the toroidal velocity ωr as defined in (43).

Figure 9
Figure 9The distribution of (M-F ) v P in spherical harmonic degrees l with the SV assimilation (Case III).

Figure 10
Figure 10 Similar to Figure , but for the distribution of (M-F ) v P in spherical harmonic orders m.

Figure 11 Figure 12
Figure 11 Similar to Figure 10, but for (M-F ) v T .
Except the differences in the data y in analysis, everything else is identical in the experiments, including the original initial state at 1900.The analysis cycle is Δt = 5 years.By this design, we can identify exactly the causes of changes in the dynamo state x: the differences between the solutions of Case I and Case II are due to assimilation of the observed field b