On long-term trends in European geomagnetic observatory biases

We investigated the European geomagnetic observatory biases over 42 years, considered as contributions of the crustal field, and generally assumed to be constant in time. To estimate these biases, we compared observatory annual means to predictions given by the continuous CM4 model, and to four other core field models for different epochs. Solar-cycle related external fields are clearly present in the residuals. Although well-known, no suitable model to minimise them exists. We found that an empirical approach, taking advantage of the homogeneity of the external influences in the European region, can minimise these influences. Their reduction is better than when the external field description included in the comprehensive CM4 model is used. At several locations clear long-term trends remain after subtraction of the core field and minimisation of external fields. We investigated whether they are due to an insufficient description of the core field secular variation by the CM4 model, or to changes in induced lithospheric fields.


Introduction
The geomagnetic field measured at any point on the Earth's surface is a combination of several magnetic contributions generated by various sources. These fields are superimposed on and interact with each other. More than 90% of the measured field is internal in origin and is generated in the Earth's outer fluid core. This part of geomagnetic field, known as core or main field, is due to electric currents sustained by a geodynamo. Also, internal in origin is the lithospheric (crustal) field generated by magnetized rocks of the crust. External in origin, the magnetospheric and ionospheric fields vary on much shorter time scales than the core field and may create magnetic disturbances as large as 10% of the geomagnetic field during magnetic storms. Other important sources are the fields induced by currents flowing with Earth's crust and upper mantle.
Magnetic observatories remain the primary source of data to estimate the temporal variation of the core field. The use of such measurements in core-field modeling, at global as well as regional scale, needs to take into account the possible contributions of the crustal components observatory data, known as crustal biases (see Langel and Hinze, 1998, and references therein). Not considering these can lead to errors of about 10% of the field for the large scales (Langel and Hinze, 1998). Two methods are generally used to minimize the crustal contributions when inverting observatory data in secular variation studies: one is to use only time derivative estimates from observatory data (Cain et al., 1983); another method is to estimate the possible contribution from the crustal fields. To determine the crustal bi-Copyright c The Society of Geomagnetism and Earth, Planetary and Space Sciences (SGEPSS); The Seismological Society of Japan; The Volcanological Society of Japan; The Geodetic Society of Japan; The Japanese Society for Planetary Sciences; TERRAPUB. ases two different methods have been used during the last years. One is to estimate the biases directly as additional unknowns fitting process when inverting observatory and satellite data for a spherical harmonic model (Langel et al., 1982). The other method is to compare the observatory annual mean values of the three field components for a given epoch, to the values predicted for that epoch by a model obtained from satellite data only (Gubbins and Bloxham, 1985;Mandea and Langlais, 2002).
The goal of our work is to study the comparability of crustal biases determined using different main field models and investigate whether temporal changes of biases are observable. Examinations of observatory biases during the MAGSAT and Ørsted missions suggest that the crustal field may have changed measurably over the 20 years at some locations (Mandea and Langlais, 2002;. We focus on the European region, because it is most densely covered by geomagnetic observatories. Observatory annual means, defined as being the average over all days of the year and all times of day, are generally assumed to mainly reflect the core field, although it is known that the external field contributions do not average out completely (Yukutake and Cain, 1987). Therefore, a part of our work was to test the significance of external field influences in annual means.
We use the approach of comparing the observatory data to existing models. An observatory annual mean value B obs at each observatory is represented as the vector sum: (1) B core is the core field contribution, which dominates magnetic field models up to around degree and order 13. Note, however, that around degrees 12 to 15 core and crustal field influences are mixed and cannot be separated. B crust is the crustal field at the location of the observatory. It may change significantly over a distance of a few kilometers. B err are other contributions, which contain influences from external field remaining in the annual means and any kind of measurement and data errors. A comparison of observatory data to B core , given by a main field model, yields the sum (B crust + B err ). As we cannot directly distinguish between these two contributions we mostly refer to the sum as observatory biases or residuals and we investigate them in detail in Section 4.
Our analysis was performed on annual mean values of the X (northward), Y (eastward) and Z (vertically downward) components from 46 European geomagnetic observatories and models based on data from the Magsat, Ørsted, CHAMP and SAC-C satellites. A general hypothesis is that the crustal contribution at any location is constant. If the crustal anomaly is purely remanent then this hypothesis is reasonable. However, as some fraction may be induced, we investigated the temporal changes of the residuals over the period from 1960 to 2001, covered as well by the European magnetic data and the CM4 model (Sabaka et al., 2004). We also studied the short-term variations observed in the residuals influenced by the external fields, linked to the solar activity, and minimize them by an empirical approach. We finally discuss possible sources of the remaining timevarying residuals.

Data
This study is based on a comparison of real observatory data to five synthetic datasets, all obtained for the European observatory locations. Synthetic data are calculated from the CM4, MAGSAT, ØRSTED, POMME and CHAOS models. These datasets are described in detail in the following.

Observatory data
The real data are the annual mean values of the X, Y and Z components from 46 European geomagnetic observatories, which are given for epochs xxxx.5, i.e. the middle of the year. We analysed data available over the time span 1960.5-2001.5, the period covered by the CM4 model, and for individual epochs according to the satellite models. Figure 1 shows the observatory spatial distribution. The list of observatories with their corresponding IAGA codes and both geographical and geomagnetic coordinates is given in Table 1. The geomagnetic coordinates refer to DGRF1980 (Definitive Geomagnetic Reference Field; http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html), this epoch being the middle of the studied time interval. The initial dataset was obtained from the World Data Center Edinburgh (http://www.geomag.bgs.ac.uk/gifs/ annual means.HTML).
Changes in location or instrumentation of observatories lead to changes in the absolute level of recordings, known as "jumps" (Jankowski and Sucksdorff, 1996). Generally, the jump values can be determined for the different components by comparative measurements. When such jumps have been reported observatory time series with the corresponding values, all data were adjusted to the level of the most recent epoch. All the time series were then carefully checked visually, by computing first differences and partly by comparing data of neighbouring locations.
We found some peculiarities in the datasets from PEG, MOS, TRO and SFS observatories. The time series for PEG observatory clearly show an unreported jump in 1987.5, both X and Y components. Moreover, there is a gap of 5 years between 1994.5 and 2000.5 and the remaining data, at least in X, look quite scattered. For these reasons we decided not to consider the PEG values after 1987.5. The time series of all three components at MOS observatory are more scattered than in other observatories, suggesting problems with technical noise or instrumentation at that site. However, as this scattering is not a systematic error, we kept all these data, being aware of their lower quality. At TRO observatory, we noticed a feature resembling a jump in the X component between 1971.5 and 1973.5. However, a similar behaviour with smaller magnitude can be found at the neighbouring KIR and SOD observatories. Although this feature is very sharp for a real magnetic field effect, the common occurrence at all three observatories suggests a natural source, and again we retained the data for our study (there is no way to determine whether this feature at TRO is amplified by overlapping of a natural effect with an artificial jump). For SFS, data prior to 1991.5 were not available at the WDC. Apart from the cases mentioned above, all available data from the other observatories are considered. However, not all time series cover the full time interval. For our analyses only observatories providing continuous data over a minimum of 9 years between 1960.5 and 2001.5 were kept. This amounts to 46 observatories shown in Fig. 1.
Some of the field models described in the following are centered on epochs xxxx.0 instead of xxxx.5 like the avail-able observatory annual means. In these cases we assumed linear secular variation to determine the observatory annual mean centered on epoch xxxx.0 from the surrounding years. Data to compute annual means for xxxx.0 epochs are not always easily available and the improvement in accuracy is considered to be too small to justify the effort.

Synthetic data
The datasets used to estimate the core field contributions at observatory locations were obtained from field models based on data from the satellites MAGSAT, ØRSTED, CHAMP and SAC-C, partly in combination with observatory data. As described below, several models are expanded to spherical harmonic degrees describing parts of the longwavelength lithospheric field and a few models contain ex-ternal field descriptions. We need the core field contributions for our study and used the internal field descriptions of each model up to spherical harmonic degree and order 14.
MAGSAT model This model (Cain, 1989) is based on the available MAGSAT vector data complemented by observatory secular variation results from September 1979 to June 1980 and is centered on 1980.0. This was the first attempt to include significant parts of the lithospheric field in global spherical harmonic models and it has been expanded to degree 63 with secular variation estimates up to degree and order 10. We used this model for its center epoch, 1980.0.
OSVM The Ørsted secular variation model, OSVM , is based on Ørsted scalar and vector data from March 1999 to September 2001 and observatory secular variation values. Satellite data are selected and corrected to minimise the influence of external fields. The model is expanded to spherical harmonic degree and order 29 with a secular variation description up to degree and order 13 and includes a description of the magnetospheric ring current. 2000.0 is the epoch on which the model is centered and for which we used it. POMME. We used version 3.0 of the POtsdam Magnetic Model of the Earth, POMME . It is based on CHAMP satellite vector and scalar data from 2000.6 to 2005.7, centered on 2003.0 with secular variation and acceleration described by a Taylor series expansion of the core field coefficients up to degree and order 16. The accuracy of the internal field description is improved compared to the earlier versions by the larger amount of available data and improved data selection and correction for external fields. The static field is expanded to spherical harmonic degree and order 60 and a magnetospheric field description is part of the model. We used this model for the central epoch, CHAOS model. This recently developed model  is based on CHAMP, Ørsted and SAC-C data measured between March 1999 and December 2005. CHAOS describes the core and crustal field up to degree and order 50 with a continuous spline representation of the coefficients up to degree and order 14 and a linear secular variation estimate for degrees 15 to 18. It is the first continuous model based on satellite data and we used it for epoch 2002.5.
CM4 The fourth version of the continuous 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 spherical harmonic degree and order 60. The model also includes descriptions of various other field contributions originating in ionosphere and magnetosphere and from induction Earth by the external source fields. We want to use this model to investigate any temporal changes in observatory biases and calculated synthetic time series for the X, Y and Z components at each observatory location for the whole time span.

Crustal Biases as Determined from Different Models
First we compared the crustal biases at all observatories that result from the subtraction of the five different core field models described above. When using CM4 model, the average of the bias for each observatory over the 42 years is used for this comparison, which is displayed in Fig. 2 (left side). In the X component, the biases are generally small and the only ones larger than ±200 nT are found at KIR, related to the well-known Kiruna magnetic anomaly, and at LRV, NUR, MNK. In the Y component, four European observatories have biases larger than ±200 nT, namely LRV, TRO, ODE and again KIR with the largest average bias. There are also a number of observatories with very low Y biases. For the Z component larger biases are more common. The maximum residual in this component appears at SOD, much higher than at KIR. The panels at the right side of Fig. 2 show the deviations of all satellite based biases from the CM4 ones. In general, biases obtained using different models for all three components are in a reasonable agreement with typical differences less than 30 nT. However, biases with low amplitudes sometimes do not agree in sign. The differences between biases based on different models are partly due to the capability of a particular model to represent the core field, and can be influenced twofold by the fact that they are obtained for different times: the observatory data are influenced differently by external field contamination and measurement errors at different years, and the lithospheric field may indeed change due to induction effects. Taking this into account, we would expect close agreement between the MAGSAT and CM4 based biases, centered closely to 1980 by the averaging, and between the results based on the other three satellite models. For the results obtained for epochs about 20 years apart we might expect larger differences. This cannot be said to be true, except possibly in the Z component.
For further comparison we plotted the observatory locations on the map of long-wavelength crustal anomalies of all three components at the Earth's surface obtained from the MF4x model (Lesur and Maus, 2006) in Fig. 3. This lithospheric field model is derived from almost 5 years of CHAMP measurements (2000)(2001)(2002)(2003)(2004)(2005) and spherical harmonic expansion up to degree 90 at low latitudes, but only degree 60 polar regions. Although this satellite altitude model is more suitable for downward-continuation to the Earth's surface than MF4 , it still has to be regarded with some caution, as any noise gets strongly amplified by the downward continuation. We also do not expect a perfect agreement because local, short-wavelength anomalies are not detected by magnetic satellite data. Unfortunately, a compilation of detailed anomalies from aeromagnetic surveys is not yet available for the whole area of interest.
The agreement between the averaged observatory biases and the anomalies on the maps was found to be low. The relative amplitudes compared to each other mostly do not agree and even the signs of the computed crustal biases and of the anomaly, as shown on the map, do not agree for several locations. Clearly small-scale anomalies are dominating the observatory biases.   based on the CM4 model (left) and differences between these and biases determined using four other magnetic core field models for different epochs (right), see text. The observatories are ordered by geomagnetic latitude, see Table 1. In the differences, zero in general indicates that no observatory annual mean data was available for the respective epoch.

Temporal Evolution of Biases and External Field Influences
For a closer analysis of the biases we calculated the residuals between observatory annual means and CM4 model predictions for each observatory and each year from 1960 to 2001. Figure 4 shows some representative cases. As an example of the most common behaviour, we show the time series for WNG observatory (Fig. 4(a)). All components biases are nearly constant, however showing a quite regular variation, within some 20 nT. In Fig. 4(b), HRB observatory is shown, representing a case in which the residuals are neither constant, nor showing a clear trend in X. Figure 4(c) shows residuals evolution at MOS where large changes in the biases are observed. For example, the X component residuals are nearly constant until 1983, thereafter show a change of about 100 nT. The most outstanding trend is found in Z, for which during 42 years, the bias changed by about 150 nT. In Fig. 4(d), TRO observatory, an apparent jump X component around 1971 (see Section 2), can be seen. In Z we again notice some trend changing the bias by more than 100 nT from 1960 to 1977 with a sim-ilar jump-like feature of smaller amplitude, around 1971. Finally, Fig. 4(e) shows KIR observatory biases with the strongest constant bias in X and Y, with slight but noticeable trends in both components. Note also, that several of the small variations seen in Fig. 4(a) can also be identified in Figs. 4(b)-(e) in addition to the already described strong effects.
Several contributions can be invoked to explain observed differences: external field influences in the observatory annual means, induction effects in the lithosphere, insufficient representation of secular variation in the CM4 model, and finally instrument drifts and data errors. To identify common features and find possible regional effects we display the residuals in form of a colour-coded matrix for each component in Fig. 5. The observatory residuals (lines of the matrices) are ordered by geomagnetic coordinates, with north at the top. The averaged residuals, which are the assumed constant crustal contributions (see Fig. 2) were subtracted from each time series to study any time dependent contributions on a comparable colour-scale. Even after subtracting the corresponding averaged residuals, a few observatories remain with outstanding biases due to strong, continuous trends (e.g. MOS, TRO). In order to underline any patterns on short time scales, the colour-scale was limited to ±20 nT, the stronger trends appearing as changes from dark blue to dark red or vice versa, without information about the real magnitude. Significant differences are seen between the components, on both short and long time scales.

External contributions
The plots shown in Fig. 5 reveal quite a similar pattern for most of the short term variations of the residuals at all observatories, which we suppose to be due to external fields still present in the annual means. Most of these variations are also superimposed on the long-term changes in the observatory crustal biases, as shown in Fig. 4. They cannot be noticed in the colour-coded matrices due to the colourscale, chosen to clearly see the short-term scale patterns. A variation pattern is most prominent in the X component with vertical stripes of maxima and minima correlated with the solar cycle (see Fig. 5(d)), a well-known influence in annual means (Alldredge, 1976;Courtillot and Le Mouël, 1976). Generally, the years immediately after the solar maximum are characterised by a decreasing number of geomagnetic storms. For these specific years, the annual mean values of X are lower than the average, as each storm reduces the field values of this component for a few hours to a couple of days at mid-latitudes. This effect is obvious in 1981-1982 and 1989-1991. For the Y component a similar, but much weaker pattern with opposite sign is observed (i.e. in the years after solar maxima the residuals increase). The Z component seems to be more influenced by effects of even shorter time scale. Although the recognition is hampered by the high number of time series with long-term trends we can detect some vertical stripes, reflecting common features to all observatories. These relative maxima and minima in time often only are of annual or bi-annual duration and are not clearly linked to the solar cycle. For Y and Z, no dependence of the variations on geomagnetic latitude was found. In X, however, the maxima and minima amplitudes are mostly observed in the southern European observatories.

Long-term trends
The long-term trends observed in Fig. 5 could be explained by an insufficient description of regional secular variation by CM4, changes in the induced crustal magnetisation, or long-term induction effects in conductivity anomalies in the lithosphere. Long-term trends are mainly seen in Z (ISK, SPT, PAG, COI, THY, HRB, MOS, MNK, MAB, TRO), and only a few in X (THY, LVV, MOS, TRO, NCK), and Y (MOS, COI, PEG). Most of the observed trends in the Z component represent a decreasing temporal bias, with only one exception, SPT. In X and Y there are no clear preferences for decreasing or increasing trends. Before discussing whether the changing biases can be linked to the observed lithospheric anomalies we try to separate the trends from the overlying external field influences.

Reduction of the External Field
In order to estimate and reduce the remaining external influence in the observatory annual means, our first approach was to take into account external and induced contributions provided by CM4, modulated by storm-time-disturbance (Dst) and Solar flux (F10.7) indices. Figure 6 shows the obtained residuals at NGK after removing different contributions (core field, core and external fields, and finally core, external and induced fields). The model predictions includ-a) 1960 1965 1970 1975 1980 1985 1990 1995  ing external fields were averaged over all hours and days of the year. In the X and Y components a clear reduction of external influence is achieved, but in the Z a significant level of variability remains. The limited performance of CM4 for this component could be attributed to the very simple conductivity model which does not take into account regional or local induction effects. The induced part of the external field variations is strongest in the vertical component. In X and Y it hardly makes a difference whether the induced field contributions described by the CM4 model are considered or not. The remaining variability data series after subtraction of all the signal described by CM4 is still order of 8 nT in X and Z, and 4 nT in Y. Similar orders of magnitudes are obtained at other observatory locations. Additionally, we analysed whether quiet day annual means better describe the internal field than all-days annual means. The observatory quiet day annual means are obtained as the average of the five most quiet days of each month over a year. One could assume that they are significantly less influenced by external fields. The dashed pink line in Fig. 6, however, shows that after subtracting the CM4 core field, a similar variation pattern as all day annual means remains, with only slightly lower amplitude remains.
To minimise the persisting external and induced fields, we tested a very simple approach based on the clearly homogeneous variation pattern for the whole region.  Table 1. d) Annually averaged Sun spot numbers.
We took advantage of the relative homogeneity of the external field influences in a small region such as Europe, reflected in the observed similarity of the short period signal in all observatories. First, it was necessary to detect and  Fig. 6. X, Y and Z residuals at NGK after subtracting: core field (green); core and external fields (red); core, external and induced fields (blue) as given by the CM4 model. Th black lines shows the residuals after subtraction of the empirical template (see Section 5). The dashed lines represent quiet-days residuals after subtracting core field. exclude those observatories with temporal patterns significantly different from the generally observed. Our approach consists of the following steps: i) We computed the standard deviations of the residuals for each component and each observatory (SD x , SD y , SD z ) as a measure for scattering of the individual time series. ii) We estimated the cumulative distribution function for SD x , SD y , SD z . Like the example in Fig. 7, all three graphs show that a large number of observatories have small scattering values compared to a few others. iii) We estimated the level at which the cumulative function is constant for at least two values of SD x , SD y , SD z . The percentage and number of observatories below those levels, as well as corresponding standard deviations are listed in Table 2. Note that column 4 does not include the same observatories in all components. Choosing those common to at least two components reduces the number to 36 observatories. These 36 observatories were used to construct templates, for each component separately, which are removed from the residuals. The templates were obtained by a median averaging of the 36 observatory residuals in each year. We subtracted the obtained template from all observatory residuals. The resulting colour-coded matrices are presented in Fig. 8. When comparing this figure with the matrices in Fig. 5, the largest change can be noticed in the X component. Here, it is already seen that the prominent stripes linked to the solar cycle are not present any more with the exception of a few southern observatories, where still some influence of the external field remains. The matrix related to Y changed slightly because the pattern in the original matrix was already smoother. In the Z matrix it is obvious that the previously observed annual and bi-annual stripes have disappeared.
To check the success of our procedure, the residuals time series were plotted before and after removing the templates, for each component. The resulting corrected time series for NGK is included as the black line in Fig. 6. It becomes clear that a better reduction of the external field variations in the annual means is achieved than by the CM4 external and induced field description. As further examples, in Fig. 9, the time series together with the template for the observatories ISK and WNG are shown, the first being an example for the southern and the second one for the northern observatories. In the original data for both ISK and WNG, the solar minima and maxima are identified in the X component. Subtracting the template, we clearly reduced the solar cycle contributions. It is obvious that a better reduction was achieved for the northern observatory than for the southern one. As already noted, applying the above procedure has    Fig. 9. Residuals in X, Y and Z components prior to (green) and after (red) applying the reduction procedure for the external field contributions, along with the corresponding templates (blue) at ISK (first column) WNG (second column) observatories (see text for details). Note the different scales.
we were not able to reach the same improvement due to the data quality: only since 1965 a proton magnetometer has been used for F measurements. Moreover, since 1969 a vector proton magnetometer for H , Z and F replaced the old H determination with the Gauss method. These two instrumental changes dramatically improved the data quality at WNG '70s. The WNG example is not unique, and for further studies we should keep in mind that new instruments and methodology in measuring the magnetic field continuously improved the data quality.

Discussion and Conclusions
In this study we analysed the observatory crustal biases of 46 European observatories as obtained using five different geomagnetic core field models. We focussed particularly on any temporal changes in the residuals between real annual mean data and CM4 model predictions. Apart from con-stant offsets between data and models, which represent permanent crustal field signature, short-period variations of external origin in the order of ±10 nT are present in all the annual means. For a more suitable use of observatory annual means for core and lithospheric field studies, we utilised an empirical procedure to minimise the external influences in these data. Although their presence is well-known, no dedicated model to remove them exists. We found that for a limited region, such as Europe, the subtractions of a simple template, making use of the homogeneity of the external influences, gives better results than the external field description included in the CM4 model. This is not surprising because even with the modulation by two indices the CM4 model is mainly a quiet time field description. Solar cycle related influence is more obvious in observatory residuals at the lower geomagnetic latitudes. That is not exactly what could be expected, as strong variations caused by geomag-netic storms characterise higher latitudes. The observed latitudinal differences could be linked to ring current effects.
One of the most striking results of this study is that longterm trends were detected in many biases series. These trends often are in the order of 2 to 5 nT/year. They are unlikely to originate from instrumental drift or data errors, being systematic over such long time intervals. The trends in biases are not grouped in a specific region, but distributed all over the investigated area. Only in a few cases do neighbouring observatories show the same crustal biases trends. Moreover, as long-term changes in biases are more often found in Z than the other components, it seems likely they are linked to changes in crustal magnetisation or induction effects in lithospheric conductivity anomalies, over some decades. Lesur and Gubbins (2000) showed that substantial induced magnetisation exists around almost all the 20 European observatories considered by them. However, they concluded that a satisfactory separation between induced and remanent parts is difficult to be obtained and our study provides first indications where the induced parts may be particularly strong.
The comparison between the anomaly maps obtained from the MF4x model (Lesur and Maus, 2006) and averaged observatory biases revealed that in most cases small scale anomalies cause the observed biases. Differences of a few nT up to 50 nT occur when different main field models are used to determine the biases. This is due to differences core field descriptions of the models and external field influences data. Some of the observatories with temporal trends in the Z biases are located on rather high anomalies (e.g. MNK, LVV), while others (e.g. COI, SPT) are not. The same is true for the other components. On the other hand, some observatories without temporal change in the biases are located on strong anomalies. Looking again at the comparison to large-scale lithospheric anomalies in Fig. 3, we get a similar result: the temporal changes in the biases are not correlated to the strength of the corresponding regional anomaly field. On the contrary, there are cases where only one of two observatories located on the same large-scale anomaly shows some temporal change, as for the Z component of MOS and BOX.
The time-varying biases could provide information about the induced contribution of magnetic anomalies located under the respective observatories. A detailed description of the small-scale magnetic anomalies for the whole of Europe does not yet exist, although efforts to achieve this are going on (World Digital Magnetic Anomaly Map, WDMAM; see http://projects.gtk.fi/WDMAM/index.html). Better understanding of the distribution of these time-varying features, not described by the CM4 field model, might be gained by regional field modelling, and by considering repeat station data to determine their spatial structure.