Skip to main content
  • Research Article
  • Open access
  • Published:

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

Abstract

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.

Background

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 et al. (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 model 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 variation. 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 secular variation (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. Fournier et al. (2011, 2013) also tested their approach with a geomagnetic field model. Aubert (2013, 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. (2011, 2014) have been continuing their effort on a new combined forward and adjoint system towards a full geodynamo model. 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 model 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 particular importance to the spin-up processes of the numerical models which is characterized by the difference between the observation and the forecast, often called (\({\mathcal {O}}\)-\({\mathcal {F}}\)). These, 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. 2004, 2015; Olsen et al. 2006, 2014) show that the typical time scales τ l of the degree l components (Stacey 1992; Hulot and Le Mouël 1994; Olsen et al. 2006)

$$ \tau_{l} = \left[\frac{\sum_{m} \left({g_{l}^{m}}\right)^{2}+\left({h_{l}^{m}}\right)^{2}} {\sum_{m} \left({\dot g}_{l}^{m}\right)^{2}+\left({\dot h}_{l}^{m}\right)^{2}}\right]^{1/2} $$
((1))

varies from over 1000 years for the dipole (l=1) to less than 100 years for higher degrees. In (1), \(\left ({g_{l}^{m}}, {h_{l}^{m}}\right)\) are the Gauss coefficients of the field, and \(\left ({\dot g}_{l}^{m}, {\dot h}_{l}^{m}\right)\) 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, 2015; Olsen et al. 2006, 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. OSSEs results also suggest long spin-up time in geomagnetic data assimilation. For example, Liu et al. (2007) showed, via their OSSEs with synthetic magnetic data and a fully nonlinear dynamo model, that the difference between the forecast and the truth reach the minimum in approximately 40 % of the magnetic free decay time, or 8000 years, if only the poloidal field (of the first eight spherical harmonic degree coefficients) is assimilated (see Figure four in Liu et al. 2007). Therefore, 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 for 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 spin 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. (2011, 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 \(\left ({\dot g}_{l}^{m}, {\dot h}_{l}^{m}\right)\) of the Gauss coefficients since, as we will describe in the next section, they provide additional constraints on the core flow beneath the core-mantle boundary (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,, which transforms the model state space to the observations space. Obviously, 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 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).

There is also a long history of attempting to invert core state from surface observations. In addition to core flow inversion initiated by Roberts and Scott (1965), Zatman and Bloxham et al. (1997) extended further the inversion of the poloidal magnetic field in the outer core. More recently, Aubert (2013, 2014) attempted a new approach to invert core dynamical state with the observed SV and the dynamo models. These inversions, however, are not data assimilation. However, the inversion results could benefit geomagnetic data assimilation. For example, the inversion by Aubert (2013, 2014) could be utilized for making “analysis” in a geomagnetic assimilation system.

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.

Methods

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 vector x

$$ \mathbf{x} = \left(\mathbf{v}, \mathbf{B}, \delta\varrho\right)^{T} $$
((2))

includes the velocity field v and the density anomaly δ ϱ in the outer core r i rr c (r i and r c are the mean radii of the inner core boundary (ICB) and CMB, respectively); and the magnetic field B in the outer core, the electrically conducting inner core rr i , and the D -layer r c rr 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 rr 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

$${} \begin{aligned} \left(\mathbf{v}, \mathbf{B}\right)^{T} &= \nabla\times \left[\left(T_{v}, T_{b}\right)^{T}\mathbf{1}_{r}\right]+ \nabla\times\nabla\times \left[\left(P_{v}, P_{b}\right)^{T}\mathbf{1}_{r}\right], \end{aligned} $$
((3))
$${} \begin{aligned} \left(P_{v}, T_{v}, P_{b}, T_{b}, \delta\varrho\right)^{T} &= \sum_{0\le m\le l}^{L_{M}}\left({{v_{l}^{m}}},{{\omega_{l}^{m}}},{{b_{l}^{m}}}, {j_{l}^{m}}, {\vartheta_{l}^{m}}\right)^{T}\\ &\quad\times{Y_{l}^{m}}(\theta,\phi) + C.C., \end{aligned} $$
((4))

where 1 r is the unit radial vector, θ is the co-latitude, ϕ is the longitude, \({Y_{l}^{m}}\) 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

$$ \mathbf{x} \,\ = \left(\mathbf{x}_{v},\mathbf{x}_{\omega}, \mathbf{x}_{b},\mathbf{x}_{j},\mathbf{x}_{\rho}\right)^{T}, $$
((5))

where the subsets are defined with the relevant spectral coefficients in (4), e.g.,

$$ \mathbf{x}_{b} = \left\{{{b_{l}^{m}}}(r_{k}) \mid 0 \le r_{k} \le r_{d};\ 0\le m\le l\le L_{M} \right\}^{T} $$
((6))

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

$$\begin{array}{@{}rcl@{}} \mathbf{B}^{o} & =& -\nabla\Psi, \end{array} $$
((7))
$$\begin{array}{@{}rcl@{}} \Psi \, \ & =& r_{s}\! \sum_{0\le m\le l}^{L_{o}} \left(\frac{r_{s}}{r}\right)^{l+1}\!\left({{g_{l}^{m}}}\cos m\phi + {{h_{l}^{m}}} \sin m\phi\right) {P_{l}^{m}}(\theta) \end{array} $$
((8))

where \({P_{l}^{m}}\) is the Schmidt normalized associate Legendre polynomial of degree l and order m, \(\left ({{g_{l}^{m}}}, {{h_{l}^{m}}}\right)\) 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 \(\left ({{g_{l}^{m}}}, {{h_{l}^{m}}}\right)\) are provided by different field models over the past 10,000 years (e.g., Jackson et al. 2000; Korte et al. 2005, 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 \(\left ({{g_{l}^{m}}}, {{h_{l}^{m}}}\right)\) in (8) and the observed \({b_{l}^{m(o)}}\) in the form of (4) via the radial component B r of the magnetic field B

$$ \begin{aligned} B_{r}^{o} & = - \frac{\partial\Psi}{\partial r} = \sum_{0\le m\le l}^{L_{o}} (l+1)\left(\frac{r_{s}}{r}\right)^{l+2}\\&\quad\times\left({{g_{l}^{m}}}\cos m\phi + {{h_{l}^{m}}}\sin m\phi\right) {P_{l}^{m}} (\theta)\\ & = -\frac{{\hat L}}{r^{2}}P_{b}^{{o}} = \sum_{0\le m\le l}^{L_{o}}\frac{l(l+1)}{r^{2}}{b_{l}^{m(o)}} {Y_{l}^{m}} + C.C. \end{aligned} $$
((9))

where \(\hat L\) is the angular momentum operator. With the definitions of \({Y_{l}^{m}}\) and \({P_{l}^{m}}\), (9) requires that

$$ \begin{aligned} {b_{l}^{m(o)}}(r) &= \frac{{r_{s}^{2}}}{l}\left(\frac{r_{s}}{r}\right)^{l}G_{m} \left({{g_{l}^{m}}}-i{{h_{l}^{m}}}\right),\\ G_{m} &= \left[\frac{2\pi(1+\delta_{m0})}{2l+1}\right]^{1/2} \end{aligned} $$
((10))

for r d rr s . The spectral coefficients of the SV are the time derivatives of (10):

$${} {\dot b}_{l}^{m(o)}(r) = \frac{{r_{s}^{2}}}{l}\left(\frac{r_{s}}{r}\right)^{l}G_{m}\left({\dot g}_{l}^{m}-i{\dot h}_{l}^{m}\right) \quad\text{for}\quad r_{d}\le r\le r_{s}, $$
((11))

where \(\left (\dot \ \right)\) means the time derivative.

SV and core state

Geomagnetic observations only provide the time series of \(\left ({{g_{l}^{m}}}, {{h_{l}^{m}}}\right)\). The SV coefficients \(\left ({\dot g}_{l}^{m}, {\dot h}_{l}^{m}\right)\) 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:

$${\dot g}_{l}^{m}(t) \approx \frac{{{g_{l}^{m}}}(t+\delta t)-{{g_{l}^{m}}}(t-\delta t)}{2\delta t}. $$

Then, the relative numerical error can be estimated as

$$\frac{{{g_{l}^{m}}}(t+\delta t)-{{g_{l}^{m}}}(t-\delta t)}{2\delta t} = {\dot g}_{l}^{m}(t) + {\mathcal{O}}\left[\frac{\partial^{3} {{g_{l}^{m}}}}{\partial t^{3}}\left(\delta t\right)^{2}\right] $$

Since \(\partial ^{3} {{g_{l}^{m}}}/\partial t^{3} \approx {\dot g}_{l}^{m}/{\tau _{l}^{2}}\) by (1) and δ tt o (the typical time intervals of data series), we have

$$\begin{aligned} \frac{{{g_{l}^{m}}}(t+\delta t)-{{g_{l}^{m}}}(t-\delta t)}{2\delta t} &= {\dot g}_{l}^{m}(t)\left[1+ {\mathcal{O}}\left(\tau_{o}/\tau_{l}\right)^{2}\right]\\ &\equiv {\dot g}_{l}^{m}(t)\left(1+\epsilon_{n}\right) \end{aligned} $$

In general, τ o ≤1 month in the field models using modern observatory and satellite data (e.g., Sabaka et al. 2004, 2015; Olsen et al. 2006, 2014), while τ l can be as short as 10 years (Christensen et al. 2012) for high-degree coefficients. 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 nanotesla at the Earth’s surface (Sabaka et al. 2015) 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}^{-}\))

$${} {\dot b}_{l}^{m} = -\frac{r^{2}}{l(l+1)}\left[\nabla_{h}\cdot\left(\mathbf{v}_{h} B_{r}\right)\right]_{l}^{m} + \eta \left[\frac{\partial^{2}}{\partial r^{2}}-\frac{l(l+1)}{r^{2}}\right] {{b_{l}^{m}}}, $$
((12))

and in the D -layer

$$ {\dot b}_{l}^{m} = \eta_{d} \left[\frac{\partial^{2}}{\partial r^{2}}-\frac{l(l+1)}{r^{2}}\right] {{b_{l}^{m}}}. $$
((13))

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 \({\dot b}_{l}^{m(o)}\) will impose the constraint on v and on the non-potential part of the poloidal field.

The latter, i.e., (13), implies that, at the top of the D -layer (r=r d ), a potential poloidal field \(b_{l}^{m(p)}\) can fully recover the observed field \({b_{l}^{m(o)}}\). However, it cannot recover the observed SV \({\dot b}_{l}^{m(o)}\) since

$$\frac{\partial^{2}b_{l}^{m(p)}}{\partial r^{2}}-\frac{l(l+1)}{r^{2}} b_{l}^{m(p)} = 0. $$

In other words, the observed SV provides the information on the non-potential part of the field that is missing in the observed field \({b_{l}^{m(o)}}\). Therefore, SV assimilation is not redundant to the field assimilation.

Indeed, our assimilation results (in Fig. 2) demonstrate clearly that assimilation of \({b_{l}^{m(o)}}\) could not reduce the differences between the forecast SV \({\dot b}_{l}^{m(f)}\) and the observed SV \({\dot b}_{l}^{m(o)}\), called (\({\mathcal {O}}\)-\({\mathcal {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.

Fig. 1
figure 1

The rms \(({\mathcal {O}}\text {-}{\mathcal {F}})_{B}\) of the magnetic field in case II (dashed line) and case III (solid line). In both cases, \(({\mathcal {O}}\text {-}{\mathcal {F}})_{B} \ll 1\) and decays monotonically after the first three analysis cycles and then levels off in the last 20 years. This shows the continuing improvement in the forecast accuracies. In addition, the \(({\mathcal {O}}\text {-}{\mathcal {F}})\) results in case III (with the assimilation of \({b_{l}^{m(o)}}\) and \({\dot b}_{l}^{m(o)}\)) are in general more than 20 % smaller than in case II (with only the assimilation of \({b_{l}^{m(o)}}\)), showing a clear improvement in forecast accuracies

Fig. 2
figure 2

Similar to Fig. 1, but for \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) of the SV. In case II (dashed line), \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B} = {\mathcal {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 \({\dot b}_{l}^{m(f)}\) and the observed SV \({\dot b}_{l}^{m(o)}\). But its magnitude is much smaller in case III (solid line), and it decays monotonically in time, indicating that the SV assimilation accelerates the spin-up process

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:

$$ \frac{\partial \mathbf{x}^{f}}{\partial t} = \mathbf{M}\left(\mathbf{x}^{f}\right), \qquad \mathbf{x}^{f}(t_{a}) = \mathbf{x}^{a}. $$
((14))

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

$$\begin{array}{@{}rcl@{}} && \mathbf{x}^{a} = \mathbf{x}^{f}+\mathbf{K}\left(\mathbf{y}-\mathbf{H}\mathbf{x}^{f}\right) \end{array} $$
((15))
$$\begin{array}{@{}rcl@{}} && \mathbf{K} \, = \mathbf{P}^{f}\mathbf{H}^{T}\left(\mathbf{H}\mathbf{P}^{f}\mathbf{H}^{T}+\mathbf{R}\right)^{-1} \end{array} $$
((16))

where K is called the gain matrix, P f and R are the error covariance matrices of the forecast x f and of the observation y, respectively. (15) is obtained to minimize the error |H·(x tx a)|2 between the analysis x a and the truth x t.

In our assimilation system MoSST DAS, the error covariance P f is calculated via three different approaches: an ensemble-based covariance analysis (Sun et al. 2007; Sun and Kuang 2015), an empirical covariance based on forecast solution properties (Tangborn and Kuang 2015), and an optimal interpolation (OI) scheme with a predefined time-invariant covariance (Kuang et al. 2009). The first approach is computationally very expensive, in which an ensemble N initial states are created x i =x f+ε i (ε i are random white noise perturbations). Their free running model solutions at a later time (often a fraction of the magnetic free decay time) are then used to calculate the covariance \(\mathbf {P}^{f} = \langle \left (\mathbf {x}_{i}-{\mathbf {\overline {x}}}\right) \left (\mathbf {x}_{i}-{\mathbf {\overline {x}}}\right)^{T}\rangle \) (\(\mathbf {\overline {x}} = \sum \mathbf {x}_{i}/N\) is the mean). The empirical covariance is assumed diagonal (no cross covariance between different degrees l and orders m) and is determined by the forecast error standard deviations and an exponentially decaying spatial correlation function. This approach can be updated at analysis time without much computational effort. The OI-type error covariance is similar to the empirical one, except that the forecast error standard deviations are assumed time-invariant (i.e., constant throughout the assimilation) and that spatial correlation does not decay exponentially with the distance. We simply use the OI-type error covariance in this study.

If only the observed field is assimilated, then

$$ \mathbf{y} = \left\{{b_{l}^{m(o)}}(r_{d})\Big\vert 0\le m\le l \le L_{o}\right\}^{T} \equiv \mathbf{y}_{b}. $$
((17))

By (2) and (5), H is linear and very simple

$$ \mathbf{H} = \left({\mathbf 0}, {\mathbf 0}, \mathbf{H}_{b}, {\mathbf 0}, {\mathbf 0}\right)^{T}, $$
((18))

where H b corresponds to the subset x b and has only non-zero entries for \({{b_{l}^{m}}}(r_{d})\) with 0≤mlL o . If the observed SV is also assimilated, then

$$\begin{array}{@{}rcl@{}} && \mathbf{y} \ = \left(\mathbf{y}_{b},\mathbf{y}_{\dot b}\right)^{T} \end{array} $$
((19))
$$\begin{array}{@{}rcl@{}} && \mathbf{y}_{\dot b} \equiv \left\{{\dot b}_{l}^{m(o)}(r_{d})\mid 0\le m\le l \le L_{o}\right\}^{T}, \end{array} $$
((20))

However, by (12) and (13), transformation between \(\mathbf {y}_{\dot b}\) and x f is a differential-functional projection and is denoted as (x f). One could of course construct an independent projection system which evaluates (x f) directly (see, e.g., Remarks 5.3.1 in Kalnay 2003). Alternatively, a linearization approximation (x f)≈H·x f could be made so that (15) can still be used.

There are different means to linearize (x f). In our current study, we create an effective observed field \({{\widetilde {b}}_{l}^{m(o)}}\) defined in the D -layer that matches both \({b_{l}^{m(o)}}\) and \({\dot b}_{l}^{m(o)}\). In this approach, \({{\widetilde {b}}_{l}^{m(o)}}\) comprises of a potential field that accounts for \({b_{l}^{m(o)}}\) and a non-potential field for \({\dot b}_{l}^{m(o)}\):

$$ \begin{aligned} {{\widetilde{b}}_{l}^{m(o)}}(r) &= \left(\frac{r_{d}}{r}\right)^{l}{b_{l}^{m(o)}}(r_{d})+ \frac{1}{2\eta_{d}}\\ &\quad\times\left(r-r_{d}\right)^{2}{\dot b}_{l}^{m(o)}(r_{d}) \quad \text{for}\quad r_{c}\le r \le r_{d}. \end{aligned} $$
((21))

Obviously, at the top of the D -layer r=r d ,

$$\begin{aligned} {{\widetilde{b}}_{l}^{m(o)}} &= {b_{l}^{m(o)}},\quad \frac{\partial{{\widetilde{b}}_{l}^{m(o)}}}{\partial r} = \frac{\partial{b_{l}^{m(o)}}}{\partial r}\\ &= - \frac{l}{r_{d}}{b_{l}^{m(o)}},\quad{{\dot {\widetilde{b}}}_{l}^{m(o)}} = {\dot b}_{l}^{m(o)}. \end{aligned} $$

The coefficients \({{b_{l}^{m}}}(r)\) in the D -layer can be approximated via Taylor series expansion

$$\begin{aligned} {{b_{l}^{m}}}(r) &= {{b_{l}^{m}}}(r_{d}) + {{b_{l}^{m}}}^{\prime}(r_{d}) \left(r-r_{d}\right) + \frac{1}{2}{{b_{l}^{m}}}^{\prime\prime}\\ &\quad\times(r_{d}) \left(r-r_{d}\right)^{2} + {\mathcal{O}}\left[{{b_{l}^{m}}}^{\prime\prime\prime}\left(r-r_{d}\right)^{3}\right], \end{aligned} $$

where implies the radial derivative. Since the approximation (21) satisfies the first three terms in the above expansion, the errors are given by the last term. Since \({{b_{l}^{m}}}^{\prime \prime \prime }\sim {{b_{l}^{m}}}/r^{3}\) in D -layer, the relative errors are of the 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

$$ \mathbf{y} = \left\{{{\widetilde{b}}_{l}^{m(o)}}(r) \mid 0\le m\le l \le L_{o};\ r_{c}\le r \le r_{d}\right\}^{T}\equiv \mathbf{y}_{{\widetilde b}}. $$
((22))

The observation projection is again linear:

$$ {\boldsymbol{\cal{H}}}(\mathbf{x}^{f})=\mathbf{H}\cdot\mathbf{x}^{f} $$
((23))

with H defined in (18). However, H b now includes non-zero entries on all grid points in the D -layer r c rr d .

We can use this approach to further construct an effective observed velocity field \({\mathbf {\widetilde {v}}}^{o}\) beneath the CMB \(\left (r = r_{c}^{-}\right)\). Since \({\dot b}_{l}^{m}\) is continuous across the CMB, by (12), (13) and (21), we have

$${} \begin{aligned} -\frac{{r_{c}^{2}}}{l(l+1)}&\left[\nabla_{h}\cdot\left({\mathbf{\widetilde{v}}}_{h}^{o} {\widetilde B}_{r}^{o}\right)\right]_{l}^{m} +\eta \left[\frac{\partial^{2}}{\partial r^{2}}-\frac{l(l+1)}{{r_{c}^{2}}}\right] {{\widetilde{b}}_{l}^{m(o)}} \\&= {{\dot {\widetilde{b}}}_{l}^{m(o)}}(r_{c}) = {\dot b}_{l}^{m(o)}(r_{d})\!\!\left[\!1-\frac{l(l+1)}{2{r_{c}^{2}}}\!\left(r_{c}-r_{d}\right)^{2}\!\right] \end{aligned} $$
((24))

For those geodynamo models without an electrically conducting D -layer, one can simply replace \({{\dot {\widetilde {b}}}_{l}^{m(o)}}(r_{c})\) in the above equation with \({\dot b}_{l}^{m(o)}(r_{c})\) continued downward from the surface observation.

Obviously, (24) is an underdetermined system, since both \({{\widetilde {b}}_{l}^{m(o)}}\) and \({\mathbf {\widetilde {v}}}^{o}\) are unknown at \(r_{c}^{-}\). But one can find the “best-fit” \({\mathbf {\widetilde {v}}}^{o}\) and \({{\widetilde {b}}_{l}^{m(o)}}\) via minimizing the following difference

$$ \begin{aligned} &\min_{{\mathbf{\widetilde{v}}}^{o}, {{\widetilde{b}}_{l}^{m(o)}}}\left|{{\dot {\widetilde{b}}}_{l}^{m(o)}}(r_{c}) + \frac{{r_{c}^{2}}}{l(l+1)}\left[\nabla_{h}\cdot\left(\mathbf{v}_{h} B_{r}\right)\right]_{l}^{m}\right.\\ &\left.\qquad\qquad\qquad\,\,\,\,- \eta \left[\frac{\partial^{2}}{\partial r^{2}}-\frac{l(l+1)}{{r_{c}^{2}}}\right] {{b_{l}^{m}}}\right|^{2}. \end{aligned} $$
((25))

We provide here only a sketch of the minimization procedure: denoting by x b and x v the vectors of the spectral coefficients of B r and of v h in (25), respectively, and then (25) is equivalent to minimize the (inner) product

$$ \left[\dot{\mathbf{x}}_{b} + {\mathbf{A}}\cdot\mathbf{x}_{b}+ {\mathbf{B}}(\mathbf{x}_{b})\cdot\mathbf{x}_{v}\right]^{T}\cdot \left[\dot{\mathbf{x}}_{b} + {\mathbf{A}}\cdot\mathbf{x}_{b}+ {\mathbf{B}}(\mathbf{x}_{b})\cdot\mathbf{x}_{v}\right] $$
((26))

with respect to x v . Note that A·x b and B(x b x v in (26) describe the diffusion and the advection terms in (25), respectively. Minimizing the product (26) leads to the solution

$$ \left({\mathbf{B}}^{T}\cdot{\mathbf{B}}\right) \cdot {\tilde {\mathbf{x}}}_{v} = -{\mathbf{B}}^{T}\cdot\left[\dot{\mathbf{x}}_{b} + {\mathbf{A}}\cdot\mathbf{x}_{b}\right] $$
((27))

which provides the spectral coefficients of the effective observed velocity field \({\mathbf {\widetilde {v}}}^{o}(r_{c}^{-})\) beneath the CMB. If \({\mathbf {\widetilde {v}}}^{o}(r_{c}^{-})\) is included, then the observation vector y is

$$ \mathbf{y} = \left(\mathbf{y}_{{\widetilde v}},\, \mathbf{y}_{{\widetilde \omega}},\, \mathbf{y}_{{\widetilde b}}\right)^{T} $$
((28))

where \(\mathbf {y}_{{\widetilde \omega }}\) includes, as shown in (3) and (4), the spectral coefficients \({\widetilde \omega }_{l}^{m(o)}\) of \({\mathbf {\widetilde {v}}}^{o}\) at \(r_{c}^{-}\). Again, the linearized observation projection (23) is achieved. However, H includes additional subsets:

$$ \mathbf{H} = \left(\mathbf{H}_{v}, \mathbf{H}_{\omega}, \mathbf{H}_{b}, {\mathbf 0}, {\mathbf 0}\right)^{T}, $$
((29))

where H v and H ω include only non-zero entries for \({\widetilde v}_{l}^{m(o)}\) and \({\widetilde \omega }_{l}^{m(o)}\) at \(r = r_{c}^{-}\), respectively.

Effective observation error covariance

Since the gain matrix K in (16) depends on the observation error covariance R, we need to determine the effective error covariance \({\mathbf {\widetilde {R}}}\) for \({{\widetilde {b}}_{l}^{m(o)}}\) which can be calculated from those of the Gauss coefficients \({{g_{l}^{m}}}\) and \({{h_{l}^{m}}}\). 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 et al. 2005; Olsen et al. 2006; Gillet et al. 2013), the Gauss coefficients, e.g. \({{g_{l}^{m}}}\), can be described in general as

$$ {{g_{l}^{m}}} = {\mathbf S}^{T}(t)\cdot{\boldsymbol\alpha}_{lm}, $$
((30))

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 includes 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:

$$ {\boldsymbol{\alpha}}_{lm} = {\boldsymbol{\alpha}}_{lm}^{t} + {\boldsymbol{\epsilon}}_{\alpha},\quad \langle {\boldsymbol{\epsilon}}_{\alpha}\rangle = {\mathbf 0},\quad \left\langle\boldsymbol{\epsilon}_{\alpha}{\boldsymbol{\epsilon}_{\alpha}^{T}}\right\rangle = {\mathbf{C}}_{\alpha}, $$
((31))

where \({\boldsymbol {\alpha }}_{\textit {lm}}^{t}\equiv \left \langle {\boldsymbol \alpha }_{\textit {lm}}\right \rangle \) is the truth (expectation) and C α is the observation error covariance matrix of α lm . Thus, by (30),

$$ \begin{array}{l} {{g_{l}^{m}}} = g_{l}^{m(t)} + \epsilon_{g},\quad g_{l}^{m(t)} \ = {\mathbf S}^{T}\cdot{\boldsymbol{\alpha}}_{lm}^{t},\quad \epsilon_{g} = {\mathbf S}^{T}\cdot{\boldsymbol{\epsilon}}_{\alpha}, \\ \left\langle{\epsilon_{g}^{2}}\right\rangle= {\mathbf S}^{T}\cdot{\mathbf C}_{\alpha}\cdot{\mathbf S} \equiv R_{g}^{lm}. \end{array} $$
((32))

Similar formulation applies to \({{h_{l}^{m}}}\), but with the error ε h . By (10) and (32), we have

$$\begin{array}{@{}rcl@{}} {b_{l}^{m(o)}}(r) & = & b_{l}^{m(t)}(r) + \epsilon_{b}(r), \end{array} $$
((33))
$$\begin{array}{@{}rcl@{}} b_{l}^{m(t)}(r) & =& \frac{{r_{s}^{2}}}{l}\left(\frac{r_{s}}{r}\right)^{l}G_{m}\left(g_{l}^{m(t)}-ih_{l}^{m(t)}\right), \end{array} $$
((34))
$$\begin{array}{@{}rcl@{}} \epsilon_{b}(r)\ \ \ \ & =& \frac{{r_{s}^{2}}}{l}\left(\frac{r_{s}}{r}\right)^{l}G_{m} \left(\epsilon_{g}-i\epsilon_{h}\right) \end{array} $$
((35))

This leads to

$$ \left\langle\epsilon_{b}\epsilon_{b}^{*}\right\rangle = \left(\frac{{r_{s}^{2}}}{l}\right)^{2}\left(\frac{r_{s}}{r}\right)^{2l}{G_{m}^{2}} \left[\left(R_{g}^{lm}\right)^{2}+\left(R_{h}^{lm}\right)^{2}\right] $$
((36))

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 (32) is replaced by \(\dot {\mathbf S}\), then we can obtain the covariance \(R_{\dot g}^{lm}\) of the SV

$$R_{\dot g}^{lm} = {\dot {\mathbf S}}^{T}\cdot{\mathbf C}_{\alpha}\cdot{\dot {\mathbf S}}, $$

and therefore the variance of \({\dot b}_{l}^{m(o)}\)

$$\begin{array}{@{}rcl@{}} {\dot b}_{l}^{m(o)}(r) & =& {\dot b}_{l}^{m(t)}(r) + \epsilon_{\dot b}(r), \end{array} $$
((37))
$$\begin{array}{@{}rcl@{}} \left\langle\epsilon_{\dot b}\epsilon_{\dot b}^{*}\right\rangle(r) & = & \left(\frac{{r_{s}^{2}}}{l}\right)^{2}\left(\frac{r_{s}}{r}\right)^{2l}{G_{m}^{2}} \left[\left(R_{\dot g}^{lm}\right)^{2}+\left(R_{\dot h}^{lm}\right)^{2}\right]. \end{array} $$
((38))

The full error covariance of \({{\widetilde {b}}_{l}^{m(o)}}(r)\) can then be determined from (21), (36), and (38).

Results

In this study, we focus only on (22), i.e., assimilation of the effective observed field \({{\widetilde {b}}_{l}^{m(o)}}\) which matches both the observed field \({b_{l}^{m(o)}}\) and the observed SV \({\dot b}_{l}^{m(o)}\) 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 (\({\mathcal {O}}\)-\({\mathcal {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 \(\mathbf {\widetilde {v}}\) in (25) and thus for implementation of the more comprehensive observation (28).

The baseline geodynamo model is the MoSST core dynamics model (Kuang and Chao 2003; Jiang and Kuang 2008) for the thermal convection of a Boussinesq electrically conducting fluid in the (rapidly rotating) outer core, confined between the electrically conducting inner core and the D -layer. The non-dimensional parameters include the Ekman number E (for viscosity), the magnetic Rossby number R o (for magnetic diffusivity), the modified Prandtl number q κ (for thermal conductivity), and the modified Rayleigh number R th (for buoyancy force). In our assimilation, the ICB and the CMB are assumed impenetrable, stress-free, and fixed heat fluxes. We also select the following parameter values:

$$ E = R_{o} = 1.25\times 10^{-6},\quad R_{o} = 1.0\, \quad R_{th} = 15 R_{th}^{c}, $$
((39))

where \(R_{\text {th}}^{c}\) is the critical Rayleigh number for purely thermal convection. In our assimilation, the numerical truncation order is L M =96. There are 20 grid points in D -layer (which is 20 km thick with η d =20 η), 80 grid points in the outer core, and 40 grid points in the inner core. With this set of the parameter values, the mean time scale of the dipole field is approximately 0.7 % of the magnetic diffusive time τ η (which is used for the time scaling of the dynamo model) and those for non-dipole components are more than an order of magnitude shorter. These are very similar to those derived from satellite magnetic data. Therefore, in our assimilation, we choose τ η =200,000 years to convert non-dimensional dynamo time to years (with this conversion, the mean time scale of the dipole field is approximately 1400 years).

We consider only the observations for the time period 1900−2000 simply because modern observatory and satellite data provide very high quality \(\left ({{g_{l}^{m}}}, {{h_{l}^{m}}}\right)\) and \(\left ({\dot g}_{l}^{m}, {\dot h}_{l}^{m}\right)\). These coefficients are from gufm1 (Jackson et al. 2000) for 1900–1962 and CM4 (Sabaka et al. 2004) for 1962–2000. We also set L o =8, lower than the highest degrees of the two models. For our research purposes, we carry out three distinct experiments:

  • Case I: Free-running model (no assimilation)

  • Case II: Assimilation of \({b_{l}^{m(o)}}\) with (17)

  • Case III: Assimilation of \({{\widetilde {b}}_{l}^{m(o)}}\) with (22)

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_{l}^{m(o)}}\) and the differences between the solutions of case II and case III are due to the assimilation of the observed SV \({\dot b}_{l}^{m(o)}\). 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

$${} R_{lm} = \left\vert\epsilon_{R}(l) {{b_{l}^{m}}}\right\vert^{2},\ \ \epsilon_{R}(l) = \epsilon_{0}(t) + \left[\epsilon_{1}(t)-\epsilon_{0}(t)\right]\frac{l-1}{L_{o}-1}, $$
((40))

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 (40) 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 (40) 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 (\({\mathcal {O}}\)-\({\mathcal {F}}\)) of the radial magnetic field B r and its SV \({\dot B}_{r}\). Instead of using traditional (\({\mathcal {O}}\)-\({\mathcal {F}}\)), we prefer the following modified definition

$$ \begin{aligned} ({\mathcal{O}}-{\mathcal{F}})_{B}^{2} &= \sum_{1\le l}^{L_{o}}\left\{\left[ \sum_{0\le m\le l}\left|\frac{{b_{l}^{m(o)}}}{b_{1}^{0(o)}}-\frac{{b_{l}^{m(f)}}}{b_{1}^{0(f)}} \right|^{2}\right]\right.\\ &\quad\times\left.\left[ \sum_{0\le m\le l}\left|\frac{{b_{l}^{m(o)}}}{b_{1}^{0(o)}}\right|^{2}\right]^{-1}\right\} \end{aligned} $$
((41))

at r=r d , i.e., the misfit normalized by the observed field strength for a given degree l. Replacing \({{b_{l}^{m}}}\) by \({\dot b}_{l}^{m}\) in (41), we have \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) of the SV. This modified \(({\mathcal {O}}\text {-}{\mathcal {F}})\) can tell us more accurately how close is the forecast to observation because it eliminates the effect of changes in the magnitude of the individual spectral coefficients. Our modified \(({\mathcal {O}}\text {-}{\mathcal {F}})\) in (41) is different from other quantities used for measuring the difference (or misfit) between the observation y and the forecast H·x f, such as the misfit normalized by the error covariances (e.g., Aubert 2014). Since, in our assimilation, the modeled error covariances are proportional in magnitude to the observed field strength; (41) is actually very similar to the normalized misfit.

Figure 1 is the \(({\mathcal {O}}\text {-}{\mathcal {F}})_{B}\) 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). Similar improvement can be also observed from \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) shown in Fig. 2.

The SV assimilation also helps accelerate the dynamo model spin-up process. For example, we can observe from Fig. 1 that the time variations of \(({\mathcal {O}}\text {-}{\mathcal {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 \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\), as shown in Fig. 2, 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, \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) 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_{l}^{m(f)}}\) and \({\dot b}_{l}^{m(f)}\) respond to the observations y b in (17) and \(\mathbf {y}_{\widetilde b}\) in (22), we examine first the \(({\mathcal {O}}\text {-}{\mathcal {F}})\) for individual degrees. In Fig. 3 are \(({\mathcal {O}}\text {-}{\mathcal {F}})_{B}\) for the degrees l≤6. Improvements are clearly shown in all six degrees, as all values are smaller in Case III than those in Case II. But we can also observe different patterns among these degrees. For example, \(({\mathcal {O}}\text {-}{\mathcal {F}})_{B}\) for the odd degrees (l=1,3,5) increase in magnitude again from around 1980 to the end of the assimilation period. But there is no such clear reversing trend in those for the even degrees (l=2,4,6): it decreases monotonically for l=2; for l=4, it oscillates with a damping amplitude; and the \(({\mathcal {O}}\text {-}{\mathcal {F}})\) for l=4 remains nearly constant after the rapid decay in the first two analysis cycles.

Fig. 3
figure 3

The \(({\mathcal {O}}\text {-}{\mathcal {F}})_{B}\) of the first six spherical harmonic degrees in case II (dashed lines) and case III (solid lines)

As shown in Fig. 4, the difference between the odd and even degrees of \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) 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.

Fig. 4
figure 4

Similar to Fig. 3, but for \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\)

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 \(\partial {{v_{l}^{m}}}/\partial r\) and \({{\omega _{l}^{m}}}\) at r=r c . Therefore, it is very convenient to examine the following two variables beneath the CMB:

$${} {\small{\begin{aligned} v_{r}' & \equiv \frac{\partial v_{r}}{\partial r} = \sum_{0\le m\le L}^{L_{M}}\frac{l(l+1)}{{r_{c}^{2}}}\frac{\partial {{v_{l}^{m}}}}{\partial r}{Y_{l}^{m}}\left(\theta,\phi\right) + C.C. \end{aligned}}} $$
((42))
$${} {\small{\begin{aligned} \omega_{r} & \equiv \left(\nabla\times\mathbf{v}\right)_{r} = \sum_{0\le m\le L}^{L_{M}}\frac{l(l+1)}{{r_{c}^{2}}}{{\omega_{l}^{m}}} {Y_{l}^{m}}\left(\theta,\phi\right) + C.C., \end{aligned}}} $$
((43))

where v r′ is poloidal and describes the up-and-down welling and ω r is toroidal and describes the differential rotation. The rms differences (\({\mathcal {M}}\text {-}{\mathcal {F}}\)) between the two variables of the forecast x f (cases II and III) and of the free-running model x M (case I) can be used to quantify the responses of the core flow to the assimilation of surface observations:

$${} \begin{aligned} (\mathcal{M}\text{-}{\mathcal{F}})_{v_{P}} & \equiv \|v_{r}^{\prime M} - v_{r}^{\prime f}\|_{2}\\ &= \left[ \sum_{0\le m\le l}^{L_{M}}\frac{l^{2}(l+1)^{2}}{{r_{c}^{4}}}\left| \frac{\partial v_{l}^{m(M)}}{\partial r}- \frac{\partial v_{l}^{m(f)}}{\partial r}\right|^{2}\right]^{1/2} \end{aligned} $$
((44))
$${} \begin{aligned} ({\mathcal{M}}\text{-}{\mathcal{F}})_{v_{T}} & \equiv \|{\omega_{r}^{M}} - {\omega_{r}^{f}}\|_{2}\\ &= \left[ \sum_{0\le m\le l}^{L_{M}}\frac{l^{2}(l+1)^{2}}{{r_{c}^{4}}}\left| \omega_{l}^{m(M)}-\omega_{l}^{m(f)} \right|^{2}\right]^{1/2} \end{aligned} $$
((45))

In the above equations, ·2 is the L 2−norm (or rms) over the CMB.

In Fig. 5 are the non-dimensional (with the scaling factor 5×10−6 year−1 for dimensional values) v r2 (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 \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) (shown in Fig. 6) and \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\) (shown in Fig. 7) increase in time, i.e., a growing divergence between the forecast state x f and the free-running model state x M.

Fig. 5
figure 5

The non-dimensional v r2 (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

Fig. 6
figure 6

The \(({\mathcal {M}}\text {-}{\mathcal {F}})\) of the poloidal velocity field v r′ as defined in (44). The dashed lines are the results without SV assimilation (case II) and the solid lines are those with the SV assimilation (case III)

Fig. 7
figure 7

Similar to Fig. 6, but for The \(({\mathcal {M}}\text {-}{\mathcal {F}})\) of the toroidal velocity ω r as defined in (45)

From Figs. 6 and 7, we can also observe that \(({\mathcal {M}}\text {-}{\mathcal {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 \({\mathbf {\widetilde {v}}}^{o}\) to increase further \(({\mathcal {M}}\text {-}{\mathcal {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 \(({\mathcal {M}}\text {-}{\mathcal {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 (44-45) with 0≤ml for a given degree l, and their distributions in m, i.e., the summation of the terms in (44-45) with mlL M for a given order m. Since, as shown in Figs. 6 and 7, the differences between the two cases are very small, we can focus only on case III without loss of generality.

In Fig. 8 is the distribution of \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) in the degree l and in Fig. 9 is its distribution in the order m. From the figures, we can find that \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) varies substantially in the spectral spaces. As shown in Fig. 8, 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 Fig. 9, the distribution in m is more broadband: the differences for 5≤m≤35 increase rapidly in time and reach comparable values in magnitude at the end of the assimilation period. However, \(({\mathcal {M}}\text {-}{\mathcal {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).

Fig. 8
figure 8

The distribution of \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) in spherical harmonic degrees l with the SV assimilation (case III)

Fig. 9
figure 9

Similar to Fig. 8, but for the distribution of \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) in spherical harmonic orders m

The distribution of \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\) of the toroidal velocity, as shown in Figs. 10 and 11, displays both similar and distinct characteristics. Its distribution in l is very similar to that of \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\), except that it peaks at a higher degree l=30. But its distribution in m (Fig. 11) 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 \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) (in Fig. 9), \(({\mathcal {M}}\text {-}{\mathcal {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.

Fig. 10
figure 10

Similar to Fig. 8, but for \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\)

Fig. 11
figure 11

Similar to Fig. 9, but for \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\)

Discussions and 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 \({{\widetilde {b}}_{l}^{m(o)}}\) (21) in the D -layer, which is then used in the sequential assimilation algorithm (15). We designed the three experiments 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 \(({\mathcal {O}}\text {-}{\mathcal {F}})\) of the field and SV, defined in (41), at the top of the D -layer are used to measure the forecast accuracies; the (\({\mathcal {M}}\text {-}{\mathcal {F}}\)) of the poloidal velocity field (44) and of the toroidal velocity field (45) 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 Figs. 1 and 2, both \(({\mathcal {O}}\text {-}{\mathcal {F}})_{B}\) and \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) in case III are more than 20 % smaller than those in case II. In particular, the improvements occur to all degrees, as shown in Figs. 3 and 4. The nearly monotonic decay in time of \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) in case III (Fig. 2) 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 \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) and \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\) between the free-running model state and those of the assimilations. As shown in Figs. 6 and 7, these differences grow rapidly in time, showing an accelerated departure of the core state with assimilation from the free-running 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 \(({\mathcal {O}}\text {-}{\mathcal {F}})\) of the magnetic field (Figs. 1 and 2).

Our results have further implications. First, even with the help of (21), the dynamo model is still not fully spun up. For example, though \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) decreases monotonically in time, the SV forecast is still very far away from the observations, as \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B} \approx {\mathcal {O}}(1)\) for all degrees (see Fig. 4). This can be shown further by the continuously growing differences \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) and \({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\) between the forecast velocity field and that of the free-running model (see Figs. 6 and 7). In addition, the differences at the end of the assimilation period are still very small, approximately 10 % in the magnitude of the velocity field of the free-running model (Fig. 5).

These suggest that much larger velocity differences \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) and \(({\mathcal {M}}\text {-}{\mathcal {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 \(\mathbf {v}_{h}^{(o)}\) and \({{\widetilde {b}}_{l}^{m(o)}}\) (more specifically, \(\partial ^{2}{{\widetilde {b}}_{l}^{m(o)}}/\partial r^{2}\)) are unknown beneath the CMB. Thus, the responses of v h to the SV assimilation, e.g., \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) (in Figs. 8 and 9) and \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\) (in Figs. 10 and 11) are needed to determine \(\mathbf {v}_{h}^{(o)}\). For example, as shown in the two figures, the non-axisymmetric (m>0) poloidal velocity \(v_{l}^{m\prime }\) around the degree l=20 and the toroidal velocity \({\omega _{l}^{m}}\) 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 of Aubert (2013, 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 \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\). As shown in Fig. 4, \(({\mathcal {O}}\text {-}{\mathcal {F}})_{\dot B}\) 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 \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) and \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{T}}\) between the forecast velocity field \({\mathbf {v}_{h}^{f}}\) and the \(\mathbf {v}^{M}_{h}\) of the free-running model (see Figs. 8, 9, 10, and 11), from our experiments could also help the inversion of core flow from the observed SV. For example, the different characteristics in \(({\mathcal {M}}\text {-}{\mathcal {F}})_{v_{P}}\) and \(({\mathcal {M}}\text {-}{\mathcal {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. Of course, the core flow response may also suggest that the free-running model solutions do not provide a good description of the small-scale dynamical processes in the core. And assimilation of geomagnetic observations could substantially improve our understanding of these small-scale processes. Regardless, our results suggest that unobservable small-scale dynamical processes are very important in interpretation and prediction of fast SV observable at the Earth’s surface.

Again, we should point out that our results could be improved in several areas: assimilation with other geodynamo models different from (39) and with more sophisticated assimilation algorithms and field models with more accurate error statistics. The former will build up an ensemble of assimilation results to provide more accurate statistics on the characteristics of the core flow to the observed SV. Improvements from the latter are also obvious. For example, we anticipate more accurate estimation of \(({\mathcal {O}}\)-\({\mathcal {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 (40) 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.

Abbreviations

CMB:

core-mantle boundary

GDAS:

geomagnetic data assimilation

ICB:

inner core boundary

SV:

secular variation

References

  • Aubert, J (2013) Flow through the Earth’s core inverted from geomagnetic observations and numerical dynamo models. Geophys. J. Int.192: 537–56.

    Article  Google Scholar 

  • Aubert, J (2014) Earth’s core internal dynamics 1840–2010 imaged by inverse geodynamo modeling. Geophys J Int.197: 1321–34.

    Article  Google Scholar 

  • Aubert, J, Fournier A (2011) Inferring internal properties of Earth’s core dynamics and their evolution from surface observations and a numerical geodynamo. Nonlin Processes Geophys.18: 657–74.

    Article  Google Scholar 

  • Backus, GE (1968) Kinematics of geomagnetic secular variation in a perfectly conducting core. Phil Trans R Soc LondonA263: 239–66.

    Article  Google Scholar 

  • Bloxham, J, Zatman S, Dumberry M (2002) The origin of geomagnetic jerks. Nature420: 65–8.

    Article  Google Scholar 

  • Christensen, U. R, Aubert J, Hulot G (2010) Conditions for Earth-like geodynamo models. Earth Planet Sci Lett.296: 487–96.

    Article  Google Scholar 

  • Christensen, U. R, Wardingsky I, Aubert J, Lesur V (2012) Timescales of geomagnetic secular acceleration in satellite field models and geodynamo models. Geophys J Int.190: 243–54.

    Article  Google Scholar 

  • Evensen, G (1994) Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J Geophys Res.99: 10143–62.

    Article  Google Scholar 

  • Fournier, A, Eymin C, Alboussiere T (2007) A case for variational geomagnetic data assimilation: insights from a one-dimensional, nonlinear, and sparsely observed MHD system. Nonlinear Process Geophys14: 163–80.

    Article  Google Scholar 

  • Fournier, A, Hulot G, Jault D, Kuang W, Tangborn A, Gillet N, et al (2010) An introduction to data assimilation and predictability in geomagnetism. Space Sci Rev. doi:10.1007/s11214-010-9669-4.

  • Fournier, A, Aubert J, Thébault E (2011) Inference on core surface flow from observations and 3-D dynamo modeling. Geophys J Int.186: 118–36.

    Article  Google Scholar 

  • Fournier, A, Nerger L, Aubert J (2013) An ensemble Kalman filter for the time-dependent analysis of the geomagnetic field. Geochem Geophys Geosys. doi:10.1002/ggge.20252.

  • Gillet, N, Jault D, Finlay CC, Olsen N (2013) Stochastic modeling of the Earth’s magnetic field: inversion for covariances over the observatory era. Geochem Geophys Geosys. doi:10.1002/ggge.20041.

  • Glatzmaier, GA, Roberts PH (1995) A three-dimensional self-consistent computer simulation of a geomagnetic field reversal. Nature377: 203–6.

    Article  Google Scholar 

  • Holme, R (2007) Large-scale flow in the core. In: Olson P (ed)Core dynamics, treaties on geophysics.. Elsevier, Amstedam.

    Google Scholar 

  • Hou, AY, Ledvina DV, da Silva AM, Zhang SQ, Joiner J, Atlas RM, et al (2000) Assimilation of SSM/I-derived surface rainfall and total precipitable water for improving the GEOS analysis for climate studies. Mon Weath Rev.128: 509–37.

    Article  Google Scholar 

  • Hulot, G, LeMouël JL (1994) A statistical approach to the Earth’s main magnetic field. Phys Earth Planet Inter.82: 167–83.

    Article  Google Scholar 

  • Jackson, A, Jonkers AR, Walker MR (2000) Our centuries of geomagnetic secular variation from historical records. Phil Trans R Soc Lond.A358: 957–90.

    Article  Google Scholar 

  • Kageyama, A, Sato T (1997) Generation mechanism of a dipole field by a magnetohydrodynamic dynamo. Phys Rev E.55: 4617–26.

    Article  Google Scholar 

  • Kalnay, E (2003) Atmospheric modeling, data assimilation and predictability. Cambridge University Press, Cambridge, U.K.

    Google Scholar 

  • Korte, M, Genevey A, Constable CG, Frank U, Schnepp E (2005) Continuous geomagnetic field models for the past 7 millennia: 1. A new data compilation. Geochem Geophys Geosys. doi:10.1029/2004GC000800.

  • Korte, M, Constable CG, Donadini F, Holme R (2011) Reconstructing the Holocene geomagnetic field. Earth Planet Sci Lett.312: 497–505.

    Article  Google Scholar 

  • Kuang, W, Bloxham J (1997) An Earth-like numerical dynamo model. Nature389: 371–4.

    Article  Google Scholar 

  • Kuang, W, Tangborn A (2011) Interpretation of core field models. In: Mandea M Korte M (eds)Geomagnetic observations and models.. Springer, Heidelberg.

    Google Scholar 

  • Kuang, W, Tangborn A, Jiang W, Liu D, Sun Z, Bloxham J, Wei Z (2008) MoSST DAS: The first generation geomagnetic data assimilation framework. Comm Comp Phys.3: 85–108.

    Google Scholar 

  • Kuang, W, Tangborn A, Wei Z, Sabaka TJ (2009) Constraining a numerical geodynamo model with 100 years of surface observations. Geophys J Int.179: 1458–68.

    Article  Google Scholar 

  • Kuang, W, Wei Z, Holme R, Tangborn A (2010) Prediction of geomagnetic field with data assimilation: a candidate secular variation model for IGRF-11. Earth Planet Space.62: 775–85.

    Article  Google Scholar 

  • Larmor, J (1919) How could a rotating body such as the Sun become a magnet?Rep Br Assoc87: 159–60.

    Google Scholar 

  • Li, D, Jackson A, Livermore PW (2011) Variational data assimilation for the initial value dynamo problem. Phys Rev E. doi:10.1103/PhysRevE.84.056321.

  • Li, D, Jackson A, Livermore PW (2014) Variational data assimilation for a forced, inertia-free magnetohydrodynamic dynamo model. Geophys J Int.199: 1662–76.

    Article  Google Scholar 

  • Liu, D, Tangborn A, Kuang W (2007) Observing system simulation experiments in geomagnetic data assimilation. J Geophys Res.doi:10.1029/2006JB004691.

  • Nilsson, A, Holme R, Korte M, Suttie N, Hill M (2014) Reconstructing Holocene geomagnetic field variation: new methods, models and implications. Geophys J Int. doi:10.1093/gji/ggu120.

  • Olsen, N, Mandea M (2008) Rapidly changing flows in the Earth’s core. Nature Geosci.1: 390–4.

    Article  Google Scholar 

  • Olsen, N, L uhr H, J S. T, Mandea M, Rother M, Tøffner-Clausen L, Choi S (2006) CHAOS-a model of the Earth’s magnetic field derived from CHAMP, Ørsted, and SAC-C magnetic satellite data. Geophys J Int.166: 67–75.

    Article  Google Scholar 

  • Olsen, N, Lühr H, Finlay C. C, Sabaka T. J, Michaelis I, Rauberg J, Tøffner-Clausen L (2014) The CHAOS-4 geomagnetic field model. Geophys J Int.197: 815–27.

    Article  Google Scholar 

  • Roberts, PH (1992) Geomagnetism. Encyclopedia Earth Syst Sci.2: 277–94.

    Google Scholar 

  • Roberts, PH, Scott S (1965) On analysis of the secular variation, 1, A hydromagnetic constraint: theory. J Geomag Geoelec.17: 137–51.

    Article  Google Scholar 

  • Sabaka, TJ, Olsen N, Purucker M (2004) Extending comprehensive models of the Earth’s magnetic field with Ørsted and champ data. Geophys J Int.159: 521–47.

    Article  Google Scholar 

  • Sabaka, T. J, Olsen N, Tyler R. H, Kuvshinov A (2015) CM5, a pre-Swarm comprehensive geomagnetic field model derived from over 12 year of CHAMP, Ørsted, SAC-C and observatory data. Geophys J Int.200: 1596–626.

    Article  Google Scholar 

  • Stacey, FD (1992) Physics of the Earth. Brookfield Press, Kenmore, Australia.

    Google Scholar 

  • Sun, Z, Kuang W (2015) An ensemble algorithm based component for geomagnetic data assimilation. Terr Atmos Ocean Sci. doi:10.3319/TAO.2014.08.19.05.

  • Sun, Z, Tangborn A, Kuang W (2007) Data assimilation in a sparsely observed one-dimensional modeled MHD system. Nonlin Proc Geophys.14: 181–92.

    Article  Google Scholar 

  • Tangborn, A, Kuang W (2015) Geodynamo model and error parameter estimation using geomagnetic data assimilation. Geophys J Int. doi:10.1093/gji/ggu409.

  • Zatman, S, Bloxham J (1997) Torsional oscillations and the magnetic field within the Earth’s core. Nature388: 760–3.

    Article  Google Scholar 

Download references

Acknowledgements

This research is funded by grants from NASA Earth Surface and Interior Program (NNZ09AK70G) and from NSF (EAR-0757880). We also thank NASA NAS for support on high-end scientific computation.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Weijia Kuang.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

WK proposed and designed the study and carried out the assimilation experiments. AT provided error covariance models and mathematical derivations for the sequential assimilation algorithm used in the experiments. Both authors read and approved the final manuscript.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kuang, W., Tangborn, A. Dynamic responses of the Earth’s outer core to assimilation of observed geomagnetic secular variation. Prog. in Earth and Planet. Sci. 2, 40 (2015). https://doi.org/10.1186/s40645-015-0071-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40645-015-0071-4

Keywords