International Geomagnetic Reference Field: the thirteenth generation

In December 2019, the International Association of Geomagnetism and Aeronomy (IAGA) Division V Working Group (V-MOD) adopted the thirteenth generation of the International Geomagnetic Reference Field (IGRF). This IGRF updates the previous generation with a definitive main field model for epoch 2015.0, a main field model for epoch 2020.0, and a predictive linear secular variation for 2020.0 to 2025.0. This letter provides the equations defining the IGRF, the spherical harmonic coefficients for this thirteenth generation model, maps of magnetic declination, inclination and total field intensity for the epoch 2020.0, and maps of their predicted rate of change for the 2020.0 to 2025.0 time period.


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_model s.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.
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 B(r, θ, φ, 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, B = −∇V , and the potential function V (r, θ, φ, t) is represented as a finite series expansion in terms of spherical harmonic coefficients, g m n , h m n , also known as the Gauss coefficients: Here, r, θ , φ refer to coordinates in a geocentric spherical coordinate system, with r being radial distance from the center of the Earth, and θ , φ representing geocentric co-latitude and longitude, respectively. A reference radius a = 6371.2 km is chosen to approximate the mean Earth radius. The P m n (cos θ) 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 m n (t), h m n (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 where g m n (T t ), h m n (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 ≤ t < T t + 5 . For T t < 2020 , the parameters (1) V (r, θ, φ, t) = a N n=1 n m=0 a r n+1 g n m (t) cos mφ + h n m (t) sin mφ P n m (cos θ) .

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: 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) -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 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  and was earlier reported by Terra-Nova et al. (2019). 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).

Table 3 Magnetic observatories contributing data used in the construction of IGRF-13
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 highlatitudes. 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   Table 4. 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).
At Earth's surface in 2020, the contribution from the dipole terms g 0 1 , 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:   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π a 3 |g 0 1 (t)|/µ 0 and is shown in blue in the figure.