Bulgarian Geomagnetic Reference Field ( BulGRF ) for 2015 . 0 and secular variation prediction model up to 2020

The Bulgarian Geomagnetic Reference Field (BulGRF) for 2015.0 epoch and its secular variation model prediction up to 2020.0 is produced and presented in this paper. The main field model is based on the well-known polynomial approximation in latitude and longitude of the geomagnetic field elements. The challenge in our modelling strategy was to update the absolute field geomagnetic data from 1980.0 up to 2015.0 using secular measurements unevenly distributed in time and space. As a result, our model gives a set of six coefficients for the horizontal H , vertical Z, total field F , and declinationD elements of the geomagnetic field. The extrapolation of BulGRF to 2020 is based on an autoregressive forecasting of the Panagyurishte observatory annual means. Comparison of the field values predicted by the model with Panagyurishte (PAG) observatory annual mean data and two vector field measurements performed in 2015 shows a close match with IGRF-12 values and some difference with the real (measured) values, which is probably due to the influence of crustal sources. BulGRF proves to be a reliable alternative to the global geomagnetic field models which together with its simplicity makes it a useful tool for reducing magnetic surveys to a common epoch carried out over the Bulgarian territory up to 2020.


Introduction
In the era of satellite observations, a variety of global models have been developed to represent the Earth's magnetic field. Because of the high altitude, they have spatial resolutions of hundreds of kilometres and focus mostly on the core field and large-scale lithospheric signals (e.g. CHAOS model series: Olsen et al., 2006;Finlay et al., 2016;GRIMM: Lesur et al., 2008). Study of intermediate-wavelength anomalies requires a combination of satellite, and/or airborne and ground measurements, as, for example, applied on the global scale for the World Digital Magnetic Anomaly Map (WDMAM) (Korhonen et al., 2007;Dyment et al., 2015) and of course the most recent comprehensive model (CM5) of the geomagnetic field  which co-estimate the major field sources using many different datasets. It uses joint inversion to describe field contributions from core, lithospheric, and external fields, along with associated Earth-induced signals. Regional models, on the other hand, are produced from data measured inside the area of interest thus being able to represent more accurately the field behaviour over that area. Unlike the global models, they are capable of modelling wavelengths in the kilometric range, and thus providing a better spatial resolution.
There are different techniques suitable for modelling the field in regions covering from a few squared geographical degrees to continental scales. Models for some countries are developed by means of revised spherical cap harmonic analysis (R-SCHA) (Thébault et al., 2006;Korte and Thébault, 2007;Qamili et al., 2010). Others successfully use simpler procedures such as second-degree polynomial fitting (e.g IT- GRF: De Santis et al., 2003;Dominici et al., 2007;Žagar and Radovan, 2012) to represent the space and temporal behaviour of the field better than the global models. Kovács et al. (2015) compared the advantages and disadvantages of polynomial and adjusted spherical harmonic analysis (ASHA) (De Santis, 1992) models over the territories of Croatia and Hungary. They concluded that because of its physical justification the obtained ASHA model is superior to the polynomial models, although the latter resulted in lower residuals. In the case of Bulgaria we chose to apply the procedure of second-degree polynomial fitting for three reasons: (1) we aim at developing a simple reference model up to 2020 which will be used when carrying out magnetic surveys, (2) our repeat stations are too limited in space and too scattered in time in order to obtain a good fit for a SCHA model, and (3) the territory of Bulgaria is small enough to neglect the curvature of the Earth in the investigated area. In this paper we present an up-to-date version of the Bulgarian geomagnetic reference model that is valid for epoch 2015 with a secular variation model until 2020. The model, characterized by a set of six coefficients for each element of the field was produced by least-squares fitting of second-degree polynomial of geographical coordinates to 473 points over the territory of Bulgaria.
Comparison of the modelled values with real measurements is a method for model validation. Of course, a full match cannot be expected because some contribution also comes from the magnetized rocks of the Earth's crust which turns the modelling into a tool for crustal anomaly investigation. We use our model also to evaluate the magnitude of the well-delineated magnetic anomaly in Panagyurishte (PAG) observatory region which was the main obstacle when in 1935-1936 experts discussed the location for building a magnetic observatory in Bulgaria. The anomaly is caused by two andesitic veins located inside the crust south of Panagyurishte which contain high levels of ferromagnetic minerals. Nevertheless, after consultation with geologists who discussed the possible viscous nature of the rock magnetization and/or shifting of the masses relative to the felsic granite above which the observatory will be built it was decided that this anomaly could be accepted as a constant.

Data
The repeat station network of Bulgaria was established in 1934 (Kostov and Nozharov, 1987). The eight points se-lected were then supplemented by seven more points in 1964 ( Fig. 1). Up to 1980 repeat stations were measured every 3 years and then, because of the small secular variation, every 5 years. The last absolute geomagnetic survey was performed in the period 1978-1980. The geomagnetic elements D, H , and F were measured at 473 points. Magnetic theodolite "Schulze-545", three quartz horizontal magnetometer (QHM) and two proton magnetometers PMP-2A were used. Reduction to the 1980.0 epoch of the observations was made according to the Panagyurishte (PAG) observatory. The last complete geomagnetic measurement on the 15 secular stations in Bulgaria was carried out in 1990. Because of lack of finance, further measurements were covering only part of the territory (usually one-quarter of the area per year). The most recent campaign lasted from 2007 to 2012 (Cholakov and Mihovski, 2010) and covered the territory of Bulgaria consecutively in four parts starting from northwest and going clockwise to southwest.
The big challenge in our modelling strategy was to update the absolute field geomagnetic data from 1980.0 to 2015.0 using measurements that are sparse and unevenly distributed in time and space. The main problem was the data extrapolation from the last measured value to the aimed 2015.0 epoch. To solve this problem we investigated the correlation between the secular variation trends of every single component in each repeat station and the secular variation trends in Panagyurishte (PAG) observatory for the same period. Thus, we obtained coefficients allowing us to calculate the secular variations of the geomagnetic field elements in repeat stations for the extended time intervals without real measurements, using data from Panagyurishte observatory and the equation: where E ϕ,λ (t 2 ) is the predicted element's secular variation in the repeat station (in nT yr −1 ) for period t 2 (in years) without field measurements, 1 the E ϕ,λ (t 1 ) is the known element's secular variation in the repeat station (in nT yr −1 ) for period t 1 (in years), and ϕ and λ are geographical coordinates (latitude and longitude). α = E PAG (t 2 ) / E PAG (t 1 ) is the ratio between the secular variations of the element for the two periods as recorded in PAG observatory. Using the above relation, secular variation values for every year (between 1980 and 2015) for each component for all repeat stations were calculated. Example plots of predicted D, F , and H secular variation values for four repeat stations (Slavotin, Bekleme, Popovo, and Michurin) can be seen in Fig. 2, together with PAG secular variations for the chosen period.
As a result of all calculations, for each geomagnetic element (H , D, Z and F ) a catalogue was obtained containing data for the Bulgarian repeat stations reduced to epochs up to 2015.0. Using those data and annual mean values from neighbouring observatories (SUR, GCK, IZN, PEG) a series of isoporic maps of the accumulated change of the main field elements for a 35-year period were prepared (Fig. 3). The respective values for the 473 field points measured in 1980 (Kostov et al., 1991) were extracted from the maps. Further, those data were used for reduction of the last absolute geomagnetic survey of Bulgarian territory performed in 1980 to the 2015.0 epoch. Accordingly, four input databases (H , D, Z, and F ) were generated and used as regressand in the least-squares regression technique.

Model generation
Although the surface polynomials were the first analytical method used to produce regional models (Haines, 1990)   In the present research we fit a second degree polynomial of the form (Buchvarov and Kostov, 1981) E ( ϕ, λ) = a 1 + a 2 ϕ + a 3 λ + a 4 ϕ 2 + a 5 ϕ λ + a 6 λ 2 , where independent variables (regressors) are 1, ϕ, λ, ϕ 2 , ϕ λ, λ 2 , and dependent variables (regressand) are geomagnetic field elements D, H , Z, F observed at 473 points over the Bulgarian territory. Latitude and longitude, expressed in decimal degrees, are referenced to a central point (25.0 • E, 42.5 • N) of the concerned area.
We use the well-known technique of least squares (Wilks, 1967) according to which if we have a function where E is the observation matrix of order N (i = 1, 2, . . ., N is the number of observations), A is the coefficient matrix with elements a j (j = 1, 2, . . ., 6) and X is the rectangular matrix of normal equations with dimension 6 × N. Then, the empirically determined statistical estimate of the coefficients having minimum variance is (Wilks, 1967) where (X T X) = C is the covariance matrix.
Other important parameters of our analysis (see for example Buchvarov, 1977) are as follows: the statistical estimation variance: the limits of the confidence intervals of the coefficients: where t N−k,γ is Student's coefficient corresponding to confidence interval γ with N − 6 of freedom; the confidence intervals of the calculated variables: whereσ E = X T C −1 Xσ 2 ; the correlation coefficient (for assessing the significance of the regression): During the calculations the following procedure was applied to remove anomalous points or possible errors in measurements (Buchvarov and Kostov, 1981). Using data from all observational points, a statistical estimate of the regression coefficientsÂ was determined for the geomagnetic field elements D, H , Z, and F . Then, the normal valueÊ ( ϕ, λ) for all observational points was determined using the obtained model coefficients, and the difference E i ( ϕ i , λ i ) between the normal value and the real observation E ( ϕ, λ) was found. If E i was bigger than 3σ , the respective observation point was removed from the analysis. After removing all "anomalous" points, the regression coefficientsÂ were determined again. Following this procedure, 17 outliers were found in D, 24 in H , 48 in Z, and 9 in F .
As mentioned above, the input data for producing the 2015.0 model are the last absolute geomagnetic field measurements (1980) updated with the prepared isoporic maps from secular measurements and neighbouring observatories. Annual mean values from 2015.0 to 2020.0 are computed by extrapolation of Panagyurishte observatory data, assuming an autoregressive secular variation to define the coefficients a 1 up to epoch 2020 as proposed by De Santis et al. (2003). We impose that the PAG observatory annual mean value of each element E t (PAG) taken at time t (in years) depends on the p previous annual means. We use a stepwise least-squares estimation of an autoregressive (AR) model of order p (De Santis et al., 2003): where p lies between p min and p max and is chosen as the optimizer of Schwarz's Bayesian criterion (Schneider and Neumaier, 2001). Calculation of the four AR models (for geomagnetic field elements D, H , Z, F ) returned p = 4 as optimal regression order, as was obtained by De Santis et al. (2003) as well.

Model results and analysis
Values of the BulGRF 2015.0 coefficients a 1 , a 2 , . . ., a 6 , using the least-squares regression method with the associated dispersionσ and correlation coefficient R are shown in Table 1. Plots of the confidence intervals of coefficients a j for the geomagnetic elements are shown in Fig. 4. The a 1 coefficients (Table 2) for each field element are computed for every year up to 2020 using Eqs. (2) and (9).
Since the 2015.0 BulGRF model was developed using independent polynomials for each element, it was necessary to check the self-consistency of the model. We checked if the geometrical constraint (Haines, 1990) is satisfied using the equation  (Finlay et al., 2016) and the latest edition of IGRF model (Thébault et al., 2015), which provides a reference field model for the same epoch. In Fig. 6 for each of the three "test" points D, H , Z, F values are displayed: (1) from real data, reduced to 2015.0 epoch, (2) from IGRF-12 model, (3) from CHAOS-6 model, and (4) predicted by the BulGRF model (together with the confidence interval). As can be expected and is seen from the plots, values of the IGRF-12 are very close to the present model except for the case of Balchik (H , Z, F elements) but all are within the confidence interval of BulGRF modelled elements. Balchik airport is situated in the northeastern end of the territory and thus suffers from boundary effects, resulting in greater in-  accuracy and wider confidence interval. Values obtained by the latest CHAOS model nearly coincide with the regionally predicted D and F elements for the three points and have opposite deviations for H and Z elements, resulting in underestimation of H and slight overestimation of Z.
Very well expressed however is the difference between the three models and the real data in PAG observatory. It is clear that if someone measures the magnetic field at a point on the Earth's surface, they cannot expect to get the value predicted by our model or the IGRF model. The reason for this is well known: there is a significant contribution which comes from the magnetized rocks of the Earth's crust -typically 200-300 nT, but it can reach bigger values in the anomalous regions. They are "shared" among the components but mostly affect the Z component (as it is in PAG and one of the airports). There are also signals, both man-made (traffic, electric facilities, etc.) and natural fields (from electric currents in the ionosphere), and the associated induced fields in the Earth but with smaller magnitude than the crustal field.
Concerning the case with PAG observatory difference, it can be clearly seen in Fig. 6 that calculated annual means are strongly below the values of BulGRF model and both global models (IGRF-12 and CHAOS-6) as well. It was expected because of the already mentioned anomaly found in the region before the observatory construction. This anomaly is also observed on the detailed exploratory magnetic anomaly map of Bulgaria prepared by Pchelarov (2000), where this region falls in a well-delineated crustal anomaly with 150 nT magnitude of Z (see for example Trifonova et al., 2012). Characteristics of the anomaly obtained using the BulGRF model are declination D shifted to the west by 20 and lower intensities in H -160 nT, Z -250 nT and F -297 nT.
The main point in the future will be to check whether this anomaly has a constant residual value in time or whether it varies along with the main geomagnetic field thus influencing the secular variation.

Conclusions
The advantages of the proposed 2015.0 BulGRF model with extension to 2020 can be summarized in three points: (1) a second-degree polynomial approximation of the geomagnetic field elements of Bulgaria is elaborated which is closer to the base level of the geomagnetic field over a limited region like the Bulgarian territory; (2) because of the long periods without secular measurements in Bulgaria there was a need of recent model compilation which now is fulfilled; and (3) BulGRF is an accurate model of the regional geomagnetic field which together with its simplicity makes it a useful tool for reducing to a common epoch the magnetic surveys carried out over the Bulgarian territory up to 2020. Data availability. Panagyurishte observatory data and secular measurements data are available from the World Data Center Edinburg (http://www.wdc.bgs.ac.uk/). The data from the 1979-1980 geomagnetic survey on the territory of Bulgaria is only available as a hard copy and it is possible to be used in the National Institute of Geophysics, Geodesy and Geography.