Evidence-based uncertainty estimates for the International Geomagnetic Reference Field

The International Geomagnetic Reference Field (IGRF) is a multi-institute model of the Earth’s magnetic field, compactly described by sets of up to 195 spherical harmonic (Gauss) coefficients to degree and order 13, which allows the continuous evaluation of the field at any location and time on or above the surface. It is developed from satellite and ground-based magnetometer data and describes the large-scale variation of the magnetic field in space and time under quiet conditions. While much effort has been made on improving the forecast of the secular variation of the field over the 5-year intervals between release and renewal, less emphasis has been placed on understanding the spatial errors from a user point of view. In this study, we estimate the large-scale time-invariant spatial uncertainty of the IGRF based on the globally averaged misfit of the model to ground-based measurements at repeat stations and observatories between 1980 and 2021. As the ground measurements are reduced to quiet-time values, the external field is minimized for the purposes of this study. We find the 68.3% confidence interval is 87 nT in the North (X) component, 73 nT in the East (Y) component and 114 nT in vertical (Z) component. Due to the Laplacian distribution of the residuals, the standard deviations are larger at 144, 136 and 293 nT, respectively.


Introduction
The International Geomagnetic Reference Field (IGRF) is a multi-institute large-scale model of the Earth's magnetic field, compactly represented by sets of up to 195 spherical harmonic (Gauss) coefficients to degree and order 13, which allows the continuous evaluation of the field at any location in time on or above the surface. It is intended to reproduce the quiet-time average magnetic field generated primarily by the core. As the dominant contribution to the field at the surface is from the core, the IGRF typically describes around 95-98% of the magnitude of a given measurement, though this depends strongly on time and location (e.g. Finlay et al. 2015).
The first IGRF was introduced by the International Association of Geomagnetism and Aeronomy (IAGA) in 1968 in response to the need for a standard or reference spherical harmonic representation of the Earth's main field (e.g., Barraclough et al. 1978). The model is updated at 5-yearly intervals to better capture the non-linear secular variation of the main field. It is produced and released by IAGA Working Group V-MOD. Over the decades the model has improved in accuracy and precision as more and better magnetic field data become available (Macmillan and Maus 2005;Thébault et al. 2015). The model presently consists of Gauss coefficients at 5-year intervals between which the field is assumed to vary linearly with time and the latest model includes a prediction for 5  years into the future from the date of release. The current model is 13th generation (IGRF-13) which is valid until 2025.0 (Alken et al. 2021a). Its intention is to describe the internal large-scale (> 3000 km) variation of the magnetic field in space and time from the year 1900 to present day. The most recent coefficients in IGRF-13 were developed using up to 14 candidate models, derived from satellite and ground-based magnetometer data (Alken et al. 2021b). While much effort has been made recently on improving the forecast of the secular variation of the field over the 5-year intervals between release and renewal (see Fournier et al. (2021b) for a review), less emphasis has been placed on understanding the spatial uncertainties from a user point of view. The call for the IGRF-13 candidate models requested uncertainties in spectral terms for the Gauss coefficients which a few teams provided (e.g., Fournier et al. 2021a). Lowes (2000) made a formal estimate of the uncertainties in a root-mean-square sense based on the degree and order of each generation of the IGRF and effect of numerical truncation of the coefficients but this is a theoretical rather than measured uncertainty.

Graphical Abstract
It is also to be noted that the more rapidly time-varying external field is explicitly not considered in this study as the IGRF does include it, though it will be present in most measurements. As the effect of the ionosphere and magnetosphere is latitude-dependent (e.g., Beggan et al. 2021;Chulliat et al. 2020) quantifying the external field effects is not difficult. However, as these uncertainties are time-dependent (daily, seasonal, solar cycle) providing an 'average' error for these sources is difficult though possible, as long as it is clear to the user what the uncertainties represent (Beggan et al. 2018). We focus on the internal parts of the IGRF that are missing from the model. The unmodelled SV (e.g., jerks), tidal and induced fields are time-dependent too, but should average to a small value over the several decades considered. However, the crustal field dominates the spatial spectrum from around degrees 13 and above, so most of the missing (omitted) part of the IGRF is from the crustal field.
Major scientific users include space physicists (for coordinate systems), near-surface geophysicists (e.g., for leveling scalar surveys) and paleomagnetists (for orienting cores). Non-scientific users are often interested in directional information for navigation, particularly declination and inclination. In this brief study, we look at approximately estimating the large-scale spatial error of the IGRF to provide general uncertainty values to users. The errors are based on the residuals of the IGRF-13 model values to ground-based measurements (e.g., at repeat stations and observatories) between 1980 and 2021 and taken together represent a globally averaged uncertainty. We describe the data and methodology in the "Data and method" section. The results are presented in the "Results" section, and we interpret the meaning and appropriate use in the "Discussion and conclusion" section.

Data and method
We use land-based repeat stations and annual means of observatory data collected between 1980 and 2021 from the World Data Centre for Geomagnetism, Edinburgh.
The starting year is chosen as it was the first IGRF model to benefit from dedicated satellite data from the MAG-SAT mission (Cain et al. 1989). It is also the period when digital instrumentation began to be rolled out at magnetic observatories.
The observatory data will have been used in the development of the IGRF models so they are not a truly independent test. However, observatory data form around 32% (6800/21,500 points) of the initial raw dataset. Repeat station measurements, which are the majority of the data, are generally not used in the construction of main field models so can be considered independent. Observatory monthly mean or annual time-series are used in main field models as they have a fixed reliable cadence that can be easily ingested. Using repeat station data is eminently possible and main field models that describe the past field incorporate such data (e.g., Jackson et al. 2000). However, IGRF model candidates tend to mainly use observatory and satellite data where possible as the time coverage (and secular variation) is uniform.
Repeat station measurements are taken very carefully by skilled observers and are checked for quality and consistency before deposition. These measurements are reduced to quiet-time local night-time values using a local observatory as a reference (Newitt et al. 1996). This removes a large proportion of any external field effects on the measurement. Similarly, annual mean values are computed using quiet-time values to represent the best estimate of the internal field signal (core and crust) for the location (Jankowski and Sucksdorff 1996). These reductions allow as fair as possible a direct comparison to the IGRF-13 model values. Figure 1 shows the locations of the repeat station and observatories used. The upper plots (Panel a) show the spatial distribution for each 5-year period (e.g., 1980 includes measurements from 1980 to 1984.99 and so on). The distribution of data is unsurprisingly concentrated in the land areas particularly Europe, South Africa, Australia and the Americas though there are many data in the east Asia region including Indonesia. The repeat stations are red circles while observatories are blue. There are relatively fewer observatory data in the period from 1980 to 1990s but more repeat stations; by 2010 there are fewer repeat stations as satellite data reduce the need for such surveys.
There are large regions of the Pacific and Atlantic Oceans and polar regions where there are few to no data available. Over the ocean regions, observatories are often placed on volcanic islands which have large crustal anomaly magnitudes (Lesur et al. 2016); indeed, Scott Base Antarctica (SBA) has one of the largest vertical residuals at -3750 nT. Volcanic island contributions to the global uncertainty may produce a bias as they are not representative of the generally lower magnetization of continental sedimentary rocks. In marine regions typically only total field (F) data are available as few vector surveys exist and are usually confined to small regions (e.g., parts of the North Sea). Offshore aeromagnetic and marine datasets typically pre-date the 1980s, so are not directly useful for this particular study as the IGRF models are themselves poorer prior to 1980 in any case. The lower plot (b) shows the position in quasi-dipole coordinates of the measurement locations which vary slowly over time due to secular variation of the field (Emmert et al. 2010); we will return to this later.
The residuals between the IGRF-13 model (which is definitive between 1980 and 2015) and ground measurements should be representative of the average uncertainties expected on the land surface. We first compute the residuals in all seven components: declination (D), inclination (I), horizontal intensity (H), North component (X), East component (Y), vertical component (Z) and total field intensity (F). In each component, we identify and remove any residuals that are greater than 3 standard deviations from the mean. This step removes any very large outliers in the dataset, but will also remove some signal from sites located over very strong anomalies too. Depending on the component, between 16330 and 19310 residual values remain. The absolute values of the residuals are sorted from smallest to largest and the confidences levels at 68.3%, 95.4% and 99.7% are computed by finding the relevant percentile value. This accounts for the Laplacian distribution of residuals in magnetic data . We also recompute the mean and standard deviation of the outlier-free residuals.

Results
The residuals are plotted by component as histograms in Fig. 2 and spatially in Fig. 3. Table 1 provides the global statistics of the residuals in each component. As hoped, the mean of the residuals are close to zero for the inclination and declination. The mean of the H, X and Z residuals are around 18-21 nT away from zero. This may be related to the near-constant magnetospheric ring current signal in the ground measurements (and deliberately not captured or corrected for in the IGRF model), in which case these non-zero means also probably reflect the bias of measurements in the northern hemisphere. Table 2 shows how the 68.3% confidence interval (CI) varies between the eight 5-year epochs (e.g., 1980-1985) and the approximate number of data per epoch. As can be observed in general, the fewer data there are in any particular epoch the larger the CI values are. A similar trend occurs in the other CI (95.4% and 99.7%) and the standard deviations too. This reflects the smaller number of the sites surveyed.
The histograms in Fig. 2 show the residuals have a Laplacian distribution rather than Gaussian, with sharp peaks around zero and long tails. This distribution arises from the geophysical processes that create magnetic fields and the fact that the power spectrum of the crustal field at the surface tends to be red rather than white (e.g., Olsen et al. 2017). The Laplacian distribution is also reflected in the difference between the standard deviations of the component residuals and the ranked confidence intervals in Table 1; for example, the 68.3% confidence interval (the 1 σ-equivalent value) is 87 nT in the North (X) component, 73 nT in the East (Y) component and 114 nT in vertical (Z) component, while the standard deviations are larger at 144, 136 and 293 nT. In declination the standard deviation is 0.39 • while the 68.3% CI is 0.18 • . In all components the standard deviation is larger than the 68.3% CI.
The spatial distribution of the absolute value of the residuals in Fig. 3 are plotted as colored circles to indicate their difference from zero. The majority of the points globally are of a blue color indicating a difference of less than 0.15 • in angle or 50 nT in magnetic strength. There are no discernible spatial patterns that are obviously related to the core field magnitude or the auroral zone influence. This suggests the majority of the difference in the residuals is related to unmodelled local crustal fields. When plotted as positive or negative residuals (not shown), a similar conclusion can be drawn as there are no hemispherical or polar-only anomaly patterns visible.
Finally, Fig. 4 shows the residuals plotted against quasidipole (QD) latitude (between ± 80 • ) for the X, Y and Z components. The residuals (gray dots) are concentrated in the northern hemisphere particularly between 20 • to 60 • QD latitude. There is no obvious increase in the residuals with high latitude as might be expected if external fields were still present. The residuals were further analyzed by dividing them into bins of 2.5 • latitude and the sorting the absolute values in each bin into 68.3%, 95.4% and 99.7% confidence intervals. The colored bars indicate the 68.3% (green), 95.4% (blue) and 99.7% (red) intervals. As can be observed, again, there is no clear pattern of increase or decrease with latitude, suggesting the residuals are from crustal anomalies rather than external. A set of plots (not shown) for declination, inclination and total field show similar patterns.

Discussion and conclusion
The IGRF aims to provide a model of the main field which is accessible to a variety of users. It is meant to give a reasonable approximation, at, near and above the Earth's surface, to that part of the Earth's magnetic field which has its origin inside the surface. However, it clearly cannot capture all internal sources due to small-scale features of the crustal field and will also not capture the non-linear features of secular variation such as geomagnetic 'jerks' (Brown et al. 2016). This is the first IGRF-based analysis of the spatial error when compared directly to semi-independent ground measurements (though see Chulliat et al. (2015) Table 12 for the World Magnetic Model analysis). The measurements are themselves biased by being land-based and only in certain regions of the world. In addition, we have pruned the dataset to remove the very largest outliers which will also reduce the overall uncertainty values. It is nonetheless pleasing to note that the IGRF model still fits measured ground data very well, as shown in Table 1 and Figs. 3 and 4.
The histograms have means of approximately zero in the angular components (Fig. 2). The deviations of the means of the horizontal, north and vertical components from zero can be accounted for by the missing contribution of the steady ring current signal which is generally removed from the IGRF candidate models but not the repeat stations or observatory annual means.
The very largest values (< 2% of the initial dataset) tend to be gross outliers which are introduced for example from poor data processing, mistakes in digitization or file formatting, though there may be some true values rejected too. We also checked the effect of choosing to set the rejection level at the 4 σ rather than 3 σ in the initial dataset. This gives the largest change in the 99.7% confidence interval, but hardly affects the 63.3% or 95.4% values. Choosing to reject at the 4 σ level also affects the average and standard deviations as all data contribute to these values, but only very slightly. The mean of the Y changes from − 2 to − 7 nT and X and Z means change by − 2 nT (not shown). Figures 3 and 4 suggest there is not a strong latitudinal control on the residuals. This is primarily due to the quiet-time reduction to remove external fields, but may also be a function of the sparsity of the ground data available too. As the external field effects are not the first-order cause of the residual differences between measurements and models, most differences are of geological origin; for example over known large anomalies like Bangui in central Africa or Ross Island in Antarctica.
In terms of useful nominal global uncertainties associated with the IGRF-13 model between 1980 and 2020, the standard deviation of the residuals in North (X) component is 144 nT, 136 nT in the East (Y) component and 293 nT in vertical (Z) component. For declination  it is 0.39 • , in inclination it is 0.29 • , and for total field intensity it is 178 nT. More carefully, as the distribution of the residuals is Laplacian, the uncertainties of X, Y and Z can be estimated to be 87, 73 and 114 nT at 68.3% confidence interval. This analysis of the comparison between ground measurements and the IGRF-13 model provides a direct, evidence-based, estimate of the global average uncertainties attributable to the missing internal field sources. The results are a useful first-order approximation of the expected difference between a field measurement and the IGRF. However, as noted in the official IGRF 'health warning' (https:// www. ngdc. noaa. gov/ IAGA/ vmod/ igrfhw. html) written by Frank Lowes in 2010: 'If you measure the magnetic field at a point on the Earth's surface, do not expect to get the value predicted by the IGRF!' .