Recovering external contribution from the monthly mean series of a given geomagnetic observatory

The differences between monthly mean values of the observed geomagnetic field and monthly values predicted by different models of the internal geomagnetic field (named “model biases”) for the time period 2000-2015 at several geomagnetic observatories are analyzed. We notice that increasing the maximum degree of the model is not always followed by the decrease of such “model bias”. The time series of these “model biases” reduced by their average resulted to be approximately the same for all models and should represent the external (non-modeled) contribution to the observed geomagnetic field. These time series for different observatories (close or away to each other) are compared and their power spectra are analyzed. Such spectra have common features like the annual and semi-annual variation with some possible sporadic cases of


Introduction
The magnetic field measured at any geomagnetic observatory is the result of several magnetic contributions generated by various sources.More than 90% of the measured field is of internal origin and is mostly generated in the Earth's outer fluid core.Also, of internal origin is the magnetic field generated by the remanent magnetization of the crust and by the induced magnetization of the crust produced from the core and external fields.The known contributions of external origin are magnetospheric and ionospheric fields, fields of ring currents and aligned field currents.
The global models of the geomagnetic field and its secular variation (SV), like IGRF, WMM, EMM, POMME, CM, GUFM, MF, etc., are based on spherical harmonic analysis (SHA).The magnetic field potential (W) is expressed as the sum of spherical harmonic functions.The internal sources contribution to the radial variation of these functions is quite different from the external sources contribution: respectively (a/r) l+1 and (r/a) l , where a is the reference radius of the Earth, and l the spherical harmonic degree.The internal source magnetic field is presented by B r , B i , B { -orthogonal components in the geocentric reference frame or X (north), Y (east), Z (downwards) as: Using different data sets of ground or/and satellite measurements in disposal, different models define the Gauss coefficients (g m l ; h m l ) of spherical harmonic expansion up to a maximum number of degrees L and their time derivative up to a maximum number of degrees L' < L.
Generally, the denser and numerous the data used by the model, the greater is the maximum degree of harmonics (L) and the shorter is the minimum wavelength of harmonics i.e. more detailed is the model description of the geomagnetic field, according to the formula [Backus et al. 1996]: (2) The separation of internal contributions into the core and the crustal (actually lithospheric) origin is based on the spherical harmonics spatial power spectrum [Lowes 1974, Langel andEstes 1982].It is known that the spherical harmonics up to maximum degree L = 13 represent the core field, while the spherical harmonics of degree greater than 13 represent the Lithosphere field and there is an overlap of contributions in between 12 and 14 [Cain et al. 1984].In fact, lithosphere sources have some small contributions to the spherical harmonics of degrees <13 and core sources have also some small contribution to the spherical harmonics of degrees >13.Some models aim to describe the Earth's magnetic field with high spatial and temporal resolution, for example CHAOS-5 model robustly determined spherical harmonics up to degree L = 90 for the lithosphere field, up to degree L = 20 for the core field and up to degree L'= 16 for the time-varying core field [Olsen et al. 2014].Temporal changes of the field of internal origin up to degree 16 are usually attributed to changes in the core field itself and temporal changes of the field of internal origin produced by the induced part of the lithosphere magnetization could dominate the core field signal beyond degree 22-24 [Hulot et al. 2009].The time-varying lithosphere field dominates and conceals the time-varying core field beyond that critical degree (22)(23)(24), just as the permanent lithosphere field dominates and conceals the permanent of the core field beyond degree 14 [Hulot et al. 2009].
It is more difficult to separate the crustal contribution into parts coming from remanent magnetization and from induced magnetization.[Lesur and Gubbins 2000] developed a new method for such separation, but its application for several observatories was unsuccessful.
The so named "crustal bias" [Mandea and Langlais 2002] in an observatory is estimated by comparing the magnetic components measured in the observatory with those predicted by a geomagnetic model truncated to its nuclear part (i.e., up to degree and order of 13).These models are based only on the satellite data sets [Bloxham and Gubbins 1986], where the crustal sources contributions are negligible, therefore the differences between observed values and model values are considered as the signature of the crustal field of induced or remanent magnetization origin.Such "crustal biases" of annual means over 42 years period are well investigated [Verbanac et al. 2007a] for 46 European geomagnetic observatories by using models based on Magsat, Ørsted, CHAMP and SAC-C satellites.Also, the temporal evolution of the observatory monthly means (over 9 years period) "crustal biases" are analyzed by Verbanac et al. [2015] using the G-field model [Lesur et al. 2015].Even in this analysis, "crustal biases" are considered the differences of the magnetic components measured at the observatory with the values predicted by a model obtained from satellite data only.Verbanac et al. [2007b] attempted to separate, interpret and explain the external field signals obtained by subtracting from the annual means the core field predicted by models truncated at the maximum degree 14.From the obtained external signal they removed magnetospheric contributions of the external field predicted by the POMME-2.5 model and a Sq variation predicted by the CM4 model.Averaging these residuals over representative European observatories, they found interesting results regarding the correlations of the residuals with external sources such as solar cycles.
Here we will study a kind of "bias", which is estimated as the difference of the monthly mean values of the geomagnetic field components registered at a given observatory with the respective values predicted by the geomagnetic field models based on ground and satellite data.The models are truncated at different maximum degrees providing different "biases" that have mostly signatures of non-modeled local sources of the crust and of external sources.Different models provide different "biases" (named "model bias") at the same observatory.We will compare the results of three different models for different observatories that are geographically close to each other (WIK: latitude 48.265°N, longitude 16.318°E, altitude 0.4 km; NGK: latitude 52.072°N, longitude 12.675°E, altitude 0.078 km; FUR: latitude 48.17°N, longitude 11.28°E, altitude 0.572 km) and two other observatories located far away from them (HER: latitude 34.43°S, longitude 19.23°E, altitude 0.026 m) (KAK: 36.23°N,longitude 140.18°E, altitude 0.036 km).
We expect that the difference ("model bias") between the observed value and predicted value given by a model should decrease as the maximum degree of spherical harmonics of the model increases.But this does not happen always for all components at every place on Earth, as it will be shown in Section 2.
The global geomagnetic models supply Gauss coefficients and their secular variation for every 0.5 year, 1 year, 2.5 years, etc.Then the geomagnetic field at a given time is calculated from such models by interpolating the Gauss coefficient values in between time intervals, i.e. each model has limitations on time distinction.When the monthly values of the internal geomagnetic field are predicted by a model, the time distinction of one month is considered.On the other side, the monthly values series of observatories are calculated by averaging the observatory minute registrations for all days of a month.Therefore the geomagnetic field variations with period less than a month (pulsations, quite daily variations Sq, substorms/storms effects) are principally removed.
By removing the averages of the differences between the predicted monthly series and the respective monthly mean of observed value series, we can recover the un-modeled external source contributions to the geomagnetic field at the given observatory (see Section 3).
The comprehensive model, CM4 [Sabaka et al. 2004], covers the whole time interval from 1960 to 2002.It has been derived from quiet-time POGO, MAGSAT, Ørsted and CHAMP satellite data in combination with observatory hourly means.The internal field is expanded to different number of maximum degree of spherical harmonics up to 65, resolving magnetic anomalies down to 611 km.The model also calculates primary and induced magnetospheric field, primary and induced ionospheric field.It offers calculation of the geomagnetic field by two steps.In each step, the user inputs the first and second maximum number of degrees.Using codes from (http:// core2.gsfc.nasa.gov/CM/codes.html), the monthly values series of core field (L = 14) and lithospheric field (with different values of maximum degree, from L =15 to L = 65) at different observatory locations are generated.
POMME-9 (http://geomag.org/models/POMME-9.html) is an internal field model representing the geomagnetic field in the region from the Earth's surface to an altitude of a couple of thousand kilometers [Maus et al. 2004[Maus et al. , 2006[Maus et al. , 2010]].The time variations of the internal field are given by a piece-wise linear representation of the spherical harmonic (Gauss) coefficients of the magnetic potential.
POMME-9 was produced from CHAMP satellite vector magnetic measurements from July 2000 up to September 2010, Ørsted satellite total field measurements from January 2010 to June 2014 and Swarm satellite vector magnetic measurements from December 2013 to January 2015.It has the same parameterization of the magnetospheric field as POMME-6, POMME-7 and POMME-8.POMME-9 extends the maximum degree of spherical harmonics to 133, resolving magnetic anomalies down to 300 km.The Gauss coefficients of every year (from 2000 to 2015) are embedded into the header (.h) files for the core and crust field.The program calculates the core geomagnetic field every day from the starting year 2000, while the crust field coefficients are static from degree 16 to 133.
Enhanced magnetic model (EMM2015) based on the ellipsoidal harmonic representation of Earth's lithospheric magnetic field [Maus 2010] represents the main field with maximum degree L =15 and the lithospheric field up to L = 720.It was compiled using a vast amount of data from satellite, marine, aeromagnetic and ground magnetic surveys.It also includes data from the European Space Agency's Swarm satellite mission.This model (https://www.ngdc.noaa.gov/geomag/EMM/EMMSurveySPH.shtml)covers the 2010-2020 period resolving magnetic anomalies down to 56 km.We tried to find a common time interval for all considered models.All considered models have a short common time interval 2000-2002.We firstly compare the monthly mean series of the NGK geomagnetic observatory with respective monthly values series predicted by these models for that short time period (Figure 1).One can see that series predicted by CM4, CM5, EMM2015 represent a smooth variation, while the series predicted by POMME-9 model have some irregular variations much like the series of observations.Even the monthly values series received by averaging the daily values series generated by POMME-9 show the same behavior.It seems that this model fits better to the observed monthly values series than other models.The CM5 model represents a greater baseline jumps from the observed monthly values series.

NON-MODELED EXTERNAL CONTRIBUTION
In Figure 2 Comparing the results of different models, one can see that for the X component, the POMME-9 model fits the observations better (see also Table 1, where the values of averaged "bias" and their errors for different models and different maximum degree are presented).While for the Y and Z components the EMM2015 model (respective averaged "bias": −1.9 nT and −7.2 nT) fits better with the observations than the POMME-9 model (respective averaged "bias": 4.9 nT and −56.1 nT).The CM5 model fits the observations worse than the two other models for all three components.It seems that CM5 model gives worse results than CM4 model.See for example Figure 4, where the results of CM5 at WIK observatory are presented, and Table 2, where the averaged "biases" of CM4 model at NGK for a long period  are presented.
The same comparison for two observatories, Furstenfeldbruck (FUR) and Wien Kobenzl (WIK) observatories are presented respectively in Figure 3 and Figure 4.As it is seen in these figures, the POMME-9 and EMM2015 models fit the observed data better than CM5 model.Also all models provide the same time variations.
In order to study the impact of the maximum degree number of the model on the fitting quality between  the observed and predicted series, the monthly values series are generated by EMM2015, POMME-9, CM4 and CM5 models with different numbers of maximum degrees at different observatories.Some results are presented in the plots of Figures 5, 6, 7, and in Table 1.
In Figure 5, one can see that decreasing the maximum number of degrees from L = 133 to L = 95, or L = 65 for the POMME-9 model worsened the fitting quality for X component, but not for Z component.
In Figure 6 (EMM2015 model), one can see that changing the maximum number of degrees from 720 to 520 gives better fitting for the X component (see also Table 1), while for Y and Z components the fitting quality is poorer.When the maximum number of degrees is decreased to 120 we see the fitting is worsening for all components.
In Figure 7       X and Z components are considerable increased when L changes from 100 to 85 then shift a little when L changes from 85 to 60 and from 60 to 40.
Regarding the CM4 model, as its validity period is 1960-2002 we have applied this model for the period 1987-2002 to WIK observatory and for the period 1960-2002 to NGK observatory.As can be seen in Figure 8 and in Table 2, there is a tendency of decreasing of the average "bias" by increasing of the maximum number of degrees from 25 to 65 for the Z component at NGK observatory, that is not so clear for X and Y components Almost the same happened for WIK observatory (the results are not shown here).

Recovering external contribution
We have considered the averaged "biases" as contributions of the crustal sources: remanent magnetization and induced magnetization from the internal sources that are not modeled by the given model.Therefore, the residuals of subtracting such average from the model "bias" should represent the contribution of non modeled external sources at a given place (observatory).
The time series of differences between the geomagnetic field values observed at an observatory (that are the superposition of magnetic fields from internal and external sources) and respective values predicted by a model of the internal geomagnetic field, are: where i = 1, 2 …N indexes the series of months, k indexes the different models.The averages of these series are: where the first term is the averaged of the internal part not modeled by the k-th model and the second one is the average of external contribution to the observations.Subtracting these averages from the series of differences  (2) (3) 2015.While the averages x -obs (ext) provide the averaged "bias" at NGK observatory that is not modeled by the given model.Their values are presented in Table 1.
The same calculations are carried out for different observatories.When plotted for the same observatory, the series of residuals calculated by different models should collapse together.Here we present the results for FUR nearby NGK observatory (see Figure 10) and two other observatories far away from NGK: HER (see Figure 11) and KAK (see Figure 12).In all figures one can see that after removing the averages, the "model bias" is the same for all models.Again, we noticed that averaged "bias" (non modeled static contribution) for CM5 model is greater than for other models (see Table 3, where the averaged "bias" for different models and different components are presented).

Comparison of external contributions from different geomagnetic observatories
Plotting the external contributions, calculated as described in Section 3, at three nearby observatories NGK, FUR, and WIK, one can see (Figure 13) that the external contributions in these observatories are almost the same.This can be seen in all components.The external contribution to the eastward component (Y), which is supposed to be the least affected by the external fields [Mandea et al. 2010], looks like noise.
In Figure 14, the external contributions at three observatories (NGK, KAK and HER), located far away each other, are plotted.It can be seen that the differences between such contributions, at the three observatories that are more than 100° in latitude or longitude apart, are small ones, especially for the X component.While the time series for the Z components are similar, but the series of north observatories (NGK, KAK) are inverted in phase from the series of the south observatory (HER).This reflects the sign change of the Z component from North to South hemisphere.Regarding the Y component, one can notice the expressed noise-like character of these series and the phase inversion can be noticed more between eastern observatories (KAK) and less eastern observatories (NGK and HER).

Spectral analyses of the external contribution
The external contributions to the monthly mean series of the geomagnetic field observed at different observatories should contain signatures from the same external sources.In order to investigate any common contributors in the series plotted in Figures 9 to 14, we have applied the same spectral analysis method (FFT) on all of these series.The frequency contents (until the Nyquist frequency) of power spectrum of   In these plots, one can see several peaks that are common for all series, lightly shifted for different components, different observatories and different models.The first one being in the frequency interval 0.0039÷0.0055month -1 (period interval 21.37÷15.15years), should belong to the time series length (16 years).The second peak is in frequency interval 0.02÷0.025month-1 (period interval 4.16÷3.2years) and the third one is in the frequency interval 0.049÷0.051month -1 (period interval 1.7÷1.6 years).These two last peaks are not so clear.However there are two distinct and clear peaks in the intervals: 0.082÷0.086month -1 (period interval 1.01÷ 0.97 years) and 0.165÷0.168month -1 (period interval 6.06÷5.96months).
The amplitudes of these peaks are different for different components and different observatories.In case of POMME-9 model, the greatest peak is reached at KAK observatory for X and Y components and at NGK observatory for the Z component at the first frequency and the feeble maximum is reached at HER observatory for the Z component.The reason is that the absolute values of X and Y components at KAK observatory are greater than those at other observatories.
It can be noticed that the frequency content of the Y component for the POMME-9 model (only for this model) represents a series of small peaks that correspond to the frequencies: 0.3555 month -1 (period 2.82 months);     0.4375 month -1 (period 2.3 months); 0.4766 month -1 (period 2 months).Maybe these frequencies can indicate some external contributions with seasonal origin.
We think that such distinction of the Y component of POMME-9 model can be explained by the way the time series of this model is generated.The model is more accurate if the present state of the magnetosphere is provided as an input in the form of magnetic indices [Lühr and Maus 2010].Unfortunately, these indices are presently not available in real-time, so during real-time evaluation of the model we have to set these indices to zero.Then the time series (especially of Y component) of this model should bear some additional magnetosphere contribution.
Interestingly none of the contributions with seasonal origin are visible in the power spectrum of the Y component in both EMM2015 and CM5 models.However the Y component in the case of the CM5 model has a spiky power spectrum, especially in the high frequency part, although there are not any distinguishable peaks.This fact indicates the contributions from external sources of different time scales.
A major problem in the spectra of finite length time series is that they suffer from "spectral leakage".While the infinite-length signal has its power concentrated exactly at the discrete frequencies, the truncated signal has a continuum of power "leaked" around these discrete frequencies.This "spectral leakage" is more evident when data series are short as our monthly time series are (192 month long).Thompson's multitaper method (MTM) provides an improved PSD estimate [Thomson 1982] by using a bank of optimal bandpass filters.These optimal filters, i.e. with minimal leakage, are derived from a set of sequences known as discrete prolate spheroidal sequences (DPSSs, also known as Slepian sequences) [Percival et al. 1993].Another prob-lem when working with truncated time series is variance.It is shown that while reducing spectral leakage one can increase variance and viceversa.Thus an important issue is the appropriate choice of parameters so that a trade-off between spectral leakage and variance is achieved [Percival et al. 1993].Wardinski and Mandea [2006] used this method in order to investigate the geographic dependence of the annual and semi-annual variations.They use the notation of previous workers, defining K, the number of Slepian tapers, and p which is basically the number of points in the bandwidth W of the time series in the frequency domain.The relation between the two parameters is K = 2p − 1.It is proven that only 2p − 1 tapers are good enough for reducing spectral leakage, while averaging their results [Percival et al. 1993] also helps reducing variance.For further details we refer the reader to formulas (1)-(3) of Wardinski and Mandea [2006].
The MTM method with parameters described above is provided by "pmtm" Matlab function.For large values of time series length N, one can construct the sequence of Slepian tapers for all the values of p usually used (1, 2, 2.5, 3, 3.5, 4).In our case (N = 192), we have used p = 3.5, that means K = 6.So we use 6 Slepian tapers which is considered accurate enough for short time series like what we use in this study [Park et al. 1987].Apart from the power spectrum density function of frequency: PSD ( f ), the "pmtm" Matlab function calculates the confidence interval.According to our calculations, for a 95% confidence level, the true spectral density function lies within the interval from 0.52 to 2.3 of calculated PSD ( f ) function.Decreasing the confidence level one can reduce this interval, but the form of PSD ( f ) function does not change.
We present here the results for the time series of three different observatories (POMME-9 model; see Figure 18).Similar plots are obtained for EMM and CM5 models that are not shown here.Beside the greatest peak corresponding to the series length (192 months), there are two clear and wide peaks centered in the frequencies: 0.08 and 0.16 month -1 corresponding to annual (12.5 months) and semiannual (6.25 months) variations.In case of the Y component PSD (only for POMME-9 model) there appear to be three peaks centered in the frequencies: 0.36, 0. 44 and 0.48 month -1 , corresponding respectively to the periods: 2.78, 2.27 and 2.08 months.

Conclusions
Different models of the geomagnetic field of internal origin use different maximum degrees of the spherical harmonic expansion of the geomagnetic field components.Comparing the differences between monthly mean values of the geomagnetic field at a given place on the Earth and the respective values predicted by different models, we noticed that these differences have almost the same time dependency for the period considered (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).These time dependencies are shifted from each other and this shift can be removed by subtracting from these differences their averages.The resulting residuals representing the external contributions for the different models (EMM2015, POMME-9, CM5) at different observatories (WIK, NGK, KAK, HER) show some common features.They are almost in the same amplitude at all observatories (order of 20 nT for the X and Z components, order of 5 nT for the Y component) are almost the same at the nearly observatories (NGK, FUR, WIK).In particular the X component varies almost the same in all observatories.The time series of the Z component at observatories located in different hemispheres show inversion of the phase differences, whilst the time series of the Y component show much more noisy-like behavior.
The spectral analyses of these time series for different models, different observatories and different components evidenced some common frequency contents, with maximum of power spectrum for periods of 16 years, 1 year, and 0.5 year.These latest variations show the presence of an annual non-ionospheric and a seasonal modulation of Sq [Wardinski and Mandea 2006].
By varying the maximum degree for a given model of the internal field, the changes on the predictions at given place are seen.We noticed that increasing the maximum degree of the model does not go along with the best fitting of predicted values to the observed values of any geomagnetic field component at any place, i.e. the un-modeled internal or external contributions are not always decreased by increasing the maximum degree.However, we noticed that the shape of the plotted time series is not affected by the change of maxi-mum degree as it is always shifted a little bit up or down according to the considered component, model and maximum degree.By removing this shift, we can easily detect the external contribution in a given place that is not modeled by the models of internal field with different maximum degree.It is easy to notice that the time dependence is practically the same for all series.
We have not studied the correlation of these unmodeled external signals with DST-index or for other indices, such as F10.7, suggested by Wardinski and Holme [2011], that would help in specifying their external sources.Removal of these signals enhances the resolution of fine-scale detail in secular variation; this is useful in considering the phenomenology of geomagnetic jerks [Wardinski and Holme 2011].
, the "model biases" are presented, i.e. the differences of monthly means of X, Y, Z components observed at Niemegk observatory (NGK) during the common time interval 2000-2015 with respective values predicted by: a) EMM2015 model of maximum degree of 720 (dashed line), b) POMME-9 model of maximum degree of 133 (solid line) and c) CM5 model of maximum degree of 85 (dotted line).

Figure 2 .
Figure 2. Differences between monthly mean values observed at Niemegk observatory and respective values predicted by EMM2015 model (dashed line), POMME-9 model (solid line) and CM5 model (dotted line).

Figure 3 .
Figure 3. Differences between monthly mean values observed at FUR observatory and respective values predicted by EMM2015 model (dashed line), POMME-9 model (solid line) and CM5 model (dotted line).

Figure 4 .
Figure 4. Differences between monthly mean values observed at WIK observatory and respective values predicted by EMM2015 model (dashed line), POMME-9 model (solid line) and CM5 model (dotted line).
Dx i , we have: Considering the un-modeled part of internal contribution is almost unchanging during the considered time interval (16 years), then: and These series of residuals do not depend on the kindex, i.e. the series are principally the same for all models and should represent the non-constant external contribution to the given observatory.These series are calculated as residuals of average subtraction from differences between monthly mean values and respective values predicted by different models.These series for three different models EMM2015 (L = 720), POMME-9 (L = 133) and CM5 (L = 85), at NGK observatory are plotted in Figure 9.The graphs are almost the same for all models and they coincide for each component X, Y, Z.That means such series have the time variation contributions of the external sources at the NGK observatory for the period 2000-

Figure 9 .
Figure 9.The external contribution to the NGK observatory, calculated as shown in the text for three models: POMME-9 (solid line), EMM2015 (dashed line) and CM5 (dotted line).

Figure 10 .
Figure 10.The external contribution to the FUR observatory, calculated as shown in the text for three models: POMME-9 (solid line), EMM2015 (dashed line) and CM5 (dotted line).

Figure 11 .
Figure 11.The external contribution to the HER observatory, calculated as shown in the text for three models: POMME-9 (solid line), EMM2015 (dashed line) and CM5 (dotted line).

Figure 12 .
Figure 12.The external contribution to the KAK observatory, calculated as shown in the text for three models: POMME-9 (solid line), EMM2015 (dashed line) and CM5 (dotted line).