Crustal field recovery and secular variation from regional and global models over Albania

GM101 ABSTRACT The static geomagnetic field of crustal origin is optionally calculated by the recent global geomagnetic field models. However, their description in global scale tends to miss some local characteristics. The same can be inferred for the rate of the geomagnetic field changes i.e. secular variation (SV). In order to depict some particularity of crustal field in local scale for a small region like Albania, two regional models are constructed: one based on the Legendre’s polynomials and the other based on a linear approximation. Both models use data from different measurement campaigns in the Albanian repeat station network and a few data from a neighbouring country. The residuals produced by these models and by the recent global models: EMM2015, POMME-9 and CHAOS-5 are calculated and compared. The SV from regional and global models are also calculated. Substantial differences between SV calculated by global models and regional models are observed.


Introduction
The geomagnetic field is the superposition of magnetic fields that are generated by different sources located inside or outside the Earth. A small part of the geomagnetic field, generated in the magnetosphere and ionosphere, depends on the solar activity and is widely studied [Jacobs 1991, Voight 1991, Olsen and Stolle 2016. Because these sources are located in the upper atmosphere and above, their respective fields are collectively named as outer geomagnetic field.
The essential contribution, named core field, is the magnetic field generated by the liquid outer core. There is a smaller part, the crustal or lithosphere field, whose sources are located in the crust or lithosphere. Because the sources of these magnetic fields are located inside the Earth, these contributions are known as internal geomagnetic field [Backus et al. 1996].
The improvement of computational facilities, especially the rise of supercomputers, has enabled workers in this discipline to model with high accuracy the main field using a huge amount of data (ground and satellite). Consequently there have been conceived global models for the main and crustal field. The site https://geomag.colorado.edu/geomagnetic-and-electric-field-models.html provides information for several magnetic field models. Regarding this paper, there is information for the models EMM2015 and POMME-9. Detailed information for the CHAOS-5 model can be found in the site http://www.spacecenter.dk/files/magnetic-models/CHAOS-5/.
Most of the crustal field originates from the permanent magnetization of ferromagnetic mineral deposits under their Curie temperature located in the crust. This field normally is static and depends on local geological structure. Remanent magnetisation is strong near the surface due to temperatures below the Curie point and dominates the short wavelength crustal field [Maus and Haak 2002]. There are also other minor sources and details can be found elsewhere [Lowrie 2007]. There is some part of the crustal magnetisation that is induced, i.e. being (roughly) parallel and proportional to the local actual field. Due to higher temperatures in the lower crust (approximate or above the Curie point) the remanent magnetisation is unstable and the induced magnetisation seems to dominate the long wavelength crustal field. This variable field normally depends on the main field and has its temporal variation [Maus and Haak 2002]. However, the static compound of the crustal field is more important from a geological perspective. Furthermore, it is a common assumption in global magnetic field modelling that the crustal field does not change in time [Macmillan and Thomson 2003].
In this paper we aim to study the crustal field for small regions like Albania by using data from repeat station measurements in Albania and the neighbouring country of FYRO Macedonia. In order to recover the static crustal field from the observed data we devise a regional model available for small region of Earth's surface. Regional models are normally based on a denser data grid than global models and for this reason produce more accurate results over the same region than the global ones [Haines 1990]. One of our regional models consists on the expansion of the geomagnetic field components in Legendre's polynomials. The use of Legendre's polynomials means that the basis' functions of the model are orthogonal in the area over which the model is applied. This means that the model has much less dependence on the abundance of available data [Burdelnaya et al. 1999] thus producing more accurate results.
Usually the global models are based on the fact the geomagnetic field in the lower atmosphere is curl-free, so it can be expressed in terms of a scalar potential. This potential is expanded in spherical harmonics of different degrees. Models are distinguished from each other by the degree of expansion of the core field L c and the maximum degree used L max [Duka et al. 2016]. The remaining part, from L c to L max , represents the crustal field also considered to be static. In this paper we have used several models like the model EMM2015 [Maus 2010], the POMME-9 model and CHAOS-5 [Finlay et al. 2015]. Other details can be found in the corresponding websites (see above). All these global models do not include the year 1995 in their validity interval and therefore there are no results for this epoch.
The geomagnetic field is a complex system that is subject to temporal variation where are observed frequencies that span a very wide range of orders from milliseconds to millions of years. These variations are indeed the manifestation of many physical phenomena. It is widely accepted that the high frequency variations have an outer origin. Whilst the variations on the time scale of years, named secular variation (SV) have mostly internal origin. It is considered to rise from the details of the flow of molten iron in the outer core [Merrill and Mcfadden 1999]. The SV is obtained for the Albanian region from both the global models and regional models.
The structure of the paper is as follows: in section 2 we give a short account on the data we have used in this study. In section 3 we describe the regional mod-els. In section 4 we report the results recovered by both regional and global models. In section 5 we discuss SV and the respective results and finally in section 6 we give some conclusions and discuss some problems worth to be investigated in the future.

Data
In the territory of Albania there is no geomagnetic observatory. The same is valid for the parts of the neighbouring countries close to the borders with Albania. Consequently the models (explained in section 3) are built using only data from the campaigns of measurements in geomagnetic repeat station. Some historical campaigns were conducted by non-domestic institutes before 1990 [Bolz 1963, Bock 1948, Rössiger 1941].
There are reported data of non-vector geomagnetic field measurements (total field F) over Albanian territory during 1990-1994 [Kane et al. 2005]. There are PEQINI ET AL.
λ λ λ λ also some sparse data for the Z (downward) component of the geomagnetic field and F for epoch 1961 or F and declination for epoch 1943 [Duka and Bushati 1991]. All these data suffer poor temporal coverage and do not provide information about the X and Y components of the magnetic field. Thus we have considered those data not useful for our purposes. Therefore, for the Albanian territory, we have used only the data from the campaigns of the epochs 1995 [Chiappini et al. 1997[Dominici et al. 2007], 2010 [Dominici et al. 2012 Figure  1). These data are obtained from the British Geological Survey, BGS (http://www.geomag.bgs.ac.uk/ data_service/data/surveydata.shtml).

Regional models
Two regional models were built. The first model of normal field is built up using the Legendre's (N-leg) polynomials. A similar but more extended model that covers a much larger region like the easternmost part of Far East Asia is described elsewhere [Burdelnaya et al. 1999]. N-leg (as well as the other model) is built in terms of latitude φ and longitude λ. This model as well as the simple linear model (see below) cover a rectangular shape wholly including Albania, ranging from 39.5o to 43.0o N and 19.0o to 21.5o E. There is chosen as the arbitrary point the central point of the selected area with coordinates 41.250o N and 20.250o E. The components of the field as well as the total field are approximated by linear expansions in normalised variables xi and yi (the Legendre's polynomials) given by (1) The set of these functions compose the basis in the mathematical space of functions for the model and are all orthogonal to each other into the selected area. In our models only the first terms of approximation are considered. Calculations (respective results are not shown) indicate that the inclusion of terms of higher order does not bring any considerable improvement to the actual results.
Mathematically the equation of the model reads (2) where B i represents each of the components or the total field at the point x i , y i of observation (φ i , λ i ) and a, b, c and d are the model's coefficients. If the coefficients do not depend on time, then the model is simply spatial. However we consider them to be linear functions of time such as a = a 0 +a 1 t, so we can write (3) The linear approximation is accurate enough for the very short period of time considered. We can use matrix formulation consequently writing (4) where y is the vector of data, M is the matrix containing the functions of the basis in the observation points and c is the vector of coefficients that are to be determined. Here "leg" refers to the N-leg model. The vectors and the matrix in this equation are of dimensions N×1, N×8 and 8×1, where N is the number of data.
The second model uses a linear approximation (N-lin) where a given component of the field as well as the total field are considered to be linear functions with latitude and longitude. Analogous models can be built for inclination and declination whenever is possible. The equation of the model reads (5) where B i represents each of the components or the total field at the observation point (φ i , λ i ) and A, B and C are the model's coefficients. With "c" subscript are denoted the latitude and longitude of the central point of the selected zone over Albania. Considering these coefficients as linear functions of time we can write (6) Physically the coefficients B 0 and C 0 are the North -South and West -East gradients of the magnetic field while B i and C i represent the rate of change of these gradients. In matrix notation equations 6 can be written as (7) where "lin" refers to the N-lin model. The vectors and the matrix in this equation are of dimensions N×1, N×6 and 6×1, where N is the number of data.
The coefficients in both models are calculated using a standard least squares procedure and are obtained for each of the components and for the total field. The standard deviations for these coefficients are computed in a classic way described in the literature [Richter 1995]. We initially calculate the correlation matrix where the diagonal elements are the squares of the standard deviations of the coefficients. This tech- nique yields the same standard deviation for the analogous coefficients of the individual components and total field.

Results of regional models
The data at our disposal are initially reduced at a height of 0 km and the coefficients of the model have values that pertain to this specific height. The reduction is done by the classic downward reduction considering that the magnetic field attenuates as a magnetic dipole field. Reduction at other heights (results are not reported here) does not produce any substantial difference hence we limit our study to the height 0 km only.
In Tables 1 and 2 are shown the coefficients for each component and the total field as well as the standard deviation for both regional models. The first thing highlighted by results is that the magnitudes of the respective coefficients for both models are almost the same. This means that the results for both models are very similar. This result is expectable because for such small areas like that over Albania both regional models are practically dominated by linear terms. Hence we can limit our analysis to the N-leg models only and similar results can be drawn for N-lin model also. The regional models used in this study produce consistent results with prior studies when the coefficients are considered to be static, i.e. by fixing the time. For the epoch 1994.75 [Chiappini et al. 1997], we find that the coefficients of our models are almost identical with the coefficients recovered by these authors.
The coefficient of the free term is the largest and very approximate to the values of the respective components or total field. The coefficients of the x and y term, that represent the north -south and west -east gradients of the geomagnetic field, are two -three orders of magnitude smaller than the free term coefficient. Practically these coefficients have magnitudes of the same order as the analogue coefficients of the simple spatial normal field model for the 1990.0 epoch [Duka and Bushati 1991]. The coefficients of time -dependent terms are too small yielding a slowly varying field, describing a typical temporal variation of the geomagnetic field known as secular variation, SV [Backus et al. 1996]. Generally, the standard deviation is from one to several orders of magnitude smaller than the respective values of the coefficients, showing the stability of the models. The only exception is the value of d 0 for the total field F which is almost twice of the respective standard deviation.
We considered the residuals as the differences be-tween the observed values at the locations of the repeat station and the expected values under the models for the same locations. In principle such residuals should represent the non-modelled field of internal contributions, mainly from local sources, i.e. crustal static field. Also the repeat stations are not located in the anomalous zones, i.e. they do not drastically influence the models. Figures 2 and 3 represent the results for the regional models. The maps for the years 2005 and 2010 are not shown here but they have very similar results.
In figure 2 are shown the maps of normal field (dotted lines) and residuals (solid lines) for each component or total field. The Y component shows a localized anomaly between the Dibra's region and the south-eastern (not shown here) region which are well-known magnetic anomalies. The other well-known Mirdita's region anomaly is not recovered. However, the other maps show a smoother spatial variance. In figure 3 are shown the analogous maps for the N-lin model, which are almost identical to the maps of the previous model.

Results of Global model
The global models used in this study have various characteristic degrees like: EMM2015 with L max = 720 and Lc = 15, POMME-9 with L max = 133 and L c = 16 and CHAOS-5 with L max = 90 and L c = 20 (this variant 1.52×10 -2 9.70×10 -2 8.24×10 -2 8.05×10 -2 5.99×10 -5 b 1 7.95×10 -4 -1.51×10 -2 -2.64×10 -3 -1.57×10 -3 1.70×10 -4 c 1 -3.74×10 -3 -1.08×10 -2 -3.36×10 -3 -5.20×10 -3 1.64×10 -4 d 1 5.53×10 -3 -2.91×10 -2 1.32×10 -3 2.03×10 -3 4.51×10 -4    is used in the present study), where L c and L max are the maximal degree of expansion of the core field and the maximum degree used by the global model. The field is calculated for each global model at the same geographical points (with the same reduced altitude) where the measurements are performed and on the points of the grid over the selected zone which is used to compute the normal field and the SV (see section 5 for SV). We obtain the residuals from the global models in two distinct ways. The first one is to calculate the residuals between the measured values and the values of the inner field (core field + crustal field) predicted by each of the global models (Figures 4, 7). The second one consists in calculating the residuals between the measured values and only the values of the core field given by each of the global models (Figures 5, 8). Additionally, we show the crustal field values predicted by each global model (Figures 6, 10 and 11).
In figure 4 are shown the maps of the residuals obtained by the first way of using the EMM2015 model (solid lines) and the region field predicted by this model (dotted lines). In figure 5 are shown the maps of residuals calculated as the difference between measured values and the values of the core field obtained by the same global model. There are a limited number of points of measurements outside Albania and this is the reason why the isolines are centred in Albanian territory.
The residuals calculated in both ways are very similar despite of some minor displacements of the isolines because there is some contribution of the outer magnetic field as well as the non-modelled crustal field which is not captured by the global models used in this study [Duka et al. 2016]. Their removal is not possible due to the lack of temporal coverage with appropriate resolution. However, analysing these maps one can safely state that both methods are very useful in isolating the crustal field signal within a given region and that they produce very similar results.
In Figure 6 are shown the corresponding maps of components of total magnetic field of the crustal field as given by EMM2015. The map is extended over the whole region. However, we should focus in the territory of Albania because only here we can be make a comparison. The maps of different components of crustal field show differences, but in the Y component maps there are common and visible local anomalies like as in the central Mirdita region and in the eastern Dibra region. However, in general there are considerable differences between maps of crustal field obtained by residuals techniques and different global models (compare Figures 4, 5 to Figures 7, 8). These differences are not only in the magnitude of the residuals but also in the shape of the isolines.
In the Figures 7 and 8 there are shown the corresponding maps of residuals (observed values minus core field values given by a global model) for the POMME-9 and the crustal field of this model. Similarly, the results obtained by calculating the residuals of CHAOS-5 model are resented in the fFigure 9. The residuals technique applied to the CHAOS-5 model yields for the X and Z components and consequently for the total field F, a systematic difference from other global models. There is a shift in the residuals magnitude even for the Y component though much smaller than for the other components. We think that such systematic difference is due to the fact that the core field for CHAOS-5 extends to degree Lc = 20 while the core field for other models extends to degree 14 or 15. Therefore when we calculate the residuals we subtract from the repeat station data a larger core field.
Despite of the systematic shift in magnitude, the qualitative similarity of the CHAOS-5 model with other models is not affected. There are slight changes in the configuration of the isolines for this model; however the "big picture" of the field remains unchanged. Recall that each of the global models is built by considering the crustal field to be static hence these fields are the same for all epochs the models are valid.
The crustal field provided directly by each global model (with spherical harmonics from Lc+ 1 to Lmax) seems to be different. The crustal fields provided by both POMME-9 and CHAOS-5 (Figures 10 and 11) are much smoother than the respective field given by EMM2015 ( Figure 6). These differences are not a surprise because the EMM2015 has a maximal degree of 720 hence the characteristic wavelength, in other words the resolution of reconstruction of the crustal field anomalies, is of the order of 50 km. Important conclusions can be drawn from the comparison of Figures 2, 6, 10 and 11. The application of the residuals technique on the normal field models yields a crustal field that is more compatible with the geological structure of the territory than that of the global models.

Secular Variation (SV)
Unfortunately, we do not have a lot of temporal data over the Albanian region to provide a detailed description of SV. However, we have calculated the SV for any component (X, Y, Z) and the total field for the epochs 2005, 2010 and 2015 (only for the models EMM2015 and POMME-9). The year 1995 was not included because none of the global models is valid for that year. We have calculated SV for these years as the CRUSTAL FIELD OF ALBANIA FROM REGIONAL AND GLOBAL MODELS       The data at our disposal are temporally reduced to the middle of a given year. This means that the external contributions of period less than 1 year are almost removed from the SV. Also, according to equation 8, it is clear that the crustal contribution due to remanent magnetization -is practically removed. Therefore, we can conclude that SV calculated by this approach represents the rate of change of the magnetic field of core origin. To eliminate the problem of the sparse network of repeat station, we apply equation 8 for synthetic fields generated by all models over a grid built upon the area in study with a step of 0.25° in both latitude and longitude. The exception from this procedure is the EMM2015 model which provides SV directly.
In the figures 12 -15 are shown the maps of SV for the year 2010 (2015 is out of the validity period for CHAOS-5). With the exception of POMME-9, the other global models offer a smooth picture of SV where it is almost the same for the whole territory. On the other hand, the POMME-9 offers a rather complex picture of SV. This result could be an overestimation because the SV for such a small area should have a smoother configuration of isolines. One can see that the normal (regional) model offers a better alternative. Not very smooth (EMM2015 and CHAOS-5) and not very complex (POMME-9).

Discussions and conclusions
The present paper discusses the methods to isolate the static crustal field from available data. In principle there are some ways how this can be achieved that are explained in the result sections. The calculation of residuals with both regional and global model yields similar results. This similarity is extended to all the years in study (except 1995 for being out of the validity period of EMM 2015) promising to be a very accurate method for recovering the static crustal field.
It is noticeable that the plots of the crustal field estimated directly by EMM2015 model (spherical harmonics from Lc = 15 to Lmax = 720) are different from the plots of the crustal field estimated by the residuals (differences between observed values and predicted by the model values). Almost the same is noticed for the two other global models. We point out that the residuals are calculated in the 11 points of Albanian repeat stations, while the crustal field of the models is calculated in a regular and denser network of points (see section 3). It is obvious that the crustal field estimated directly by the global models varies smoothly over small regions like Albania, and this variation would be smoother as Lmax gets smaller (720 for EMM2015, 133 for POMME-9 and 95 for CHAOS-5). In principle the EMM2015 model should provide a more detailed spatial variation of the crustal field. However the POMME-9 model provides a more realistic crustal field. May be this is due to fact that this model used more data from the considered region than other models.
On the other hand the crustal field calculated by the residuals technique applied to the normal field models have very small differences when are compared to each other. Actually there are observed some small differences between the values of residuals for the both models. There are two additional terms in the N-leg model that account for these differences.
The comparison of the crustal fields obtained by the residuals technique for the regional normal field models and the global models show quantitative and qualitative similarities, apart the CHAOS-5 model where the similarity is only qualitative due to its large Lc. These similarities are persistent through all epochs indicating that the obtained field by the residuals technique is the essential part of the static crustal field.
Comparing different plots of crustal field estimations by different models and different methods some anomalies are identified in the north-eastern and eastern Albania like in the Dibra's region as well as southeastern regions of Albania that are supported by geological evidences (Kane et al., 2005). However, we could not identify consistently by all models and methods the anomaly in the Mirdita zone. May be this is due to the fact the Albanian region presents a sharp spatial variation in the geological properties of the different formations and there is needed a denser repeat station network to detect such anomalies. A closer look at the Y component maps produced by calculating the global models residual shows clearly the anomaly mentioned above but it is not evident in the maps of the X and Z components.
The model EMM2015 produces a more distinctive crustal map than the other global models. The crustal field seems to be quite complex in the case of this model, whilst is quite smooth in the case of the other global models. However, these results are not conformant with geological and geophysical results [Kane et al. 2005]. These authors report a crustal field that does not fit with the crustal fields given by the global models. The same is valid for the regional mod-    We must consider that these authors had much more magnetic data and their reconstruction is thus more accurate. According to Duka et al. [2016], in the geomagnetic data observed at ground observatories there is a weak contribution from outer sources that is not modelled by the global models. Due to the small area under study these contributions exhibit the same time variation on every repeat station. These contributions can be removed in the case when there is a Geomagnetic Observatory in the region with at least several years of continuous registrations. Unfortunately this is not the case of Albania. Therefore, in the residual maps calculated here by various global models there are some weak outer signals included. These additional signals may explain the small deviations of the residual plots of the fields produced by the different global models.
Unfortunately, our results have not explored details of the geological structure of the upper Crust of Albanian region, like as the Kane et al., 2005 work about the ophiolitic belt. The reason is they have used a very dense data set of total field. Apart of this, in our work we could not consider the relief contribution to the data that is not negligible in the mountainous part of Albania that constitutes a major part of the area under study. But, we can say that the crustal contribution recovered here should represent mostly the remaining magnetization contribution of the lithosphere (deep Crust).
The SV given by the regional models seems to be more reliable than the SV given by different global models. The SV cannot be subject to large changes for small areas like Albania (as given by POMME-9), but also cannot be very smooth (as given by EMM2015 and CHAOS-5). However, the description of SV would be more accurate if there were a better time coverage.
As we mentioned in the introduction section, the amount of data at our disposal is limited and with low temporal resolution and this affects the quality of the regional models. In order to enlarge the quantity of data, synthetic data produced by global models can be used to fill the parts of Albania not covered with measurements. Something similar is described in Burdelnaya et al., 1999. This approach and its consequences remains to be tested in future studies.