A geomagnetic reference model for Albania , Southern Italy and the Ionian Sea from 1990 to 2005

Taking advantage of the measurements undertaken during the Albanian and Italian magnetic repeat station networks since 1990, as well as of a selected set of Ørsted satellite total field measurements, a magnetic reference model for the region comprising the Albanian territory, the southern part of the Italian Peninsula, and the Ionian Sea is presented. The model, designed to model the components of the main geomagnetic field for epochs between 1990 and 2005, has been developed by means of spherical cap harmonic analysis applied to a cap of semiangle 8°, larger than that investigated to take into account the appropriate spatial wavelength content of the main geomagnetic field over the region. The goodness of the fit to the real data suggests that the model can be used as a reference model to reduce magnetic surveys developed in the area during the time of validity of the model. Mailing address: Dr. Luis R. Gaya-Piqué, Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143 Roma, Italy; e-mail: pique@ingv.it


Introduction
The importance of a magnetic repeat station network, where the three-component magnetic field can be measured at intervals of some years, is enhanced in those situations in which no magnetic observatories are operating in the near neighbourhood.The values of the magnetic components, and especially the secular variation field deduced from these measurements, can be used to define regional models to produce magnetic charts or reduce magnetic surveys to a given epoch.The Balkan Peninsula is one of the European regions with the poorest coverage of magnetic observatories, so the need for a regular monitoring of the magnetic field at some given points is fundamental.
A scalar survey (measurements of total intensity F only) was undertaken in epoch 1990.0 by the ex-Geophysical Enterprise of Tirana that covered all Albanian territory.The normal field model derived from those data was presented by Duka and Bushati (1991).In the framework of a joint project between the Center of Geochemistry and Geophysics of Tirana, the Physics Department of Tirana University, and the National Institute for Geophysics and Volcanology of Italy, a vector magnetic survey (F as well as declination D and inclination I were measured: Chiappini et al., 1997) covering the Albanian territory was carried out during September 1994September (epoch 1994.75).75).More recently, in August 2003(epoch 2003.6) a restricted scalar survey (6 points at repeat stations in the Albanian geomagnetic network and two new points) was carried out during the short visit of one of the authors (ADS) to Albania, funded by the Italian Foreign Ministry.
In this paper we take advantage of all these datasets, together with contemporary measurements from the Italian repeat station network and with a selected total intensity data set from the Ørsted satellite mission to develop a regional reference model for the considered area, i.e.Albania, the southeastern part of the Italian Peninsula, and the nearby Ionian Sea.The model, developed by means of Spherical Cap Harmonic Analysis (SCHA) with polynomial time dependency, represents the values of the magnetic field better (with smaller scale) than global models such as the International Geomagnetic Reference Field (IGRF), and is, therefore, of utility in reducing magnetic surveys carried out in the same region.

Data used
Scalar values at 1990.0 and 2003.6 from 31 and 8 Albanian stations, respectively, together with vector and scalar measurements reduced to epoch 1995.0 from a total of ten Albanian re- peat stations were used to develop the reference model.Details on most of these stations can be found in Chiappini et al. (1997).
In order to define a model not only for the Albanian territory but also for the Ionian Sea and the southeastern part of Italy, measurements at seven locations from the Italian magnetic repeat station network (Coticchia et al., 2001) at epochs 1990.0,1995.0, and 2000.0 were also taken into account.A comparison of the real data with the values predicted by the IGRF showed that two repeat stations, one Italian (Madonna di Servigliano) and one Albanian (Rubik), had such high bias values with respect to the reference model, that we finally decided not to include them in the inversion to produce the model.Figure 1a shows the spatial distribution of the Italian and Albanian stations used.
To obtain proper temporal behaviour of the model and improve the stability of the inversion, considering the lack of vector data surveys at different epochs in the Albanian territory, we decided to synthesize X, Y, and Z components at 1990.0 and 2000.0 for the same 9 Albanian stations reducing the real vector measurements of 1995.0 to those two epochs, using the secular variation predicted by the Italian geomagnetic reference model (ITGRF;De Santis et al., 2003), a model that has been demonstrated to predict the temporal change of the internal magnetic field better than the IGRF.
A selected subset of total intensity magnetic field measurements from the Ørsted satellite was also included to develop the model.The inclusion contributes to a homogeneous coverage of the studied region, especially in sea areas.Moreover, the satellite data act as a boundary to develop a three-dimensional model, valid not only at sea level but also at any distance between the surface and the satellite height.A total of 30 scalar values measured between 1999.5 and 2002.5 were selected according to restricted criteria to reduce the presence of external magnetic fields in the data.In this way, only nightside data were selected (in order to reduce the effect of ionospheric fields) for periods in which the 3-hourly index Kp was lower than 1+.The absolute value of the Dst index was also limited to a maximum of 10 nT.The data chosen are distributed in a height range between 650 and 850 km above the Earth's surface.The spatial distribution of these satellite data is shown in fig.1b.
The ground and satellite data were weighted differently according to the reciprocal of the variance σ 2 .This variance represents the contributions to the measurements of the crustal field σc and the measurement error σm.The ground data (vector and scalar) were assumed to have σc =50 nT and σm=10 nT, whereas no crustal contribution was assigned to the satellite data, and σm was set equal to that of ground data.

Model
The reference model was obtained applying Spherical Cap Harmonic Analysis (SCHA; Haines, 1985).This choice represents an improvement with respect to the previously presented polynomial models for the region (Chiappini et al., 1997), due to the fact that a SCHA model allows for the computation of the field component values through expressions that satisfy Laplace's equation.Moreover, the radial variation of the magnetic field is implicitly described in the model, without need to assume a dipolar continuation of the field like as with polynomial models.Although a similar approach was used by Chiappini et al. (1999), the field was analysed for a fixed epoch only, without secular variation modelling.
The solution of Laplace's equation for the magnetic potential due to internal sources over a spherical cap in spherical coordinates (r, θ, φ) can be written as where the polynomial time dependency is included.The spherical cap harmonic coefficients g m k,q and h m k,q are those which determine the model.The number of coefficients depends on the maximum spatial and temporal indices of the expansion, K and Q respectively.The associated Legendre functions P m n k (m) (cosθ) that satisfy the boundary conditions (a zero of the potential or its derivative with respect to colatitude at the border of the cap; Haines, 1985) have integer order m but a generally non-integer degree nk(m), where k is the index used to order the different roots for a given order m.Legendre functions were computed using a procedure proposed by Olver and Smith (1983), since it seems to provide more reliable values than the original approach suggested by Haines (1985) when the cap is rather small (Thébault et al., 2002).
The magnetic components are obtained as the derivatives of eq.(3.1) in spherical coordinates, since the potential is non-observable.The fact that vector measurements are combined with total field measurements introduces a nonlinearity in the equations involved to obtain the coefficients of the model.To avoid this problem, a first order Taylor expansion of the total magnetic field intensity, as a square root function of the X, Y, and Z components, was used (Haines and Newitt, 1997).
After many tests, the parameters that defined the best model in terms of fit to the input data and spatial and temporal behaviour corresponded to K = 2 and Q = 1, covering the period between 1990.0 and 2005.0.The coefficients were obtained through a least squares regression procedure.Haines and Torta (1994) pro- posed a statistical estimation to include coefficients in the regression or remove them according to an F level.In our case this F level was fixed to 0 due to the small number of coefficients involved in the regression, so no coefficients were rejected.
The cap over which the model was originally defined had a semiangle of 3°.Nevertheless some problems arose when such small cap was considered, as typically happens when SCHA is applied to small caps.Figure 2 shows the behaviour of the magnetic field components at 1990.0 over the studied region.It is clear the presence of fictitious oscillations in the charts, especially for the Y component.This is due to the lack of presence of lower harmonics in the spherical cap harmonic expansion (3.1): the smaller the dimensions of the cap, the higher the degree of the Legendre functions (for a 3°cap, the minimum degree nk(m) involved was roughly 45).To avoid these prob-  lems, we decided to enlarge the cap over which to impose the boundary conditions up to 8°s emiangle, in order to cover the most significant harmonics (the minimum degree is approximately equal to 12 for this cap for a maximum spatial index K = 2).
The final model has a total of 18 coefficients (table I).The IGRF2000 model computed at 2000.0 was previously removed from the input data to avoid mathematical problems related to the combination of small and large values during the inversion procedure when determining the coefficients, and to be used as the first value in the iterative process to obtain the model.The value of the IGRF2000 computed for a given location at 2000.0 has therefore to be added to the values provided by the coefficients in table I to obtain the final field value.

Results and conclusions
Table II shows the goodness of the SCHA model in terms of fit to the input data.The spherical cap model improves the fit to the observatory X component and to the satellite total field with respect to IGRF by about 50%, whereas the improvement is around 30% for the Y, Z, and F ground elements.Figure 3 shows the difference    between the SCHA and the IGRF models for epoch 2000.0.This difference is between 20 and − 60 nT for all the components in the area from which data were available, being a bit higher as we move towards areas without data.
Once the validity of the model has been demonstrated, we can recommend its use as reference model over the investigated region, with particular application to the production of regional charts as those shown in fig. 4 for epoch 2000.0 at sea level, to reduce magnetic surveys developed in the area, or to study phenomena related to the Earth's internal magnetic field.For example, the temporal evolution of the isolines in the magnetic charts for the Y component shows a clear westward drift of the east component of the magnetic field (fig.5).The value of the drift deduced from these charts is about 0.6°per year, which is comparable to the values observed elsewhere in the Atlantic hemisphere (Barraclough and Malin, 1999).Other kinds of studies can be developed through the inspection of the secular variation of the magnetic field.The secular variation for all the magnetic elements for epoch 1995.0 (fig.6), obtained from the SCHA model as the field differences between epochs 1995.5 and 1994.5, confirms that the region under study presents low values for the temporal variation of the geomagnetic field for this period (see for example Gubbins, 1990).
All these results will be confirmed when new real vector data are available for the region.It is planned to make a vector survey over the Italian region as well as in Albania during the period 2004-2005, and the model presented in this paper will be a great help to know a priori the values one can expect to find during the surveys.

Fig
Fig. 1a,b.a) Location of the Albanian and Italian repeat stations from which vector and scalar (stars) and only scalar (squares) magnetic field data were used to develop the model; b) location of the satellite data used.

Fig. 2 .
Fig. 2. Maps (in nT) for X (top left), Y (top right), Z (bottom left), and F (bottom right) elements for epoch 1990.0 at sea level obtained from a SCHA model developed on a 3°semiangle cap.Clear large oscillations appear as symptom of lack of the proper spatial spectral content.

Fig. 3 .
Fig. 3. Difference (in nT) for X (top left), Y (top right), Z (bottom left), and F (bottom right) elements between SCHA and IGRF models at sea level for epoch 2000.0.

Fig. 4 .
Fig. 4. Maps (in nT) for X (top left), Y (top right), Z (bottom left), and F (bottom right) elements for epoch 2000.0 at sea level obtained from a SCHA model developed on an 8°semiangle cap.

Fig. 6 .
Fig. 6.Maps for the secular variation of X (top left), Y (top right), Z (bottom left), and F (bottom right) magnetic elements for epoch 1995.0 at sea level obtained from the SCHA model.Contour lines are shown at 1 nT/year interval.

Table I .
Coefficients of the Albanian GeomagneticReference Model developed by SCHA.

Table II .
RMS fit to input data of SCHA and IGRF models (in nT).