IZMIRAN candidate field model for IGRF-13

The International Geomagnetic Reference Field (IGRF) model is a combination of the several models developed by independent groups of scientists using different approaches for the selection of input data and methods for calculating harmonic coefficients. This approach allows for mutual comparison of individual models and for their combination to obtain the most reliable values of the harmonic coefficients. This letter provides a brief description of methods for building the IZMIRAN Earth’s main magnetic field model, submitted to the IAGA Working Group V-MOD for creating IGRF-13. Special efforts were made to obtain as uniform coverage of the entire Earth’s surface as possible with observations. The surface was divided into a grid of approximately equal cells. Then the data for geomagnetically quiet intervals were placed in the corresponding cells and a median filter was applied to select the data in each cell. Spherical harmonic coefficients up to degree 13 were calculated for the interval 2014-Jan to 2019-Aug with a time step of 10 days and extrapolated to 01.01.2010.


Introduction
During recent decades the most important sources of vector magnetic field data for modeling the Earth's main magnetic field were dedicated satellites designed for the magnetic surveys (MAGSAT, CHAMP, and Swarm missions). Now Swarm project (Friis-Christensen et al. 2006) provides the scientific community with a uniform data set for the entire IGRF 2015-2020 modeling period that allows us to derive candidate model based only on Swarm data.
The most common way to describe the spatial distribution of the Earth's main magnetic field is to fit the observed magnetic field with a set of spherical harmonics (Barton 1997). The general approach to deriving global models of the Earth's magnetic field is well-known for a long time (Chapman and Bartels 1940), and we will skip its detailed overview here. An extended list of references is given in Alken et al. (2020, in press).
However, there are several problems that should be taken into account when building a model and it can be done in various ways.
The method of calculating the spherical harmonic coefficients allows us to separate the sources of the magnetic field lying inside and outside the surface on which the magnetic field is measured by a satellite. The ionospheric currents are external sources of the Earth's total magnetic field. But they are located inside the orbits of satellite and are included in the spherical harmonic coefficients describing the internal part of the Earth's magnetic field.
During quiet geomagnetic conditions when currents induced in conducting Earth are small, the external sources of the magnetic field (such as magnetopause, ring, and field-aligned currents), should not theoretically affect the coefficients of the internal part, but they would be completely excluded if we had simultaneous measurements on the entire surface or they were constant during the whole period of data collection. In reality, to get as uniform coverage of the entire surface as possible, data for several days are required and external sources may change during this time and perfect separation is not possible.
Usually, to find spherical harmonic coefficients, the method of minimizing the difference between measured magnetic field values and calculated (modeled) values is used. If the surface around the Earth where magnetic Petrov and Bondar Earth, Planets and Space (2021) 73:46 field measurements were carried out is unevenly covered by measurement points, the method results in a better fit of the data in areas where the density of these points is greater, and it results in a worse fit in areas where there are fewer points. This is especially evident for data from the high-latitude and polar orbit satellites when there are significantly more measurements per unit of surface area at the high latitudes compared to the low latitudes. Usually, such latitudinal uneven distribution is taken into account by introducing weight coefficients proportional to the cosine of latitude into the minimization equation.

Data selection
To calculate the model, we used Swarm 1-s vector data from all three satellites for January 2014-August 2019. The data were checked for errors, and all outliers and erroneous data were deleted. Observatories and magnetic survey repeat stations data were not used. The main reason for the rejection of these data was uneven spatial space distribution. A common method for reducing the influence of the ionospheric and magnetospheric currents is a special selection of time intervals when these currents are minimal. Usually, the selection is based on various indicators of the geomagnetic and solar activity, local time, and the sun's illumination of the ionosphere. Each research group used its own set of indexes which can coincide partially, see for example Olsen et al. (2000), Finlay et al. (2015).
In our model we used the next set of indices.
Planetary Kp index: we selected only 3-h intervals when Kp was less than 1.5 and Kp for the previous 3-h interval was less than 2.5.
ASY/SYM index: then, from the selected by Kp intervals, intervals when SYMH variation during this 3-h interval was more than 20 nT were rejected. AE-index: for years 2018-2019 when final AE was not available, real-time (QuickLook) AE plots were digitized. From each selected data only data when AE was less than 100 nT and the absolute value of SYMH was less than 20 nT were used.
We suppose that the used criteria allowed us to select intervals with the quiet and stable magnetosphere and ionosphere when the effects of the ionospheric and magnetospheric currents were minimal. Such a set of selection criteria is a compromise between the acceptable level of the external magnetic field variations and the spatial data coverage.

Deriving of the model
To obtain a homogeneous distribution of points in space, the Earth's surface was divided into 6 degrees (~ 660 km) latitudinal belts, and each belt was divided into about 660 km width cells by the longitude, such that the number of the cells in every belt was an integer. The lower latitudinal belt had 60 cells, the high-latitude one had only three, but the area of each cell was approximately the same. Then all selected data were placed into the appropriate cell in accordance with the position of the satellite. Usually, one whole day (without the choice of quiet geomagnetic conditions) of Swarm data fills 70% of cells, 2 days-95%. Taking into account the selection of data only for quiet intervals, we used the 10-day interval for calculating the individual set of the spherical harmonic coefficients. To get the best temporal resolution, we derived models for every fifth day, using 5 days before and 5 days after the date. After filling the cells the median filter was used to select only one Earth's magnetic field value for each component for each cell and each satellite, so in most cases, nine data points in each cell were used for the spherical harmonics calculation. Then we used only 10-day intervals for which more than 90% of all cells were filled to obtain homogeneous spatial data distribution. Extended paper to describe the model algorithm in more detail is in preparation. The coefficients of spherical harmonic were calculated to a maximum spherical harmonic degree and order of 13 and plotted as a function of date. Figure 1 shows an example of such graphs. Most of the spikes on the plots result from intervals with the worst filling of cells.
To check the effects of the daylight ionosphere, the same calculations were done with an additional condition-the Sun elevation angle at the ionosphere height should be greater than 95°, so only nighttime data were used. Smoothed values of the harmonic coefficients do Fig. 1 Variations of spherical harmonic coefficients g10, g11, h11 over the period 2014-2020. Red lines show parabolic fit to data not change significantly, but due to a decrease in the number of data points in each cell and the number of cells with data, scattering of harmonic coefficients values significantly increased. So, the IZMIRAN candidate model was derived taking into account all nighttime and daytime data.
Parabolic fits to the data show that for g10 and h11 fit is virtually linear, but for g11 non-linearity is apparent and fit should be a parabolic. Many other harmonic coefficients also show visible non-linearity, sometimes very strong (see g22, Fig. 2). Some spherical harmonic coefficients even could not be described by a parabolic curve and require cubic fit (see h33 Fig. 2).

IGRF candidate model for 2020
The method described above allowed us to construct variations of the spherical harmonic coefficients until the end of August 2019, and it became necessary to extrapolate them to 01.01.2020 to submit for IGRF-13. If the harmonic coefficients time variation is linear, this extrapolation does not cause problem. However, if there is a noticeable non-linearity, different results can be obtained for such coefficients, depending on the used extrapolation model and the interval selected for extrapolation. For example, for the spherical harmonic h33 shown in Fig. 2 various extrapolation methods give extrapolated value from -546 to -542 nT. To get the IZMIRAN model for 01.01.2020 we used data for the last 2 years. The extrapolated spherical harmonic coefficients were calculated as parabolic extrapolated data for the interval of 01.09.2017-31.08.2019 on 01.01.2020. For each spherical harmonic a three-pass parabolic approximation was used. After the first fit, points with a deviation of more than two standard deviations from the fit were discarded, and the next fit was made based on the remaining data. This procedure was repeated two times. The new (final) IZMIRAN model for IGRF-2015 was calculated in the same way from fits for the interval of 01.01. 2014-31.12.2015 and fitted values on 01.01.2015 were presented. Harmonics coefficient errors were calculated as mean square deviation from fit divided by square root of used point numbers.

Secular variation (SV) model
It is much more difficult to estimate the secular variation in the case of a noticeable non-linear harmonic coefficients time variation. For some harmonics (for example h33, h22, g32) approximations for 2015-2019 and 2018-2019 were significantly different (SV significantly change during 5 years interval) and extrapolation of 2-year data set 3 years forward is very unreliable. So, SV_2022.5 was calculated on the base of 5 years (08.2014-08.2019) parabolic approximation. This choice leads to a larger error at the beginning of the period 2020-2025 compared to the estimation for the last 2 years, but since there are no assumptions about a possible change of the secular variation in the future, we believe that such conservative estimation is reasonable (Fig. 3).

Conclusions
Recent Swarm satellite magnetic mission served as the basis for a high-precision model of the main magnetic field for 2020. In this letter, we briefly describe the Earth's main magnetic field and SV models developed in IZMIRAN based only on Swarm 2014-2020 data. The results were presented to IAGA Working Group V to derive the IGRF-13 model. Temporal variations of spherical harmonic coefficients show that for some harmonics, even a quadratic approximation cannot describe the variations correctly and one of the reasons for the differences in models may be the approximation method itself. It is especially important for SV values variation.
Abbreviations IAGA : International Association of Geomagnetism and Aeronomy; IGRF: The International Geomagnetic Reference Field; SV: Secular variation.