A model of Earth’s magnetic field derived from 2 years of Swarm satellite constellation data

More than 2 years of magnetic field data taken by the three-satellite constellation mission Swarm are used to derive a model of Earth’s magnetic field and its time variation. This model is called SIFMplus. In addition to the magnetic field observations provided by each of the three Swarm satellites, explicit advantage is taken of the constellation aspect of Swarm by including East–West magnetic intensity and vector field gradient information from the lower satellite pair. Along-track differences of the magnetic intensity as well as of the vector components provide further information concerning the North–South gradient. The SIFMplus model provides a description of the static lithospheric field that is very similar to models determined from CHAMP data, up to at least spherical harmonic degree n=75\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=75$$\end{document}. Also the core field part of SIFMplus, with a quadratic time dependence for n≤6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n \le 6$$\end{document} and a linear time dependence for n=7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=7$$\end{document}–15, demonstrates the possibility to determine high-quality field models from only 2 years of Swarm data, thanks to the unique constellation aspect of Swarm. To account for the magnetic signature caused by ionospheric electric currents at polar latitudes we co-estimate, together with the model of the core, lithospheric and large-scale magnetospheric fields, a magnetic potential that depends on quasi-dipole latitude and magnetic local time.


Introduction
Swarm, a satellite constellation mission comprising three identical spacecraft, was launched on November 22, 2013. Two of the Swarm satellites, Swarm Alpha and Swarm Charlie, are flying almost side-by-side in nearpolar orbits of inclination 87.4 • at an altitude of about 465 km (in November 2015) above a mean radius of a = 6371.2 km. The East-West separation of their orbits is 1.4 • in longitude, corresponding to 155 km at the equator. The third satellite, Swarm Bravo, flies at a slightly higher (about 520 km altitude in November 2015) orbit of inclination 88 • .
Each of the three satellites carries an Absolute Scalar Magnetometer (ASM) measuring Earth's magnetic field intensity, a Vector Fluxgate Magnetometer (VFM) measuring the magnetic vector components and a three-head Star TRacker (STR) mounted close to the VFM to obtain the attitude needed to transform the vector measurements of the VFM magnetometer to a known coordinate frame. Time and position are obtained by on-board GPS. All Swarm data are available at http://earth.esa.int/ swarm. Quite a number of models of the recent geomagnetic field have been derived during the last few years. One class of model is based on the combined analysis of data from several satellite missions (in particular Ørsted, CHAMP and Swarm), sometimes also including ground observatory data in order to obtain an improved description of field time variations. Examples of such models are: the CHAOS series (e.g., Finlay et al. 2015Finlay et al. , 2016Olsen et al. 2006Olsen et al. , 2014, the GRIMM model series (e.g., Lesur et al. 2008Lesur et al. , 2010, the POMME models (e.g., Maus et al. 2005Maus et al. , 2006 and the Comprehensive Model (CM) series (e.g., Sabaka et al. 2002Sabaka et al. , 2004Sabaka et al. , 2015. Other recent models are based on Swarm satellite magnetic data alone. Prominent examples are some of the Open Access *Correspondence: nio@space.dtu.dk Division of Geomagnetism, DTU Space, Technical University of Denmark, Diplomvej 371, 2800 Kongens Lyngby, Denmark candidate models for IGRF 2015 that are collected in a special issue of Earth, Planets and Space (see Thebault et al. 2015, for an overview), and dedicated lithospheric field models that were derived from Swarm observations after subtracting model values of the core and large-scale magnetospheric field (e.g., Kotsiaros 2016;Thebault et al. 2016).
The present paper describes a model of the Earth's magnetic field that has been derived only from the first 28 months of Swarm data. It is an extension of the Swarm Initial Field Model (SIFM) of Olsen et al. (2015) that includes more recent data as well as vector gradient estimates. Shortly after launch a difference in the measurements taken by the ASM and VFM was observed; this so-called VFM-ASM disturbance field issue had not been solved when SIFM was derived. This effect resulted in degraded vector gradient data at that time, and therefore, only scalar intensity gradient estimates (no vector gradients) were used for SIFM. However, a procedure for correcting the magnetic field data for the "VFM-ASM disturbance field" has been found in the meantime (c.f. Lesur et al. 2015), as discussed in detail by Tøffner-Clausen et al. (2016), which allows us to now also include vector gradient data in our new model.
The goal of the investigations presented in the present paper is threefold. Firstly, we study the impact of including more data of the same kind as used for SIFM (i.e., scalar, vector and scalar gradient data) on the model results; we denote this extended SIFM model as SIFM x . Secondly, we investigate model improvement by including also vector gradient data. And thirdly we look at some systematic behavior of model residuals in the polar regions, in particular their dependence on magnetic local time (MLT) and the Interplanetary Magnetic Field (IMF), and assess their possible impacts on core and lithospheric field models. Our final model, denoted as SIFMplus, coestimates a magnetic potential that depends on quasidipole latitude and MLT. An additional model presented for comparisons that includes vector gradient data but no specific treatment of ionospheric currents in the polar region is called SIFMplus noMLT .

Data and model parameterization
We used 28 months (November 26, 2013-March 30, 2016 of magnetic data from the three Swarm satellites. Data were selected using the same criteria as for the SIFM model : In particular, we select data (vector and scalar) from dark regions only (sun at least 10 • below the horizon) for which the strength of the magnetospheric ring current, as measured by the RC index (Olsen et al. 2014), varied by at most 2 nT/h. At quasidipole (QD) latitudes (Richmond 1995) equatorward of ±55 • we require that the geomagnetic activity index Kp ≤ 2 0 , while for regions poleward of ±55 • QD latitude the weighted average over the preceding 1 h of the merging electric field at the magnetopause (e.g., Kan and Lee 1979) has to be below 0.8 mV/m. The vector components of the magnetic field were taken for non-polar latitudes (equatorward of ±55 • QD latitude), while only scalar data were used for higher latitudes.
In contrast to the selection of magnetic vector and scalar field data (Kp ≤ 2 0 , |dRC/dt| < 2 nT/h, a condition that in 2014-2015 was fulfilled for 39 % of the time) we allow for higher geomagnetic activity when selecting gradient data (Kp ≤ 3 0 , |dRC/dt| < 3 nT/h, which is fulfilled for 60 % of the time). We also use scalar and vector gradient data from the dayside but excluded low-latitude (QD latitudes < ±10 • ) dayside data to avoid contamination by the Equatorial Electrojet.
The East-West gradient is approximated by the difference δB EW = ±[B A (t 1 , r 1 , θ 1 , φ 1 ) − B C (t 2 , r 2 , θ 2 , φ 2 )] of the magnetic observations measured by Swarm Alpha (also referred to here as SW-A) and Charlie (SW-C), where B may be either the scalar intensity F or one of the three magnetic vector components. Here t i , r i , θ i , φ i , i = 1-2 are time, radius, geographic co-latitude and longitude of the two observations. The sign of the difference was chosen such that δφ = φ 1 − φ 2 > 0. For each observation B A (from SW-A) fulfilling the above selection criteria we selected the corresponding value B C (from SW-C) that was closest in co-latitude θ, with the additional requirement that the time difference |δt| = |t 1 − t 2 | between the two measurements should not exceed 50 s.
The North-South gradient is approximated by the difference δB NS = ±[B k (t k , r k , θ k , φ k ) − B k (t k + 15 s, r k + δr, θ k + δθ, φ k + δφ)] of subsequent data measured by the same satellite (k = A, B or C) 15 s later, corresponding to an along-track distance of ≈115 km (≈1 • in latitude). The sign of the difference was chosen positive if δθ > 0, otherwise negative.
Not only the data selection but also the basic model parameterization follows closely that of SIFM. The model parameters consist of spherical harmonic expansion coefficients for the magnetic scalar potential V and sets of Euler angles describing the rotation of the satellite vector measurements from the magnetometer frame to the star tracker frame.
The magnetic field vector B = −∇V is derived from the magnetic scalar potential V = V int + V ext consisting of a part, V int , describing internal (core and lithospheric) sources, and a part, V ext , describing external (mainly magnetospheric) sources and their Earth-induced counterparts. Both parts are expanded in terms of spherical harmonics.
For the internal part this yields where (r, θ , φ) are geographic coordinates, P m n are the associated Schmidt semi-normalized Legendre functions, g m n , h m n are the Gauss coefficients describing internal sources, and N int = 80 is the maximum degree and order of the internal expansion (for SIFM N int = 70 was chosen). Coefficients up to degree n = 6 include a quadratic dependence on time (i.e., a secular acceleration), while coefficients of degree n = 7-15 vary linearly in time (linear secular variation). This yields 80 × 82 + 15 × 17 + 6 × 8 = 6863 coefficients describing the internal part of Earth's magnetic field.
The parameterization of external magnetic field contributions is also similar to that of our previous models, with an expansion of near magnetospheric sources in the Solar Magnetic (SM) coordinate system (up to n = 2, with special treatment of the n = 1 terms) and of remote magnetospheric sources in Geocentric Solar Magnetospheric (GSM) coordinates (also up to n = 2, but restricted to order m = 0). We solve for an RC-baseline correction (described by SM dipole coefficients that explicitly vary in time) in bins of 5 days (for m = 0), resp. 30 days (for m = 1), which in total results in 238 parameters describing the external field part of the model. See Sect. 3 of Olsen et al. (2014) for details on this parameterization of magnetospheric field contributions.
Experience with previous modeling efforts, like SIFM, reveals that model residuals (δF = F obs − F mod ) are typically considerably larger in the polar regions compared to non-polar latitudes. This is likely due to ionospheric cross-cap currents and associated auroral electrojets, below the satellite altitude. In an attempt to describe the magnetic signatures of these currents we co-estimate a spherical harmonic expansion of an additional potential part that depends on QD latitude θ QD and magnetic local time τ: Note that there are no m = 0 (i.e., zonal) terms, which means no variation independent of MLT, since that part is handled by the SH expansion of the core and crustal field, cf. Eq. (1). This results in 420 additional coefficients g m,MLT n , h m,MLT n , which are determined by scalar field and scalar gradient data at all latitudes, and by vector field data at non-polar latitudes. To stabilize the solution at non-polar latitudes we add 10,000 synthetic zero-value data points randomly distributed in MLT and non-polar (±55 • ) QD latitudes. Adding the 420 coefficients of Eq. 2 to the 238 magnetospheric terms results in 658 parameters describing external field contributions. Finally, we co-estimate the Euler angles describing the rotation between the vector magnetometer frame and the star tracker frame in bins of 10 days (i.e., 3 × 85 sets of angles for each of the three satellites Alpha, Bravo and Charlie) which results in an additional 765 model parameters.
The 8286 model parameters are estimated from almost 15 × 10 6 observations (373,985 scalar data, 3 × 1,272,456 = 3,817,368 vector data, 4,360,011 estimates of scalar gradients and 3 × 2,048,239 = 6,144,717 estimates of vector gradients) by means of an Iteratively Reweighted Least-Squares approach using Huber weights, using SIFM as starting model. The gradient data were handled by taking the difference of the design matrices corresponding to the two positions t i , Gradient dayside data do not contribute to the core field part of the model (i.e., internal Gauss coefficients up to n = 15), whereas the remaining parts of the model are constrained by all data. No model regularization has been applied. Figure 1 presents Mauersberger-Lowes spectra and degree correlation for various models. For reference we also show results for the MF7 (Maus 2010) and CHAOS-6 (Finlay et al. 2016) field models, which were derived mainly from CHAMP satellite data. The lithospheric field spectra, presented in Fig. 1a, show very good agreement between SIFMplus and MF7 as well as CHAOS-6 at least up to degree n = 70 or 75. The degree correlation ρ n , shown in the bottom left part of the figure, is above 0.9 for n < 73. This is a considerable improvement compared to SIFM (red curve), for which ρ n drops below 0.9 already at degrees n ≥ 60. Inclusion of 9 months of additional Swarm scalar, vector and scalar gradient data leads to a major increase in correlation, as is obvious when comparing the red (SIFM) and yellow (SIFM x ) curves. Adding vector gradient data and modeling of the MLT dependence of polar cap currents results in a further improvement (model SIFMplus, green curves), although spectra and coherence of a model without such a MLT-dependent polar field (SIFMplus noMLT , magenta curves) are very similar to those of SIFMplus.

Results and discussion
Models derived from magnetic intensity data alone suffer from the Backus effect (e.g., Backus 1970;Stern and Bredekamp 1975). The high-degree lithospheric field coefficients of the models SIFM and SIFM x are primarily constrained by scalar gradient data, which do not provide information concerning small-scale features at low latitudes. Models that were derived without vector gradient data, like SIFM and SIFM x , are therefore erroneous near the magnetic equator. This finding from SIFM (cf. Fig. 3 of Olsen et al. 2015) is confirmed by the top right part of Fig. 2 where the difference in B r between SIFM x and MF7 reveals a clear pattern following the magnetic equator. Adding vector gradient data reduces the Backus effect, and indeed, the difference of the lithospheric radial magnetic field between SIFMplus and MF7, shown in the top left part of the figure, indicates better agreement between these two models.
Although there is hardly any difference between the lithospheric models with (SIFMplus) and without (SIFM noMLT ) MLT-dependent polar ionospheric field when looking at the geomagnetic spectra (the spectrum of the model difference is below 0.3 nT 2 for all degrees n), maps of B r at Earth's surface differ by up to 18 nT, as can bee seen from the bottom left part of Fig. 2. This behavior is due to the fact that the model differences are, as expected, concentrated in the polar regions, whereas the spectrum measures the global average. Accounting for a MLT-dependent ionospheric field in Swarm satellite data may therefore indeed improve lithospheric field models.
The first time derivative (SV) between the two models differs in polar regions at Earth's surface by less than 4 nT/year, which is much weaker than the SV signal (up to 200 nT/year). As expected, there are almost no model differences at non-polar latitudes.
Despite these model improvements, the rather high altitude (about 450 km as of April 2016) of the satellites means that Swarm is not yet in an optimal configuration for determination of small-scale lithospheric structures, compared to what is possible with data collected by the CHAMP satellite. CHAMP was flying at altitudes below 330 km during the last 2 years of its mission, which yields a crustal field power at degree n = 75 that is about 25 times stronger compared to Swarm at its present altitude. Lithospheric field models derived from low-altitude CHAMP data are therefore probably still superior to models determined from Swarm, but the gradient concept of Swarm discussed in this paper is nonetheless promising for the analysis of future Swarm data at lower altitude. The bottom right part of Fig. 2 shows B r at Earth's surface from the final SIFMplus model. Inclusion of gradient data not only improves the highdegree lithospheric field but also is beneficial for the determination of secular variation, as was shown by Olsen et al. (2015). However, inclusion of vector gradient data does not seem to further improve the determination of secular variation. This is evident in the right panel of Fig. 1 which shows degree power (top) and degree correlation (bottom) for the first time derivative (secular variation) of the core field. Adding 13 months more data improves the model (compare SIFM, red curve, and SIFM x , orange), but no further improvement was achieved when adding vector gradient data (SIFMplus noMLT , magenta). This is likely because the secular variation part of the model does not suffer from the Backus effect as much as the lithospheric part, since the vector field data used are sufficient to eliminate any Backus effect signature at the degrees (n ≤ 15) of the secular variation model. However, solving for a MLT-dependent polar ionospheric field (model SIFMplus, green curves) yields a slightly better agreement with CHAOS-6. Although that model does not account for a MLT-dependent polar ionospheric currents, it is based on data when IMF B z > 0, which likely reduces the disturbing effect of these currents. Table 1 lists the number of data points, together with Huber-weighted means and root-mean-squared (RMS) misfit values between the observations and the predictions of the models SIFMplus (top), SIFMplus noMLT (middle) and SIFM x (bottom). Most remarkable is the reduction of the polar scalar misfit F polar by 10 % of SIFMplus compared to SIFMplus noMLT . There is also a slight misfit reduction of 2-3 % of the polar scalar gradient data and of the non-polar vector data, resulting in a RMS misfit of the non-polar radial component of less than 2 nT.  Also remarkable is the achieved RMS misfit of the gradient data which (for non-polar dark conditions) is below 200 pT for the N-S scalar gradient, about 270 pT for the radial and North-South components δB r,NS , δB θ ,NS of the vector gradient, and slightly higher (about 350 pT) for the East-West component δB φ,NS . Note that these N-S gradients have been obtained using data from the same instrument.

SW-A SW-B SW-C SW-A -SW-C
The RMS misfit of non-polar dark E-W gradient data is slightly higher compared to the N-S gradient data: 370 pT for the scalar gradient and between 480 and 580 pT for the vector gradients, with largest value again for the East-West component δB φ,EW . The dayside RMS misfits are slightly higher due to enhanced ionospheric contributions, both for N-S and for E-W gradient data.
Since November 4, 2014, no ASM scalar data are available for Swarm Charlie due to a fatal instrument failure. As a consequence, the VFM magnetometer of Charlie is calibrated using "mapped" scalar data from the Alpha satellite, as described in Tøffner-Clausen et al. (2016). This leads to marginally degraded scalar measurements on Charlie after November 4, 2014 (since F = |B VFM | instead of F ASM have to be used), resulting in an increase of the RMS misfit of δF EW, non-polar, dark from 290 pT for the first 7 months of the constellation (April-November 2014) to 340 pT between November 2014 and January 2016. This is confirmed by the time series of the E-W gradient residuals (observed gradient estimates minus their model predictions), presented in Fig. 3 for the scalar gradient and for the three vector gradient components, considering non-polar dark data. The solid blue line shows the monthly mean value of the residuals, while the shaded blue area indicates their ±1σ estimates. The time average of these estimates of standard deviation σ corresponds to the non-polar dark RMS values listed in the last column of Table 1. As expected, the scatter of the scalar gradient residuals of δF EW, non-polar, dark is smaller than that of the vector gradient. The running mean values of �δF (top panel) and of �δB θ (third panel) increase slightly since end of 2015. We attribute this to a not yet optimal determination of the time-dependent scale values of the VFM instrument on Swarm Charlie. More data will certainly help to reduce this end effect. The enhanced scatter of �δF during periods when the satellites are in a dawn dusk orbit (the times of which are shown by the red and blue dashed vertical lines) indicates enhanced ionospheric field signature near the terminator (separating the sunlit side of the Earth from the night time areas). Figure 4 shows QD-latitude/MLT maps of the scalar field signature at 400 km altitude, synthesized from Eq. 2, for the Northern (left) and Southern (right) polar regions, respectively.

Polar cap scalar residuals
These patterns reveal many features corresponding to the well-known current systems associated with plasma convection in the polar cap ionosphere. Most investigations of this current system have been made for geomagnetic active conditions, but our analysis confirms that magnetic fields associated with these current systems are also present during the geomagnetic quiet times (|dRC/dt| < 2 nT/h, E m < 0.8 mV/m) and for "dark" conditions (sun more than 10 • below horizon) that is typically used for geomagnetic field modeling. Although we have not explicitly selected polar data according to the Kp index, the above-mentioned selection criteria result in Number N of data points, (Huber-weighted) mean, and RMS misfit (in nT) of scalar (F), vector (B r , B θ , B φ ), N-S gradient (δF NS , δB r,NS , δB θ,NS , δB φ,NS ) and E-W gradient (δF EW , δB r,EW , δB θ,EW , δB φ,EW ) data, at polar (>±55 • ) and non-polar (< ±55 • ) QD latitudes and for dark (sun at least 10 • below horizon) and sunlit conditions data periods for which Kp < 3 0 (Kp < 2 0 ) is fulfilled for 99 % (88 %) of the time, which confirms that also according to the Kp index the data set is representative of geomagnetic quiet conditions. Previous modeling efforts that attempt to describe polar ionospheric currents include the Comprehensive Model series (Sabaka et al. 2002(Sabaka et al. , 2004(Sabaka et al. , 2015 and the GRIMM model (Lesur et al. 2008). However, neither of these models show the typical convection cell pattern (Fig. 4). For GRIMM a possible reason is that the modeling was done using dipole coordinates which seems to be less optimal for describing the polar current systems that are organized with respect to the magnetic pole (i.e., QD coordinates) rather than the geomagnetic (dipole) pole.
As discussed in the previous section, modeling the polar cap ionospheric currents by co-estimating the potential V MLT of Eq. 2 reduces the scalar polar misfit by 10 % and up to 30 % in the auroral oval. But despite this reduction there is still considerably more scatter at polar latitudes compared to non-polar regions. Part of this scatter is due to the dependence of the polar ionospheric currents on the IMF, which is not accounted for in the average maps shown in Fig. 4.
We therefore divided the SIFMplus model residuals according to the direction (B y , B z ) of the IMF and determined mean residual maps for each of the four possible cases (B y , B z positive or negative). Note that these maps do not show the total QD latitude/MLT dependence but only the variability with IMF on top of the mean maps presented in Fig. 4. These resulting difference maps, shown in Fig. 5, show larger residuals for B z < 0, which is expected since Fig. 3 Residuals of differences ("gradient") between Swarm Alpha and Charlie at non-polar latitudes and for non-sunlit conditions. Blue curves show monthly means; shaded blue area indicates ±1σ. The dashed red and blue vertical lines indicate the satellites are at dusk or dawn local time negative B z often results in higher geomagnetic activity. It should, however, be noticed that our selection of geomagnetic quiet times results in fewer data points for negative B z (about 76,000 out of 313,000 data points) compared to B z > 0, and thus the residual maps for negative B z are obtained from fewer data. The number of data points for B y > 0 (about 166,000) is comparable to those for negative B y (about 147,000).
By using the potential V MLT of Eq. 2, which depends on QD latitude and MLT, we assume that the QD coordinates describe an orthogonal frame. This is, however, not strictly true. Considering coordinates as well as magnetic Northern Hemisphere Southern Hemisphere vector components in a consistent way in the QD-frame is discussed by Richmond (1995) and applied to magnetic satellite data from the CHAMP and Swarm missions by Laundal et al. (2016) in an attempt to consistently determine ionospheric and field-aligned currents in the polar regions. Encouraged by the achieved misfit reduction of 10 % that we achieve by using the approximation of Eq. 2, we plan to implement the correct approach used by Laundal et al. (2016) in a co-estimation with the core and lithospheric field.

Conclusions
We derived a new model of Earth's magnetic field from more than 2 years of Swarm satellite constellation data. The model is an extension of the Swarm Initial Field Model (SIFM) of Olsen et al. (2015) by adding more recent Swarm measurements, by including vector gradient data at non-polar latitudes and by co-estimating a polar ionospheric field that depends on magnetic local time (MLT).
The SIFMplus model provides a description of the static lithospheric field that is very similar to that of the MF7 (Maus 2010) and CHAOS-6 ) models which were determined from CHAMP data, up to at least spherical harmonic degree n = 75. Also the core field part of SIFMplus, with its quadratic time dependence for n ≤ 6 and a linear time dependence for n = 7 − 15, demonstrates the possibility to determine high-quality field models from only 2 years of Swarm data, thanks to the unique constellation aspect of Swarm.
Co-estimation of the magnetic field caused by polar cap ionospheric currents reduces the polar scalar RMS model misfit by 10 % (up to 30 % in the auroral oval) compared to a model that does not account for such currents. However, despite these improvements there is a considerably larger scatter of the model residuals in polar regions. This is partly due to the dependence of the current systems on the direction of the IMF. We plan to account for this dependency in future.
The SIFMplus model, and software to evaluate it, is available from www.spacecenter.dk/files/magnetic-models/ SIFMplus/.