Parent magnetic field models for the IGRF-12GFZ-candidates

Abstract

We propose candidate models for IGRF-12. These models were derived from parent models built from 10 months of Swarm satellite data and 1.5 years of magnetic observatory data. Using the same parameterisation, a magnetic field model was built from a slightly extended satellite data set. As a result of discrepancies between magnetic field intensity measured by the absolute scalar instrument and that calculated from the vector instrument, we re-calibrated the satellite data. For the calibration, we assumed that the discrepancies resulted from a small perturbing magnetic field carried by the satellite, with a strength and orientation dependent on the Sun’s position relative to the satellite. Scalar and vector data were reconciled using only a limited number of calibration parameters. The data selection process, followed by the joint modelling of the magnetic field and Euler angles, leads to accurate models of the main field and its secular variation around 2014.0. The obtained secular variation model is compared with models based on CHAMP satellite data. The comparison suggests that pulses of magnetic field acceleration that were observed on short time scales average-out over a decade.

Background

The International Geomagnetic Reference Field (IGRF) is a geomagnetic field model used for numerous scientific and industrial applications. It is updated every 5 years (Macmillan et al. 2003; Maus et al. 2005; Finlay et al. 2010). The model is determined by a group of geomagnetic field modellers associated with the International Association of Geomagnetism and Aeronomy (IAGA) and built by comparison of different model candidates generated by scientists affiliated to different institutions around the world. A single institution can propose only one series of models. In this short paper, we present the model candidates proposed by scientists of the GFZ German Research Centre for Geosciences, in collaboration with scientists from other institutions.

Three model candidates were provided for the IGRF-12: two main field model snapshots, one for 2010 and one for 2015, and a predictive linear secular variation (SV) covering years 2015 to 2020. The main field model for 2015 and the SV model were derived from a parent model built from a combination of Swarm satellites and observatory data. This parent model includes a complex time-dependent parameterisation of the core field, a static representation of the lithospheric field, the external fields and their induced counterparts. Weaker signals, such as the field generated by the tidal motion of the oceans, are not modelled. Hereinafter, we do not present the derivation of this parent model but a very similar model built following the same approach and using a slightly longer time series of satellite data that have been re-calibrated. In the same way, the parent model of the main field snapshot for 2010 is not described. It has been derived from observatory and CHAMP satellite data but otherwise follows the same model parameterisation as the parent model derived from Swarm data.

The Swarm constellation of satellites was launched in November 2013, but the satellites reached their survey orbits only by mid-April 2014. At this early stage of the mission, the data are not yet fully calibrated and a specific data set has been provided by the European Space Agency (ESA) to be used for modelling purposes. However, some difficulties had to be handled in order to use these data. In particular, each satellite is carrying two instruments for magnetic measurements, and a correction had to be applied to explain the observed differences between the calculated magnetic field strength from the vector fluxgate measurements (VFM) and the magnetic total intensity obtained form the absolute scalar measurements (ASM). The first part of the second section of this paper is dedicated to the description of this correction. We made the choice to present here the latest version of the different correction processes we studied and the field model associated with it. It will be shown in the last section of this paper that our IGRF candidates are very similar to the model obtained with this corrected data set. Of course, the parent model of our IGRF candidates also used corrected data, but with a slightly less robust correction process than the one presented below. Outside this re-calibration, the processing path used to obtain accurate models of the magnetic field is similar to that of previous models of the magnetic field model series GRIMM (Lesur et al. 2008, 2010; Mandea et al. 2012; Lesur et al. 2015) We note however that, unlike during the CHAMP epoch, only a very short time span of satellite data was available at the end of September 2014 when the IGRF candidates were submitted. Therefore, observatory data had to be used to obtain robust models. Furthermore, specific regularisation processes had to be introduced in order to obtain models of acceptable quality.

The next section described the methods used to obtain the field models. In the first step of the modelling effort, data are selected and processed. This is described in detail in the second sub-section. The model parameterisation is explained in the third sub-section and the data inversion process in the fourth. The results are presented in the third section and discussed.

Methods

Corrections of the ASM-VFM differences

On each Swarm satellite, two types of instruments are providing magnetic data. The scalar instrument provides ASM data. These are absolute data, possibly corrected for stray fields (e.g. the magnetic field generated by the torquers), but otherwise not sensitive to temperature changes or ageing. In contrast, the vector data provided by the fluxgates - i.e. VFM data - need full calibration. After correction of stray fields and possible temperature drift, nine parameters per satellite have to be estimated for this calibration:

• Three scaling values, one for each of the three sensors, in the directions (defined below) E ₁, E ₂ and E ₃, respectively. These three scaling parameters are s ₁, s ₂ and s ₃.

• Three offset values, one for each of the three sensors, in the directions E ₁, E ₂ and E ₃, respectively. These three offset parameters are o ₁, o ₂ and o ₃.

• Three angles, called the non-orthogonality angles, that are calculated to insure that the three magnetic field components are in orthogonal directions. These three angles are a ₁₂, a ₂₃ and a ₃₁. They describe deviations from 90° of the angles between the E ₁ E ₂, E ₂ E ₃ and E ₃ E ₁ sensor directions, respectively – i.e. if a ₁₂, a ₂₃ and a ₃₁ are zero, then the E ₁, E ₂ and E ₃ sensor directions are already along orthogonal directions. In the process of estimating these three angles, the E ₁ direction is assumed fixed and is not modified. The re-orientation of the obtained orthogonal set of directions relative to the Earth-fixed, Earth-centred coordinate system is performed at a later stage of the processing, sometimes simultaneously with the field modelling process.

The VFM sensor E ₁, E ₂, and E ₃ directions correspond roughly to the direction perpendicular to the satellite boom oriented down, the direction perpendicular to the boom oriented right relative to the satellite flying direction, and the direction along the boom oriented toward the scalar magnetometer, respectively. From the experience gained during previous satellite missions, it is known that the calibration parameters estimated on ground have to be re-estimated in flight. A description of the parameters defined for CHAMP and Ørsted satellites can be found in (Merayo et al. 2000; Olsen et al. 2003; Yin and Lühr 2011; Yin et al. 2013). In the case of Swarm satellites, the non-orthogonality angles and the offsets are not expected to change with time, whereas the scaling is likely to change slowly with time because of the ageing of the magnetometers. Due to the structure and mechanical properties of the magnetometers onboard Swarm satellites, it is expected that the rate of change of these scaling values is the same for the three orthogonal directions.

The calibration process requires that the strength of the magnetic field measured by the VFM matches the ASM scalar data. If data are collected over a full day, this requirement allows the estimation of snapshot values of these nine parameters. This is possible because the orientation of the measured magnetic field changes during the orbits and the day, relative to the sensors. Nonetheless, the estimated calibration parameters are not very robust and often only few of these nine parameters are estimated, or they are estimated on periods longer than one day. Technically, this match is obtained by least-squares, adjusting the parameters to minimise the sum of the squared differences between ASM and VFM field strengths (see Olsen et al. 2003; Yin and Lühr 2011). For a datum i, the relation between ASM and VFM data is:

$$ F_{\text{ASM}}(i) = \sqrt{B_{c1}^{2}(i) + B_{c2}^{2}(i) + B_{c3}^{2}(i)} + \epsilon(i), $$

((1))

where F _ASM are the scalar ASM data. B _c1, B _c2 and B _c3 are corrected vector VFM data in the E ₁, E ₂ E ₃ directions, respectively. ε is the difference between the ASM and corrected VFM field strengths. The relation between corrected VFM data and measured VFM data is:

$${} {\fontsize{8.8pt}{9.6pt}\selectfont{\begin{aligned} \left\{ \begin{array}{l} B_{c1} = \frac{B_{o1} - o_{1}}{s_{1}}\\ B_{c2} = \left[ \frac{B_{o2} - o_{2}}{s_{2}} - \sin a_{12} \; B_{c1} \right] \frac{1}{\cos a_{12}}\\ B_{c3} = \left[ \frac{B_{o3} - o_{3}}{s_{3}} - \cos a_{31} \, \sin a_{23} \; B_{c2} - \sin a_{31} \; B_{c1}\right] \frac{1}{\cos a_{23} \, \cos a_{31}} \end{array} \right. \end{aligned}}} $$

((2))

where B _o1, B _o2 and B _o3 are observed vector VFM data. We note that these relations are non-linear. Also, B _c2 depends on B _c1, and B _c3 depends on B _c1 and B _c2. These three relations become independent as soon as the non-orthogonality angles are zero. Slightly different, but equivalent, relations have been used for Ørsted, CHAMP and in official Swarm processing chains (Olsen et al. 2003; Yin and Lühr 2011). The inverse transform to compute the observed vector field from its corrected values is:

$${} {\fontsize{8.8pt}{9.6pt}\selectfont{\begin{aligned} \left\{\! \begin{array}{l} B_{o1} = s_{1} B_{c1} + o_{1}\\ B_{o2} = s_{2} \left[ \cos a_{12} \, B_{c2} + \sin a_{12} \, B_{c1} \right] + o_{2} \\ B_{o3} = s_{3} \left[\cos a_{31} \left[ \cos a_{23} \, B_{c3} + \sin a_{23} \, B_{c2} \right] + \sin a_{31} \, B_{c1} \right] + o_{3}\,{.} \end{array} \right. \end{aligned}}} $$

((3))

To illustrate the difficulties with Swarm data, we used the full Swarm A magnetic data set (baseline 0301/0302), where the usual processing steps that apply the calibration on a daily basis are dropped. Instead, the parameters estimated on ground are applied, and a constant scaling factor, the same for all three directions, is recalculated so that roughly acceptable values of field measurements, in nT, are obtained. The data are also selected to keep only valid data as indicated by the different flags. The derivation of this data set is described in Tøffner-Clausen 2014. The data set is very large (24,345,289 vector data values); therefore, the data are sub-sampled to one vector measurement every 20 s and rejected during obvious manoeuvres. This way we obtained a data set from 25 November 2013 to 30 September 2014, consisting of 1,188,891 vector data values. The data set is relatively clean, and only few noisy data can be identified after this process as it will become evident later. We then apply our own processing, as defined above, where we assume the offsets and angles to be constant in time, and a linear variation of the scaling with time. The calculated offset, scaling and angle values are given in Table 1. The residuals of the least-squares fit to the data are shown in Figure 1. We observed unexpected large residuals that vary rapidly in time. Along a single orbit, variations can be as large as 7 nT peak-to-peak. Such rapid variations can be explained neither by our processing methods nor with the daily estimates of the nine parameters in place in the operational Swarm magnetic processing chain.

Figure 1

Residuals of the least-squares fit to ASM data. Black dots: residual differences between ASM data and the magnetic field strength calculated from VFM data for the set of calibration parameters given in Table 1. The model does not include a Sun position dependence of the offsets. Red line: estimated local time for data points selected within 1° of the equator when the satellite is flying North.

Table 1 Offsets, scaling and angles obtained for models with and without Sun dependence

The main expected changes in the satellite environment along a single orbit are the magnetic field strength and direction, and also the temperature. The latter depends on the satellite lighting orientation. We therefore display the obtained residuals as a function of the Sun position relative to the satellite, where the Sun position is parameterised by two angles, α and β. The angle α gives the Sun direction angle relative to the E ₂ VFM direction, whereas the β angle is the angle between the −E ₁ VFM direction and the projection of the Sun direction on the E ₁ E ₃ VFM plane. The angle α is 90° and β is zero when the Sun is nearly above the satellite, shifted by 13° because of the angle between the boom and satellite body. The angles α and β are both 90° when the Sun is in the E ₃ direction. With such a parameterisation of the Sun position, the satellite local time (LT), shown in red in Figure 1, maps in the α angle. We therefore show in Figure 2 the residuals to the least-squares fit, displayed as a function of the Sun position for three different time spans, each corresponding to periods between dates where the satellite is flying in dawn-dusk orbits - i.e. 06h00, 18h00 LT.

Figure 2

Residuals to the least-squares fit displayed as a function of the Sun position. Same residuals as in Figure 1 that are now displayed as a function of the Sun position relative to the satellite and parameterised with the two angles α,β. The colour scales correspond to the residual amplitudes. The angles indicated on the Hammer projection plots are equal to 90° −α. The angles indicated on the polar plots correspond to the β values. The Sun is in the (−E ₁) direction - i.e. nearly above the satellite, when 90° −α=0 and β=0. The Sun is beneath the E ₂ E ₃ plane when β is outside the [−90:90] range. Top panel: residuals from 25 November 2013 to 25 February 2014. Central panel: residuals from 25 February 2014 to 10 July 2014. Bottom panel: residuals from 10 July 2014 to 30 September 2014. On the central panel, two obvious small-scale features are circled in black.

The three panels of Figure 2 display similar patterns, although the central panel is of the opposite sign. This is consistent with a small magnetic perturbation carried by the satellite, generated in the vicinity of the VFM sensors, and that depends on the Sun’s position. In such a scenario, the central panel where the satellite is in a descending mode on the dayside of the Earth - i.e. flying toward South on the dayside - should have roughly opposite sign anomalies compared to the two other panels where the satellite is flying North on the dayside of the Earth. We also observe that the anomalies differ dependent on the Sun being on one side of the satellite or the other - i.e. if α is smaller or larger than 90°. We note that the maximum perturbation is not observed when the Sun is just above the satellite, but rather slightly on the side. Finally, we see that there are areas, as those circled in black on the central panel in Figure 2, where the anomalies are displaying small-scale structures.

These observations suggested a simple modification of the offset parameterisation: to add a Sun position dependence to the otherwise constant offset values. This writes:

$$ \left\{ \begin{array}{l} o_{1}={o_{1}^{0}} + \sum_{l=1}^{30} \sum_{m=-l}^{m=l} o_{1}^{l,m} \; {Y_{l}^{m}}(\alpha,\beta)\\ o_{2}={o_{2}^{0}} + \sum_{l=1}^{30} \sum_{m=-l}^{m=l}o_{2}^{l,m} \; {Y_{l}^{m}} (\alpha,\beta)\\ o_{2}={o_{3}^{0}} + \sum_{l=1}^{30} \sum_{m=-l}^{m=l} o_{3}^{l,m} \; {Y_{l}^{m}} (\alpha,\beta) {\;,}\\ \end{array} \right. $$

((4))

where the ${Y_{l}^{m}} (\alpha,\beta)$ are spherical harmonics (SH). The parameters ${o_{1}^{0}}$, ${o_{2}^{0}}$ and ${o_{3}^{0}}$ are the constant offset values, and $o_{1}^{l,m}$, $o_{1}^{l,m}$, $o_{1}^{l,m}$ are the parameters for the Sun position dependence of the offsets. We choose a maximum SH degree of 30 arbitrarily. Recent experiments have shown that this maximum degree can be reduced (Lars Tøffner-Clausen, personal communication). In the results presented below, we also used a SH representation for the scaling factor up to SH degree 10, but this can probably be dropped if a sensor temperature dependence is assumed instead. Using this parameterisation of the anomaly, we proceed as before and adjust the set of parameters to minimise the differences between the ASM readings and the strength of the magnetic field observed by the VFM instruments through a least-squares fit.

The residuals of the least-squares fit are shown in Figure 3. They remain weak and do not present a structured signal, neither in time nor as a function of the Sun position. Some residuals are still large for very short periods of time, probably associated with satellite manoeuvres. The processing is therefore a success. The Sun-dependent parts of the offsets are displayed in Figure 4. The constant values obtained for the offsets, non-orthogonality angles and scaling, and the slope of the latter are given in Table 1.

Figure 3

Residuals of the fit to ASM data with a Sun position dependence of the offsets. Residual differences between ASM data and the magnetic field strength calculated from VFM data for the set of calibration parameters given in Table 1. The model does include a Sun position dependence of the offsets. The red line is the estimated local time for data points selected within 1° of the equator when the satellite is flying North.

Figure 4

Sun position dependence of the offsets. Perturbation vector strength in the E ₁, E ₂ and E ₃ directions in the VFM reference frame, as a function of the Sun position relative to the satellite. The projection and the definition of the angles are the same as in Figure 2.

The most interesting results are associated with Sun dependence of the offsets along the E ₁ direction. The Sun dependence shows two strong offset anomalies:

• A large negative anomaly when the Sun is nearly above and behind satellite A, slightly on the right side of the satellite when looking in the flight direction.

• A large positive anomaly, just before the Sun lowers below the horizontal plane of the satellite, on its left side.

The offsets in the E ₂ component also show a relatively large anomaly when the Sun is slightly behind the satellite on its right side when looking in the flight direction. These Sun-dependent offsets, which are supposed to be independent of time, correct the VFM data so that the anomalies presented in Figure 2 vanish.

Overall, the correction is successful. It can certainly be improved and stabilised, but what is described above was the best model available in December 2014. The differences between ASM and VFM values for the other two satellites are much weaker. Nonetheless, the processing described here leads also to a good fit between ASM and VFM data. With such a correction, the standard deviation of the differences between ASM values and the field strength, as estimated from VFM data, is around 0.17 nT for all three satellites.

Data selection

The magnetic field models were derived from the three Swarm L1b Baseline 0301/0302 satellite data series that have been processed and corrected as described in the previous section. The corrected magnetic field values in the VFM reference frame are rotated in the required reference frame - typically the Earth-centred, Earth-fixed, North, East, Center reference frame (NEC), when necessary. We also use observatory hourly means as prepared by Macmillan and Olsen 2013.

For all these data, the time is defined in Modified Julian Days 2000 (MJD) that counts the days since 1 January 2000 at 00:00h. The time period is limited from MJD 4749 to 5479 - i.e. 2013.0 to 2015.0. Over this time span, Swarm satellite data were available for:

• Satellite A: from MJD 5078.14 to 5369.87

• Satellite B: from MJD 5080.02 to 5369.88

• Satellite C: from MJD 5086.00 to 5369.87

The observatory data were available only up to MJD 5219.

The selection criteria used for these data are similar to those used in GRIMM series of models - e.g. Lesur et al. 2010. We recall these criteria for completeness.

The satellite and observatory vector data are selected in the solar magnetic (SM) coordinate system between ±55° magnetic latitude for magnetically quiet times according to the following criteria:

• Positive value of the Z-component of the interplanetary magnetic field (IMF B _z ),

• Sampling points are separated by 20 s at minimum,

• Data are selected at local time between 23 : 00 and 05 : 00, with the Sun below the horizon at 100 km above the Earth’s reference radius (a=6371.2 km),

• D s t values should be within ±30 nT and their time derivatives less than 100 nT/day, and

• Quality flags set to have accurate satellite positioning and two star cameras operating.

At high latitudes - i.e. polewards of ±55° magnetic latitude - the three-component vector magnetic data are used in the NEC system of coordinates. Their selection criteria differ from those listed above on the following way:

• Data are selected at all local times and independently of the Sun position.

We point here to the fact that provisional values of the D s t index are used for selection and modelling. The definitive D s t index values are likely to be different. At the time of data selection, no observatory data were available to define a more suitable selection index (e.g. the VMD index defined in Thomson and Lesur (2007)).

Model parameterisation

Away from its sources, the magnetic field can be described as the negative gradient of potentials associated with sources of internal and external origin:

$$ \begin{aligned} {} &\mathbf{B}= \mathbf{-\nabla} \{ V_{{\mathrm{i}}}(\theta,\phi,r,t) + V_{{\mathrm{e}}}(\theta,\phi,r,t) \} \\ {}&V_{{\mathrm{i}}}(\theta,\phi,r,t)= a {\sum\nolimits}_{l=1}^{L_{{\mathrm{i}}}}{\sum\nolimits}_{m=-l}^{l} \left(\frac{a}{r}\right)^{l+1} {g_{l}^{m}}(t) {Y_{l}^{m}}(\theta, \phi) \\ {}&V_{{\mathrm{e}}}(\theta,\phi,r,t)= a {\sum\nolimits}_{l=1}^{L_{{\mathrm{e}}}}{\sum\nolimits}_{m=-l}^{l} \left(\frac{r}{a}\right)^{l} {q_{l}^{m}}(t) {Y_{l}^{m}}(\theta, \phi) \end{aligned} $$

((5))

where ${Y_{l}^{m}}(\theta,\phi)$ are the Schmidt semi-normalised spherical harmonics (SH). θ,ϕ,r and a are the colatitude, longitude, satellite radial position and model reference radius, respectively, in geocentric coordinates. We use the convention that negative orders, m<0, are associated with sin(|m|ϕ) terms whereas null or positive orders, m≥0, are associated with cos(m ϕ) terms.

For the largest wavelengths of the field generated in the core and lithosphere (here, assumed up to SH degree L _i=18), the reference radius used in Equation 5 is a=3485 km. It corresponds to the radius of the core. This choice has no effect on the final field model unless regularisation or constraints are applied on the Gauss coefficients. These are parameterised in time from 2013.0 to 2015.0, using order six B-splines ${\psi ^{6}_{i}}(t)$, with half-year time intervals between spline nodes. The time dependence of the Gauss coefficients is therefore given by:

$$ {g_{l}^{m}}(t)=\sum_{j=1}^{N_{\mathrm{t}}}g_{lj}^{m} \, {\psi^{6}_{j}}(t), $$

((6))

where N _t=9. For the core and lithospheric field of SH degree greater than 18, the reference radius is set to a=6371.2 km. The maximum SH degree used for modelling the field of internal origin is 30, although a constant field, defined in Lesur et al. (2013), covering all SH degrees from 30 to 100 is subtracted from the data so that only very small contributions from the lithospheric field remain unmodelled. The remaining parts of the internal field model are the induced fields. They are modelled in a SM system of coordinates using four coefficients, for N _e=4 different 6-month time intervals, scaling the internal part of the D s t index - i.e. the I s t. These coefficients are the SH degree l=1 coefficients and the zonal SH degree l=2 coefficient. The magnetic potential for the induced field is therefore:

$${} \begin{aligned} &V_{{\text{induced}}}(\theta,\phi,r,t)= \\ &\qquad a {\sum\nolimits}_{j=1}^{N_{\mathrm{e}}} \left\{ {\sum\nolimits}_{m=-1}^{1} \left(\frac{a}{r}\right)^{2} g_{1j}^{m \, \text{Dst}} \, {Y_{1}^{m}}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}})\right.\\&\quad\qquad\qquad+\left. \left(\frac{a}{r}\right)^{3} g_{2j}^{0 \, \text{Dst}} \, {Y_{2}^{0}}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}}) \vphantom{\left\{ \sum_{m=-1}^{1} \left(\frac{a}{r}\right)^{2} g_{1j}^{m \, \text{Dst}} \, {Y_{1}^{m}}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}})\right.}\right\} \; {\mathcal{H}}_{j}(I\!st) \end{aligned} $$

((7))

where the function ${\mathcal {H}}_{j}(X)$ takes the value X in the time interval [t _j :t _j+1] and is zero otherwise. θ _S, ϕ _S are the colatitudes and longitudes in the SM reference frame. For observatory data, we also co-estimate crustal offsets.

The external field parameterisation also consists of independent parts. A slowly varying part of the external field model is parameterised over each 6-month time interval by a degree l=1 order m=0 coefficient in the geocentric solar magnetospheric (GSM) system of coordinates, and two coefficients of SH degree l=1 with orders m=0 and m=−1 in the SM system of coordinates. The rapidly varying part of the external field is controlled using the external part of the D s t index - i.e. the E s t - and the IMF B _y time series. Here again, 6-month time intervals are used. Four scaling coefficients for the E s t are introduced in each interval in the SM system of coordinates: three for SH degree l=1 and orders m=−1,0,1 and one for SH degree l=2 and order m=0. One scaling coefficient for the IMF B _y is introduced in each time interval for SH degree l=1 and order m=−1 in the SM system of coordinates.

Overall, the parameterisation of the external magnetic potential is:

$${} {\fontsize{7.9pt}{9.6pt}\selectfont{\begin{aligned} &{\mathbf V}_{\mathrm{e}}(\theta,\phi,r,t)= r {\sum\nolimits}_{j=1}^{N_{\mathrm{e}}} \left\{ \, q_{1j}^{0 \, \text{GSM}} \, {Y_{1}^{0}}(\theta_{\mathrm{G}}, \phi_{\mathrm{G}})\right\} \; {\mathcal{H}}_{j}(1) \\ &\quad + r {\sum\nolimits}_{j=1}^{N_{\mathrm{e}}} \left\{ \,q_{1j}^{0 \, \text{SM}} \, {Y_{1}^{0}}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}}) + \; q_{1j}^{-1 \, \text{SM}} \, Y_{1}^{-1}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}}) \right\} \; {\mathcal{H}}_{j}(1) \\ &\quad + r {\sum\nolimits}_{j=1}^{N_{\mathrm{e}}} \left\{\! {\sum\nolimits}_{m=-1}^{1}\! q_{1j}^{m \, \text{Dst}} \, {Y_{1}^{m}}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}}) \,+\, \left(\frac{r}{a}\right) q_{2j}^{0 \, \text{Dst}} \, {Y_{2}^{0}}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}})\!\right\} \! {\mathcal{H}}_{j}(E\!st) \\ &\quad + r {\sum\nolimits}_{j=1}^{N_{\mathrm{e}}} \left\{ \,q_{1j}^{-1\, \text{IMF}} \, Y_{1}^{-1}(\theta_{\mathrm{S}}, \phi_{\mathrm{S}}) \right\} \; {\mathcal{H}}_{j}(\text{IMF} \, {B_{y}}) \end{aligned}}} $$

((8))

θ _G, ϕ _G and θ _S, ϕ _S are the colatitudes and longitudes in GSM and SM reference frames, respectively.

We used independent external field parameterisations for the satellite and observatory data. For the latter, we impose that $q_{1j}^{0 \, \text {SM}}$ is set to zero to avoid co-linearities with the observatory crustal offsets.

Independently of the parameterisation of the magnetic field, we also want to estimate the so-called Euler angles between the VFM orthogonal set of measurements and the star camera reference frame, such that the measured vector field in the VFM reference frame can be mapped into an Earth-centered, Earth-fixed reference frame. The latter reference frame is usually the NEC system of reference. We actually only compute corrections of predefined Euler angles, for a series of 30-day windows. The parameterisation and algorithm we used are detailed in Rother et al. (2013) and are not repeated here.

Inversion process

The relationships between data and model are given by Equations 5 to 8 and lead to a linear system of equations:

$$ {\mathbf d} = {\mathbf A}_{d} \, {\mathbf g} + {\mathbf \epsilon} $$

((9))

to be solved, where d is the data vector, g is the vector of Gauss coefficients defining the model and ε is this part of the data that cannot be explained by our model. We note that for mid-latitude data, the Equations 5 have to be rotated into the SM system of coordinates because the data at mid-latitudes are selected in that system. However, for the internal part of the model, the data in the direction Z _SM are not used as they are strongly contaminated by the external fields. This approach has been used in all models of the GRIMM series. More details are provided in Lesur et al. (2008).

As we try to estimate a core field model with a relatively complex behaviour in time up to SH degree 18, the model is not fully resolved by the data and needs regularisation. The regularisation chosen minimises the integral of the squared second time derivative of the core field model radial component over the model time span. Further, at two single epochs close to the end point of the model, the integral of the squared first time derivative of the radial component of the core field model is minimised. This regularisation method defines three functionals at the core radius c:

$${} \begin{aligned} &\Phi_{{1\,\text{or}\,2}}(t)=\frac{1}{4\pi \, c^{2}} \int_{\Omega(c)} \left| \partial_{t} B_{r}(\theta, \phi, r, t) \right|^{2} \; d\omega \\ &\text{for } t=t_{1} \text{ and } t=t_{2}, \text{and } \\ &\Phi_{{3}}=\frac{1}{4\pi \, c^{2} \, (t_{2}-t_{1})}\int_{t_{1}}^{t_{2}} \int_{\Omega(c)} \left| {\partial_{t}^{2}} B_{r}(\theta, \phi, r, t) \right|^{2} \; d\omega \, dt, \end{aligned} $$

((10))

which have to be minimised for the core field model Gauss coefficients. Unlike Lesur et al. (2010), we control the first and second time derivatives of the radial field, as these time derivatives are poorly resolved by the brevity of input data. Ultimately, this leads to three systems of linear equations, one for each functional:

$$ {\mathbf 0}= \lambda_{i} \; {\mathbf L}_{Bi} \, {\mathbf g} \hspace{1cm} i=1,2,3. $$

((11))

Equation 9 and 11 should be solved simultaneously by least-squares where the λ _i are scalars that need to be adjusted so that the satellite and observatory data are fit to their expected noise level.

Since the Swarm satellite data set contains data during manoeuvre days and since the most recent observatory data are ‘only’ quasi-definitive data, we expect a higher number of outliers than usual and a general distribution of residuals deviating significantly from a Gaussian distribution. We therefore use a re-weighted least-squares algorithm where at iteration j+1, the datum value d _i in Equation 9 is associated with a weight $w_{i}^{j+1}$ where:

$${} {w}_{i}^{j+1} = \left\{ \begin{array}{ll} \frac{1}{\sigma_{i}} &\text{for } |{d}_{i}-{\mathbf A}_{d} \cdot {\mathbf g}^{j}| \le k_{i} \, \sigma_{i},\\ \frac{1}{\sigma_{i}} \left[\frac{k_{i} \, \sigma_{i}}{|{d}_{i}-{\mathbf A}_{d} \cdot {\mathbf g}^{j}|} \right]^{1-\frac{a_{i}}{2}} \;\;\;\; & \text{for } |{d}_{i}-{\mathbf A}_{d} \cdot {\mathbf g}^{j}| > k_{i} \, \sigma_{i}. \end{array} \right. $$

((12))

The g ^j are the set of Gauss coefficients obtained at iteration j. σ _i is the prior standard deviation of the noise associated with the datum d _i . The control parameters σ _i , k _i and a _i are set before starting the iterative process and are given in Table 2. They depend on the data type. The iterative process is started by setting a _i =2 for all data points - i.e. assuming a Gaussian distribution of residuals.

Table 2 Satellite and observatory data weight parameters, residual means and rms

At each iteration, we therefore minimise the functional:

$${} \Phi_{j} = \left[ {\mathbf d}-{\mathbf A}_{d} \, {\mathbf g}^{j} \right]^{t} {\mathbf W}^{j} \left[ {\mathbf d}-{\mathbf A}_{d} \, {\mathbf g}^{j} \right] + \sum_{i=1}^{3} \lambda_{i} \left[{\mathbf L}_{Bi} \, {\mathbf g}^{j} \right]^{t} \left[{\mathbf L}_{Bi} \, {\mathbf g}^{j} \right]\!, $$

((13))

where j is the iteration number, and W ^j is a diagonal matrix which elements are $({w}_{i}^{j})^{2}$, with ${w}_{i}^{j}$ defined by Equation 12.

Once an acceptable solution is obtained, a threshold filtering is applied to get rid of some remaining outliers. Then, the Euler angles are estimated and the inversion process is restarted. The final solution is reached when the Euler angles are not changing significantly from one iteration to the next.

Results and discussion

In Table 2, we present the fit to the data for all the used subsets of data, together with the control parameters σ _i , k _i and a _i . These are specified independently for each of the satellites when necessary. The residual rms for the low-latitude magnetic field components X _SM, Y _SM and Z _SM data are in the same range as for the GRIMM models based on CHAMP satellite data (see Lesur et al. 2015). The Z _SM component is included, as it enters in the Euler angles and external field parameter estimation, but it does not affect the estimation of the fields of internal origin directly. The misfits for the high-latitude data are slightly smaller than those obtained with CHAMP data for similar model parameterisation. It is not clear yet why this occurs. The magnetic data selected for the extended local time window are only used for the Euler angle estimation. Since early evening and late morning data are included, the fit to the Z _SM component is slightly degraded. Nonetheless, the fit to these data remains surprisingly good.

The misfits for the observatory data are comparable to those of the satellite data, although slightly degraded for the mid- and low-latitude SM components. The mean of the residuals is relatively high for the observatory X _HL component which is an indication of the strongly non-Gaussian residual distribution. Overall, accounting for the corrections applied to satellite data and the obligation to use preliminary D s t index values for the selection and modelling, the fit to the data is acceptable.

Figure 5 shows the histograms for three of the satellite data components. The residuals for other data types show similar distributions. The distributions, plotted in semi-logarithmic scale, clearly show the deviation from the Gaussian distribution. Our chosen weights correspond to a Gaussian distribution of residuals for small errors, and significant tails for large errors as shown by the blue dotted line. Our chosen weights and associated distributions of residuals allow prior error distributions that are compatible with the posterior distributions. It proves to be an important step for obtaining high-quality magnetic field models.

Figure 5

Residual histograms for three of the satellite data components. Histogram of the residuals for the X _SM, Y _SM and Z _HL data types, scaled by their prior σ given in Table 2. Also, the Gaussian distribution is given in green, and the distributions corresponding to the weights in Equation 12 for large residuals are shown in blue.

The power spectra of a snapshot of the model core field and secular variation are presented in Figure 6 for the epoch 2014.0. The epoch 2014.0 corresponds to a time where satellite data and observatory data coverage overlaps. The results discussed here are labelled DCO(o) for dedicated core field where ‘o’ indicates the use of observatory data for a modelling approach largely based on development made in the Swarm SCARF DCO L2 framework (Rother et al. 2013). The power spectra of the IGRF-GFZ candidates are also displayed for comparison. The IGRF main field model candidate has been built from its parent model snapshot for 2014.0, forwarded in time by one year using the averaged SV of the parent model from 2013.5 to 2014.5. The IGRF SV candidate for 2015 to 2020 is this average SV derived from the parent model. While the differences for the two main field models in Figure 6 only reflect the SV between 2014.0 and 2015.0 (the reference date for the IGRF), some significant differences between parent SV models appear above SH degree 8 (not shown in Figure 6). Besides the data updates since the IGRF candidate submission, these differences may be caused by changes in the preliminary D s t values, the corrections applied on Swarm data and of course slightly different values for the λ controlling the regularisation. These increasing differences from SH degree 8 give some indications on the SV model robustness. The poor modelling of the SV above SH degree 8 is to be expected given the limited time span of available Swarm and observatory magnetic field vector data.

Figure 6

Power spectra of the modelled magnetic fields. Power spectra of the magnetic field model and its SV, calculated for 2014.0 at the Earth’s surface. The spectra for the model presented in this paper - i.e. DCO(o) - for the GFZ-IGRF candidates and the power of the differences between the two models are also shown.

Figure 7 presents maps of the vertical down component of the core field model and its SV for epoch 2014.0. The main field model map shows the usual reverse patches in the southern hemisphere and the undulations of the magnetic equator. A close comparison with equivalent maps for 2005 shows small evidence of westward drift in the southern Atlantic (not shown). Similarly, the SV map displays the usual weakness of the SV over the Pacific and Antarctic regions. With only one year of data, it is not possible to obtain reasonably accurate models of the acceleration. We can nonetheless estimate the average acceleration since the CHAMP epoch at large scales. This is displayed in Figure 8, where the acceleration has been calculated at the Earth’s surface by scaled differences between GRIMM-4 model SV for 2005 (Lesur et al. 2015) and our SV estimate for 2014.0. Interestingly, large peaks of acceleration during the CHAMP era, first shown in Lesur et al. (2008) and as large as 30 nT /y ², average-out over several years. These acceleration pulses (as described in Chulliat et al. (2010)) can then be interpreted as magnetic field short-term disturbances over an otherwise smooth field evolution. The only region where this average acceleration displays significant values is under Eastern Asia. Yet, it is not clear if that is a transient or permanent feature of the geomagnetic field.

Figure 7

Maps of the magnetic field model and its SV. Maps of the vertical down component of the magnetic field model for year 2014.0 at the core mantle boundary. Top: snapshot of the core field. Bottom: snapshot of the SV.

Figure 8

Map of the modelled averaged acceleration. Map of the vertical down component of the magnetic field average acceleration between 2005.0 and 2014.0 at the Earth’s surface.

Figure 9 shows time series of the Euler angle corrections for 30-day segments, estimated over the satellite data time span, for all three satellites. The Euler angles are constant during each 30-day period, but otherwise, no constraints have been applied to limit their amplitudes or variations. Although these corrections never exceed 20 arc sec, the variability of the Euler angles is surprisingly high. This is very likely not due to a lack of stability of the optic bench, which rigidly links the VFM sensors to the star cameras that define the spacecraft reference frame. The apparent correlation of the angle variations after the first third and before the last third of the period - i.e. at MJD ≃5170 and MJD ≃5300 - is clearly associated with the changes in northward versus southward flight directions of the satellites at these epochs in the selected nighttime data set. As indicated in Rother et al. (2013), the Euler angles may accumulate and respond to remaining noise or signals due to various sources, as well as hidden timing errors or poor/incomplete modelling of the external fields. The local time dependence is less obvious than it was for the CHAMP data, possibly because not enough Swarm data have been accumulated yet. For CHAMP data, the corrections were also of much larger amplitudes (up to 50 arc sec).

Figure 9

Euler angle corrections. Time series of estimated Euler angle corrections for constant rotations from VFM into S/C system over 30-day segments. There are three angles for each satellite. Each symbol defines the start of a 30-day interval.

Conclusions

In the view of providing candidates for the new IGRF-12, we processed a combination of Swarm satellite and observatory data. In order to obtain robust models of the main field, we first applied a correction to the satellite data such that the strength of the magnetic field, as measured by the VFM instrument, matches the data measured by the absolute scalar instrument. Our underlying hypothesis is that the discrepancies between these two instruments are closely linked to the position of the Sun relative to the satellite. So far, the corrections obtained with such a hypothesis seem suitable. The results presented here are actually part of a longer term study, made in close collaboration with the European Space Agency and several other European institutions. The ultimate goal is to identify the original source of this perturbative signal such that, firstly, the Swarm satellite data can be corrected to remove that signal and, secondly, a new design can be adopted in future satellite missions to avoid such difficulties. As stated in the text, the model presented here is not optimal but corresponds to our best results available at the end of 2014.

Using these corrected data, we have been able to derive accurate models of the geomagnetic field around epoch 2014.0. Since less than a year of satellite data had been accumulated, the SV model - i.e. the magnetic field linear variation in time - is only accurate for the longest wavelengths (up to SH degree 8). We therefore cannot compare meaningfully maps of the SV at the CMB with those obtained at earlier epochs. Nonetheless, the average acceleration of the magnetic field - i.e. its second time derivative - between 2005 and 2014.0 is remarkably weak at the Earth’s surface compared with the acceleration values obtained during the CHAMP era. This observation strongly supports the view that the observed acceleration peaks, which have been associated with magnetic jerks in 2003, 2007 and 2010, correspond to short-term disturbances of the field over an otherwise slowly and, most of the time, smoothly evolving magnetic field.