Kalmag: a high spatio-temporal model of the geomagnetic field

We present the extension of the Kalmag model, proposed as a candidate for IGRF-13, to the twentieth century. The dataset serving its derivation has been complemented by new measurements coming from satellites, ground-based observatories and land, marine and airborne surveys. As its predecessor, this version is derived from a combination of a Kalman filter and a smoothing algorithm, providing mean models and associated uncertainties. These quantities permit a precise estimation of locations where mean solutions can be considered as reliable or not. The temporal resolution of the core field and the secular variation was set to 0.1 year over the 122 years the model is spanning. Nevertheless, it can be shown through ensembles a posteriori sampled, that this resolution can be effectively achieved only by a limited amount of spatial scales and during certain time periods. Unsurprisingly, highest accuracy in both space and time of the core field and the secular variation is achieved during the CHAMP and Swarm era. In this version of Kalmag, a particular effort was made for resolving the small-scale lithospheric field. Under specific statistical assumptions, the latter was modeled up to spherical harmonic degree and order 1000, and signal from both satellite and survey measurements contributed to its development. External and induced fields were jointly estimated with the rest of the model. We show that their large scales could be accurately extracted from direct measurements whenever the latter exhibit a sufficiently high temporal coverage. Temporally resolving these fields down to 3 hours during the CHAMP and Swarm missions, gave us access to the link between induced and magnetospheric fields. In particular, the period dependence of the driving signal on the induced one could be directly observed. The model is available through various physical and statistical quantities on a dedicated website at https://ionocovar.agnld.uni-potsdam.de/Kalmag/.


Introduction
Separating the different contributions to the Earth's magnetic field from direct measurements of it is a difficult task.The main reason making this problem complex is the wide range of spatial and temporal scales overlapping one another.The core field, which is sustained by dynamo action in the Earth's outer core, is at the Earth's surface the dominant large-scale field, and it evolves on timescale ranging from months to millennia.On the opposite, the lithospheric field is dominant at small scales.Emanating from the remnant magnetization of the rocks lying within the crust, it follows the motions of the latter and therefore varies very slowly with time.External sources, such as the magnetospheric fields or the ionospheric field are driven by thermospheric winds and solar radiations.Their direct link to solar activity make them subject to intense variations from very short up to decadal timescales.These fluctuations induce currents within the electrically conducting parts of the crust and the mantle which in return generate a secondary magnetic field.Induction processes also occurs within the oceans.The circulation or tidal motions of the latter within the ambient magnetic field create electrical currents which also produce a secondary field.
From the seventeenth century to today, geomagnetic data have been continuously accumulated.First collected during marine and land surveys, measurements of the Earth's magnetic field were quickly complemented by instrumentation installed within groundbased observatories [see Jackson and Finlay (2007)].The development of aviation in the 1950s offered another support to measure the field.But the biggest step in geomagnetic monitoring certainly comes from the rise of low-orbiting satellite missions.Starting in 1965 with the POGO mission, many spacecrafts dedicated to geomagnetic field modeling were later launched.These include non-exhaustively the Mag-Sat, the Oersted and the CHAMP spacecrafts and the Swarm constellation.
Technical constraints to build geomagnetic field models strongly depend on the type of data to be assimilated.Satellite missions provide measurements at a high frequency.The algorithms they are feeding therefore need to be adapted to treat a large amount of observations.Land, marine and airborne surveys operate at the level or slightly above the Earth's surface.As a consequence, the contribution of the small-scale lithospheric field to the Graphical abstract data they produce is important.Accounting for this field requires to model it at a very high resolution, a technical challenge.
Many models of the geomagnetic field have been proposed over the last decades [see Hulot et al. (2015)], and most of them were obtained with a regularized least square approach.This is the case for the CHAOS model series from Olsen et al. (2006) to Finlay et al. (2020), the comprehensive models by Sabaka et al. (2002Sabaka et al. ( , 2015Sabaka et al. ( , 2018Sabaka et al. ( , 2020)), the GRIMM models by Lesur et al. (2008Lesur et al. ( , 2010Lesur et al. ( , 2015)), the POMME models by Maus et al. (2005Maus et al. ( , 2010)), or the gufm1 by Jackson et al. (2000).Least square methods are very efficient numerically, but the usually considered reweighed version can only provide unique solution.On the opposite, Bayesian inversions are computationally demanding but results are expressed in terms of posterior distributions, providing therefore predictions of mean solutions together with their associated uncertainties.Bayesian inversion in the context of geomagnetic field modeling was initiated by Gillet et al. (2013).Considering ground-based observatory, survey and satellite data, they could derive the COV-OBS model spanning the 1840 − 2010 time win- dow, a period which was recently increased to 2020 by Huder et al. (2020).Similar efforts have been followed by Holschneider et al. (2016) in a study where emphasis was put on better characterizing the spatial properties of the different magnetic sources through correlation kernels.Extending this work to the time domain, and sequentializing the problem, Baerenzung et al. (2020), Ropp et al. (2020) could derive geomagnetic field models from the combination of a Kalman filter and a smoothing algorithm.This approach conserves all the advantages of the Bayesian method proposed by Gillet et al. (2013) and alleviates most of its drawbacks.In particular, the dimension of the system, the amount of observations to be assimilated, or the non linear link between certain magnetic sources, are not anymore strong limiting factors.
In this paper, we present the extension of the Kalmag model by Baerenzung et al. (2020) to the twentieth century.Deriving only from CHAMP and Swarm data, Kalmag covered the 2000.5 − 2020 time period, and was a candidate for the IGRF-13 model [see Alken et al. (2021)].The present version resulted from the assimilation of extra measurements taken by ground-based observatories, POGO, MagSat and Oersted satellites and during land, airborne and marine (L.A.M.) surveys.To assimilate the latter type of data, we introduced a statistical approximation within the Kalman filter algorithm enabling us to resolve the lithospheric field up to spherical harmonics degree and order ℓ = 1000 .Therefore, there is no need to subtract the lithospheric contribution to L.A.M. survey observations with high-resolution models such as the WDMAM by Lesur et al. (2016), the EMM model by Maus (2010) or the recent model of Thébault et al. (2021), to build the model.In addition, a small-scale lithospheric field model could be recovered without preprocessing of the data.
The article is organized as follows.In the first part, the dataset used to construct the model and the selection criteria applied are presented.In the second part, the different magnetic sources, their prior characterization and dynamical behavior are detailed.At the end of this section, the various formulations to assimilate data, update the model and sample it are provided.In "Results" section, the properties of our model for the core field, the secular variation, the lithospheric field and external and induced fields are discussed.The article ends with a discussion and some concluding remarks.

Data
The proposed model was derived from either vector field or intensity measurements of the geomagnetic field taken from 1900.0 to today by satellites, ground-based observatories and during land, airborne and marine surveys.Satellite observations from five different missions were considered.These are, the POGO (1965POGO ( -1971) (e.g., Cain and Sweeney 1973), the MagSat (1979MagSat ( -1980) (e.g., Langel and Estes 1985a), the Oersted (since 1999) (e.g., Neubert et al. 2001), the CHAMP (2000-2010) (e.g., Rother et al. 2000), and the SWARM (since 2013) (e.g., Olsen et al. 2013) missions.For ground-based observatories, hourly mean vector fields provided by the World data center for geomagnetism from 1886 (e.g., Macmillan and Olsen 2013) and selected through the procedure which is detailed in the following, were used to derive secular variation data, the latter being used only to constrain the core field evolution.These types of observations, feeding also other models such as the CHAOS series by Olsen et al. (2006), Finlay et al. (2020), the C3FM by Wardinski and Holme (2011), Wardinski et al. (2020) or the COV-OBS model by Gillet et al. (2013), Huder et al. (2020) were here obtained by first averaging vector field measurements over 0.1-year time windows.The resulting mean values b(t) were then used to derive secular variation data γ (t) The location of each observatory taken into account is displayed with black triangles in Fig. 1.For aeromagnetic, land and marine survey data, three compilations served the model derivation (e.g., Quesnel et al. 2009).The first one is provided by British Geological Survey at www. wdc.bgs.ac.uk/, the second one by the National Oceanic and Atmospheric Administration at maps.ngdc.noaa.
gov and the third one is made accessible by the U.S. geophysical survey at www. mrdata.usgs.gov.The positions of L.A.M. data are typically given through the latitude, longitude and altitude location of the measuring vessel.For airborne measurements, whenever altitude was provided by radar altimeter it was corrected above land surfaces with the ETOPO1 global relief model of Amante and Eakins (2009).
Before being assimilated, each data containing vector information, such as North, East, Down or declination, inclination and intensity, was projected in geographic spherical coordinates.The resulting dataset was then subject to selection.The main purposes of this procedure are to avoid the contribution of the dayside ionospheric field which is not modeled, to operating during low geomagnetic activity and, for satellite observations, to be weakly perturbed by the substorm auroral electrojet.The latter two criteria were fulfilled through a selection based on the values of independently derived indices, respectively, a given threshold on the Kp geomagnetic index and the required positiveness of the z-component of the interplanetary magnetic field (IMF).The Kp threshold was set to 2 − for satellite data and to 4 − for all other observations.To limit the contri- bution of the dayside ionospheric field, only night-time measurements (when the sun is below the horizon) were kept at magnetic latitude lying between ±60 • .This constraint was nevertheless relaxed for MagSat data for which the satellite followed a dawn-dusk orbit and for some land survey data which were either not dated precisely enough to determine their local solar time, or only used to derive the lithospheric field model.Note also that for CHAMP and SWARM satellites, it was also required that measurements were taken when both the vector field magnetometer and the star tracker were functioning in nominal mode.
Finally, each L.A.M. surveys and satellite dataset were subsampled.For POGO, MagSat, and Oersted satellites, a rate of 1 datum every 10s (0.1Hz) was chosen.For CHAMP satellite, the sampling rate was increased to 0.2Hz.For SWARM, only satellites Alpha and Bravo are considered with a simultaneous sampling rate of 0.1Hz.Distance criteria were applied to subsample L.A.M. surveys data.In a first selection, a minimum distance of 5 km between any data point within 1-h time windows was imposed.Every measure lying too close to the previously selected ones were removed.The resulting dataset was then split in 8 subsets in which the minimum distance was set to 40 km.Therefore, at a given epoch within the Kalman filter algorithm, data from each of these subsets were sequentially assimilated whenever they were available.In Table 1, the time period, the selection criteria and the type and total number of measurements associated with each dataset are summarized.

Magnetic sources
Seven sources compose the Kalmag model.These are a core field ( b c ), a lithospheric field ( b l ), an induced/resid- ual ionospheric field ( b ii ), a remote ( b rm ), a close ( b m ) and a fluctuating ( b fm ) magnetospheric fields and a source associated with field-aligned currents ( b fac ).Except for b fac , each of these sources b s is assumed to derive from a potential V s such as b s = −∇V s .For b fac , as in Sabaka et al. (2004) the currents themselves are assumed to derive from a potential V fac .Waters et al. (2001) has shown that under this assumption the resulting magnetic field could be expressed as b s = −r × ∇V s .
The potentials V s are then expanded in spherical har- monics (SH) such as potentials of internal and external origin are, respectively, given by: Where Y ℓ,m are Schmidt semi-normalized spherical har- monics of degree ℓ and order m considered, respectively, up to ℓ max and m max , a s is a reference radius, and g s,ℓ,m (t) (later referred as g s ) are the spherical harmonics coeffi- cients expressed at a s .Each field is projected in a given spherical coordinate system {r, θ s , φ s } as indicated in (1) Table 2.These systems can either be geographic (GEO), magnetic (MAG), solar magnetic (SM), or geocentric solar magnetospheric (GSM) (see Laundal ( 2017)).Depending on the observations which are being assimilated, the spatial resolution of the lithospheric field is varied.Whereas for CHAMP and Swarm data, the latter is expanded up to ℓ = 150 , it is only modeled up to ℓ = 100 for other satellite measurements.Since L.A.M. survey data are taken close to the Earth's surface, they contain a strong contribution of the small scale lithospheric field.To assimilate such measurements, the lithospheric field is therefore parameterized up to spherical harmonics degree ℓ = 1000 with an approximation of the associated covariance matrix between 100 < ℓ ≤ 1000 as detailed in the following.

Sequential modeling
The Kalmag model is constructed sequentially through a Kalman filter approach [see Kalman (1960)].This technique proceeds in two alternating steps, namely a forecast and an analysis.In the forecast, the model is propagated in space and time until some measurements become available.Then the analysis takes place and the model is updated accordingly to them.Because this method provides the posterior distribution of the model only given the previously assimilated data, it is complemented by a smoothing algorithm.Performing backward in time, this algorithm enables us to correct the model at any time according to the complete dataset.

Dynamical model
The spatio-temporal evolution of the various sources composing the geomagnetic field is of complex nature.
Involving nonlinear couplings, a large range of spatial and temporal scales, some regimes which are not yet numerically achievable or simply not sufficiently well characterized, the dynamics of the Earth's magnetic field cannot be directly simulated.This is why, as initiated by Gillet et al. (2013) in the context of geomagnetic modeling, we chose simplified stochastic equations, namely auto-regressive processes (or ARPs), to predict the evolution of the different fields.Mimicking dispersion and memory effects occurring within dynamical systems, such processes are computationally cheap to simulate and are formulated within a Gaussian framework as required by the Kalman filter approach.A priori, each source is characterized by its own process which is independent from the others.As shown in Appendix A, ARPs in their sequential form, can be described by the following general relation: where z s is a quantity characterizing the sth source to be propagated, F s (�t) is the parameter of the ARP and ξ i (t, �t) is a temporal Gaussian white noise spatially characterized by the distribution N 0, z s is the stationary state covariance matrix associated with z s .Except for the lithospheric field which is assumed to be static, and for the core field which evolution is prescribed by a second-order process, the dynamics of each source is controlled by a first-order ARP.In this case, z s (t) simply corresponds to the vector of SH coefficients g s (t) associated with the sth field and the parameter of the process is given by: (3) where τ s (ℓ) is a parameterized scale-dependent charac- teristic time which is specified for each source in the following.For the core field, the use of a second-order ARP induces a coupling between the field itself ( g c ) and its first time derivative ( ∂ t g c ).Therefore, z c = (g c , ∂ t g c ) T and the parameter of the process is given by: where τ c (ℓ) is also chosen to be scale dependent.Con- trary to first-order ARPs where the stationary state covariance matrices are given by ∞ z s = ∞ g s , for the core field it reads: as shown by Hulot and Le Mouël (1994).With the proposed setup, ∞ z s and τ s (ℓ) completely defines the dynamical behavior of the ARPs.The covariance matrices characterizing the stationary state of each source are assumed to derive from energy spectra E s (ℓ, a s ) expressed at given radii a s such as: where N m is the number of modeled spherical harmonics coefficients per degree ℓ , and R is given by R(ℓ) = ℓ + 1 and R(ℓ) = ℓ for internal and external sources, respec- tively.The shape of each energy spectrum is imposed.It can either be flat, such as E s (ℓ) = A 2 s or identical to the correlation kernels proposed by Holschneider et al. (2016) that we refer as of C-based type with E s = A 2 s (2ℓ + 1)R(ℓ) , where A s is the magnitude of the spectrum.For most sources, the dipole part is assumed to be independent from the rest of the spectrum such as s .Under these assumptions, the radii a s , the amplitudes A s and the dipole magnitudes D s form the free parameters of the stationary state covariance matrices ∞ z s .Characteristic timescales are parameterized by power laws such as τ s (ℓ) = M s ℓ −α s with given magni- tudes ( M s ) and slopes ( α s ) which are for some sources allowed to continuously vary from one range of spherical harmonics to the other.The ARP's parameters were estimated through a machine learning algorithm with a subsample of CHAMP and Swarm data as detailed in Baerenzung et al. (2020).The same values are used in this study.They are reported in Table 3. (5) Note that here, the energy spectrum of the lithospheric field is split into two ranges.In the first one, between ℓ = 1 and ℓ = 74 , the spectrum is of the C-based type and exhibits a characteristic radius of a l = 6287 km and a magnitude of A l = 0.16 nT.These are values obtained by Baerenzung et al. (2020).In the second range, between ℓ = 75 and ℓ = 1000 , the spectrum is flat with a l = 6367.9km and A l = 6.5 nT.In this case, the param- eters were estimated through a least square fit between ℓ = 75 and ℓ = 400 of the energy spectrum associated with the WDMAM model of Lesur et al. (2016).
The source associated with field-aligned currents, as well as the components at SH degree larger than ℓ = 1 of the fluctuating magnetospheric field, exhibit very small characteristic timescales of, respectively, τ fac (ℓ) = 1 min and τ fm (ℓ > 1) = 18 min.These timescales being smaller than the time step of the Kalman filter algorithm (here set to 30 min), the associated fields are assumed to temporally evolve as a white noise but are correlated in space and time during the analysis.Setting a priori a zero mean for both fields their covariance can be expressed as:

Filtering, smoothing, sampling
The prior statistical properties as well as the dynamics of the different magnetic sources being characterized, assimilation can be initiated.As a first step, a vector z containing the spherical harmonics coefficients of each ( 8) field is constructed.For the full model, z is composed of N M = 1002696 entries.The lithospheric field which is expanded up to ℓ = 1000 is filling more than 99.9% of the z vector.With such a model dimension, the size of the covariance matrix associated with z , namely z , should be of N M × N M ∼ 10 12 .Yet computations with such a matrix would be numerically impossible.This is why we approximate the predicted uncertainties of the small scale lithospheric field (for 101 ≤ ℓ ≤ 1000 ) by only keeping its variance information (the one associated with each of its spherical harmonics coefficients).Under such an assumption the dimension of z reduces to N M × N M ∼ 10 8 , a computationally conceivable size.This strong approximation, which induces a complete loss of the predicted spatial correlations of the lithospheric field beyond ℓ = 100 is evaluated in "Lithospheric field" section.
To forecast z , each parameter matrix F s of equation 3 are incorporated in a global matrix F .The same operation is performed for the stationary state covariance matrices ∞ which are assembled into the covariance matrix ∞ .Given F and ∞ , the covariance matrix associated with the Gaussian white noise of the full model forecast step reads ˜ = ∞ − F ∞ F T .Therefore, the evolution of the mean model and its covariance from time step k − 1 to step k is then given by: Table 3 Magnetic sources parameters as described in "Dynamical model" section and evaluated by Baerenzung et al. (2020) The prior spatial covariance matrices are derived from energy spectra expressed at some radii a s which are either flat with  2022) is applied.In their developments, they showed that at first order, the predicted intensity I k could be related to the predicted magnetic field B k through the relation , where Ĩk is the intensity derived from the mean magnetic field E[B k ] .Note that this projection is realized with the mean magnetic field prediction.Therefore, no iteration over the updated solutions is required.When d k correspond to secular varia- tion data, H k = 0 for each source except for the core field where H k projects its associated secular variation on the data.Once all data have been assimilated, the different modeled epochs are corrected through a smoothing algorithm (see Rauch et al. (1965)).Starting at the final step of the Kalman filter it performs iteratively backward in time through the following relations: where d corresponds to the full dataset.The smooth- ing algorithm only provides snapshots of the posterior model, therefore the resulting solution does not contain information about temporal correlations.Although the posterior covariance between the model at different epochs can be analytically derived, obvious storage limitations makes this option numerically inapplicable.This (11 is why we introduced a formulation to sample ensembles from the posterior model which are correlated both in space and time.Starting with an ensemble z e randomly drawn from the last state of the Kalman filter solution, the algorithm proceeds similarly to the smoothing algorithm, backward in time with: where ζ e is a random realization from the Gaussian dis- tribution characterized by a 0 mean and a covariance matrix given by: Note that to correct deviations, due to sampling errors, between the ensembles and the true posterior means, the ensembles were recentered at each epochs accordingly to the mean smoothing solutions.For this study, we used an ensemble of 1024 members.

Model construction
To construct the model, the time step of the Kalman filter algorithm was set to t = 30 min.Nevertheless, whenever the distance between two analysis windows exceeded this value, t was increased accordingly.With a dataset covering the 20 th century and the last 22 years, the direct approach would have been to start assimilating measurements in 1900.0 and to progress forward in time until today.However, we did not proceed this way.Instead, the Kalman filter simulation was initiated in 2000.5 to first assimilate ground-based observatories, CHAMP and then Swarm data until 2022.18, the last epoch at which measurements were currently available.The smoothing algorithm was then applied to update the model within this time window.In a third part, the smoothing solution in 2000.5 was used as a restart file to assimilate, backward in time, the measurements taken prior to this date.Finally, the smoothing algorithm was applied from 1900 to today with a slight modification beyond 2000.5 which is detailed in Appendix B. Two reasons motivated this choice of splitting the assimilation process.The first one is to possess a wellresolved large-scale lithospheric field before assimilating survey data.This, in order to be able to distinguish the gain of assimilating such observations on this part of the field.To this end, the lithospheric field was fully modeled up to SH degree ℓ = 150 during the CHAMP and Swarm eras.The full solution (mean and covariance) was then truncated at ℓ = 100 in 2000.5 to restart the Kalman filter between 2000.5 and 1900.0.Beyond ℓ = 100 only ( 17) the mean and variance were kept from the CHAMP and Swarm solution.This part was finally extrapolated with a zero mean and the prior variance of equation 7 between ℓ = 150 and ℓ = 1000 .The second reason for splitting the assimilation process was motivated by the fact that the older the measurements, the lower their accuracy and spatial coverage.Yet, before assimilating any measurement an outlier detection is performed.The latter process consists in checking that the measurements do not excessively deviate from their predicted values, in particular that each vector field or intensity measure lies within the 95.6% confidence interval of the model prediction.On top of this selection, the misfit of the sequentially assimilated tracks was evaluated.Whenever the misfit value exceeded the imposed threshold of 3, the corresponding track was dismissed.The algorithm to detect outliers performs better when the model accuracy is high, which occurs when the data quality and coverage are good.This is why starting with a very well constrained solution in 2000.5 and assimilating data backward in time enabled us to optimize the detection process.
Over the entire time span of the model, each source, except the one associated with field-aligned currents, is stored every 0.1 year, setting up the temporal resolution of the model to this time step.However, to better track the evolution of rapidly evolving sources, such as the close and fluctuating magnetospheric fields or the residual ionospheric/induced fields, the latter were stored every 3 hours during the CHAMP and Swarm eras and every 5 days between 1900 and 2000.5.

Main field and secular variation
The Kalman filter and smoothing algorithms provide a model in terms of mean solution and associated covariance matrix.Combining these two quantities gives a precise knowledge of locations where the solution is reliable and where it is not.As an illustration for the main field, i.e., the sum of the core field and the lithospheric expanded up to SH degree ℓ = 20 , Fig. 2 shows at different epochs the radial component of the mean field (isocontours) and its associated standard deviation (color maps).Locations where the maps are red correspond to locations where the mean solution is likely to deviate strongly from the true field.On the opposite, within blue and purple areas the model predicts that the true and the mean predicted field are close.These maps are complemented on their bottom right by a global measure of the predicted uncertainty.It corresponds to the r.m.s.standard deviation given in nT and expressed as: where σ is the standard deviation associated with the radial component of the field and is the Earth's surface.
Until the 1960s, uncertainty maps exhibited a strong dichotomy between the Northern and Southern hemispheres.Whereas in the North, the standard deviation associated with the radial component of the field does not globally excess 25 nT, it reaches and even exceeds 50 nT in the South.The difference of predicted uncertainties is particularly important between land and oceanic surfaces reflecting the lack of measurements taken over the latter, the location where the field is best resolved is Europe.This is a benefit of the high density of ground-based observatories operating at this place and during this time period.When looking at the r.m.s.standard deviation, the year 1920 slightly stands out with σ = 44 nT, whereas this value oscillates around σ ∼ 50 nT in 1910, 1930 and 1940.This phenomenon can be explained by the multiple land and marine surveys occurring at and around this epoch and which are offering a large data coverage of the globe (see Fig. 1).In 1960, the global resolution of the model is improved and the North-South dichotomy mostly disappears.Two reasons explain this gain of accuracy.The first one is the dense spatial coverage of survey data at this epoch (see Fig. 1).The second one is the time proximity of the POGO mission which started in 1965.One can also observe that observatories still play an important role to reduce the posterior variability as it is the case in and around Europe and Japan.In 1970, the jump in accuracy of the model is striking.At this period lying within the POGO era, the standard deviation associated with the Kalmag solution is strongly reduced.However, the model predicts a higher possible variability around the magnetic dip equator.This phenomenon is the transcription of the Backus effect, or more generally the "perpendicular error" effects within the model.Indeed, as first recognized by Backus (1970), to be then generalized by Lowes (1975), when constructing a geomagnetic field model with intensity measurements alone, larger errors will contaminate the model near the equator.This effect is surely affecting our mean solution, but covariance information enables us to quantify it.With MagSat observations, which cover less than a year ( 1979 − 1980 ), the model precision is equivalent to the one obtained with POGO data except around the dip equator where vector field measurements eliminate the "perpendicular error" effects induced by the assimilation of intensity data.The map in 1990 highlights the importance of low-orbiting satellites to recover the Earth's magnetic field.Lying between MagSat and Oersted missions, in the middle of almost 20 years without satellite measurements, the solution obtained at this time is strongly degraded.It presents levels of uncertainties equivalent to the 1960 ones except in Northern America and Russia where the coverage with ground-based observatories has since been increased.The situation is ameliorated with Oersted measurements and becomes even better with CHAMP and Swarm observations.With the high-quality instrumentation of CHAMP and Swarm satellites, the model is extremely precise and this is almost everywhere at the Earth's surface.It is however worth noting that the constellation of Swarm satellites permits to obtain a slightly more accurate solution than the unique CHAMP spacecraft.
When looking at the mean secular variation (SV) and its associated standard deviation as displayed at similar epochs in Fig. 3, one can observe that the dichotomy in accuracy between the North and the South is also present for this quantity.The dichotomy persists until the year 2000, but with a lower contrast after 1960.Ground-based observatory data are of particular importance to constrain the secular variation, as locations where their density is high always coincide with areas of low posterior variability.Globally, uncertainties are decreasing with time except between 1970 and 2000, where the r.m.s.standard deviation fluctuates due to the lack of persistent low-orbiting satellite missions.In addition, the distribution of uncertainties over the different spatial scales is not homogeneous.Instead, small scales typically exhibit a higher posterior variability relatively to their mean signal than large scales.This effect can be observed in Fig. 4 where time series between 1900 and 2022 of the 68.2% confidence interval associated with some selected SH coefficients are displayed in red.In this figure, it is clearly visible that the larger the degree of the coefficient (from left to right and top to bottom), the larger its posterior standard deviation relatively to its mean values.The COV-OBS.x2 model of Huder et al. (2020), exhibits a similar behavior as its predicted 68.2% confi- dence intervals (blue areas) show.Although the two models are mostly consistent with one another, small differences can nevertheless be distinguished, in particular in the predicted standard deviations.Until ∼ 1920 their level is lower for COV-OVS.x2,they become equivalent between COV-OVS.x2 and Kalmag until ∼ 1960 to be lower for Kalmag afterwards.
To precisely characterize the spatio-temporal resolution of the secular variation over the model time span, we computed the ratio C ġ (ℓ, k) between the Fourier power spectra of the mean secular variation and its associated standard deviation for 20 years time periods.This quantity, which was proposed by Gillet et al. (2015), can be expressed as: (20) where ĝc,ℓ,m (k) is the Fourier transform of the secular variation, and σ ĝc ,ℓ,m (k) is its associated standard devia- tion.To estimate the latter quantity, we used an ensemble of 1024 Fourier transform of secular variation time series.In Fig. 5, C ġ (ℓ, k) is displayed for 6 different time win- dows.The blue and red areas correspond to spatio-temporal scales which are, respectively, well resolved and not resolved.At early times, between 1900 and 1920, only some limited amount of temporal scales of the SV up to SH degree ℓ = 4 are resolved.The situation slightly improves between 1920 and 1960 where some signal up to SH ℓ = 6 can be accurately recovered, and this down to a few years for the largest spatial scales.The emergence of satellite missions and the increase of groundbased observatory and survey data helps improving the model resolution between 1960 and 2000.During this time interval some spherical harmonics coefficient up to degree ℓ = 5 are either partially or fully resolved down to time periods lower than a year.Reaching such a temporal resolution is impossible with the secular variation data derived from annual differences of observatory measurements.It can therefore only be achieved thanks to the high temporal coverage of satellite and survey data.In agreement with our previous results and with the study of Gillet (2019), the secular variation is best resolved during the CHAMP and Swarm eras, where spatial scale up ℓ = 15 can be partially resolved down to periods of approximately 5 years, and 2-year fluctuations can be very well captured up to ℓ = 10.

Lithospheric field
As previously mentioned, the lithospheric field model was built in multiple steps.During the CHAMP and Swarm eras, it was fully modeled up to SH degree ℓ = 150 .After applying the smoothing algorithm, the lithospheric in 2000.5 was divided in three parts.In the first one, between ℓ = 1 and ℓ = 100 , the full smoothing solution (mean and associated covariance matrix) was kept.In the second part, between ℓ = 101 and ℓ = 150 , only mean and variance information were considered.Finally, between ℓ = 151 and ℓ = 1000 a zero mean and the variance derived from equation 7 with parameters of Table 3 were a priori imposed.The Kalman filter algorithm was then launched backward in time with this prior lithospheric field between 2000.5 and 1900.
Keeping only variance information within the Kalman filter algorithm is a strong approximation.Before implementing it, this approximation was tested during the CHAMP and Swarm eras.For this evaluation phase, the lithospheric field was fully modeled up to ℓ = 30 and partially modeled (keeping only variance information) between ℓ = 31 and ℓ = 150 .The remaining part of the model was simulated normally and the dataset used is the one described in "Data" section.The resulting model is referred as the PR model.With this setup, comparisons with the solution obtained at full resolution (FR model) can be performed.In a first simulation, it was observed that the posterior variance associated with the approximated solution had a tendency to be underestimated.In particular, the transition between the degree variance (the sum of the variances at a given degree) at SH degree ℓ = 30 and ℓ = 31 exhibited a pronounced discontinuity.To partially correct this effect, variances beyond the transition were increased by a multiplication factor.The latter was imposed to vary linearly with the degree of the SH expansion, and forced a smooth transition as well as a level of variance at the last modeled degree corresponding to stationary state variance of equation 7.Because of the latter operation, the lithospheric field resolution was increased to ℓ = 200 , a degree at which the signal at satellite altitude becomes very low as shown by Olsen et al. (2017).
The results of this evaluation phase are displayed in Fig. 6.On the left panels, the mean downward component of the lithospheric field at the Earth's surface is shown for both the solution obtained at full resolution (top) and the one obtained at partial resolution (bottom).These two maps look very similar and most features which can be recovered by the FR model are present in the PR model.This aspect is confirmed by the map which exhibits the difference between the two mean solutions (top right).Only at the level of Antarctica, Eastern Europe and Western Russia, discrepancies become quite intense.These discrepancies coincide with relatively large-scale errors (up to ℓ = 70 ) as shown with crosses by the energy spectrum at the Earth's surface of the also highlights their proximity.The latter reaches a minimum of 0.915 at ℓ = 66 and stabilizes around the mean value of 0.979 beyond ℓ = 100 .The energy spectra asso- ciated with the standard deviations show that the model where only variance information was updated, had a tendency to underestimate the level of predicted uncertainties.Although the technique previously mentioned to rescale the variance was applied, it did not completely resolve this issue.Nevertheless, the fact that the smallscale lithospheric field was only marginally affected by the proposed modeling approximation comforted us to implement it for the complete model derivation.
The lithospheric field resulting from the assimilation of the entire dataset is first analyzed through energy spectra at the Earth's surface.In the left part of Fig. 7, the spectra of the mean, the standard deviation and the prior standard deviation of the lithospheric field are displayed with black lines.In this solution, energy populates the entire range of modeled scales.However, the mean field is predicted to be globally reliable only up to SH degree ℓ ∼ 450 , where the spectrum of the mean and the spec- trum of the standard deviation cross one another.In Although our solution is apparently closer at any degree to the LCS-1 model than to the WDMAM model, the examination of the degree correlation (right panel of Fig. 7) indicates that this aspect is only true up ℓ = 150 .Beyond this value, even if ρ ℓ is relatively low, the correla- tion between Kalmag and WDMAM (red line) is higher.Contrary to the degree correlation between LCS-1 and Kalmag which decays smoothly, the one associated with Kalmag and WDMAM presents two transitions.One of them is at SH degree ℓ = 100 , the spatial scale delimit- ing the satellite data solution ( ℓ ≤ 100 ) from the survey data solution ( ℓ > 100 ) of the WDMAM model.The other transition occurs at ℓ = 150 , the degree beyond which our model is only constrained by survey data.This second drop in ρ ℓ may be explained by the lower spatial resolution that our solution exhibits in certain areas.This phenomenon can be observed in Fig. 8 where the downward components of WDMAM (top left) and Kalmag The intense signals predicted by WDMAM in the Southern parts of the Pacific, the Atlantic and the Indian oceans, or on large portions of continental areas are mostly absent in our solution.It is however worth noting that WDMAM does not only derive from direct measurements of the geomagnetic field, but also from the combination of ocean floor age map, relative plate motions and geomagnetic polarity time scale (see Dyment et al. (2015)).Logically, the difference between the downward component of both models (top right of Fig. 8) is larger at these oceanic and land locations than anywhere else.On the opposite, discrepancies are reduced in most areas where the standard deviation associated with the large scale part of the field (up to ℓ = 100 ) is low (map on the bottom right).These uncertainty predictions which are tied to data coverage (see Fig. 1) therefore provide a good approximation of locations where the Kalmag model is likely to be well resolved.
The model being expressed in terms of posterior distributions, it can be used as a prior information to assimilate new data when some of them become available, and therefore be updated accordingly.measurements is of 18nT.Globally the predictions of the truncated model (right column) are closer to the data than the EMM predictions (third column).Of course Afghanistan is a particular location and no claim is made here that the Kalmag model would be globally more accurate than the EMM model since this is certainly not the case.However, this example shows that the method proposed in this study is well suited to construct regional high-resolution models of the lithospheric field and this even when data coverage is not optimal.

Magnetospheric and induced fields
With the proposed approach, magnetospheric and induced fields are jointly estimated with the rest of the model.A priori, the field generated by the currents flowing in the outer magnetosphere ( g rm ) is predicted to evolve slowly with time ( τ g rm = 10.3 years) in comparison to other external sources.A posteriori, such a behavior is confirmed as illustrated by the evolution of the annual mean dipole component of E[g rm ] projected in magnetic coordinates and shown in the left panel of Fig. 10 with circles.Note that prior to 1953, our model cannot correctly extract this field and the latter oscillates around 0 with a large posterior variance.However, g rm alone can- not explain decadal variations of external sources as they can be detected at the Earth's surface or at the altitude of low-orbiting satellites.The rapidly evolving magnetospheric components also exhibit long-term trends whenever the latter can be captured.This effect can be observed when comparing the annual mean dipole component of E[g rm ] to the one of E[g rm + g m + g fm ] shown with a continuous line in Fig. 10.During satellite eras, the latter is always found to be more intense than the former, meaning that the ring current can generate some persistent annual signal as already documented by Lühr and Maus (2010).With our current method, this signal can only be recovered when temporal data coverage is high enough due to the fact that E[g m ] and E[g fm ] exhibit very low memory timescales.A possible way to improve the AR processes characterizing these sources would be to consider some extra timescales accounting for the slow varying part of the field generated by the ring current.The cycle of approximately 10.5 years highlighted by Huder et al. (2020) with the COV-OBS.x2model (shown with dashed lines in Fig. 10) is also present in our solution.Although the mean solutions of both models slightly differ from one another, the COV-OBS.x2dipole always lies within the 68.7 confidence interval predicted by our model (gray areas in Fig. 10).
To evaluate the model over short periods of time and when all sources are predicted to be well separated, we now compare predictions of the azimuthal component of the model with ground based observatory measurements taken at four different locations, Hermanus, Niemegk, Canberra and Kakioka.Observatory data being only assimilated to constrain the core field secular variation, they can be considered as independent measurements for external and induced fields.In order to make visual comparisons possible and to remain within the conditions the model was built in, only hourly night-time measurements and predictions were kept to be then averaged over 10 days time periods.The results are reported in the right panel of Fig. 10 with red lines for observatory data, black lines for the full model predictions and blue lines for the predictions of the core field alone.Globally, monthly and annual variations of B θ are well captured by the model.Only during the time gap between the CHAMP and Swarm missions, when external sources are not updated anymore, predictions and observations differ strongly.One can also notice that the core field does not seem to be contaminated by external or induced fields, as its evolution does not reproduce the rapid variations observed in the data.The largest discrepancies between predictions and observations are in the magnitude of the signals.Intense excursions are not predicted by the model.The reason for this is that the model was trained on a dataset selected for very quiet magnetic conditions [see Baerenzung et al. (2020)].Therefore, the selection algorithm of the Kalman filter prevents the assimilation of data containing a too strong signal from external sources.A recalibration of the model for more general conditions would certainly solve this issue.
Finally, our model contains a source for induced/ residual ionospheric fields.The latter is a priori uncorrelated from magnetospheric fields.Yet rapid variations of external fields generate currents within the Earth's interior, which in return induce a secondary magnetic field (e.g., Schmucker 1985;Langel and Estes 1985b;Olsen et al. 2005;Finlay et al. 2020).The intensity and temporal evolution of the induced field depends on the conductivity of the crust, the mantle and the core.Under the assumption that conductivity only depends on depth, each spherical harmonics coefficient of the induced field will be linked the same coefficient of the external field through the relation: where ι is the induced field, ǫ the external fields, and Q is referred as the Q-response.In our model, ι = g ii and ǫ = g rm + g m + g fm , where g rm is projected in magnetic coordinates.
In the particular case discussed by Olsen et al. (2005), where the mantle is assumed to be insulating until a given depth d followed by a superconductor t) .Focusing on the dipole compo- nent of induced and external fields, and assuming a depth of d = 1200 km, leads to ι 1,0 (t) = Q1,0 ǫ 1,0 (t) with Q1,0 = 0.27 as estimated by Langel and Estes (1985b) with POGO data.In the left panel of Fig. 11, the evolution of, respectively, ǫ 1,0 and ι 1,0 (t)/ Q1,0 is displayed between 2019.45 and 2019.65 with, respectively, red and black lines.In order to concentrate on rapid variations only (we recall that external and induced field were stored every 3 hours during the CHAMP and Swarm eras), temporal scales larger than 15 days have been filtered out from both time series.Furthermore, we chose the [2019.45, 2019.65]time period because temporal coverage of Swarm data is optimal during this interval.The two time series in Fig. 11 follow one another quite closely and Q−1 1,0 seems appropriate to rescale the induced field.Over the current Swarm time span, induced and external fields exhibit a Pearson correlation ρ = Cov(ǫ, ι)/(σ ǫ σ ι ) , calculated here with the mean Kalmag solutions, of ρ = 0.79 .It is of ρ = 0.84 over the time interval of Fig. 11 and of ρ = 0.73 over the CHAMP era.This lower correlation value is probably caused by the uncertainty level of external and induced fields which are higher during the CHAMP mission than during the Swarm one.However, the particular 1-D conductivity model leading to Q1,0 is known to be imperfect.More complex conductivity profiles are required to better model induction processes within the Earth's interior.
We now investigate the Q-response predicted by our model when keeping the assumption that the conductivity within the Earth is only depth-dependent, but relaxing the constraint about its profile.For this evaluation, we operate in spectral space.Considering only dipole components of ι and ǫ and applying a Fourier transform to equation 22 the latter becomes: From this equation, the real and imaginary parts of Q(k) are, respectively, given by: To evaluate these two quantities we considered induced and external fields during the [2015.0,2021.0]time interval when the model reaches its peak accuracy.In the right panel of Fig. 11, Re{ Q(2π/k)} and Im{ Q(2π/k)} aver- aged at period T i = 2π/k i over [T i , 2T i ] are, respectively, displayed with red and black continuous lines.For comparisons, the real and imaginary parts of Q1,0 as well as the Q-response (referred as Q O ) estimated by Olsen et al. (2005) with a realistic conductivity model are, respectively, shown with dashed and dotted lines.The general behavior of the Q-response we recover is coherent with our prior knowledge about it.Indeed, for short periods of time the real part of Q is much more intense than its (23) ι1,0 (k) = Q1,0 (k)ǫ 1,0 (k) .
(24) imaginary part and its decay pattern is close to the one predicted by Olsen et al. (2005).However, in comparison to Re{ QO } , Re{ Q} is globally underestimated.This effect might be due to the fact that induced fields vary rapidly with time, and when no data is feeding the model, its mean value tends quickly toward 0 contrary to the remote and close magnetospheric fields which evolves slower.The behavior of the imaginary part of Q , which reflects the temporal lag of the induced field response, is on the contrary very similar to the one predicted by the direct model of Olsen et al. (2005).

Conclusion
In this study, we proposed a method to assimilate different types of geomagnetic data in order to construct a high spatio-temporal model of the Earth's magnetic field.The model being expressed in terms of posterior distribution, it reflects the quality and spatial coverage of the measurements it is derived from.At the beginning of the twentieth century, the main field and the secular variation are quite uncertain in the Southern hemisphere and more particularly in oceanic areas and in Antarctica.With the first data collected by loworbiting satellites, these two fields gain in precision and become very reliable during the CHAMP and Swarm eras.We demonstrated that the rapid dynamics of the core field could be captured by the model.However, the spatial resolution at which timescale fluctuations are recovered is not constant over time and strongly depends on the spatial scale considered.Typically, rapid variations can only be accurately modeled at large spatial scale.On the opposite, fluctuations of the secular variation at high spherical harmonics degree can only be resolved for long periods of time.The model reaches its peak accuracy both spatially and temporally during the CHAMP and Swarm eras.It is therefore mandatory that such satellite missions are perpetrated in the future to better understand the nonlinear and wave dynamics occurring within the Earth's outer core (see Aubert and Gillet (2021), Gillet et al. (2021)).
To be able to consider land, airborne and marine survey observations, which contain an intense contribution of the small-scale lithospheric field, the latter was modeled up to spherical harmonic degree ℓ = 1000 .However, this operation could not be per- formed directly, since the dimension of the associated covariance matrix would have forbidden any numerical computation.We therefore introduced, and conclusively evaluated, a statistical approximation where only mean and variance information were updated beyond ℓ = 100 .The resulting mean solution exhibits highly detailed structures on every areas where data coverage was dense enough.Furthermore, the part of the covariance which is still fully modeled (up to ℓ = 100 ) pro- vides a rough estimation of locations where the mean is likely to be well resolved.
An important aspect of the proposed approach is that whenever new observations become available, the model can be updated accordingly without restarting the entire assimilation process.The example presented with the dataset taken above Afghanistan demonstrates the flexibility of the method.
As for the core field, the accuracy of external and induced fields is not constant over the model time span.While signal of the remote magnetospheric field could be extracted from 1953 on, rapidly evolving sources such as the close and fluctuating magnetospheric fields or the induced field, could only be separated from the data when the latter exhibit a high temporal coverage.In general, optimal solution for such field was obtained during satellite eras and in particular during the CHAMP and Swarm ones.The global behavior of external fields is in agreement with previous studies of it (see Lühr and Maus (2010); Huder et al. (2020)).However, the training of the model under very quiet magnetic conditions forbids the reproduction of most intense external field variations.A recalibration of the model under more general conditions appears therefore as necessary.Although magnetospheric and induced fields were a priori assumed to be independent, their connection revealed itself a posteriori.Through the proposed approach we showed that external and induced fields could be jointly estimated from direct measurements of the geomagnetic field although the process characterizing their evolution remain quite simplistic.A refined parametrization of their dynamical behavior would certainly enhance the ability of the algorithm to extract such sources from the data.
The model will be frequently updated (at least once every 2 months), in particular with Swarm and observatory data.Furthermore, it can be accessed through different physical and statistical properties on a dedicated website at: https:// ionoc ovar.agnld.uni-potsd am.de/ Kalmag/.

Appendix A: Sequentialization
Describing the evolution of a given g quantity by continuous first and second-order auto-regressive processes can be preformed through the following relations: (A.1) where ω is a Gaussian white noise scaled by the factor σ .
Introducing z = g for first-order ARP and z = g, ∂tg T for second-order ARP, equations (A.1) and (A.2) can be written as: with A = 1/τ and ζ = σ ω for first-order ARP and: for second-order ARP.
The homogeneous solution of equation A.3 is given by: where F = exp(−A�t) is the parameter of the ARP as expressed in equations 4 and 5.The general solution of equation A.3 is simply z(t + �t) = Fz(t) + ξ , where the white noise ξ characterized by the distribution N (0, �) is chosen here to force the process to remain stationary.Under such a constraint, one can write that Therefore, calculating the spatial covariance of both sides of the solution z(t + �t) = Fz(t) + ξ and rearranging the result gives ˜ = ∞ − F ∞ F T .

Fig. 1
Fig. 1 Locations (dots) and epoch (color) of each land, airborne and marine survey measurement.Black triangles correspond to every ground-based observatories feeding the model with data

Fig. 2
Fig. 2 Standard deviation associated with the radial component at the Earth's surface of the sum of the core field and the lithospheric field expanded up to spherical harmonics degree ℓ = 20 .Each panel corresponds to a different epoch which is displayed on their bottom left.Isocontours show the mean B r = 0 solution.White triangles represent the locations of ground-based observatories available at the presented epochs.On the bottom right of each map is indicated the r.m.s standard deviation σ in nT

Fig. 3
Fig. 3 Same as Fig. 2 for the secular variation.Each quantity is expressed in nT/yr

Fig. 4 Fig. 5
Fig. 4 Time series between 1900 and 2022 of selected spherical harmonics coefficients (indicated on the top of each panel) of the secular variation at the Earth's surface.68.2% confidence interval of the Kalmag solution (red areas) and the COV-OBS.x2solution (blue areas) of Huder et al. (2020)

Fig. 6
Fig. 6 Lithospheric field at the Earth's surface expanded up to spherical harmonics degree ℓ = 150 .Left: mean downward component solution for the FR model estimated with full covariance information (top) and for the PR model estimated with variance only information from ℓ = 30 (bottom).Top right: difference between FR and PR models mean downward components.Bottom right: energy spectra of the means (solid lines), the standard deviations (dashed lines) and the difference (crosses) between the FR model (black lines) and PR model (blue lines) discontinuity in the spectrum of the mean at SH degree ℓ = 150 indicates that even up to ℓ ∼ 450 a non-negligible portion of the crustal signal remains unmodeled.Nevertheless, comparisons with the FR model previously discussed (blue lines and dots) demonstrate that the assimilation of survey data helps to better constrain the large-scale lithospheric field.Indeed, the mean signal of the final solution has gained in intensity, and its standard deviation has decreased.In the same figure, the spectra of the difference with two other lithospheric models, the WDMAM model by Lesur et al. (2016) (red dots), and the LCS-1 model by Olsen et al. (2017) (green dots), are also shown.

Fig. 7 Fig. 8 Fig. 9
Fig. 7 Left: energy spectra at the Earth's surface of the mean (continuous lines), the standard deviations (dashed lines), and the prior standard deviation (dash dotted line) for the Kalmag lithospheric model (black lines) and the FR model (blue lines).Dots represent the spectra of the difference between the Kalmag model and the WDMAM model by Lesur et al. (2016) (red), the LCS-1 model by Olsen et al. (2017) (green), and the FR model (blue).Right: degree correlation ρ ℓ between the Kalmag lithospheric field model and the WDMAM model (red lines), and the LCS-1 model (green lines) To illustrate this aspect, airborne intensity measurements taken above Afghanistan in 2006 and 2008 were put aside from the dataset serving the model derivation.They are now used to update the lithospheric field following the method detailed in Appendix C. The locations at which each measure was taken during these surveys are shown with colored dots (blue for 2006 red for 2008) in the bottom left panel of Fig. 9.The downward component of the mean prior lithospheric field, which comes from the smoothing solution taken up to ℓ = 1000 in 2006.0, is shown on the top left panel.Its resolution was increased to ℓ = 2000 before the Kalman filter simulation was launched.The result of the assimilation process is shown through the downward component of the mean posterior field in 2009.0 in the second panel of the top row of Fig. 9. On this map, it can be seen that structures which were completely invisible in the prior model appear in the posterior one.In particular, high-intensity anomalies could be detected along the Southern and Western border of Afghanistan.The field in the central part of the land is globally weaker.Such patterns are also predicted by the EMM 2017 model of Maus (2010) as shown on the third panel of the top row.They are nevertheless of lower magnitude, and less detailed due to the resolution of the model which is limited to ℓ = 790 .To make the compari- son with the EMM solution possible, the posterior mean was truncated at SH degree ℓ = 790 .The resulting down- ward component is shown in the top right of the figure.Now the two models are looking more alike.Nevertheless, discrepancies in predicted intensity still remain.In order to assess the degree of compatibility of the different models with the observations, the absolute value of the difference between a subset of the measurements and the intensities predicted by the sum of the core and the different lithospheric field solutions was computed.The results are shown on the bottom panel below each corresponding downward components.The model exhibiting the higher degree of freedom, displayed on the second column, is without surprise the model which can better explain the data.As shown on the bottom of the map, the r.m.s.difference between the model and the

Fig. 10
Fig. 10 Left: annual average of the mean dipole component in magnetic coordinates associated with the remote magnetospheric field g rm (line with circles), the sum of the remote g rm , close g m , and fluctuating g fm magnetospheric fields (black line), and the COV-OBS.x2model by Huder et al. (2020) (dashed line).The gray area represents the confidence interval predicted by E[g rm + g m + g fm ] ± σ .Right: azimuthal component of the geomagnetic field taken during night-time and averaged over 10 days time periods at the level of four ground observatories, Niemegk (top left), Kakioka (top right), Hermanus (bottom left), and Canberra (bottom right).Red lines correspond to observatory data, black lines to the full Kalmag model predictions and blue lines to the Kalmag core field predictions

Fig. 11
Fig. 11 Left: mean dipole component in magnetic coordinates and evaluated at a radius r = 6371.2km of the sum of all magnetospheric sources (red line) and of the induced field rescaled by the inverse Q-response Q−1 1,0 discussed in the manuscript (black line).Components associated with time periods larger than 15 days have been filtered out.Right: real (red) and imaginary (black) parts of the Q-response.Full lines: Kalmag solution estimated during the [2015, 2021] time interval and averaged over [T i , 2T i ] period intervals.Lines with circles: Q-response estimated by Olsen et al. (2005) with a 1-D model of the conductivity within the Earth's interior.Dashes: Q-response associated with a simplified 1-D conductivity model of the Earth's interior (see text)

Table 1
Dataset used to derive the model.Missions (first column) and their time span (second column)

surveys data Source Coordinate ℓ max m max
Baerenzung et al. (2020)nd external sources The characteristic timescales are parameterized by τ s (ℓ) = M s ℓ −αs After the forecast, whenever measurements are available, the model is updated accordingly.This operation is performed through a Bayesian inversion which reads:where R k is the covariance matrix associated with meas- urement errors, K k is the Kalman gain matrix and H k is the operator projecting the model to the observations d k at iteration k.R k is chosen to be diagonal with con- stant standard deviations of 0.1 nT for intensity data [seeQuesnel et al. (2009)] and vector field measurements, and of 4.85 nT/yr for each component of secular variation data as we estimated it with a similar algorithm used to calibrate Kalmag (seeBaerenzung et al. (2020)).When d Mauerberger et al. (2020)y measurements, the linearization approach proposed byMauerberger et al. (2020);Schanner et al. (