\documentclass{eps} %\documentclass[jgrga]{AGUTeX} %\documentclass[draft,jgrga]{AGUTeX} %\documentclass[article] %%\usepackage[utf8x]{inputenc} %\usepackage{graphics} %\usepackage{multirow} %\usepackage{graphicx} \usepackage{enumerate} \usepackage{amsmath} %\usepackage{amsfonts} %\usepackage{amsbsy} %\usepackage{epsfig} %\usepackage{color} %\usepackage{rotating} \usepackage{lineno} \usepackage{url} \usepackage{stackengine} \stackMath \setcounter{page}{1} \linenumbers*[1] \newcommand{\bR}{\textbf{R}} \newcommand{\Atau}{|\Delta t|} \newcommand{\bF}{\textbf{F}} \newcommand{\bd}{\textbf{d}} \newcommand{\bb}{\textbf{b}} \newcommand{\bx}{\textbf{x}} \newcommand{\by}{\textbf{y}} \newcommand{\bz}{\textbf{z}} \newcommand{\bH}{\textbf{H}} \newcommand{\bG}{\textbf{G}} \newcommand{\bK}{\textbf{K}} \newcommand{\bxi}{\boldsymbol{\xi}} \newcommand{\bsigma}{\boldsymbol{\Sigma}} \newcommand{\red}{\textcolor{red}} \newcommand{\bB}{{\textbf B}} \newcommand{\bg}{{\textbf \gamma}} \newcommand{\br}{B_r} \newcommand{\bu}{{\textbf u}} \newcommand{\bgamma}{{\boldsymbol \gamma}} \newcommand{\blue}{\textcolor{blue}} \newcommand{\MM}{\mathcal{M}} \newcommand{\MME}{\mathcal{M}_E} \newcommand{\MMAR}{\mathcal{M}_\Gamma} \newcommand*{\rom}[1]{\expandafter\@slowromancap\romannumeral #1@} \newcommand{\Cov}{\mathrm{Cov}} \newcommand{\E}{\mathrm{E}} \newcommand{\tbz}{\tilde{\bz}} \newcommand{\bbz}{\bar{\bz}} %\newcommand{\exp}{\mathrm{exp}} \DeclareMathOperator*{\argmax}{arg\,max} %\begin{document} \title{The Kalmag model as a candidate for IGRF-13} \author{Baerenzung Julien, Institute for Mathematics, University of Potsdam, Potsdam, Germany, baerenzung@gmx.de\\ Holschneider Matthias, Institute for Mathematics, University of Potsdam, Potsdam, Germany, matthias.holschneider@gmail.com \\ Wicht Johannes, Max Planck Institute for solar system research, G\"ottingen, Germany, wicht@mps.mpg.de\\ Lesur Vincent, Institut de Physique du Globe, Paris, France, lesur@ipgp.fr \\ Sanchez Sabrina, Max Planck Institute for solar system research, G\"ottingen, Germany, sanchezs@mps.mpg.de} \abstract{ We present a new model of the Geomagnetic field spanning the last 20 years and called Kalmag. Deriving from the assimilation of CHAMP and SWARM vector field measurements, it separates the different contributions to the observable field through parameterized prior covariance matrices. To make the inverse problem numerically feasible it has been sequentialized in time though the combination of a Kalman filter and a smoothing algorithm. The model provides reliable estimates of past, present and future mean fields and associated uncertainties. The version presented here is an update of our IGRF candidates, the amount of assimilated data has been doubled and the considered time window has been extended from $[2000.5,2019.74]$ to $[2000.5,2020.33]$. } \keywords{Geomagnetic field, secular variation, assimilation, Kalman filter, machine learning} \begin{document} \maketitle \section{Introduction} \noindent The Earth's magnetic field has different sources. Classically, we distinguish internal and external sources below and above the site where the field is measured. The three principal internal contributions are the core field, the lithospheric field and the magnetic fields induced in the ocean or within the crust and upper mantle. At large scales, the core field clearly dominates. Sustained by dynamo action in the liquid outer core, it is predominantly dipolar, and varies on timescales ranging from months to millennia. Magnetic field reversals, rare events on even longer times scales, are not considered here. On scales that correspond to spherical harmonics (SH) degree beyond about $16$, the core field is dominated by the lithospheric field coming from the magnetized rocks of the crust. Since the Earth's mantle evolves extremely slowly, the lithospheric field, can be considered as almost static. The external fields that are generated by electrical currents in the ionosphere and in the magnetosphere, on the other hand can vary extremely rapidly in time. Because the sources of the magnetospheric field (the ring current, the magnetopause and magnetotail currents) are distant from the Earth, only its large scale contributions can be detected by the low-orbiting magnetic satellites or by measurements on Earth's surface. This is not the case for the ionospheric field which is generated closer to the Earth's surface. Because the dynamical behavior of both ionospheric and magnetospheric fields is controlled by solar radiations (and also thermospheric winds for the former), they are closely tied to solar activity, and can vary on very short timescales. The external fields can therefore also induce a non negligible secondary field inside the electrically conducting parts of the mantle, crust and oceans. Other potentially important induced fields are created because the oceans move relative to the core fields. \noindent Disentangling the different field contributions is a difficult task since they often overlap in spatial scale and time scale. Many field models therefore resort to a regularization in space and time and only use selected data. Some example are the CHAOS model series by \cite{Olsen2006,Finlay2016}, the comprehensive models by \cite{Sabaka2002,Sabaka2015,Sabaka2018,Sabaka2020}, the GRIMM models by \cite{Lesur2008,Lesur2010,Lesur2015}, the POMME models by \cite{Maus2005,Maus2010}, or the gufm1 by (\cite{Jackson2000}). The COV-OBS model by (\cite{Gillet2013}) and the model recently proposed by \cite{Ropp2020} are the only models that use a Bayesian approach instead of regularization. Usually only vector field measurements taken during night time, under geomagnetic quiet conditions, and at low to mid magnetic latitude are considered in order to minimize the contribution of fields created by currents in the ionosphere or in auroral regions (field-aligned currents, DP2, auroral electrojet). Four main contributions then remain, the core, the lithospheric, the magnetospheric and the induced fields. The core and the lithospheric fields are usually treated as one internal source described by one set of spherical harmonics coefficients. %({\color{red} common only if they are seen with respect to a comon radius of expansion}). Small scale contributions beyond degree $16$ or so are supposed to be of lithospheric origin and are static in time while the larger scales represent the varying core field. Most of the models also treat the induced fields and the magnetetospheric field in a simplified way. The field induced by ocean circulation and the the fields created in the magnetotail and the magnetopause are either neglected or estimated separately. Only the magnetic field generated by the ring current and the respective induced part then remain to be modeled. Yet, when assuming a 1d electrical conductivity profile for Earth's mantle, the axisymmetric magnetospheric field and the related induced field can be parameterized by the Disturbance short-time (Dst) index proposed by \cite{Sugiuara1963}. The Dst or other similar indices are independently estimated from observations and can serve as model input. \noindent To achieve an optimal separation of the different contributions, proper temporal parametrization of the different sources is mandatory. In many models (CHAOS, GRIMM, CM, COV-OBS), the time dependency of the core field is modeled by B-splines for an a priorily fixed time step. This imposes a relatively smooth evolution of the field, excluding the rapid variations attributed to external fields. In addition, with algorithms based on regularized least square approaches, only the time derivatives of the core field are penalized and no constraints on the morphology of the field itself are imposed. With the COV-OBS model, \cite{Gillet2013} went a step further in characterizing a priori the spatio-temporal behavior of the core field. They assumed that its dynamical evolution was controlled by a specific second order auto regressive process which can reproduce the temporal statistical properties of the core field which have been characterized with both observatory measurements (see \cite{DeSantis2003,Lesur2017}) and numerical simulations of the geodynamo (see \cite{Bouligand2016}). \noindent Although a good calibration of the temporal constraints is crucial for deriving magnetic field models from observatory and satellite data, some key information can also be extracted from the morphology and the spatial correlation structure of the different fields. \cite{Holschneider2016} have shown that the use of appropriate parameterized correlation kernels, could greatly improve the separation of the different components of the Earth's magnetic field. Working with observatory data for a single epoch, they could detect the spatial signature of the core, lithospheric, magnetospheric and ionospheric fields. Their Bayesian approach even allowed the quantification of uncertainties. \noindent The Kalmag model we propose here combines such a technique with the sophisticated temporal correlation functions introduced by \cite{Gillet2013}. Since ground based observatories and satellite missions such as Oersted, SAC-C, Champ, or Swarm, have produced or are still generating a huge amount of data, a block inversion would be numerically impossible. This is why we decided to assimilate the data sequentially using a Kalman filter approach combined with a smoothing algorithm. \noindent The article is organised as follows: in the next section, the data selection criteria and the modeling strategy are detailed. We first present the different magnetic sources that are taken into account, we then show how they are a priori characterized, and how such prior information can be modeled through auto regressive processes. Based on these processes, the equations for the Kalman filter approach and the smoothing algorithm are then given. Finally the methodology used to derive the different candidate models for IGRF-13 is explained. In section {\textbf{Results and discussion}}, we present and discuss the outcomes of the model. However, since the spatio-temporal prior characterization of each modeled magnetic source is parameterized, we first show how these parameters are evaluated to be then incorporated in the model. The article ends with some concluding remarks and perspectives to improve the Kalmag model. \section{Methods}\label{modelingStrategy} \subsection{Data}\label{data} For the moment, the Kalmag model only uses the vector field measurements of the CHAMP and SWARM low orbiting satellites. We sample CHAMP data at a rate of $1$ datum every $5$ seconds and only use measurements where the vector field magnetometer (VFM) and the star tracker (STR) instruments were functioning nominally. Very early in the SWARM mission (September $2014$), the scalar magnetometers on satellite Charlie stopped operating properly. We therefore only consider data from the Alpha and Bravo satellites, using simultaneous sampling every 10 seconds. For the construction of the IGRF-13 candidate models, which we will refer to as the Kalmag candidates, a two times lower sampling rate was used. Futhermore, we only used data up to $2019.74$ for the Kalmag candidates but now extended this to $2020.33$. For latitudes between $60^\circ$ north and south, only night time data are considered. Furthermore, independently of the satellites locations, the following selection criteria are also applied: \begin{itemize} \item The $z$-component of the Interplanetary magnetic field(IMF) is positive . \item The geomagnetic activity index $Kp \leq 2^0$. \end{itemize} All in all the dataset is composed of $2\, 985\, 442$ vector field measurements for CHAMP and $4 \,606\, 159$ for SWARM. \subsection{Magnetic sources}\label{magneticSources} \noindent The different contributions to the observations are described in terms of magnetic sources of either internal or external origin. Except for the field produced by field-aligned currents ($b_{fac}$), each of these contributions $b_i$ is deriving from a potential $V_i$: %\begin{linenomath*} \begin{linenomath}\begin{equation} b_i(r,\theta_s,\phi_s,t) = -\nabla V_i(r,\theta_s,\phi_s,t)\ . \end{equation}\end{linenomath} For $b_{fac}$, we followed the study of \cite{Waters2001} and express it through the potential $V_{fac}$ as following: \begin{linenomath}\begin{equation} b_{fac}(r,\theta_s,\phi_s,t) = -{\bf r}\times \nabla V_{fac}(r,\theta_s,\phi_s,t)\ . \end{equation}\end{linenomath} Note that depending on the source, the spherical coordinate system $\{r,\theta_s,\phi_s\}$ the magnetic field is expressed in may differ. Here four types of systems are used: geographic (GEO), magnetic (MAG), solar magnetic (SM), and geocentric solar magnetospheric (GSM). %The various magnetic sources and associated coordinate system are given in table \ref{magneticSourcesTable}. %\end{linenomath*} Each potential $V_i$ is then expanded in spherical harmonics (SH), which for internal and external sources respectively read: %\begin{linenomath*} \begin{linenomath}\begin{eqnarray} V^I_i(r,\theta_s,\phi_s,t) &=& a_i \sum_{\ell\leq \ell_{max}} \sum_{m=-\tilde{m}}^{m=\tilde{m}} \left(\frac{a_i}{r}\right)^{l+1}g_{i,\ell,m}^I(t)Y_{\ell,m}(\theta_s,\phi_s) \ ,\\ V_i^E(r,\theta_s,\phi_s,t) &=& a_i \sum_{\ell\leq \ell_{max}} \sum_{m=-\tilde{m}}^{m=\tilde{m}} \left(\frac{r}{a_i}\right)^{l}g_{i,\ell,m}^E(t)Y_{\ell,m}(\theta_s,\phi_s) \ . \end{eqnarray}\end{linenomath} %\end{linenomath*} The $Y_{\ell,m}$ are Schmidt semi-normalized spherical harmonics of degree $\ell$ and order $m$, $\ell_{max}$ is the maximum of degree expansion, $a_i$ is a reference radius, and $g_{i,\ell,m}(t)$ (later referred as $g_i$) are the spherical harmonics coefficients expressed at $a_i$. $\tilde{m}$ is the maximum order considered for the spherical harmonics expansion. A complete expansion, referred as standard, requires $\tilde{m} = \ell$. However some sources, in particular external fields, are known to have a strong zonal signature (see \cite{Finlay2017}), and are therefore restricted to either zonal spherical harmonics modes with $\tilde{m} = 0$ or to an expansion we refer as zonal iso where $\tilde{m} = 1$. \begin{table*}[h] \caption{Magnetic sources considered in the model. The second column corresponds to the coordinate system each field is expressed in. GEO stands for geographic, SM for solar magnetic, MAG for magnetic and GSM for geocentric solar magnetetospheric. $\ell_{max}$ is the maximum degree of the SH expansion, for the three following types of decomposition: standard with $m=[-l,l]$, zonal with $m=0$ and zonal iso where $m=\{0,1,-1\}$. } \centering \begin{tabular}{| l | p{20mm} | p{25mm} | p{25mm} | } \hline Source & Coordinate & $\ell_{max}$ & SH decompostion \\ \hline Core $g_{c}$& GEO & $20$ & Standard \\ \hline Lithospheric $g_{l}$& GEO & $76$ & Standard \\ \hline Remote magnetospheric $g_{rm}$& GSM & $1$ & Zonal \\ \hline Close magnetospheric $g_{m}$& SM & $15$ & Zonal \\ \hline Fluctuating magnetospheric $g_{fm}$& SM & $15$ & Zonal iso\\ \hline Residual ionospheric/ induced $g_{ii}$ & MAG & $50$ & Zonal iso \\ \hline Field-aligned currents $g_{fac}$ & SM & $15$ & Zonal iso \\ \hline \end{tabular} \label{magneticSourcesTable} \end{table*} \noindent The Kalmag model is composed of $7$ sources. $3$ of them are of internal origin, the core field ($g_c$), the lithospheric field ($g_l$) and the induced / residual ionospheric field ($g_{ii}$). $g_c$ and $g_l$ are expressed in the geographic coordinate system and expanded in SH with the standard decomposition. $g_{ii}$ is expressed in the solar magnetic coordinate system and its SH decomposition is restricted to $\tilde{m} = 1$ (zonal iso). $3$ sources are used to characterize the magnetospheric field. A remote one ($g_{rm}$) in GSM which is purely dipolar and zonal, and $2$ close sources ($g_m$ and $g_{fm}$) expressed in the SM coordinate system. $g_m$ is purely zonal and it is accompanied by $g_{fm}$ a fluctuating source expanded with the zonal iso SH decomposition. Finally, the source associated with field-aligned currents is expressed in the SM coordinate system and restricted to the zonal iso SH expansion. \noindent The nature of the $7$ sources composing the Kalmag model, the coordinate system they are expressed in, and their spherical harmonics truncation level are listed in table \ref{magneticSourcesTable}. \subsection{Prior characterization of spatial and temporal correlations}\label{priorCov} \noindent To obtain an optimal separation of the various contributions to geomagnetic observations, proper prior characterization of the different magnetic sources is mandatory. Following the studies of \cite{Hulot1994,Gillet2013,Holschneider2016}, full space-time covariance matrices are used to characterize each magnetic source $g_i$. Assuming that $E[g_i]=0$ the latter read: %\begin{linenomath*} \begin{linenomath}\begin{equation} E\left[\begin{pmatrix}g_i(t)\\ g_i(t+\Delta t)\end{pmatrix} \left(g_i(t)^T g_i(t+\Delta t)^T\right)\right] = \begin{pmatrix} \Sigma_{g_i}^\infty & c_i(\Delta t)\Sigma_{g_i}^\infty \\ \Sigma_{g_i}^\infty c_i(\Delta t)^T & \Sigma_{g_i}^\infty \\ \end{pmatrix}\ , \label{priorCovariance} \end{equation}\end{linenomath}%\end{linenomath*} \noindent where the matrix $\Sigma_{g_i}^\infty$ corresponds to the stationary spatial covariance, and $c(\Delta t)$ is a temporal correlation matrix depending on the time lag $\Delta t$. $\Sigma_{g_i}^\infty$ is assumed to derive from energy spectra $E_i^\infty(\ell,a_i)$, expressed at given radii $a_i$, such as: \begin{linenomath}\begin{eqnarray} \Sigma_{g_i}^\infty(\ell,m,\ell^\prime,m^\prime, r=a_i) & =& E[g_i(\ell,m)g_i(\ell^\prime,m^\prime)] \nonumber \\ &=& \frac{E_i^\infty(\ell,a_i)}{N_m F(\ell)} \delta(\ell-\ell^\prime)\delta(m-m^\prime) \label{sigmaInf} \end{eqnarray}\end{linenomath} where $N_m$ is the number of modeled spherical harmonics coefficients per degree $\ell$, and $F$ is the pre-factor of the energy spectra given by $F(\ell) = \ell+1$ and $F(\ell) = \ell$ for internal and external sources respectively. Two types of spectra are used for the model, flat ones, with $E_i^\infty(\ell) = A_i^2$ where $A_i$ is a magnitude, and spectra of the form $E_i^\infty(\ell) = A_i^2 (2\ell+1) F(\ell)$, referred as C-based spectra, making equation \ref{sigmaInf} equivalent to the correlation kernels proposed by \cite{Holschneider2016}. Only $2$ sources are characterized by a flat spectrum, the core field and the induced / residual ionospheric field. This choice was driven by the evaluation we performed in section {\textbf{Parameter estimation}} where such a parametrization enabled us to better explain the data. \noindent Note that the covariance matrices of equation \ref{sigmaInf} are diagonal. More complex covariance structures could be used, accounting for the correlations between different magnetic modes as they appear in dynamo simulations for instance (see \cite{Sanchez2019}). However, the large amount of data available for this study are sufficient to properly constrain the model and permits such general prior assumptions. \noindent Temporal constraints are prescribed by the type of correlation functions introduced by \cite{Gillet2013,Gillet2015} in the context of geomagnetic modeling. They are deriving from the auto-regressive processes further discussed below, and read: \begin{linenomath}\begin{equation} c_i(\Delta t) = \exp\left[-\Atau/\tau_i(\ell)\right] \label{AR1corr} \end{equation}\end{linenomath} for first order processes and: \begin{linenomath}\begin{equation} c_i(\Delta t) = \left(1 + \left(\Atau/\tau_i(\ell)\right)\right)\exp\left[-\left(\Atau/\tau_i(\ell)\right)\right]\label{AR2corr}\ \ \end{equation}\end{linenomath} for second order processes. $\tau_i(\ell)$ are scale dependent characteristic timescales. For the core field, \cite{Christensen2004,Lhuillier2011} have shown that its characteristic timescales $\tau_c(\ell)$ could be approximated by a power law such as $\tau_c(\ell) = \tau_{SV} \ell^{-1}$ with $\tau_{SV}$ the secular variation timescale. In this study we decided to use such a power law description of $\tau_i$ for each magnetic source except for the lithospheric field, leading to: \begin{linenomath}\begin{equation} \tau_i(\ell) = M_i \ell^{-\alpha_i} \end{equation}\end{linenomath} with amplitudes $M_i$ and exponent $\alpha_i$. In our parametrization of the problem we therefore use four main parameters to characterize each source: 1) the amplitude $A_i$, 2) the (virtual) source radius $a_i$, 3) the time scale amplitude $M_i$ and 4) the time scale slope $\alpha_i$. In addition, because of the specific behavior of the dipole components, for the core field (see \cite{Christensen2004,Lhuillier2011}) but also for the magnetetospheric sources (see \cite{Sugiuara1963,Finlay2016}), the spatial and temporal properties of each source's dipole (except for the lithopsheric field) is treated separately from the remaining SH coefficients. These parameters are directly estimated with a subsample of the dataset following the procedure described in section {\bf Parameter estimation}. \subsection{Sequentialization}\label{priorCov} For the following developments, the parameters characterizing the prior covariance structures ($A_i$, $a_i$, $M_i$ and $\alpha_i$ ) are assumed to be known. Instead of performing a full Bayesian block inversion with the covariance matrices given by equation \ref{priorCovariance} as a prior information, we proceed in a recursive way through the Kalman filter approach proposed by \cite{Kalman1960}. To do so, dynamical equations are required to forecast the statistical properties of the different modeled sources. As previously mentioned, the covarfiance structures we wish to a priori impose are deriving from autoregressive processes. In their continuous form the latter are given for first and second orders by respectively: \begin{linenomath}\begin{eqnarray} \partial_t g_{i,\ell,m}(t) & + & \frac{1}{\tau_i(\ell)}g_{i,\ell,m}(t) = \sigma_{i_1}(\ell) \dot{\omega}_{i_1}(t) \label{AR1} \\ \partial^2_t g_{i,\ell,m}(t) & + & \frac{2}{\tau_i(\ell)}\partial_tg_{i,\ell,m}(t)+ \frac{1}{\tau_i^2(\ell)}g_{i,\ell,m}(t) = \sigma_{i_2}(\ell) \dot{\omega}_{i_2}(t) \label{AR2} \end{eqnarray}\end{linenomath} where $\dot{\omega}_{i_1}(t)$ and $\dot{\omega}_{i_2}(t)$ are Gaussian white noises scaled by the factors $\sigma_{i_1}(\ell)$ and $\sigma_{i_2}(\ell)$ respectively. These equations have explicit solutions which satisfy: %\begin{linenomath*}\be \begin{linenomath}\begin{equation} z_i(t+\Delta t) = F_i(\Delta t) z_i(t) + \xi_i(t,\Delta t) \end{equation}\end{linenomath}%\end{linenomath*} where the temporal Gaussian white noise $\xi_i$ is spatially (in terms SH coefficients) characterized by the distribution $\mathcal{N}\left({0},\Sigma_{z_i}^\infty - F_i\Sigma_{z_i}^\infty F_i^T \right)$, and where $\Delta t$ can either be positive or negative. \noindent For magnetic sources characterized by first order auto regressive processes, $z_i = g_i$ and $F_i$ is given by: \begin{linenomath}\begin{equation} F_i(\ell,\Delta t) = \exp\left[ -\Atau / \tau_i(\ell) \right] . \end{equation}\end{linenomath} The core field evolution is prescribed by a second order auto regressive process, so the field itself and the secular variation are dynamically tied together. In this case, $z_i = (g_i,\partial_t{g}_i)^T$ and: %\begin{linenomath*} \begin{linenomath}\begin{equation} F_i(\ell,\Delta t) = \begin{pmatrix} 1+\Atau / \tau_i(\ell)& \Delta t \\ -\Delta t/ \tau_i^2(\ell) & 1-\Atau/ \tau_i(\ell) \\ \end{pmatrix} \exp\left[-\Atau / \tau_i(\ell)\right] \ . \end{equation}\end{linenomath} %\end{linenomath*} %where $\zeta (\ell) = \left(\tau_i(\ell)/ \tau_i(\ell)\right)$. Note that the covariance matrix of the core field when the later is assumed to be in a stationary state is given by: \begin{linenomath}\begin{equation} \Sigma_{g_c,\partial_t{g}_c}^\infty = \begin{pmatrix} \Sigma_{g_c}^\infty & 0 \\ 0 & \Sigma_{g_c}^\infty/\tau_c^2(\ell) \\ \end{pmatrix} \ , \end{equation}\end{linenomath} as shown by \cite{Hulot1994}. \subsection{Sequential assimilation}\label{assimilation} \noindent The Kalmag model consists in a vector $\bz$ containing the spherical harmonics coefficients of every magnetic source (including the SH expansion of the secular variation). With the decomposition detailed in section {\textbf{Magnetic sources}}, $\bz$ contains $6624$ SH coefficient entries. Its evaluation is performed with a Kalman filter algorithm which proceeds sequentially in two steps. In the first step, the forecast, the evolution of the mean model $E[\bz]$ together with its associated covariance matrix $\bsigma_{\bz}$ are predicted until observations become available. In the second step, namely the analysis, the model is corrected to better reflect the data through a Bayesian inversion. \noindent To predict the simultaneous evolution of the different magnetic sources with the auto-regressive processes presented in the previous section, a matrix $\bF$ containing all the matrices $F_i$, and a matrix $\tilde{\bsigma} =\bsigma^\infty - \bF \bsigma^\infty \bF^T$ characterizing the white noise of the complete evolution model, are constructed. The evolution of the mean model and its covariance from time step $k-1$ to step $k$ is then given by the forecast: \begin{linenomath}\begin{eqnarray} E[\bz_{k|k-1}] & = & \bF_{k-1} E[\bz_{k-1}] \label{fullForecastMean}\\ \bsigma_{\bz_{k|k-1}} & =& \bF_{k-1} \bsigma_{\bz_{k-1}} \bF_{k-1}^T + \tilde{\bsigma} \label{fullForecastCov}\ . \end{eqnarray}\end{linenomath} At iteration $k$, whenever measurements are available, the model is updated with the formulations: \begin{linenomath}\begin{eqnarray} \bK_k & =& \bsigma_{\bz_{k|k-1}}\bH_k^T\left(\bH_k \bsigma_{\bz_{k|k-1}}\bH_k^T \right)^{-1} \label{KalmanGain}\\ E[\bz_{k|\bd_k}] &=& E[\bz_{k|k-1}] + \bK_k\left(d_k - \bH_k E[\bz_{k|k-1}]\right) \label{analyseMean} \\ \bsigma_{\bz_{k|\bd_k}} & =& \left(\textbf{I} - \bK_{k}\bH_k\right) \bsigma_{\bz_{k|k-1}}\label{analyseCov} \end{eqnarray}\end{linenomath} where $\bK_k$ is the Kalman gain matrix and $\bH_k$ is the operator projecting the model to the data $d_k$ at iteration $k$. \noindent Note that the time step of the algorithm has been set to $\Delta t = 30$ minutes. Within this time window most of the magnetic sources are assumed to be static. However, thespatio-temporal correlations of the FAC source as well as the non dipolar part of the close magnetetospheric field are modeled within this time window for reasons detailed in section \textbf{Parameter estimation}. \subsection{Smoothing}\label{smoothing} With the Kalman filter algorithm, one gets access to the distribution $p(\bz_k|\bd_k)$, where $\bd_k$ corresponds to all the measurements up to iteration $k$. To obtain $p(\bz_k|\bd)$ the posterior distribution of the model at iteration $k$ given the entire dataset $\bd$, one can apply a smoothing algorithm. In this study we chose the formulation of \cite{Rauch1965}. Starting at the last iteration of the Kalman filter algorithm, the smoothing algorithm performs iteratively backward in time accordingly to the following steps: \begin{linenomath}\begin{eqnarray} \bG_{k-1} & = & \bsigma_{\bz_{{k-1}|\bd_{k-1}}} \bF_k^T \bsigma_{\bz_{k|k-1}}^{-1} \\ E[\bz_{{k-1}|\bd}] & = & E[\bz_{{k-1}|\bd_{k-1}}] + \bG_k\left( E[\bz_{k|\bd}] - E[\bz_{k|{k-1}}] \right) \\ \bsigma_{\bz_{{k-1}|\bd}} & = & \bsigma_{\bz_{{k-1}|\bd_{k-1}}} + \bG_{k-1} \left( \bsigma_{\bz_{k|\bd}} - \bsigma_{\bz_{k|{k-1}}}\right) \bG_{k-1}^T \ . \end{eqnarray}\end{linenomath} \noindent The combination Kalman filter-smoothing algorithm was also chosen by \cite{Ropp2020} for their IGRF-13 main field candidate. However, their approach differs from ours in many aspects. In particular, their core field evolution is prescribed by an Euler scheme and they estimate the secular variation through the fluctuation of the field within $3$ month time windows. Here the secular variation is tied to the core field evolution through the AR2 process. Its evaluation is therefore achieved by through dynamical link and correlation with the core field. \subsection{Candidate models}\label{candidate} The models that we proposed as candidates for the IGRF-13 in $2020.0$ are the Kalman filter solutions after the last analysis step in $2019.74$ (September the $27^{th}$). This solution was forwarded in time until $2020.0$ using the forecast of equations \ref{fullForecastMean} and \ref{fullForecastCov} with propagators $\bF$ and noise covariance $\tilde{\bsigma}$ for a time step of $\Delta t=0.26 yr$. The secular variation candidate is the mean secular variation estimation in $2020.0$. The associated uncertainties were obtained by taking the square root of the diagonal elements of the covariance matrix $\Sigma_{\partial_t{g}_c}(t = 2020)$, providing the standard deviation corresponding to each SH coefficients of $\partial_t{g}_c$. Our internal field candidate model in $2020.0$ contains the sum of the mean core field and the mean lithospheric field at this epoch $\left(E[g_c] + E[g_l]\right)$. The uncertainties estimates were derived from the covariance matrix $\Sigma_{g_c+g_l} = \Sigma_{g_c} + \Sigma_{g_l} + \Sigma_{g_{cl}} + \Sigma_{g_{cl}}^T$ in $2020.0$, where $\Sigma_{g_{cl}}$ is the cross covariance between the core field and the lithospheric field. The square root of each diagonal elements of $\Sigma_{g_c+g_l}$ provides the standard deviation associated with $E[g_c] + E[g_l]$. Finally, our candidate for the DGRF $2015.0$ model was constructed as our $2020.0$ internal field model, except that the core and the lithospheric fields were taken from the smoothing solution. The Kalmag model presented below uses additional data from September $27$ $2019$ until April $2020$. \section{Results and discussion}\label{Results} \subsection{Parameter estimation}\label{parameterEstimation} In this section the parameters characterizing the different magnetic sources, the spectra amplitudes $A_i$ and radius $a_i$, the characteristic timescales magnitudes $M_i$ and slopes $\alpha_i$ are evaluated. The spectral resolution as well as the spherical harmonics expansion chosen to model the different fields are a priori imposed (see table \ref{magneticSourcesTable}). The spherical harmonics coefficients of the different sources at the different models time are calculated with our Kalman filter scheme with a subsample of the data between $2001.0$ and $2018.0$. In order to avoid measurements taken by CHAMP or SWARM during strongly magnetically disturbed epochs, and such as no permanent bias due to the static part of the magnetic field generated by the ring current remain in the data, measurements deviating by $60$ nT in intensity from the CHAOS-6 internal field model and a yearly estimation of a degree $1$ external field expressed in the SM coordinate system are removed from the set. After this operation, a sample of $N_{est} = 247\, 453$ vector field measurements regularly spaced in time is kept and used to estimate the parameters for the different sources. \noindent This estimation procedure is initialized with a first guess for each parameter. For internal (or external) sources, radii lower (or larger) than the Earth's radius are chosen. For the external sources we assume a characteristic time scale of one day and set the slopes associated with $\tau_i(\ell)$ to $\alpha_i=0$. The same is used for the induced and the ionispheric field. The lithospheric field, on the other hand, is assumed to be static. As mentioned above, several authors report that the core field time scales are inversely proportional to the spherical harmonics degree (except for the axial dipole), implying $\alpha_c=1$. We start our estimation with an initial guess of $\alpha_c=0$ and $M_c=30$ years and check whether we nevertheless recover the results suggested by the other authors. Given a set of parameters, we perform the Kalman filter assimilation described above with the data subset. Before each analysis step we can calculate how well the model predict the data with the relation: % \begin{linenomath}\begin{eqnarray} % F_k^{pred} & = & - \log{ \left|\bH_k \bsigma_{\bz_{k|k-1}}\bH_k^T \right|} \nonumber \\ % & -& \left(d_k - \bH_k E[\bz_{k|k-1}]\right)^T\left(\bH_k \bsigma_{\bz_{k|k-1}}\bH_k^T \right)^{-1} \nonumber \\ % & \times& \left(d_k - \bH_k E[\bz_{k|k-1}]\right) \ . % \end{eqnarray}\end{linenomath} \begin{linenomath}\begin{equation} F_k^{pred} = - \log{ \left|\bH_k \bsigma_{\bz_{k|k-1}}\bH_k^T \right|} \ - \left(d_k - \bH_k E[\bz_{k|k-1}]\right)^T\left(\bH_k \bsigma_{\bz_{k|k-1}}\bH_k^T \right)^{-1} \left(d_k - \bH_k E[\bz_{k|k-1}]\right) \ . \end{equation}\end{linenomath} Summing $F_k^{pred}$ over all $k$ iterations provides the measure for the model compatibility with the data. We randomly explore the multi-dimensional parameter space, seeking to maximise $\sum_k F_k^{pred}$. The final values from this parameter search are given in table \ref{parameterTable}. Remember that for each source with the exception of the lithopsheric field, we distinguish between the dipole spatial and time scales and the spatial and time scales of the other harmonics. \begin{table*} \caption{Magnetic sources parameters as described in section {\bf magnetic sources}. The prior spatial covariance matrices are deriving from energy spectra expressed at some radii $a_i$ which are are either flat with $E_i^\infty(\ell) = A_i^2$ or of the C-based type with the form $E_i^\infty(\ell) = A_i^2 (2\ell+1) F(\ell)$ where $F(\ell)=\ell+1$ and $F(\ell)=\ell$ for respectively internal and external sources. The characteristic timescales of equations \ref{priorCovariance}, \ref{AR1corr} and \ref{AR2corr} are parameterized by $\tau_i(\ell) = M_i \ell^{-\alpha_i}$.} \centering \begin{tabular}{| l | p{20mm} | p{20mm} | p{20mm} | p{35mm} | p{10mm} | } \hline Field & Spectrum &radius $a$(km) & $A$ (nT) & $M$ & $\alpha$ \\ \hline Core & Flat & $ 3456 $ & D: $ 1.12 \times 10^5$ \newline $ 9.74 \times 10^4$ & $\tau_c(1)$: $ 935 $ yrs\newline $ M(\ell\geq2) = 514 $ yrs& ${}$\newline $1.06 $ \\ \hline Lithospheric & C-Based &$ 6287 $ & $ 0.16 $ & $ \infty $ & $ 0 $ \\ \hline Close magnetospheric & C-Based & $ 12524 $ &D: $ 9.16 $ \newline $ 1.88 $ & $\tau_m(1)$: $ 1.54 $ days \newline $ M(\ell\geq2) = 18 $ min &${}$\newline $ 0$ \\ \hline Remote magnetospheric &C-Based & $ 235570 $ & $ 7.3$ & $ 10.31 $ yrs & $ 0 $ \\ \hline Fluctuating magnetospheric & C-Based &$ 13028 $ & D: $ 3$ \newline $ 4.56 $ & $\tau_{fm}(1)$: $0.36 $ day \newline $\tau_{fm}(2)$: $0.55$ days \newline $ M(\ell\geq3) = 4$ days & ${}$\newline ${}$\newline$ 1.15$ \\ \hline Residual ionospheric/ induced & Flat &$ 6324 $ & D: $ 5.48 $ \newline $ 4.39 $ & $\tau_i(1)$: $ 0.71$ day \newline $ M(\ell\geq2) = 1.76 $ day & ${}$\newline $ 0.93 $ \\ \hline Field-aligned currents & C-Based &$ 7917 $ & D: $ 0 $ \newline $ 1.22 $ & $\tau_{fac}(1)$: $ 0 $\newline $ M(\ell\geq2) = 1 $ min & ${}$\newline$ 0 $ \\ \hline \end{tabular}\label{parameterTable} \end{table*} Figure 1 shows the static energy spectra projected at the Earth's surface that define the spatial covariance structure of equation \ref{sigmaInf} for the optimal parameters. A comparison with the CHAOS-6.9 core field model of \cite{Finlay2016} (black circles) and the LCS-1 lithospheric field model of \cite{Olsen2017} demonstrates the close agreement. The other internal source taken into account in our model is the residual ionospheric/ induced field ($g_{ii}$). It is dipole dominated and exhibits an almost flat spectrum at the Earth's surface as illustrated by the blue line in figure 1. Without the restriction to magnetically quiet data, this source would be much more energetic. This is also the case for the magnetic field generated by field-aligned currents ($g_{fac}$), which reaches a similar amplitude as $g_{ii}$ (dashed line in figure 1). The last sources are the external magnetospheric fields, which we model with a close ($g_m$), a remote ($g_{rm}$) and a fluctuation components ($g_{fm}$). Together they exhibit a strong dipole and their energy spectra are rapidly decaying. The different sources cover a large variety of timescales, ranging from minutes to centuries. For the core field, the characteristic time associated with its non dipolar part reads $\tau_c(\ell) = 514 \ell^{-1.06}$. This power law is close to the slower estimate of \cite{Lhuillier2011} given by $\tau(\ell) = 470 \ell^{-1}$, but suggests about $10\%$ longer time scales. For the core dipole, the estimation algorithm yields $\tau_c(1) = 935$ years. Since we only consider data over a $17$ year period, this estimate is likely not very precise but nevertheless illustrates that the dipole evolves much slower than the other harmonics. The residual ionospheric and induced fields vary very rapidly in comparison, with time scales between one hour and one day, independent of the length scale. Sometimes assumed to be static (see \cite{Olsen2014,Finlay2016}), the remote magnetospheric field has a time scale of $\tau_{rm}\sim10.3$ years in our study, a value close to the solar cycle. The part of the external fields typically associated with the ring current are the degree $\ell=1$ contribution of the close and fluctuating magnetetospheric fields. Whereas the purely zonal part $g_{m}$ exhibits a characteristic timescale of $\tau_m(1) = 1.5$ days, the fluctuating part $g_{fm}$, assumed to have SH order $m=\{0,1,-1\}$ here, varies faster with $\tau_{fm}(1) \sim 8$ hours. For the small scales magnetospheric field, the time scales of the zonal contributions are shorter than those of the degree one contributions. Note that for $g_m$, $\tau_m(\ell \geq 2) = 18$ minutes, a characteristic time lower than the $30$ minutes time step of the Kalman filter algorithm. In such a case, where the estimation of $\tau$ was leading to lower values than the algorithm time step, the slope of the parameterized timescale was set to zero, and the source was only characterized by a spatio-temporal covariance structure of the form given by equation \ref{priorCovariance}. Its evolution from one time step to the other was also treated as a temporal white noise. Finally, the fastest varying source is the field-aligned currents with $\tau_{fac}(\ell) = 1$ minute. Together with the non dipolar part of $g_m$, $g_{fac}$, with its zonal structure and short memory of its past, strongly resembles the observed disturbance along satellite tracks discussed in \cite{Finlay2017}, and generally affecting the construction of small scale lithopsheric field models (see \cite{Thebault2017}). \subsection{Model results}\label{model} The optimal model parameters described in the previous section are fixed in the sequential Kalman filter assimilation. The model seeks to describe the data with the spherical harmonics source coefficients $g_i(t)$, which we call the Kalmag geomagnetic field model. Because the parameters were derived from data at low geomagnetic activity, we also have to restrict the final model data. This is done on the fly by testing how much a forecast differs from the data. Whenever the difference lies outside the $95.4\%$ confidence interval predicted by the slow varying sources (the ones exhibiting some characteristic time larger or equal than a day at some degree $\ell$) the associated data points are dismissed. All in all, $28.6\%$ of the originally selected data (see section {\textbf{Data}}) were dismissed. We recall that the forecast time step is set to $30$ minutes. The entire model (mean and covariance of $\bz$) is stored every $0.25$ year, although outputs could be saved down to every time steps. Figure 2 shows the Kalmag energy spectra at the Earth's surface for the core and the lithospheric field for the epoch $2015.0$. For degree $\ell<=15$, the standard deviations (SD) of both fields are comparable and exceed the mean value of the lithospheric field. Moreover, the SD of the combined field is smaller than the SD of the individual fields. This illustrates that we cannot separate core and lithospheric contribution at these large scales. It also indicates that the prior level of variance of the lithospheric field, as estimated in the previous section, is therefore simply the extrapolation of the small scale stationary spectrum towards the larges sales. Figure 3 compares energy spectra for three types of solutions for the main field (left) and the secular variation (right) in $2015.0$. The Kalman filter solution (thin gray lines and symbols), the solution after the smoothing algorithm (thick black lines and symbols), and a third solution for a $5$ year forecast from $2010.0$ (thin black lines and symbols). Continuous lines show the mean, dashed lines the standard deviation, triangles the differences to the DGRF-13 final field model, and circles the difference to the CHAOS-6.9 SV model. Not surprisingly, the forecast yields the largest uncertainties. The smallest uncertainties are achieved in the smoothed solution, since the smoothing process allows to take information from the future into account. For the field itself, the fact that the differences to the DGRF-13 final model are similar to the model uncertainties, indicates that these uncertainties are reliably estimated. For the secular variation, the predicted uncertainty levels seems to be slightly overestimated, at least for the Kalman filter and the smoothing solutions. The maximum resolution achieved for the SV is $\ell=16$ for the smoothing solution, beyond this value the SD becomes larger that the mean signal. On figure 4 are displayed various estimations of the radial (left), azimuthal (middle) and longitudinal (right) secular variation at the level of several ground based observatories over the period $2000.0-2025.0$. Blue dots correspond to SV estimations deriving from ground based observatory measurements. They are obtained by taking annual differences of the measured magnetic field averaged over $0.1$ years. The black lines are evaluations of the SV through the CHAOS-6.9 model. The blue and yellow lines are respectively the IGRF-13 secular variation and the Kalmag candidate SV. The red area is the Kalmag mean secular variation plus and minus $2$ standard deviation ($\sigma$). Between $2000.6$ and $2020.33$ the outcomes of the smoothing solution are shown whereas outside this time window the secular variation is estimated with the forecast step the Kalman filter. Finally, the red dashed lines are the mean SV $\pm 2\sigma$ coming from $5$ year forecast simulations. The first observation one can make is that whenever the secular variation deriving from observatory data exhibits a smooth evolution, the latter is well reproduced by the Kalmag model. We can also notice that at least until $2019.0$, the CHAOS-6.9 SV is always lying within $95.4\%$ confidence interval ($E\left[ \partial_t B\right] \pm 2 \sigma$) predicted by our model. Because the Kalmag model is only deriving from the CHAMP and SWARM measurements, data are missing between $2010.7$ and $2013.8$. This translates into a global increase of uncertainty predictions as it can clearly be witnessed for the longitudinal component of the SV in Mawson or Tuntungan. However, with the combination of the Kalman filter with the smoothing algorithm, the data gap does not lead to any particular issue to connect the two satellite eras since such an approach enables us to account for any space time correlations. As already shown through the energy spectra of figure 3, the forecast algorithm is quite accurate to predict the future states of the secular variation. The three hindcast simulations covering the periods $2005-2010$, $2010-2015$ and $2015-2020$ are a confirming it. Nevertheless, in particular locations where the SV exhibits rapid variations as in M'Bour, the simple auto regressive dynamics propagating the core field, fails to not only reproduce but also bound the real evolution of the SV. This calls for using more complex forecast models, able for example to account for the nonlinear interactions between the core field and a time dependent outer core flow as in \cite{Barrois2017,Baerenzung2018,Sanchez2019}. For the incoming $5$ years, both the IGRF-13 or our candidate SV models (which are everywhere quite close to one another) are lying well between the $ \pm 2 \sigma$ predicted error bars of the updated Kalmag model. In M'Bour however, recent observations tend to show a rapid increase of the azimuthal component of the SV. If it is not followed by a decrease, core field predictions at this location using the IGRF model may rapidly deviate from reality. The last result analyzed in this study, is a comparison of the different candidates for the IGRF-13 main field in $2020.0$, and the field as it can be evaluated with measurements taken after $2020.0$. As shown with figure 3 and discussed previously, the accuracy of the model deriving from the smoothing algorithm, which takes into account knowledge beyond the epoch of evaluation, is higher than the Kalman filter solution where the model derivation only accounts for previously assimilated data. We could also observe that the longer the forecasts the lower the accuracy of the model. For the construction of the IGRF-13 model, measurements were only available up to maximum $2019.75$, so the difference between the updated Kalmag model (which derives from data assimilated up to $2020.33$) and the various candidates can be considered as the errors of the candidates predictions. These errors are displayed in figure 5 through their energy spectra evaluated at the Earth's surface. Whereas the error spectrum of the Kalmag candidate is drawn with a thick black line, the ones associated with the other candidates are shown with thin gray lines. Because the Kalmag model may exhibit a permanent bias, the difference between the model and the candidate may represent an erroneous evaluation of the error. Therefore the candidates errors were also computed using another model taking recent data into account (up to March 2020), the CHAOS-7.2 model of \cite{Finlay2020} which is also, in its first version, the parent model of the DTU candidate for IGRF-13. The spectra of these error evaluations are shown with dashed lines on figure 5. When compared to the Kalmag model, the Kalmag candidate appears to be the most accurate prediction of the main field in $2020.0$ with an error level lower than every other candidate at any SH degree $\ell$. When the comparison is performed with the CHAOS-7.2 model, the Kalmag candidate globally remains the most precise estimation of the $2020.0$ field up to $\ell=8$. However, at smaller scales the DTU candidate is closer to its parent model CHAOS-7.2, but the level of approximated error of the Kalmag candidate remains extremely low. \section{Conclusion}\label{Conclusion} We presented in this study a new approach to derive a Geomagnetic field model from direct measurements of the Earth's magnetic field. Performing sequentially in time, the Kalmag model, which is the combination of a Kalman filter and a smoothing algorithm, enables us to consider complex prior covariance structure to characterize both spatially and temporally the different magnetic sources composing the observable field. The evaluation of the parameters controlling the statistical properties of each modeled source reveals the large variety of spatial and timescales populating the Earth's magnetic field, and reinforces the idea of treating the assimilation of geomagnetic data sequentially in time. By allowing the presence of a large scale lithospheric field independent from the core field, we could show that with the prior characterization we chose, the two sources could not be separated. Furthermore, although the sum of the two fields can be very accurately estimated, the level of uncertainty associated with each individual source is directly linked to the prior variance of the lithospheric field. This implies a maximum resolution for the core field of spherical harmonics degree $\ell=\sim15$. Its time derivative however can be accurately estimated up to $\ell=16$. Globally, the model provides reliable uncertainty quantification for whether past, present or future field estimates. It also permits, through the spatio temporal correlations a priori imposed, to consistently connect the CHAMP and the SWARM satellite eras. For short term forecasts, as the derivation of the IGRF model requires it, we could observe that our approach can be more accurate than other existing methods. This is certainly due to the fact that the secular variation is estimated through its dynamical correlation with the core field and is not a fit to the past evolution. There is nevertheless still some room for improvement. Considering more physically based dynamical equations to constrain the evolution of the various fields, such as dynamo simulations for the core field, would certainly improve the separation of the different sources, and provide more accurate predictions of future states. The temporal window covered by the model could also be extended by taking data from previous satellite missions but also ground based observatories or magnetic surveys. \section{List of abbreviations} \begin{itemize} \item SH: Spherical harmonics. \item SV: Secular variation. \item SD: Standard deviation. \end{itemize} \section{Funding} This work has been funded by the German Research Foundation (DFG) within the Priority Program SPP1788 ``Dynamic Earth''. \section{Availability of data and materials} Champ data can be downloaded at https://isdc.gfz-potsdam.de/champ-isdc/access-to-the-champ-data/ Swarm data can be downloaded at ftp://swarm-diss.eo.esa.int/Level1b/Entire$\_$mission$\_$data/MAGx$\_$LR/ The $Kp$ index can be downloaded at ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/ The IMF indices can be downloaded at https://spdf.gsfc.nasa.gov/pub/data/omni/low$\_$res$\_$omni/ The model presented here is available upon request. \section{Competing interests} The authors declare that they have no competing interests. % \bibliographystyle{agufull08} % \bibliography{biblio} \section{Author's contributions} Baerenzung Julien produced the Kalmag model. Holschneider Matthias contributed to the theoretical developments. Lesur Vincent provided his expertise on satellite data, and the algorithms to calculate the different coordinate transforms required for the model. Wicht Johannes and Sanchez Sabrina participated to the elaboration of the model requirements, and to the redaction of the manuscript. \acknowledgments This work has been supported by the German Research Foundation (DFG) within the Priority Program SPP1788 ``Dynamic Earth''. \begin{thebibliography}{32} \providecommand{\natexlab}[1]{#1} \expandafter\ifx\csname urlstyle\endcsname\relax \providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else \providecommand{\doi}{doi:\discretionary{}{}{}\begingroup \urlstyle{rm}\Url}\fi \bibitem[{\textit{{B{\"a}renzung} et~al.}(2018)\textit{{B{\"a}renzung}, {Holschneider}, {Wicht}, {Sanchez}, and {Lesur}}}]{Baerenzung2018} {B{\"a}renzung}, J., M.~{Holschneider}, J.~{Wicht}, S.~{Sanchez}, and V.~{Lesur} (2018), {Modeling and Predicting the Short-Term Evolution of the Geomagnetic Field}, \textit{Journal of Geophysical Research (Solid Earth)}, \textit{123}(6), 4539--4560, \doi{10.1029/2017JB015115}. \bibitem[{\textit{Barrois et~al.}(2017)\textit{Barrois, Gillet, and Aubert}}]{Barrois2017} Barrois, O., N.~Gillet, and J.~Aubert (2017), Contributions to the geomagnetic secular variation from a reanalysis of core surface dynamics, \textit{Geophysical Journal International}, \textit{211}(1), 50--68, \doi{10.1093/gji/ggx280}. \bibitem[{\textit{{Bouligand} et~al.}(2016)\textit{{Bouligand}, {Gillet}, {Jault}, {Schaeffer}, {Fournier}, and {Aubert}}}]{Bouligand2016} {Bouligand}, C., N.~{Gillet}, D.~{Jault}, N.~{Schaeffer}, A.~{Fournier}, and J.~{Aubert} (2016), {Frequency spectrum of the geomagnetic field harmonic coefficients from dynamo simulations}, \textit{Geophysical Journal International}, \textit{207}, 1142--1157, \doi{10.1093/gji/ggw326}. \bibitem[{\textit{{Christensen} and {Tilgner}}(2004)}]{Christensen2004} {Christensen}, U.~R., and A.~{Tilgner} (2004), {Power requirement of the geodynamo from ohmic losses in numerical and laboratory dynamos}, \textit{nature}, \textit{429}, 169--171, \doi{10.1038/nature02508}. \bibitem[{\textit{{De Santis} et~al.}(2003)\textit{{De Santis}, {Barraclough}, and {Tozzi}}}]{DeSantis2003} {De Santis}, A., D.~R. {Barraclough}, and R.~{Tozzi} (2003), {Spatial and temporal spectra of the geomagnetic field and their scaling properties}, \textit{Physics of the Earth and Planetary Interiors}, \textit{135}, 125--134, \doi{10.1016/S0031-9201(02)00211-X}. \bibitem[{\textit{{Finlay} et~al.}(2016)\textit{{Finlay}, {Olsen}, {Kotsiaros}, {Gillet}, and {T{\o}ffner-Clausen}}}]{Finlay2016} {Finlay}, C.~C., N.~{Olsen}, S.~{Kotsiaros}, N.~{Gillet}, and L.~{T{\o}ffner-Clausen} (2016), {Recent geomagnetic secular variation from Swarm and ground observatories as estimated in the CHAOS-6 geomagnetic field model}, \textit{Earth, Planets, and Space}, \textit{68}, 112, \doi{10.1186/s40623-016-0486-1}. \bibitem[{\textit{{Finlay} et~al.}(2017)\textit{{Finlay}, {Lesur}, {Th{\'e}bault}, {Vervelidou}, {Morschhauser}, and {Shore}}}]{Finlay2017} {Finlay}, C.~C., V.~{Lesur}, E.~{Th{\'e}bault}, F.~{Vervelidou}, A.~{Morschhauser}, and R.~{Shore} (2017), {Challenges Handling Magnetospheric and Ionospheric Signals in Internal Geomagnetic Field Modelling}, \textit{Space Science Reviews}, \textit{206}(1-4), 157--189, \doi{10.1007/s11214-016-0285-9}. \bibitem[{\textit{{Finlay} et~al.}(2020)\textit{{Finlay}, {Kloss}, {Olsen}, {Hammer}, and {T{\o}ffner-Clausen}}}]{Finlay2020} {Finlay}, C.~C., C.~{Kloss}, N.~{Olsen}, M.~{Hammer}, and L.~{T{\o}ffner-Clausen} (2020), {DTU candidate models for IGRF-13}, \textit{Earth, Planets, and Space}. \bibitem[{\textit{{Gillet} et~al.}(2013)\textit{{Gillet}, {Jault}, {Finlay}, and {Olsen}}}]{Gillet2013} {Gillet}, N., D.~{Jault}, C.~C. {Finlay}, and N.~{Olsen} (2013), {Stochastic modeling of the Earth's magnetic field: Inversion for covariances over the observatory era}, \textit{Geochemistry, Geophysics, Geosystems}, \textit{14}, 766--786, \doi{10.1002/ggge.20041}. \bibitem[{\textit{{Gillet} et~al.}(2015)\textit{{Gillet}, {Jault}, and {Finlay}}}]{Gillet2015} {Gillet}, N., D.~{Jault}, and C.~C. {Finlay} (2015), {Planetary gyre, time-dependent eddies, torsional waves, and equatorial jets at the Earth's core surface}, \textit{Journal of Geophysical Research (Solid Earth)}, \textit{120}, 3991--4013, \doi{10.1002/2014JB011786}. \bibitem[{\textit{{Holschneider} et~al.}(2016)\textit{{Holschneider}, {Lesur}, {Mauerberger}, and {Baerenzung}}}]{Holschneider2016} {Holschneider}, M., V.~{Lesur}, S.~{Mauerberger}, and J.~{Baerenzung} (2016), {Correlation-based modeling and separation of geomagnetic field components}, \textit{Journal of Geophysical Research (Solid Earth)}, \textit{121}, 3142--3160, \doi{10.1002/2015JB012629}. \bibitem[{\textit{{Hulot} and {Le Mou{\"e}l}}(1994)}]{Hulot1994} {Hulot}, G., and J.~L. {Le Mou{\"e}l} (1994), {A statistical approach to the Earth's main magnetic field}, \textit{Physics of the Earth and Planetary Interiors}, \textit{82}(3-4), 167--183, \doi{10.1016/0031-9201(94)90070-1}. \bibitem[{\textit{{Jackson} et~al.}(2000)\textit{{Jackson}, {Jonkers}, and {Walker}}}]{Jackson2000} {Jackson}, A., A.~R.~T. {Jonkers}, and M.~R. {Walker} (2000), {Four centuries of geomagnetic secular variation from historical records}, in \textit{Astronomy, physics and chemistry of H$^{+}$$_{3}$}, \textit{Philosophical Transactions of the Royal Society of London Series A}, vol. 358, p. 957, \doi{10.1098/rsta.2000.0569}. \bibitem[{\textit{{Kalman}}(1960)}]{Kalman1960} {Kalman}, R.~E. (1960), { A New Approach to Linear Filtering and Prediction Problems}, \textit{Journal of Basic Engineering}, \textit{82}, 35--45, \doi{10.1115/1.3662552}. \bibitem[{\textit{{Lesur} et~al.}(2008)\textit{{Lesur}, {Wardinski}, {Rother}, and {Mandea}}}]{Lesur2008} {Lesur}, V., I.~{Wardinski}, M.~{Rother}, and M.~{Mandea} (2008), {GRIMM: the GFZ Reference Internal Magnetic Model based on vector satellite and observatory data}, \textit{Geophysical Journal International}, \textit{173}, 382--394, \doi{10.1111/j.1365-246X.2008.03724.x}. \bibitem[{\textit{{Lesur} et~al.}(2010)\textit{{Lesur}, {Wardinski}, {Hamoudi}, and {Rother}}}]{Lesur2010} {Lesur}, V., I.~{Wardinski}, M.~{Hamoudi}, and M.~{Rother} (2010), {The second generation of the GFZ Reference Internal Magnetic Model: GRIMM-2}, \textit{Earth, Planets, and Space}, \textit{62}, 765--773, \doi{10.5047/eps.2010.07.007}. \bibitem[{\textit{{Lesur} et~al.}(2015)\textit{{Lesur}, {Whaler}, and {Wardinski}}}]{Lesur2015} {Lesur}, V., K.~{Whaler}, and I.~{Wardinski} (2015), {Are geomagnetic data consistent with stably stratified flow at the core-mantle boundary?}, \textit{Geophysical Journal International}, \textit{201}, 929--946, \doi{10.1093/gji/ggv031}. \bibitem[{\textit{{Lesur} et~al.}(2017)\textit{{Lesur}, {Wardinski}, {Baerenzung}, and {Holschneider}}}]{Lesur2017} {Lesur}, V., I.~{Wardinski}, J.~{Baerenzung}, and {Holschneider} (2017), {On the frequency spectra of the core magnetic field Gauss coefficients}, \textit{Physics of the Earth and Planetary Interiors}, \doi{https://doi.org/10.1016/j.pepi.2017.05.017}. \bibitem[{\textit{{Lhuillier} et~al.}(2011)\textit{{Lhuillier}, {Aubert}, and {Hulot}}}]{Lhuillier2011} {Lhuillier}, F., J.~{Aubert}, and G.~{Hulot} (2011), {Earth's dynamo limit of predictability controlled by magnetic dissipation}, \textit{Geophysical Journal International}, \textit{186}, 492--508, \doi{10.1111/j.1365-246X.2011.05081.x}. \bibitem[{\textit{{Maus} et~al.}(2005)\textit{{Maus}, {L{\"u}hr}, {Balasis}, {Rother}, and {Mandea}}}]{Maus2005} {Maus}, S., H.~{L{\"u}hr}, G.~{Balasis}, M.~{Rother}, and M.~{Mandea} (2005), \textit{{Introducing POMME, the POtsdam Magnetic Model of the Earth}}, p. 293, \doi{10.1007/3-540-26800-6_46}. \bibitem[{\textit{{Maus} et~al.}(2010)\textit{{Maus}, {Manoj}, {Rauberg}, {Michaelis}, and {L{\"u}hr}}}]{Maus2010} {Maus}, S., C.~{Manoj}, J.~{Rauberg}, I.~{Michaelis}, and H.~{L{\"u}hr} (2010), {NOAA/NGDC candidate models for the 11th generation International Geomagnetic Reference Field and the concurrent release of the 6th generation Pomme magnetic model}, \textit{Earth, Planets, and Space}, \textit{62}, 729--735, \doi{10.5047/eps.2010.07.006}. \bibitem[{\textit{{Olsen} et~al.}(2006)\textit{{Olsen}, {L{\"u}hr}, {Sabaka}, {Mandea}, {Rother}, {Toeffner-Clausen}, and {Choi}}}]{Olsen2006} {Olsen}, N., H.~{L{\"u}hr}, T.~J. {Sabaka}, M.~{Mandea}, M.~{Rother}, L.~{Toeffner-Clausen}, and S.~{Choi} (2006), {CHAOS-a model of the Earth's magnetic field derived from CHAMP, Oersted, and SAC-C magnetic satellite data}, \textit{Geophysical Journal International}, \textit{166}, 67--75, \doi{10.1111/j.1365-246X.2006.02959.x}. \bibitem[{\textit{{Olsen} et~al.}(2014)\textit{{Olsen}, {L{\"u}hr}, {Finlay}, {Sabaka}, {Michaelis}, {Rauberg}, and {T{\o}ffner-Clausen}}}]{Olsen2014} {Olsen}, N., H.~{L{\"u}hr}, C.~C. {Finlay}, T.~J. {Sabaka}, I.~{Michaelis}, J.~{Rauberg}, and L.~{T{\o}ffner-Clausen} (2014), {The CHAOS-4 geomagnetic field model}, \textit{Geophysical Journal International}, \textit{197}, 815--827, \doi{10.1093/gji/ggu033}. \bibitem[{\textit{{Olsen} et~al.}(2017)\textit{{Olsen}, {Ravat}, {Finlay}, and {Kother}}}]{Olsen2017} {Olsen}, N., D.~{Ravat}, C.~C. {Finlay}, and L.~K. {Kother} (2017), {LCS-1: a high-resolution global model of the lithospheric magnetic field derived from CHAMP and Swarm satellite observations}, \textit{Geophysical Journal International}, \textit{211}(3), 1461--1477, \doi{10.1093/gji/ggx381}. \bibitem[{\textit{{Rauch} et~al.}(1965)\textit{{Rauch}, {Striebel}, and {Tung}}}]{Rauch1965} {Rauch}, H.~E., C.~T. {Striebel}, and F.~{Tung} (1965), {Maximum likelihood estimates of linear dynamic systems}, \textit{AIAA Journal}, \textit{3}(8), 1445--1450, \doi{10.2514/3.3166}. \bibitem[{\textit{{Ropp} et~al.}(2020)\textit{{Ropp}, {Lesur}, {Baerenzung}, and {Holschneider}}}]{Ropp2020} {Ropp}, G., V.~{Lesur}, J.~{Baerenzung}, and M.~{Holschneider} (2020), {Sequential modelling of the Earth’s core magnetic field}, \textit{Earth, Planets, and Space}. \bibitem[{\textit{{Sabaka} et~al.}(2002)\textit{{Sabaka}, {Olsen}, and {Langel}}}]{Sabaka2002} {Sabaka}, T.~J., N.~{Olsen}, and R.~A. {Langel} (2002), {A comprehensive model of the quiet-time, near-Earth magnetic field: phase 3}, \textit{Geophysical Journal International}, \textit{151}, 32--68, \doi{10.1046/j.1365-246X.2002.01774.x}. \bibitem[{\textit{{Sabaka} et~al.}(2015)\textit{{Sabaka}, {Olsen}, {Tyler}, and {Kuvshinov}}}]{Sabaka2015} {Sabaka}, T.~J., N.~{Olsen}, R.~H. {Tyler}, and A.~{Kuvshinov} (2015), {CM5, a pre-Swarm comprehensive geomagnetic field model derived from over 12 yr of CHAMP, Oersted, SAC-C and observatory data}, \textit{Geophysical Journal International}, \textit{200}, 1596--1626, \doi{10.1093/gji/ggu493}. \bibitem[{\textit{{Sabaka} et~al.}(2018)\textit{{Sabaka}, {T{\o}ffner-Clausen}, {Olsen}, and {Finlay}}}]{Sabaka2018} {Sabaka}, T.~J., L.~{T{\o}ffner-Clausen}, N.~{Olsen}, and C.~C. {Finlay} (2018), {A comprehensive model of Earth's magnetic field determined from 4 years of Swarm satellite observations}, \textit{Earth, Planets, and Space}, \textit{70}, 130, \doi{10.1186/s40623-018-0896-3}. \bibitem[{\textit{{Sabaka} et~al.}(2020)\textit{{Sabaka}, {T{\o}ffner-Clausen}, {Olsen}, and {Finlay}}}]{Sabaka2020} {Sabaka}, T.~J., L.~{T{\o}ffner-Clausen}, N.~{Olsen}, and C.~C. {Finlay} (2020), {CM6: a comprehensive geomagnetic field model derived from both CHAMP and Swarm satellite observations}, \textit{Earth, Planets, and Space}, \textit{72}(1), 80, \doi{10.1186/s40623-020-01210-5}. \bibitem[{\textit{{Sanchez} et~al.}(2019)\textit{{Sanchez}, {Wicht}, {B{\"a}renzung}, and {Holschneider}}}]{Sanchez2019} {Sanchez}, S., J.~{Wicht}, J.~{B{\"a}renzung}, and M.~{Holschneider} (2019), {Sequential assimilation of geomagnetic observations: perspectives for the reconstruction and prediction of core dynamics}, \textit{Geophysical Journal International}, \textit{217}(2), 1434--1450, \doi{10.1093/gji/ggz090}. \bibitem[{\textit{{Sugiuara}}(1963)}]{Sugiuara1963} {Sugiuara}, M. (1963), \textit{{Hourly values of equatorial Dst for the IGY}}, Greenbelt, Md, NASA, Goddard Space Flight Center. \bibitem[{\textit{{Th{\'e}bault} et~al.}(2017)\textit{{Th{\'e}bault}, {Lesur}, {Kauristie}, and {Shore}}}]{Thebault2017} {Th{\'e}bault}, E., V.~{Lesur}, K.~{Kauristie}, and R.~{Shore} (2017), {Magnetic Field Data Correction in Space for Modelling the Lithospheric Magnetic Field}, \textit{Space Science Reviews}, \textit{206}(1-4), 191--223, \doi{10.1007/s11214-016-0309-5}. \bibitem[{\textit{{Waters} et~al.}(2001)\textit{{Waters}, {Anderson}, and {Liou}}}]{Waters2001} {Waters}, C.~L., B.~J. {Anderson}, and K.~{Liou} (2001), {Estimation of global field aligned currents using the iridium{\textregistered} System magnetometer data}, \textit{Geophysical Research Letters}, \textit{28}(11), 2165--2168, \doi{10.1029/2000GL012725}. \end{thebibliography} %\section{Preparing illustrations and figures} % % \begin{figure}[h!] % \begin{center} % \includegraphics[width=0.8\linewidth]{./figure1.eps} % \caption{Stationary state energy spectra at the Earth's surface of the different magnetic sources given in table \ref{magneticSourcesTable}. % Black and red circles are the energy spectra of respectively the CHAOS-6.9 core field model of \cite{Finlay2016} % in $2015.0$ and the LCS-1 lithopsheric field model of \cite{Olsen2017}. % } % \label{priorSpectra} % \end{center} % \end{figure} % % \begin{figure}[!ht] % \begin{center} % \includegraphics[width=0.95\linewidth]{./figure2.eps} % \caption{Energy spectra at the Earth's surface in $2015.0$ of the mean (continuous lines) standard deviation (dashed lines) % and initial prior standard deviation (circles) associated with the core field (thin black lines and symbols), the lithopsheric field % (gray lines and symbols) and the sum of the core and lithopsheric fields (thick black lines). % }\label{internalSpectra} % \end{center} % \end{figure} % % \begin{figure}[!ht] % \begin{center} % \includegraphics[width=0.95\linewidth]{./figure3.eps} % \caption{ Energy spectra at the Earth's surface in $2015.0$ associated with the main field (left) and the secular variation (right). Thin gray lines and symbols and thick black lines and symbols, correspond to the solutions of the Kalman filter and the smoothing algorithms respectively. The thin black lines and symbols are associated with the outcomes of a $2010.0-2015.0$ forecast simulation. Continuous and dashed lines are assigned to the energy spectra of respectively the different mean solutions and their associated standard deviation. Triangles and circles represent the spectra of the difference between the Kalmag mean models and the DGRF model in $2015.0$ for the main field and the CHAOS-6.9 model for the secular variation respectively. }\label{spectraComparison} % \end{center} % \end{figure} % % \begin{figure}[!ht] % \begin{center} % \includegraphics[width=0.95\linewidth]{./figure4.eps} % \caption{Time series between $2000.0$ and $2025.0$ of the radial (left) azimuthal (middle) and longitudinal (right) secular variation evaluated at the level of the following ground based observatories: Qaanaaq (THL), Niemegk (NGK), M'Bour (MBO), Tuntungan (TUN), Hermanus (HER) and Mawson (MAW).}\label{SVobs} % \end{center} % \end{figure} % % \begin{figure}[!ht] % \begin{center} % \includegraphics[width=0.95\linewidth]{./figure5.eps} % \caption{Energy spectra at the Earth's surface of the difference between every IGRF-13 candidate model, and the latest version of the Kalmag model (continuous lines) as well as the CHAOS-7.2 model of \cite{Finlay2020} (dashed lines) in $2020.0$. The spectra associated with the Kalmag candidate are highlighted by thick black lines.}\label{IGRFcomparisons} % \end{center} % \end{figure} \end{document}