NOAA/NCEI and University of Colorado candidate models for IGRF-13

The International Geomagnetic Reference Field (IGRF) is a set of parameters representing the large-scale internal part of Earth’s magnetic field. The 13th generation IGRF requested candidate models for a definitive main field for 2015.0, a provisional main field for 2020.0, and a predictive secular variation covering the period 2020.0–2025.0. The University of Colorado (CU) and the National Centers for Environmental Information (NCEI), part of the National Oceanic and Atmospheric Administration (NOAA), have produced these three candidate models for consideration in IGRF-13. In this paper, we present the methodology used to derive our candidate models. Our candidates were built primarily from Swarm satellite data, and also relied on geomagnetic indices derived from the ground observatory network. The ground observatories played a crucial role as independent data in validating our candidates. This paper also provides a retrospective assessment of the CU/NCEI candidate model to the previous IGRF (IGRF-12) and discusses the impact of differences between candidate and final IGRF models on global model errors.


Introduction
The International Geomagnetic Reference Field (IGRF) is a mathematical representation of the large-scale timevarying portion of Earth's internal magnetic field. It is regularly updated (typically every 5 years) by an international collaboration of research scientists under the leadership of the International Association of Geomagnetism and Aeronomy (IAGA, http://www.iaga-aiga.org) Division V Working Group on geomagnetic field modeling (V-MOD, https ://www.ngdc.noaa.gov/IAGA/vmod). The IGRF is used routinely in geomagnetic research as well as operational navigation and pointing systems. In 2019, the V-MOD Working Group issued an international call soliciting candidate models for the 13th generation IGRF (hereafter IGRF-13). The call requested candidates for the main field for epochs 2015.0 and 2020.0, and a predictive secular variation for the period 2020.0-2025.0. The 2015.0 main field candidates were combined into a final model which has now replaced the older 2015.0 model from IGRF-12 (Thébault et al. 2015) and is called a Definitive Geomagnetic Reference Field (DGRF-2015) as it is not expected to be updated again. The 2020.0 main field candidates were also combined into a final model, which will likely be updated in a future generation of IGRF, and so it is named IGRF-2020. A summary of the final IGRF-13 model may be found in Alken et al. (2020c) and the procedures used to evaluate all submitted candidates and combine them into a final IGRF-13 are detailed in Alken et al. (2020b).
In this paper, we describe the three candidate models developed by the National Centers for Environmental Information (NCEI) and the University of Colorado (CU) for IGRF-13. Our modeling methodology follows the approach used for our contributions to previous IGRF generations Maus et al. 2005bMaus et al. , 2010 with some differences which will be described below. Our three candidate models were derived from two parent models, which included additional parameters (beyond the requested IGRF internal field parameters) in order to account for additional geomagnetic field sources present in near-Earth measurements. In "Methodology" section, we discuss the parameterization of our parent models and the data sources used to determine the model parameters. In "Results" section, we discuss the results of our model inversions. In "Validation" section, we present our validation of our modeling results, which is followed by a retrospective assessment of the CU/NCEI candidate model to the previous IGRF (IGRF-12) in "Retrospective assessment" section.

Field modeling
We start from the magnetic potential due to an internal source, where r, θ, φ are radial distance, co-latitude, and longitude in a geocentric spherical coordinate system, t is time, a = 6371.2 km is the approximate mean Earth radius, g m n (t), h m n (t) are the time-dependent Gauss coefficients representing the internal field, N is the series truncation limit, for which we choose a value of 50, and P m n (cos θ) are the Schmidt semi-normalized associated Legendre functions. We represent the time dependence of the Gauss coefficients up to spherical harmonic degree 15 using basis splines (B-splines), where N j (t) are cubic B-splines with no interior knots over the model interval. Since we do not place interior knots, this effectively makes each g m n (t), h m n (t) a quadratic polynomial over the model interval. For degrees 16 through N, we use static coefficients g m n (t) = g m n and h m n (t) = h m n to represent the time-independent lithospheric field.
When fitting satellite measurements, it is also important to account for a low SH degree external field, B ext (r, t) representing magnetospheric ring and tail current sources. We do not co-estimate an external field for our IGRF-13 candidates, instead using the CHAOS-6 external field model (Finlay et al. 2016;Olsen et al. 2014) as a priori information, which is parameterized by the RC index, which tracks the symmetric part of the magnetospheric ring current and is derived from a network of low and mid-latitude observatories. The model which we then fit to the satellite observations is given by where B int (r, t) = −∇V int (r, t) represents the internal field sources parameterized by the Gauss coefficients as discussed above.
The Swarm vector measurements are recorded by a fluxgate magnetometer co-located with three star cameras on an optical bench system situated on a boom. Since the Swarm attitude parameters are referenced to the star camera frame, also known as the common reference frame (CRF), the measurements recorded in the vector fluxgate magnetometer (VFM) frame need to be rotated to CRF before completing the rotation to an Earth-fixed frame with the attitude rotation. Therefore, the vector measurement in a North-East-Center (NEC) Earth-fixed frame is given by where B VFM i is the VFM measurement made at (r i , t i ) , R NEC CRF represents the attitude rotation from CRF to NEC and is provided in the Swarm Level 1b data files as a time series of quaternions, and R CRF VFM is a rotation from VFM to CRF, parameterized by three time-dependent Euler angles, α(t), β(t), γ (t) . We use a 1-2-3 Euler angle convention, and our rotation matrix R CRF VFM can be found in Alken et al (2020a, Eqs. 23−26). We use B-splines for the time dependence of the Euler angles, where α(t) = (α(t), β(t), γ (t)) is the vector time series of Euler angles and similarly for the coefficients α j . Here, we choose the N j (t) to be cubic B-splines with 30-day separation between knots. This allows a slow temporal change in the alignment parameters which could be caused by thermal variations on the satellite. The Gauss coefficients parameterizing the internal field are then computed by minimizing the cost function where (4) and F i is the Swarm scalar field measurement made at (r i , t i ) . The weights w i are designed to achieve a uniform weighting over all latitudes and longitudes, and are detailed in Alken et al (2020a, p. 15). The minimization problem is nonlinear, and so we use a Levenberg-Marquardt iterative method to determine the final set of Gauss coefficients and Euler angles.

Data selection and preprocessing
We first partition the Swarm A and B data into halforbital tracks, which span from one pole to the other. For each track, we compute a root-mean-square (rms) difference with respect to spherical harmonic degrees 1-15 of the CHAOS-6-x9 core field model (Finlay et al. 2016) for each vector component as well as the total field component. Tracks whose rms difference exceeds 200 nT in any component are discarded from the analysis. This procedure aims to remove anomalous data recorded during satellite maneuvers. CHAOS-6-x9, which is built from data through April 2019, will contain secular variation forecast errors when applied to measurements in late 2019; however, the chosen threshold of 200 nT far exceeds these errors, so the procedure will only discard measurements contaminated by non-geophysical noise. Next, we select data for geomagnetically quiet periods using the Kp and RC indices (Olsen et al. 2014) with the following criteria, • Kp index does not exceed 2 from a time period starting 2 h before the track start time to the track end time. • The temporal change of the RC index |dRC/dt| does not exceed 3 nT/h from 3 h prior to the track start time to the track end time.
Next, we remove ionospheric field contributions on the dayside with the following criteria, • For quasi-dipole (QD) (Richmond 1995) latitudes equatorward of ±55 • , the track's local time of ascending or descending node is between midnight and 05:00. • For QD latitudes poleward of ±55 • , the sun is at least 10 • below the horizon, according to the zenith angle. • For QD latitudes poleward of ±55 • , we use only scalar field measurements from Swarm to minimize perturbations caused by the high-latitude field-aligned current systems. Below ±55 • , we use all components (vector and scalar).

Parent models
We built two parent models from two different time periods, one for the DGRF-2015 main field candidate, and the other for the IGRF-2020 main field and secular variation candidates. For the DGRF-2015 parent model, we used a symmetric time window centered on 2015.0. Since the first Swarm data is available on 26 November 2013, our window included data until 5 February 2016. Our final DGRF-2015 coefficients were determined by evaluating the Gauss coefficients g m n (t), h m n (t) at t = 2015.0 and truncating them to SH degree and order 13. For the IGRF-2020 parent model, we used a 3-year time window from 16 September 2016 to 16 September 2019. In our previous experience, at least 3 years of data are needed to allow the core field to change enough in order to resolve the smaller high SH degree and order Gauss coefficients and their secular variation and acceleration. Our secular variation model for 2020.0-2025.0 was taken to be the first time derivative of the Gauss coefficients evaluated at the timestamp of the last data point available in our time window (16 September 2019). Denoting this time as t f , the secular variation coefficients are then ġ m n (t f ),ḣ m n (t f ) , truncated to SH degree and order 8, which was requested by the IGRF-13 call. The main field coefficients for 2020.0 were computed by performing a linear extrapolation of our Gauss coefficients from t f to 2020.0 using the secular variation at t f , These IGRF-2020 coefficients were then truncated to SH degree and order 13 for our final IGRF-13 candidate model. Figure 1 shows the model residuals ǫ i , f i of both parent models as spatial maps, combining Swarm A and B together. We gridded the residual data into 2 • × 2 • latitude/longitude bins and plot the median of each bin. For the DGRF-2015 parent model, we find systematically larger residuals in the Atlantic region in all components. We do not have a full explanation for this, but Alken et al. (2020b, Fig. 4) shows similar features in the B z component for all DGRF-2015 candidate models, which they attribute to rapid field variations around 2015 caused by the widely reported geomagnetic jerk occurring in late 2014 (e.g., Torta et al. 2015). Our parameterization of the Gauss coefficients with quadratic polynomials over 2 years and

Results
2 months may be insufficient to capture the full spatiotemporal behavior of this jerk event. We also note larger residuals in the B y component for both the DGRF-2015 and IGRF-2020 parent models, which is due to the highlatitude ionospheric currents which are not fully removed from the Swarm data. Finally, the total field residuals |B| exhibit structures at low latitudes, which we attribute to ring current signals which are not fully removed from the data by the CHAOS-6 magnetospheric model. We present the residual statistics of both parent models in Table 1.

Comparison with observatory network
We compared the prediction of our SV model with independent data from the global geomagnetic observatory network, using the hourly mean database compiled by the British Geological Survey (Macmillan and Olsen 2013). We limited the observatories to those located within ±55 • geomagnetic dipole latitudes to exclude signals originating in the high-latitude ionosphere, which are difficult to model and remove from ground-based data. Observatories with baseline errors and large data gaps were omitted from the analysis. A total of 64 geomagnetic observatories distributed across the world were used for the analysis, as shown in Fig. 2. The median date of the latest data available is 30 June 2019. We select data for geomagnetically quiet conditions (Ap index ≤ 10) and discard data outside a local time window of 00:00-05:00, to minimize dayside ionospheric signals in the data. We then fitted cubic splines with knots separated by 1 year separately to the B x , B y , and B z components. The secular variation of the geomagnetic field at an observatory for a  where S i (t) is the vector spline fitted to the B x , B y , B z data for observatory i. We linearly extrapolate the splines to 2020.5 to obtain SV estimates at the model epoch 2020.0. An example of the spline fitting for the Honolulu (HON) observatory is given in Fig. 3. HON is located in a region which has experienced strong secular acceleration over the past 15 years, which makes it challenging to predict future magnetic field variations. The right panels show, for the B x , B y and B z components, the raw data (blue), filtered data (green), and spline fit S i (t) (red). The left panels show, for the same components, the annual differences of filtered hourly means (black), spline secular variation Ṡ i (t) (red curves), and the SV prediction of our IGRF-13 candidate model at 2020.0 (red circles). The vertical green line indicates the time of the last available SV estimate which did not rely on spline extrapolation for this station. Table 2 presents the means and standard deviations of the differences between our SV model and spline fits for epoch 2020.0 at the 64 geomagnetic observatories shown in Fig. 2. The mean and standard deviation values reflect errors both in our SV model and in the extrapolation procedure used for the observatory splines. Swarm data were available until 16 September 2019, while observatory data were only available until 30 June 2019. Therefore, the extrapolation of the observatory splines will likely be less accurate than extrapolation of the IGRF-2020 model to the epoch 2020.0. Nevertheless, the observatory analysis provides an independent (11) S i = S i (t + 6 months) − S i (t − 6 months), verification of the candidate SV model presented in this paper.

Comparison with selected observatory data
In another validation step, following a procedure already used in the previous IGRF preparation , we compared our IGRF-13 candidate SV with the observed SV at selected INTERMAGNET observatories where quasi-definitive data (see Chulliat et al. 2017, and references therein) are available. When validating the model during the IGRF-13 preparation, we used data until November 2019. We extended this analysis with data until July 2020 while writing this paper. Observatory 1-min data were downloaded from INTERMAGNET (ftp://ftp.seismo.nrcan.gc.ca/intermagnet) and pre-processed as follows. At low and mid-latitudes (quasi-dipole latitudes equatorward of ±55 • ), data for solar elevations greater than −10 • were removed. At all latitudes, data such that Kp > 2o , |Dst| > 20 nT and |dDst/dt| > 2 nT/h were removed. We then calculated daily means and monthly means from the selected night-time and geomagnetically quiet time minute data. The SV was obtained by taking differences between monthly means at +6 months and −6 months. Examples of comparisons between the CU/NCEI candidate SV and data collected at the Guam (GUA) and Belsk (BEL) observatories are shown in Figs. 4, 5. These Figures also show model predictions from recent (predictive) IGRF SV models (from IGRF-10 to IGRF-13), the most recent (retrospective) DGRF (released as part of IGRF-13) and the (predictive) CU/NCEI SV candidate to IGRF-12.

Fig. 2 Observatories used in the model validation
The SV at Guam (Fig. 4) underwent a strong and sudden acceleration after mid-2017 on the B z component. A much more moderate acceleration occurred on B x and there was no acceleration at all on B y . As a result, the predictions from IGRF-12 and the CU/NCEI candidate to IGRF-12 were much more accurate on B x and B y than on B z . Due to the strong acceleration, even the DGRF is a poor approximation of the SV at Guam over 2015-2020. The CU/NCEI candidate to IGRF-13 is in good agreement with the SV on B x and B z at the time the model was calculated, i.e., at the end of 2019; the agreement is slightly less good for B y . The IGRF-13 SV is closer to SV values at the end of 2018 for all components. Recent data at this observatory (until and Fig. 3 An example of SV determination at Honolulu (HON). The red lines on the left panels represent the SV spline model with annual differences of filtered observatory data in black. Right panels show the raw data (blue), filtered data (green), and spline fits (red). The circles in the left panel show our predicted SV at 2020.0. The vertical green lines indicate the timestamp of the latest available SV estimate which did not require spline extrapolation including July 2020) suggest that the acceleration on the B z component remained strong and almost constant until recently, as slightly increased on B x . As a result, our IGRF-13 SV candidate currently provides a more accurate SV estimate on B z and B x at Guam than IGRF-13. The opposite situation seems to be happening at Belsk (Fig. 5). There, the SV changed at a slower pace over the past 15 years. Over the past 6 months, the acceleration changed sign on the B y component, and perhaps also on B x and B z (but this remains to be confirmed). As a result, on B x and B y , our candidate SV to IGRF-13, which by design was close to the SV at the end of 2019, is currently less accurate than the IGRF-13 SV, which was closer to the SV at the end of 2018. On the B z component, both models overshoot the current SV for an unknown reason.
The observatory comparisons shown here (and the ones not shown) usefully complement the quantitative comparisons presented in the previous subsection. Specifically, we checked that our candidate SV model didn't provide predictions too far from the SV calculated from the latest available INTERMAGNET data, and from the SV that could be expected if current trends (acceleration) were to continue in the near future.

Retrospective assessment
A retrospective analysis of all IGRF-12 candidate secular variation models was performed by the IGRF-13 task force and is reported in Alken et al. (2020b). In this analysis, the secular variation forecast by each candidate model was compared with the secular variation derived from actual measurements at ground-based observatories over 2015-2020. Results (see Table 6 in Alken et al. 2020b) show that the CU/NCEI candidate to IGRF-12  performed slightly worse than the final IGRF-12 SV model: root-mean-square (RMS) differences were 12.10 nT/year (vs. 11.47 nT/year) on the B x component, 11.22 nT/year (vs. 10.75 nT/year) on B y , and 19.2 nT/year (vs. 18.50 nT/year) on B z . We attribute this relative under-performance to the 2014 geomagnetic jerk (Torta et al. 2015) and the subsequent pulse of secular acceleration around 2016−2017 (Alken et al. 2020a). Like in this cycle's model (cf. "Parent models" section above), no attempt was made to forecast the future SV in the CU/NCEI candidate to IGRF-12. The candidate SV was taken to be the truncated parent model SV at epoch 2014.3 (which was the parent model epoch). As a result, the change in SV associated with the 2016-2017 pulse, which started around 2014, led to larger SV errors in the CU/NCEI candidate model than in some models which attempted some forecast and happened to get the acceleration partially right. An example of this can be seen at Belsk on the B z component between 2015 and 2020 (Fig. 5). As long as there exists no consistently superior SV forecasting method, there is value in averaging multiple models based on different modeling philosophies when building the IGRF. When using the IGRF as a geomagnetic reference model (e.g., for navigation or orientation), what matters is the total model error, i.e., the differences between model predictions and the actual geomagnetic field at a given time and location. The total error is the sum of (a) omission errors due to un-modeled fields (mostly crustal and external) and (b) commission errors due to the imperfect modeling of the core field and its secular variation. The cumulative effect of the SV error discussed above is the largest contribution to the commission error; it generally grows over time, unless the secular acceleration changes sign during the 5-year cycle. When the IGRF is updated, the SVs for past cycles (generally the last two cycles) are updated using all the data available, which has the effect of reducing not only the SV error but also the total error. Figure 6 shows this effect for the last three IGRF cycles. We plotted the declination root-mean-square (RMS) error at high latitudes in the Northern Hemisphere (latitudes greater than 55 deg.) for IGRF-10 (forecast only, 2005, Maus et al. 2005a, IGRF-11 (forecast only, 2010, Finlay et al. 2010, IGRF-12 (forecast only, 2015-2020, Thébault et al. 2015 and IGRF-13 (past and current epochs only, 2005-2020, Alken et al. 2020c). We used the omission error value calculated as part of the error assessment of the World Magnetic Model (WMM), a degree and order 12 spherical harmonic model of the main field and its secular variation developed by NCEI and the British Geological Survey (Chulliat et al. 2020). The difference in omission error caused by slightly different truncation degrees in WMM and IGRF is assumed to be small, as most of the omission error comes from higher degree and order crustal fields as well as external fields. The commission error was calculated by comparing each model with a retrospective, satellite-based, time-varying core field model over 2005-2020 developed using a similar approach as the one used to develop parent models for this IGRF candidate. As the declination is generally not reliable in the vicinity of magnetic dip poles (Chulliat et al. 2020, and references therein), the area where the horizontal component is smaller than 2000 nT was excluded from the error calculation. The total declination error in the Northern Polar Cap was largest towards the end of the 2015-2020 cycle, reaching almost 1 degree just before the release of IGRF-13. We attribute this larger error to the 2016 secular acceleration pulse following the worldwide 2014 jerk (Torta et al. 2015). This pulse occurred at the very beginning of the 2015-2020 interval and had a larger cumulative effect over 5 years than the 2006, 2009 and 2012 pulses  which occurred during previous cycles. The declination error was about 25% smaller at the end of the 2005-2010 and 2010-2015 cycles. However, each IGRF update was successful in reducing the commission error to a level much smaller than that of the omission error, and as a result the total error of IGRF-13 over 2005-2020 is very close to its floor, which is the omission error. The slightly larger mid-cycle errors are likely due to the non-linear SV (i.e., the acceleration), which cannot be accounted for by the IGRF model as long as this model is linear. Figure 6 also shows the time evolution of the total declination error in the Northern Polar Cap for the CU/ NCEI IGRF-12 candidate model. There is very little difference (less than 0.01 deg.) between the total error when using this candidate model and the final IGRF-12. We found similar results for other components (I, F, and D at mid-latitudes and in the Southern Polar Cap, not shown).

Fig. 6
Estimated declination root-mean-square error in the Northern Polar Cap (latitudes greater than 55 deg., except the area where the horizontal component is smaller than 2000 nT) for recent IGRF SV models (IGRF10, IGRF11 and IGRF12), the most recent DGRF (IGRF13) and the CU/NCEI SV candidate model to IGRF12 We conclude that the differences in SV error between IGRF-12 and the CU/NCEI IGRF-12 candidate model (and likely most other candidate models) are insignificant from a practical standpoint.