Introduction

The International Geomagnetic Reference Field (IGRF) is a set of spherical harmonic coefficients which can be input into a mathematical model in order to describe the large-scale, time-varying portion of Earth’s internal magnetic field between epochs 1900 A.D. and the present. The IGRF is produced and maintained by an international task force of scientists under the auspices of the International Association of Geomagnetism and Aeronomy (IAGA) Working Group V-MOD. This thirteenth generation IGRF has been derived from observations recorded by satellites, ground observatories, and magnetic surveys (see Appendix 1 for a list of World Data System data centers and services). IGRF is routinely used by the scientific community to study Earth’s core field, space weather, electromagnetic induction, and local magnetic anomalies in the lithosphere. It is also widely used in satellite attitude determination and control systems and other applications requiring orientation information.

Earth’s core field changes continuously and unpredictably on timescales ranging from months to millions of years. In order to account for temporal changes on timescales of a few years, the IGRF is regularly revised, typically every 5 years. Table 1 summarizes the current and past generations of IGRF. Each generation is composed of a set of model coefficients representing the internal time-varying geomagnetic field, which are provided in 5-year intervals. The years for which coefficients are provided are called model epochs. The coefficients of a certain epoch represent a snapshot of the geomagnetic field at that time, and can be labeled either as a Definitive Geomagnetic Reference Model (DGRF) or as an IGRF. DGRF models are so labeled because they have been built from the best available data sources of that time period and therefore are unlikely to be improved in future IGRF revisions. Models labeled as IGRF are non-definitive, and will likely be revised in the future as more data are collected. DGRF models have been built only starting in 1945. Details of the history of IGRF can be found in Barton (1997) and Macmillan and Finlay (2011). Past generations of IGRF models are archived at https://www.ngdc.noaa.gov/IAGA/vmod/igrf_old_models.html. Since later IGRFs can revise model parameters for past epochs, it is important to record which generation of IGRF was used to process a particular dataset, so that the original data can be recovered and reprocessed with the latest generation of IGRF if needed.

Table 1 Summary of IGRF generations, validity periods, release years, and references
Full size table

In this paper, we focus on the thirteenth generation of IGRF, known hereafter as IGRF-13. IGRF-13 provides a DGRF model for epoch 2015.0, an IGRF model for epoch 2020.0, and a predictive IGRF secular variation model for the 5-year time interval 2020.0 to 2025.0. For epochs 1900.0 to 2010.0, the IGRF-13 model coefficients are unchanged from IGRF-12. IGRF-13 was finalized in December 2019 by a task force of IAGA Working Group V-MOD. In the following sections, we will describe the IGRF model, provide the final set of IGRF-13 coefficients, and briefly discuss large-scale features of the geomagnetic field at Earth’s surface as revealed by the updated model.

Mathematical formulation of the IGRF model

The IGRF describes the main geomagnetic field \({\mathbf {B}}(r,\theta ,\phi ,t)\) which is produced by internal sources primarily inside Earth’s core. The IGRF is valid on and above Earth’s surface, where the main geomagnetic field can be described as the gradient of a scalar potential, \({\mathbf {B}} = -\nabla V\), and the potential function \(V(r,\theta ,\phi ,t)\) is represented as a finite series expansion in terms of spherical harmonic coefficients, \(g_n^m, h_n^m\), also known as the Gauss coefficients:

$$\begin{aligned} V(r,\theta ,\phi ,t) = a \sum _{n=1}^N \sum _{m=0}^n \left( \frac{a}{r} \right)^{n+1} \\\quad \left[ {{g_n}^m(t)} \cos {m\phi } + {{h_n}^m(t)} \sin {m\phi } \right] {P_n}^m (\cos {\theta }) \end{aligned}.$$
(1)

Here, \(r,\theta ,\phi\) refer to coordinates in a geocentric spherical coordinate system, with r being radial distance from the center of the Earth, and \(\theta ,\phi\) representing geocentric co-latitude and longitude, respectively. A reference radius \(a = 6371.2\) km is chosen to approximate the mean Earth radius. The \(P_n^m(\cos {\theta })\) are Schmidt semi-normalized associated Legendre functions of degree n and order m (Winch et al. 2005). The parameter N specifies the maximum spherical harmonic degree of expansion, and was chosen to be 10 up to and including epoch 1995, after which it increases to 13 to account for the smaller scale internal signals which can be captured by high-resolution satellite missions such as Ørsted, CHAMP and Swarm. The Gauss coefficients \(g_n^m(t), h_n^m(t)\) change in time and are provided in units of nanoTesla (nT) in IGRF-13 at 5-year epoch intervals. The time dependence of these parameters is modeled as piecewise linear, and is given by

$$\begin{aligned} g_n^m(t)&= g_n^m(T_t) + (t - T_t) {\dot{g}}_n^m(T_t), \end{aligned}$$
(2)
$$\begin{aligned} h_n^m(t)&= h_n^m(T_t) + (t - T_t) {\dot{h}}_n^m(T_t), \end{aligned}$$
(3)

where \(g_n^m(T_t), h_n^m(T_t)\) are the Gauss coefficients at epoch \(T_t\), which immediately precedes time t. The model epochs in IGRF-13 are provided in exact multiples of 5 years starting in 1900 and ending in 2020 (see Table 2), so that \(T_t \le t < T_t + 5\). For \(T_t < 2020\), the parameters \({\dot{g}}_n^m(T_t), {\dot{h}}_n^m(T_t)\) represent the linear approximation to the change in the Gauss coefficients over the 5-year interval spanning \([T_t,T_t+5]\). They may be computed in units of nanoTesla per year (nT/year) as

$$\begin{aligned} {\dot{g}}_n^m(T_t)&= \frac{1}{5} \left( g_n^m(T_t + 5) - g_n^m(T_t) \right) , \end{aligned}$$
(4)
$$\begin{aligned} {\dot{h}}_n^m(T_t)&= \frac{1}{5} \left( h_n^m(T_t + 5) - h_n^m(T_t) \right) . \end{aligned}$$
(5)

The main field coefficients are not yet known for \(T_t = 2025\), and so for the final 5 years of model validity (2020 to 2025 for IGRF-13), the coefficients \({\dot{g}}_n^m(2020), {\dot{h}}_n^m(2020)\) are explicitly provided (see last column of Table 2) in units of nT/year. Details on the individual candidate secular variation forecasts and the procedure used to combine them into a final set of \({\dot{g}}_n^m(2020),{\dot{h}}_n^m(2020)\) may be found in Alken et al. (2020b) and references therein.

Table 2 13th generation International Geomagnetic Reference Field
Full size table

The 13th generation IGRF

In August 2017, during an IAGA V-MOD Working Group meeting held in Cape Town, South Africa, a task force of volunteer geomagnetic modelers was assembled to oversee the call for IGRF-13 candidate models and their evaluation. In March 2019, the task force issued an international call for three candidates:

  • A DGRF main field model for the epoch 2015.0

  • An IGRF main field model for the epoch 2020.0

  • An IGRF linear secular variation model for the time period 2020.0 to 2025.0.

Fifteen teams representing over 30 international institutes responded to the call. The number of teams and institutions who participated in IGRF-13 exceeded that of any previous generation. The task force received 11 DGRF main field candidates for epoch 2015.0, 12 IGRF main field candidates for 2020.0, and 14 IGRF secular variation candidates for 2020.0-2025.0. Following recent IGRF conventions, the main field candidates for IGRF-13 describe the spatial variation of the field to a maximum spherical harmonic degree and order of 13, while the secular variation candidates extend to a maximum degree and order of 8. Each of the 15 teams was managed by a team leader from the lead institution, and many teams also included personnel from supporting institutions. The 15 lead institutions for IGRF-13, including references to their candidate model papers, are: (1) British Geological Survey (UK) (Brown et al. 2020); (2) Institute of Crustal Dynamics, China Earthquake Administration (China) (Yang et al. 2020); (3) Universidad Complutense de Madrid (Spain) (Pavón-Carrasco et al. 2020); (4) University of Colorado Boulder (USA) (Alken et al. 2020a); (5) Technical University of Denmark (Denmark) (Finlay et al. 2020); (6) GFZ German Research Centre for Geosciences (Germany) (Rother et al. 2020); (7) Institut de physique du globe de Paris (France) (Fournier et al. 2020; Vigneron et al. 2020; Ropp et al. 2020); (8) Institut des Sciences de la Terre (France) (Huder et al. 2020); (9) Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation (Russia) (Petrov and Bondar 2020); (10) Kyoto University (Japan) (Minami et al. 2020); (11) University of Leeds (UK) (Metman et al. 2020); (12) Max Planck Institute for Solar System Research (Germany) (Sanchez et al. 2020); (13) NASA Goddard Space Flight Center (USA) (Sabaka et al. 2020; Tangborn et al. 2020); (14) University of Potsdam (Germany) (Baerenzung et al. 2020), and (15) Université de Strasbourg (France) (Wardinski et al. 2020). Some of the lead institutes listed above also acted as supporting institutions to other teams. The supporting institutes which are not listed above include: Geoscience Institute (Spain), Hebei GEO University (China), Institute of Earthquake Forecasting, China Earthquake Administration (China), Institute of Geophysics, China Earthquake Administration (China), Kyushu University (Japan), Nagoya University (Japan), National Space Science Center, Chinese Academy of Sciences (China), Observatori de l’Ebre (Spain), Observatorio Geofísico de Toledo (Spain), Real Observatorio Geofísico de la Armada (Spain), Space Research Institute of the Austrian Academy of Sciences (Austria), The Institute of Statistical Mathematics (Japan), Tokyo Institute of Technology (Japan), Université de Nantes (France), and University of Tokyo (Japan).

Data recorded by the Swarm satellite mission (Friis-Christensen et al. 2006) and the ground observatory network (see Table 3) played a crucial role in the development of many of the IGRF-13 candidate models. Data from the Ørsted (Olsen et al. 2000), CHAMP (Reigber et al. 2002), SAC-C (Colomb et al. 2004), Cryosat-2, and CSES (Shen et al. 2018) missions were also used by some of the teams. The IGRF-13 task force voted to calculate the final main field models for epochs 2015.0 and 2020.0 as the medians of the Gauss coefficients of all the candidate models. The task force voted to use a robust Huber weighting in space to determine the final secular variation model for 2020.0 to 2025.0. Further details of the candidate models, the evaluation process, and the final model determination are provided in Alken et al. (2020b).

Table 3 Magnetic observatories contributing data used in the construction of IGRF-13
Full size table

IGRF-13 model coefficients and maps

Table 2 lists the IGRF-13 spherical harmonic Gauss coefficients, which can be used with Eq. (1) to determine the geomagnetic potential (and vector geomagnetic field) anywhere on or above Earth’s surface. This table serves as a published record of IGRF-13, which should allow users to ensure they use the correct model coefficients for a particular epoch compared with previous generations. The main field coefficients are given in units of nT, and the predictive secular variation coefficients (last column) are given in units of nT/year. These coefficients are available in digital form from https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html along with software to compute magnetic field components at different times and spatial locations, in both geocentric and geodetic coordinate systems.

Figure 1 shows global maps of the IGRF-13 declination (D), inclination (I), and total field magnitude (F) on Earth’s surface at 2020 in Miller cylindrical projection. Taken together, these three quantities fully describe the vector magnetic field at Earth’s surface. The green contour lines represent zero. For the declination component (top panel), these are the agonic lines on which a magnetic compass needle would point to true geographic north. For the inclination map (middle panel), the green contour line of zero inclination shows the magnetic dip equator, which approximately aligns with the geographic equator except for a large, well-known southward deviation over South America. The F map (bottom panel) shows that the largest field intensities occur in Siberia in the northern hemisphere and in the Southern Ocean between Australia and Antarctic in the southern hemisphere. We also see a region of significantly weaker field (compared to an idealized dipole), centered over South America, which is known as the South Atlantic Anomaly. In this region, the inner Van Allen radiation belt comes closest to Earth’s surface, which has important consequences for satellite instrumentation and human safety in low Earth orbit. Interestingly, a new second minimum is becoming more pronounced over the southern Atlantic. This feature is described in more detail in Rother et al. (2020) and Finlay et al. (2020) and was earlier reported by Terra-Nova et al. (2019).

Fig. 1
figure 1

Maps of declination (top), inclination (middle) and total field (bottom) at the WGS84 ellipsoid surface for epoch 2020. The zero contour is shown in green, positive contours in red, and negative contours in blue. White asterisks indicate locations of the magnetic dip poles. Projection is Miller cylindrical

Full size image

Figure 2 shows the predicted average change of the D, I, and F components on Earth’s surface during the 2020 to 2025 interval from IGRF-13. At low and middle latitudes, the map of dD/dt (top panel) predicts the largest declination changes in the South Atlantic Anomaly region and also in the polar regions, with northern polar declination changing more than in the southern polar region. The dI/dt map (middle panel) predicts the largest changes over Brazil, where the magnetic dip equator has moved relatively rapidly over the past few decades. The features seen in dF/dt (bottom panel) near South America predict a deepening and westward movement of the South Atlantic Anomaly, continuing a trend observed over the past century (Finlay et al. 2010a, Fig. 3).

Fig. 2
figure 2

Maps of predicted annual secular variation in declination (top), inclination (middle) and total field (bottom) at the WGS84 ellipsoid surface averaged over 2020 to 2025. The zero contour is shown in green, positive contours in red, and negative contours in blue. White asterisks indicate locations of the magnetic dip poles. Projection is Miller cylindrical

Full size image
Fig. 3
figure 3

Motion of the magnetic dip pole (red) and geomagnetic pole (blue) since 1900 from IGRF-13 in the northern hemisphere (left) and southern hemisphere (right). The scale provides an indication of distance on the WGS84 ellipsoid that is correct along lines of constant longitude and also along the middle lines of latitude shown. Note the left and right panels use different longitude ranges. The maps use stereographic projection. International and provincial boundaries are drawn in the left panel

Full size image

Figure 3 presents the positions of the geomagnetic poles and dip poles as given by IGRF-13 for 1900 to 2020, and the predicted positions in 2025. The geomagnetic poles are calculated from the three dipole (\(n=1\)) Gauss coefficients and correspond to where the magnetic dipole axis intersects a sphere of mean Earth radius 6371.2 km. These poles are antipodal and are also known as centered dipole poles (Laundal and Richmond 2017, Eq. 14). The geomagnetic poles can be used to specify the relative orientation of Earth’s magnetic field with respect to the Sun, and they are often used in magnetospheric studies for this purpose. The magnetic dip poles are defined as the locations where the main magnetic field as a whole is normal to Earth’s surface, represented by the WGS84 reference ellipsoid. Equivalently, they can be defined as the locations where the magnetic field component tangent to the ellipsoid vanishes. Here, we use the full set of IGRF-13 coefficients to spherical harmonic degree N. Magnetic dip poles provide a key reference for local orientation when navigating on or close to Earth’s surface at high-latitudes. For a perfect dipole field, the geomagnetic and dip poles would nearly coincide, but not exactly since the geomagnetic poles are defined with respect to a sphere of mean Earth radius, while the dip poles are defined with respect to the WGS84 ellipsoid. However, as can be seen in the figure, there are significant differences between the two due to the non-dipolar structure of Earth’s magnetic field. The geomagnetic and dip pole locations are provided in Table 4.

Table 4 Magnetic pole position since 1900 determined from IGRF-13 in units of degrees. Latitudes are provided in the WGS84 geodetic system
Full size table

Figure 4 shows the speed of the two magnetic dip poles. The north magnetic dip pole experienced a strong acceleration from about 1960 to 2000, but has seen a modest deceleration over the past 20 years, peaking at 55.8 km/year in 2002.5 and slowing slightly to 50.6 km/year in 2017.5. IGRF-13 forecasts a speed of 39.8 km/year in 2022.5, however we caution that past IGRF forecasts contained significant errors (Finlay et al. 2010b). As an example, IGRF-12 predicted a north dip pole speed of 42.6 km/year for 2017.5 (Thébault et al. 2015), compared with the IGRF-13 value of 50.6 km/year. Uncertainties present in IGRF models are further discussed by Lowes (2000).

Fig. 4
figure 4

Average speed of the magnetic dip poles over each 5-year epoch, plotted at the midpoint between epochs (i.e., the speed over 2015–2020 is shown at 2017.5). The value for 2020–2025 is a forecast

Full size image

At Earth’s surface in 2020, the contribution from the dipole terms \(g_1^0,g_1^1,h_1^1\) accounts for over 93% of the power in the main geomagnetic field. It is therefore instructive to monitor the temporal change in the dipole moment, which is defined as:

$$\begin{aligned} M(t) = \frac{4 \pi }{\mu _0} a^3 \sqrt{g_1^0(t)^2 + g_1^1(t)^2 + h_1^1(t)^2}. \end{aligned}$$
(6)

Figure 5 presents the change in the dipole moment of the geomagnetic field since 1900 as predicted by IGRF-13 (red). We see a clear downward trend in the dipole strength since the beginning of the last century, which is continued in 2020 and also in the forecast for 2025. This steady downward trend extends back at least as far as 1600 (Merrill et al. 1996; Constable and Korte 2015), although archeomagnetic and paleomagnetic records have revealed much lower dipole moments thousands of years in the past (Panovska et al. 2019). Due to sparsity of data, archeomagnetic and paleomagnetic studies often estimate the dipole strength along the rotation axis, ignoring the off-axis terms \(g_1^1,h_1^1\). This so-called axial dipole moment is defined as \(M_A(t) = 4\pi a^3 |g_1^0(t)| / \mu _0\) and is shown in blue in the figure.

Fig. 5
figure 5

Dipole and axial dipole moment time series derived from IGRF-13. The value for 2025 is a forecast

Full size image

IGRF-13 online data products

Further general information about IGRF: https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html

Coefficients of IGRF-13 in ASCII format: https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf13coeffs.txt

Fortran software to compute magnetic field components from coefficients: https://www.ngdc.noaa.gov/IAGA/vmod/igrf13.f

Linux C software to compute magnetic field components from coefficients: https://www.ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz

Windows C software to compute magnetic field components from coefficients: https://www.ngdc.noaa.gov/IAGA/vmod/geomag70_windows.zip

Python software to compute magnetic field components from coefficients: https://www.ngdc.noaa.gov/IAGA/vmod/pyIGRF.zip

Online calculation of magnetic field components for IGRF-13: https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml and http://geomag.bgs.ac.uk/data_service/models_compass/igrf_calc.html and http://wdc.kugi.kyoto-u.ac.jp/igrf/point/index.html

Archive of previous generations of IGRF: https://www.ngdc.noaa.gov/IAGA/vmod/igrf_old_models.html

Candidate models contributing to IGRF-13 and task force evaluation reports: https://www.ngdc.noaa.gov/IAGA/vmod/IGRF13/