Time variation of the European gravity ﬁeld from superconducting gravimeters

SUMMARY Ground-based gravity observations have the potential to add useful information to the interpretation of data from the new satellite gravity missions (CHAMP, GRACE, GOCE). We examine 4.5 yr of data from eight superconducting gravimeters (SGs) associated with the Global Geodynamics Project (GGP), from 1997 to 2001, and simulate the time variations that might be seen by a satellite. Signals that are removed from the gravity data before spatial averaging are the solid Earth and ocean tides, a global atmospheric loading using a vertical perfect gas law for the atmosphere, International Earth Rotation Service (IERS) polar motion and instrument drift. The 1-d gravity residuals form the basis of an interpolated minimum curvature grid that we spatially average and analyse using both surface polynomials and empirical orthogonal functions (EOFs). A clear annual component is present that, if truly regional, should be easily detectable by a satellite such as GRACE. The signal is consistent with expected continental water storage, which provides some interest for the future comparison of ground and satellite data.


INTROD U C T I O N
The framework of this study lies within the Global Geodynamics Project (GGP) that collects data from a network of approximately 20 superconducting gravimeters (SGs) distributed worldwide . Within the western part of Europe, eight gravity stations were operating between 1997 July 1 and 2001 December 31, though not all simultaneously. Most of the stations are clustered in the central region, but Metsahovi, Finland, lies some distance to the north (Fig. 1). Stations in the GGP have been providing data for a number of years for two distinct research areas: one for local studies of gravity and the other for gravity changes associated with global phenomena. The former have targeted phenomena such as tidal variations (gravimetric delta factors), ocean loading (tidal and non-tidal), atmospheric loading (local and global) and hydrology (local effects only). Global gravity changes are associated for example with polar motion, the nearly diurnal resonance effect of the free core nutation, and the search for the Slichter triplet and core modes (see e.g. Hinderer & Crossley 2000;Meurers 2001).
Here, we consider the spatial combination of various gravity timeseries. Of the approximately 20 GGP stations, only two groups, those in Europe and in Japan, have enough stations to form a reasonable spatial average. Preliminary results of shorter data series (1-2 yr) from Europe have already been presented (Crossley & Hinderer 1999, 2002Crossley et al. 2003), the difference here is that we extend the former analyses to a data set of 4.5 yr and introduce new hydrology computations. Like the previous studies, this paper concentrates on the ground gravity field; we have yet to make a direct comparison between ground and satellite data like the one initiated by Neumeyer et al. (2003) between monthly samples of CHAMP and superconducting gravimeter (SG) data during one year (2001).

DATA PROC E S S I N G
The stations and their locations are given in Table 1 and Fig. 1. The time period is from the start of the GGP (1997 July 1) up to the end of 2001. Several station changes were made during this time so the coverage is not entirely homogenous (Crossley et al. 2003). All stations have been regularly reporting data to the GGP, except for Medicina (MC) whose data is available through the work of Zerbini et al. (2001). We report first that a histogram of the distance between station pairs (not shown) confirms that the distance range of 200-1000 km is well covered and that station Metsahovi (ME) extends the largest dimension to 2000 km. Our database is the GGP uncorrected 1-min data, available (except MC) through the International Centre for Earth Tides (ICET, http://www.oma.be/KSB-ORB/ICET/index.html).
We start by subtracting well-modeled components from each station, beginning with a theoretical local tide (including ocean tide loading) using local tidal gravimetric factors (δ, κ) obtained from independent tidal analyses. All waves are fitted up to monthly periods and nominal tidal factors (1.16, 0 • ) are used for waves of  semi-annual and longer periods. We also subtract a nominal atmospheric effect, using an admittance of −0.3 µGal mbar −1 , and geodetic polar motion (including annual and Chandler terms) effects from the IERS web site. These nominal models are perfectly adequate for the 15-d averages of the gravity field we will eventually interpret. A critical step is to remove the instrument drift and various offsets that occur in such data. Depending on the station and which version of the instrument is present, offsets can either be frequent in some instruments (5-10 yr −1 ) or relatively rare (as a result of severe lightening strikes, for example). Helium refills are easily identified and removed because their timing is well known. More problematic are gaps in the data of several days or more, the danger being that an offset could be ignored or wrongly corrected. Fortunately, most of the sites are monitored with absolute gravimeters that help greatly in identifying major offsets. Additionally, they also provide the most common route to amplitude calibration (see e.g. Imanishi et al. 2002). Phase calibration is now usually measured, rather than estimated (Richter & Wenzel 1991;van Camp et al. 2000); again the nominal calibration errors (0.5 per cent in amplitude and approximately 5 s in time delay) are sufficient for this study.
For each series, we simultaneously remove offsets and a linear drift function, because these effects interact with each other. As the secular trend of the gravity field is not as important as the temporal variations, we specify a zero mean value. To include such secular terms in longer studies, the use of absolute gravimeters would be essential to unravel the geophysical drift of the stations. Most instrument drifts are between 1 and 4 µGal yr −1 (Table 1), an impressively small amount for a relative gravimeter. The final residuals are shown in Fig. 2, in which the breaks in the series for stations Brussels (BE), Moxa (MO) and Wettzell (WE) are clear. Station WE is the only one occupied by two different instruments, the older one with a large negative drift (Crossley et al. 2003).

ATMO S P H E R I C L OA D I N G
Of the two most important environmental effects on a gravimeter, the atmosphere and hydrology, the former has received more widespread attention and has reached a more advanced stage of analysis. The local atmospheric loading and attraction, already corrected in the residuals using a single admittance, can be improved by considering three increasingly better treatments of loading and attraction (Boy et al. 2002): (i) global surface pressure; (ii) global surface pressure and temperature, with vertical structure modelled using the perfect gas law; and (iii) fully 3-D modeling using pressure, temperature and humidity data with altitude. In this study, we use method (ii), with National Centers for Environmental Predictions (NCEP) reanalysis data (Kalnay et al. 1996) and a simple vertical structure as used by Merriam (1992). For the processing of GRACE data, Boy & Chao (2002 have shown that method (iii) requires a numerically intensive computation using more recent and precise atmospheric data sets such as NCEP or European Centre for Medium-range Weather Forecasts (ECMWF) that may introduce differences in the seasonal zonal coefficients of approximately 10 per cent. This therefore will clearly be the preferred method for future treatments of regional and global gravity data.
We show in Fig. 3 the differences between local pressure loading and the differences introduced using methods (i) and (ii) above for station ME, which has the largest atmospheric loading of all the stations. Although at times the differences can reach 1-2 µGal over short intervals of several weeks, there are no long-term differences that would lead to a seasonal (particularly annual) signal. We expect this not to be true for case (iii), so eventually we need to implement this method for the GGP analysis. Fig. 4 shows the residuals from Fig. 2 combined together after the use of method (ii) and filtered to 1-d averages. It can be seen that there is considerable variability and any regional correlation between stations is masked by local effects.

SPATIA L AV E R A G I N G
We need to spatially average the data to bring out the large-scale regional variations. With only eight stations, it is not possible to estimate reliable global spherical harmonics of the gravity field, so we proceed more intuitively. First, we apply a robust interpolation method to grid the data to a rectangular latitude/longitude grid (ignoring the surface curvature), using a minimum curvature algorithm. This provides the smoothest possible surface going through all the data points and enables contour maps to be produced. Such  maps do not do any spatial averaging of the field and neighbouring stations with conflicting series (e.g. BE and MB) still show up as inconsistencies. To do spatial averaging, we proceed two ways, the first is a simple polynomial approach consistent with previous work and the second uses empirical orthogonal functions (EOFs).

Simple approach
We fit a 2-D polynomial function to the minimum curvature grid and take this surface as the smoothed daily average. Though   wavelength of approximately 500 km, which is suitable for the analysis here. We then resample the data at the station locations on the surface and decimate the results to 15-d averages (Fig. 5). It is seen that station ME is now anomalous, as expected, as a result of its distance from the other stations and the lack of neighbouring stations as a control. Otherwise, as anticipated, the remaining stations now exhibit more coherent fluctuations, not just over short periods, but over the whole record.
To further characterize these fluctuations, we take a simple arithmetic mean of all stations except ME, for each 15-d sample, and plot the mean field with its standard deviation, as in Fig. 6. We interpret this as a rough estimate of the variation of the mean central European gravity field; a formal estimate of the error varies between 0.5 and 1.5 µGal with an average of 1.1 µGal.

The EOF approach
We later recognized that the spatial averaging can be better treated using the empirical orthogonal function (EOF) approach, as is routinely done for oceanography and atmospheric data. The idea is to analyse the whole spatial-temporal data set at once, by seeking the singular values of the data. When ordered from largest to smallest, these singular values and their associated eigenvectors represent the principal components, or most significant modes in the data. Numerous references to the basic idea can be found; one of the most popular is Wilkes (1995). We were fortunate to find the useful web site of Pierce (2003) that contains a fully functional FORTRAN code that can be downloaded to do the analysis.
We here sketch only the outline of the method. The 1-d station data is decimated to 15-d samples and, for each time sample (110 in all), we construct the minimum curvature surface, as above. We then fill each column of a system matrix with these values, packed by rows (longitude), so that each 15-d map becomes a new column. The single value decomposition (SVD) of this matrix yields the singular values, each of which is associated with a 2-D eigenvector (a spatial map) and its time variation, the principal value (a 1-D series). Fig. 7 shows the results for the first eigenmode, which is to be interpreted as the dominant spatial pattern of the gravity field. The central stations are seen to be reasonably coherent, whereas stations BE and MB more often than not contradict each other. This we explain by noting that station BE is the least reliable station of the network, being one of the founding SG sites with an early model of the GWR full-size gravimeter and also being subject to high cultural noise as a result of its central Brussels location. The nearby station MB has a newer compact lower-noise model that is better sited in a mine. Station MC also has a distinct signature, i.e. is less coherent with the other stations, perhaps as a result of its location on the southern side of the Alps. It is important to note that the first mode explains 47 per cent of the total variance of the data and the first four modes together account for 97 per cent of the variance. Fig. 8(a) shows the first principal component (the time variation of the eigenmode in Fig. 7), compared to the polynomial averaged field from the simple polynomial approach. Both series show similar peaks during the winter months. The agreement is encouraging and demonstrates that the annual signal first obtained using the polynomial surface method is not an artefact of this method of processing.

HYDROL O G Y
We have not allowed for local hydrology, on the grounds that effects confined to one station should in principle be smoothed out when averaging over distances of several 100 km. This is not a very good assumption unless we have sufficient stations to do a proper average, but in any case a regional gravity field should respond better to coherent continental-size changes than local effects. It is important to emphasize that the common practice of using local water table changes to derive a hydrology admittance (e.g. Crossley et al. 1998;Virtanen 2001) will remove both a local effect and a continental-size effect at the same time. Unless the porosity is precisely known, which it never is, using an admittance based on the correlation between gravity and a local hydrology observation (e.g. the water table) can therefore remove exactly the continental signal we seek! This is analogous to removing long-period tides where one must use nominal tidal gravimetric factors and not fitted ones, at semi-annual periods and longer. A satellite will naturally do the ideal spatial averaging of all hydrological effects (see Rodell & Famiglietti 1999).

DISC U S S I O N A N D C O N C L U S I O N S
Returning to Fig. 8(a), we see that a clear annual signal is present with peak gravity values in the winter and low values in the midsummer (except for 2001). Before looking at other data, we speculated this fluctuation might be the result of the annual loss of soil moisture and thinning of the water table in the summer. To check this, we computed the climate-induced hydrology variations at the nine GGP stations using the model of Milly & Shmakin (2002a,b), adding together the effect of soil moisture and snow cover. The mean variation of the eight central stations is shown in Fig. 8(b) along with the first principal component, confirming that both signals are similar in phase, departing from each other only towards the end of the data. The amplitude of the predicted hydrology effect is ±3 µGal, compared (indirectly) with the mean gravity of approximately 1.5 µGal from the polynomial averaging in Fig. 8(a). Our hydrology signal is similar to that predicted earlier by Van Dam et al. (2001) who found annual changes of 2-3 µGal for the European GGP stations. At some of the other GGP stations (e.g. Bandung, Indonesia) the annual signals from continental hydrology can be as large as 10 µGal. The agreement in Fig. 8(b) is suggestive, but does not prove a causal relationship because there are other annual signals we have yet to consider. One is the atmospheric signal resulting from seasonal changes in air pressure and temperature [referred to as method (iii) above; see also Simon 2002]. Another is the changes in station elevation (as a result of a variety of causes) that alters indeed gravity and can be observed directly by GPS. This latter signal needs to be   removed before comparison with satellite data because a satellite does not respond to elevation changes, only to the total gravitational potential. We know of only one station (MC), where a sufficiently comprehensive geodetic analysis of gravity, GPS and hydrology has been made (Zerbini et al. 2001(Zerbini et al. , 2002Romagnoli et al. 2003). In that case, the annual GPS vertical motion led to a smaller gravity signal (using the free air gradient conversion factor) than that observed, so a significant annual term is still to be explained. We have shown that the variation of the ground gravity signal over Europe is not inconsistent with the expectations of regional hydrology, at least for the 4.5 yr of our study. Put another way, if the hydrology predictions in Fig. 8(b) are accurate, then such a signal must appear somewhere in local gravity. Ground SG stations, when carefully processed and combined, can yield a regional gravity field whose mean value can be estimated to perhaps 1 µGal or so, over approximately 500 km. This approaches the best GRACE estimates using 5 yr of data for the extreme case of a large hydrological signal (Wahr et al. 1998). In future, it remains to be seen how well the GGP data compares with the satellite determination of the European gravity field.

A C K N O W L E D G M E N T S
We are profoundly grateful to the various European SG station operators for sending their data to ICET for the GGP. Thanks to David Pierce for his web site and the EOF program. This research was supported by the French Research Ministry (Action Concertie Incitative, ACI, Observation de la Terre).