Continuous GPS measurements of postglacial adjustment in Fennoscandia 1 . Geodetic results

[1] Project BIFROST (Baseline Inferences for Fennoscandian Rebound Observations, Sea-level, and Tectonics) combines networks of continuously operating GPS receivers in Sweden and Finland to measure ongoing crustal deformation due to glacial isostatic adjustment (GIA). We present an analysis of data collected between August 1993 and May 2000. We compare the GPS determinations of three-dimensional crustal motion to predictions calculated using the high-resolution Fennoscandian deglaciation model recently proposed by Lambeck et al. [1998a, 1998b]. We find that the maximum observed uplift rate ( 10 mm yr ) and the maximum predicted uplift rate agree to better than 1 mm yr . The patterns of uplift also agree quite well, although significant systematic differences are evident. The root-mean-square residual rate for a linear error model yields estimates of rate accuracy of 0.4 mm yr 1 for east, 0.3 mm yr 1 for north, and 1.3 mm yr 1 for up; these figures incorporate model errors, however. We have also compared the values for the observed radial deformation rates to those based on sea level rates from Baltic tide gauges. The observational error for the vertical GPS rates required to give a reduced c of unity is 0.8 mm yr . The time series do exhibit temporal variations at seasonal frequencies, as well as apparent low-frequency noise. An empirical orthogonal function analysis indicates that the temporal variations are highly correlated among the sites. The correlation appears to be regional and falls off only slightly with distance. Some of this correlated noise is associated with snow accumulation on the antennas or, for those antennas with radomes, on the radomes. This problem has caused us to modify the radomes used several times, leading to one of our more significant sources of uncertainty.


Introduction
[2] The last 800 kyr of the current ice age have been characterized by a series of ''glacial cycles,'' each with a period of $100 kyr [e.g., Broeker and van Donk, 1971].
Within each cycle a relatively slow glaciation phase, culminating in massive ice complexes over most of the highlatitude continental regions, was followed by a much more rapid deglaciation event.For example, during the Last Glacial Maximum, which occurred just $20 kyr ago, the ice sheets reached thicknesses of 2 -3 km or more in Fennoscandia, Canada, Antarctica, Greenland, Siberia, and Arctic Canada [Denton and Hughes, 1981].Remarkably, in a matter of only 15 kyr a large proportion of ice disintegrated, leaving most of these regions ice-free and raising worldwide ocean levels by over 100 m [e.g., Chappell and Shackleton, 1986].The redistribution of surface ice water mass implied by these glaciation/deglaciation episodes has induced an appreciable and ongoing isostatic adjustment of the planet.
[3] Glacial isostatic adjustment (GIA) is manifested in a wide variety of past and present-day geophysical observables that have previously been used to study this process, including time series of ancient sea level elevations (relative to present-day sea level), the modern tide gauge record, gravity anomalies, present-day secular variations in the global gravity field, and present-day secular variations in the Earth's rotational state (variations in length of day and motion of the rotation pole) [e.g., Milne, 1998].These observations provide an indirect inference of present-day ongoing crustal deformation.Direct high-accuracy measurements of crustal deformation, even in a particular region, were not possible, however, before the advent of space geodetic techniques, and even these techniques have only quite recently been able to achieve the required accuracy.Very long baseline interferometry (VLBI) has been available as a high-accuracy geodetic technique for over 20 years; several studies have now addressed the effects of GIA on VLBI determinations of site velocities [e.g., James and Lambert, 1993;Mitrovica et al., 1993].Unfortunately, the global VLBI network is extremely sparse, and only 2 -3 sites exhibit much sensitivity to GIA.A small error in the velocity determination of these sites will unduly influence conclusions regarding mantle viscosity and ice-sheet histories [Mitrovica et al., 1993[Mitrovica et al., , 1994b].
[4] Geodesy with the Global Positioning System (GPS) affords us several advantages relative to geodesy with VLBI related to the relatively low cost of the GPS receivers ($$15,000 or less).As a result of this low cost, it is feasible to deploy a dense network of receivers across a region.The detailed pattern of deformation may thus be inferred.The low cost of GPS receivers moreover makes it financially feasible to dedicate a group of GPS receivers to permanent sites within a region in order to acquire continuous measurements.In principle, given sources of error that are sufficiently steady state, it should be possible to reduce the noise and thus determine estimates of velocity much more quickly and accurately than with conventional campaign GPS measurements.The number of GPS receivers in permanent GPS networks is quickly outpacing the number used for conventional campaign high-accuracy geodetic measurements [e.g., Segall and Davis, 1997].[5] Geodesy with GPS has been steadily and significantly improving in precision and accuracy over the past 10 years.This improvement has mainly been due to advances in the GPS satellite constellation, GPS receiver design, and analysis techniques.The demonstrated repeatability of horizontal position estimates obtained from GPS is currently at the few-millimeter level on regional and local scales and at the 10-mm level on global scales [e.g., Blewitt, 1993].The repeatability in the vertical baseline component is typically 3-5 times worse.The level of accuracy achievable in a single day with the GPS technique is thus, in principle, equal to that achievable with VLBI.
[6] In August 1993 we established a network of permanently operating GPS receivers in Sweden [BIFROST Project, 1996].A number of sites in Finland were also temporarily occupied.In 1994In -1996, permanent GPS sites were established in Finland.The GPS sites within these networks, along with several already existing sites of the International GPS Network for Geodynamics (IGS) sites in Norway operated by the Norwegian Mapping Authority, make up a dense regional Fennoscandian GPS network of an intersite spacing of $100 km and a total area of $2000 Â 2000 km, situated within the area covered by ice at the Last Glacial Maximum.Investigators from three Nordic and two North American institutions formed Project BIFROST (Baseline Inferences for Fennoscandian Rebound, Sea level, and Tectonics) [BIFROST Project, 1996].One of the primary goals of BIFROST is to use the three-dimensional velocity vectors from the BIFROST GPS network to provide a new GIA observable for the determination of Earth structure and Fennoscandian ice history.In this paper, we report on the first results from this effort.We will present the BIFROST GPS networks and data sets and describe the analysis of the GPS data.We include a discussion of errors since the GIA signals that we are attempting to measure are quite small (subcentimeter per year).

The BIFROST GPS Networks
[7] The BIFROST GPS networks (Figure 1) are composed of the permanent GPS networks of Sweden (SWE-POS 2 ) and Finland (FinnRef 2 ).Table 1 describes the histories of the BIFROST sites.Table 1 contains the dates and configurations of the original installations, as well as  [8] Next, we discuss the individual aspects of both BIFROST networks.

The SWEPOS GPS Network
[9] The Swedish nationwide multipurpose network of 21 permanent GPS stations, SWEPOS, was established in 1993.On 1 July 1998, SWEPOS attained full operational capability for real-time positioning at the meter accuracy and for postprocessing applications with centimeter accuracy.Real-time positioning at the centimeter/decimeter level is planned for 2002.The National Land Survey of Sweden (Lantma ¨teriverket, or LMV) is responsible for the maintenance and the operation of the SWEPOS network.
[10] The SWEPOS network (Figure 1) currently consists of 21 continuously operating GPS stations.The ''standard'' SWEPOS monument (designed by the LMV and denoted as ''SWEPOS'' in Table 1) consists of a 3-m-tall concrete circular pillar atop a concrete platform.At five sites (Kiruna, Lovo ¨, Ma ˚rtsbo, Norko ¨ping, and Skelleftea ˚) a second pillar is available to serve as an alternate and as a platform for test measurements.The pillars are built on bedrock with clear lines of sight from the tops of the pillars down to elevation angles of 10°and often lower.Each pillar is supported by four internal steel rods set 1 m into the underlying rock.Heating coils are wound helically around each concrete pillar.Insulating material consisting of helically wound corrugated plastic sheet and rock wool sur-rounds the wire-wrapped pillar.A temperature sensor is fit into a small cavity inside the pillar and is connected to a thermostat unit in the instrument cabin.The thermostat maintains the temperature of the sensor above 15°C.On the top of each pillar is a plate for the attachment of the GPS antenna, tribrach, and adaptor.
[11] The pillars at each site are surrounded by a network of steel pins, driven into the rock so that their tops protrude a few centimeters above the surface.This local network, covering an area of $15 m Â 15 m, is used to monitor the stability of the concrete pillars.The GPS antenna is removed from the pillar and replaced with a theodolite, which is used to measure the horizontal and vertical angles to the steel pins.Through resection, the position of the pillar can be calculated.In this manner, the local position and orientation of the pillar may be monitored to better than 1 mm.The first such measurements were obtained during summer 1993 and repeated annually except for the monument in Leksand where measurements are carried out monthly.The results of these measurements are described in section 5.4.
[12] Several sites have slight variations in monumentation.The Onsala site has a different monument due to its earlier construction as an IGS site.The Onsala monument consists of a 1-m-tall pillar with a square cross section and without heating control or insulating material.The Jo ¨nko ¨ping pillar is 1 m shorter than the standard pillar (for air traffic safety) and also is not heated.The Lovo ¨and Ma ˚rtsbo monuments are SWEPOS monuments built over preexisting pillars of rectangular cross section.In addition, the multipath environment at these two stations might be worse than at others due to preexisting construction.
[13] Each SWEPOS site has a 3 m Â 2 m hut housing the GPS receivers, backup batteries, computer, and Internet connection.All SWEPOS sites are equipped with two or more GPS receiver systems, AOA SNR-8000 and an Ashtech Z-XII, connected to a single Dorne-Margolin-type antenna.At four stations, dual-frequency GPS/Global Navigation Satellite System (GLONASS) receivers are also installed.The stations are equipped with a power backup system, which can run the station for 48 hours if the main power fails.All stations are connected to the control center via leased 64 kbit lines and a redundant 19.2 kbit X.25 line.

The FinnRef GPS Network
[14] Planning for the FinnRef network (Figure 1) started at the Finnish Geodetic Institute (FGI) at the end of 1992, when it was decided that a network of 12 stations would be established.Candidate sites were chosen with several criteria in mind: (1) the network should cover the country so that the maximum land uplift differences could be sampled; (2) the stations should be built on bedrock and there should be open sky above an elevation angle of 15°; (3) absolute gravity has been or can be measured on the spot; and (4) stations should easily be connected to the precise leveling network and to the telephone and electricity networks.Criterion 2 has generally, but not universally, been met, and in most cases the horizon is 10°or lower.Planning, construction, and use of the network are described in more detail elsewhere [e.g., Koivula et al., 1998Koivula et al., , 2002]].
[15] At the Joensuu, Kuusamo, Vaasa, Virolahti, Olkiluoto, Kivetty, and Romuvaara stations we constructed heated wooden cabins of area 1.5-2 m Â 3 m to house the GPS and other electronics.Existing buildings were used at all other sites.Some stations are located close to other institutions where on-site personnel can assist in case of minor problems.Kevo is on the premises of the Subarctic Research Center of the University of Turku, Oulu is at the Aarne Karjalainen Observatory of the University of Oulu; Sodankyla ¨is visited weekly by local staff of the Sodankyla ¨Geophysical Observatory; and Tuorla is at the Astronomical Observatory of the University of Turku.At Kivetty, Olkiluoto, and Romuvaara, there are also contact individuals who can check the stations.Metsa ¨hovi is at the Space Geodetic Observatory of the FGI.
[16] Three different types of antenna platforms are used for FinnRef.The standard configuration is a 2.5-m-high steel grid mast, which is used at Joensuu, Kuusamo, Sodankyla ¨, Tuorla, Vaasa, and Virolahti.A similar mast, but 5 m high, is used at Kevo.In the case of the 2.5-m mast the thermal expansion effects amount to a height variation of <±1 mm during the annual temperature cycle.This variation is considered to be acceptable.Around the antenna masts, there are three reference bench marks, and connections to the first-order leveling network have also been established.
[17] Two stations have higher masts.There is an anchored 25-m-highsteel grid mast at Metsa ¨hovi, and at Oulu, there is a cylindrical 8-m steel mast.In both cases the height of the GPS antenna is stabilized with an Invar rod.The antenna is isolated from the mast with an attachment piece and a spring system, which is anchored to the bedrock with an Invar rod or wire [Paunonen, 1993].
The system for Oulu was adapted from that at Metsa ¨hovi.Three stations, Olkiluoto, Kivetty, and Romuvaara, were built in cooperation with Posiva Oy, a company which isresponsible for locating sites for disposal of nuclear waste.Local networks around these sites are remeasured semiannually inorder to locate possible deformations.For this reason more stable concrete pillars were chosen at these sites [Chen and Kakkuri, 1994].
[18] All stations are equipped with Ashtech Z-XII GPS receivers, Dorne-Margolin-type antennas, modems, and power supplies.The exception is Metsa ¨hovi, where an AOA SNR-8100 receiver is in use.At Metsa ¨hovi, there is also an external H maser; at all other stations the receiver's internal oscillator is used.Except Metsa ¨hovi and Tuorla, all antennas are equipped with a radome.These have proven less than satisfactory from a radio propagation viewpoint, but it was decided not to change the antenna mount further because of the experience with SWEPOS.
[19] Data are collected using a sampling interval of 30 s and a 5°elevation angle cutoff.During the 1998/1999 winter, CB00 software was installed into all receivers, and the sites were equipped with Vaisala PTU 220 meteosensors.The station histories are summarized in Table 1.

Data Analysis and Geodetic Results
[20] The dual-frequency GPS phase and pseudorange data were processed using the second release of GIPSY software developed at Jet Propulsion Laboratory (JPL) [e.g., Webb and Zumberge, 1993].Dual-frequency phase and pseudorange data from a single 30-hour period acquired from all the sites in the network are analyzed simultaneously.(There is a 3-hour overlap at each end of each observing session.)The GPS data are decimated to achieve an effective sample rate of 300 s; decimation is performed to maintain a manageable level of utilized disk space.For each 30-hour data set we estimated the usual set of parameters, including oscillator (''clock'') corrections, site positions, atmospheric zenith delay parameters, and ambiguity parameters.Satellite orbit parameters were highly constrained to the values distributed by the IGS based on a solution involving a global network of GPS sites.Temporal variations in the clock and atmosphere parameters are modeled as independent random walks [Webb and Zumberge, 1993].We adopted a minimum elevation angle of 15°for all stations.Corrections for the motion associated with ocean loading and solid Earth tides were incorporated in the model [Webb and Zumberge, 1993].
[21] For the SWEPOS sites, which now have multiple antennas, the SNR-8000 data are used in the solution up until 1 August 1998 and the Z-XII data thereafter.For a period following this date we performed a number of solutions with both GPS receivers and determined differences at the 1-mm level or less.
[22] We adopted the value of 10 mm for the uncertainties in the phase measurements at each frequency.The instrumental uncertainties for such measurements are much smaller, perhaps 1 -3 mm [Spilker, 1996].However, experience within the GPS community has shown that the scatter of the time series is greater than the theoretical value based on instrumental noise only.The increase in the scatter above the predicted value can, of course, be attributed to unmodeled phase variations, which may or may not have a white noise (or nearly white noise) nature.In section 5 we discuss a number of errors which might contribute to this increased scatter, and we investigate the spectral characteristics of the site position variations.It is important to remember throughout the paper, though, that the uncertainties for the estimated parameters, including site position and therefore velocity and geophysical parameters, are approximations.
[23] Data processing utilizes a ''no-fiducial'' technique described by Heflin et al. [1992] wherein station coordinates have weak a priori constraints.The results presented in section 3.2 were achieved without fixing the estimated phase biases to integer values, since the software cannot automatically handle such an extensive data set.Independent tests with a smaller data set (SWEPOS stations only) have shown that bias fixing (i.e., ambiguity fixing) leads to a decrease in the formal errors of $20-30%.
[24] All geodetic positions obtained in the GIPSY analysis are finally transformed into a terrestrial reference frame obtained from M. Heflin (personal communication, 2000).This transformation creates a slight inconsistency since the satellite orbits are referred to various releases of International Terrestrial Reference Frame (ITRF) [e.g., Sillard et al., 1998].The problem that concerns us here is that the frames have slightly different net rotations and translations.Using a set of core stations constituted by those that are jointly present in pairwise successive reference frames and applying weights as given by the velocity uncertainties, the rotations and translations derived by least squares adjustment can amount to 1 mm yr À1 at the Earth's surface in any component.We have therefore estimated and applied the interframe rotation and translation parameters to correct the time series of station positions for these biases.The results are thus determined in a rigid frame that is comoving with the Heflin frame, which for our sites is dominated by the realization of the Eurasian plate motion.We, however, are interested in deformations relative to the Eurasian plate motion.
[25] Removing the motion of the Eurasian plate requires a decision regarding the type of motion the comoving Eurasian (i.e., Eurasia-fixed) frame should be allowed.The requirement to observe deformation from a rigid cotravelling frame implies the need to suppress a deformation mode conveyed by a scale rate parameter.Although any residual rotation in one corner of the region could easily be corrected for by adding small rigid rotations, the case is slightly more intricate when one simultaneously considers translations.Considering horizontal motion, the virtue of the GPS data in application to GIA lies in the ability to resolve intersite motions; subtracting the motion of a rigid frame will have no influence on relative deformation.In the case of radial motions, GPS data will not only be used to study relative deformation but also offers the prospect of studying the vertical motion of the crust, for example, incomparison with the sea surface.An ''absolute'' frame is therefore desired.
[26] Allowing for translations in the comoving frame will affect the estimates of vertical components of site motion.One could argue that the European stations are well established and stable, so that the site motions in the Heflin frame, after correction for vertical motion from one or a set of models, would provide a proper regional vertical reference.The associated comoving frame would be constructed by estimating rotation and translation rate parameters using the subset of velocity estimates at those European IGS stations that also were used in the projection stage.
[27] On the other hand, the large number of stations worldwide included in the Heflin frame should provide a stable constraint.In the larger, global set, local vertical motion will appear as less correlated.Accepting this argument, the consequence is to not allow relative vertical motion between the comoving frame and the Heflin frame, i.e., use only the horizontal projection of the site motion vectors in the Heflin frame.This method brings about an advantage, namely, that the observed data do not have to be ''corrected'' with a model, thereby avoiding problems of circular arguments at the stage where the network rates are interpreted.The associated comoving frame is simply constructed by solving for rotation rates only.We have therefore adopted this method.
[28] The analysis described above yields a time series for each station of three-dimensional position in the ''Eurasia comoving'' reference frame.Time series for the BIFROST stations and Tromsø, a nearby IGS site, are shown in Figure 2).

Antenna-Related Issues
[29] Several of the time series in Figure 2 for sites of the SWEPOS network exhibit one or more ''jumps.''We do not believe that these jumps represent motions of the GPS antenna.The jumps are associated with (1) removal and replacement of the GPS antennas, (2) changes in antenna radomes, and (3) rapid changes in snow accumulation.Snow accumulation is discussed in section 5.3.
[30] In the removal and replacement of the GPS antennas to perform the local site surveys mentioned in sections 2.1 and 2.2, the GPS antenna is positioned on the monument by means of a threaded hole and a standard 5/8-inch surveyor's bolt attached to a metal plate that has been permanently set into the concrete at the top of the pillar.The GPS antenna is screwed onto the bolt until it refuses to rotate.When the antennas are removed and replaced, the orientation of the antenna is checked to insure that it is the same as before removal.Given that the surveyor's bolt has 5 threads per inch, a rather large orientation error of 45°w ould lead to a vertical displacement of only 0.6 mm and to no horizontal displacement.The ''jumps'' in Figure 2 associated with these surveys, on the other hand, can be at the 10-mm level.
[31] A likely explanation for these jumps is that very small differences in antenna orientation lead to changes in phase errors because of electromagnetic coupling [Elo ´segui et al., 1995;Jaldehag et al., 1996b] and antenna phase center variations [Schupler et al., 1994].Both these sources of error are potentially elevation and azimuth angle dependent, and in the case of the former the position relative to the pillar and metal plate is critical.If the phase errors induced by these two phenomena are represented as a series of spherical harmonics, with angular arguments of azimuth and elevation angles, then the contribution from the l = 1 term is indistinguishable from the contribution to phase variations from a site position offset.
[32] Other apparent jumps occur when there were changes in the antenna radomes.The original radomes installed on the SWEPOS sites were designed at Delft University of Technology.During the winter of the first year of the experiment, snow accumulated significantly on these radomes, and our observations led us to believe that this accumulation could be reduced by a redesigned radome having no horizontal surfaces.Redesigned radomes (''type A'') were emplaced in the winter and spring of 1995.We later discovered that the paint process used on these radomes was defective.These radomes were thus removed in the spring and summer of 1996 and later that year were replaced with improved radomes (''type B'').Each of these changes appears to produce offsets in [33] As an ad hoc treatment for these errors, we have simply estimated three-dimensional offsets in position on the epochs at which radomes were changed, the GPS antennas removed and replaced, or the antenna rotated.These changes are summarized in Table 1.The site velocity was assumed to be constant for the entire experiment.This ad hoc procedure is not very satisfying, since the existence of the offsets is an indication of an error source which could conceivably have a temporal variation and therefore could effect the estimate of the rate.

Determination of Station Velocities
[34] In this section we report and compare several methods for determining the station velocities.As a result of the analysis described above, the time series of station positions are in a consistent reference frame; it is, in principle, simply a matter of fitting a straight-line component by component and site by site to the time series.This method does not yield determinations of the correlations of the errors in the rate estimates, but these are formally very small since the orbit parameters were highly constrained in the original solutions.In order to gain a quantitative understanding of the effects of errors that are difficult or impossible to model, we present several different analyses for the rates.Each analysis uses different models for the variation of site position with time as well as different editing criteria.
[35] In the following, the standard deviations we report are the so-called ''standard errors.''These standard errors are based on the phase uncertainties used in the daily least squares analysis, described above, propagated through that analysis to yield standard deviations for daily determinations.(The standard errors are the uncertainties shown in Figure 2.) The errors for estimates obtained on different days are assumed to be uncorrelated (notwithstanding the 3-hour overlap).These standard deviations are then propagated through the analyses described in sections 3.2.1 -3.2.3 to yield the standard errors for the rate parameters.In general, the reduced c 2 postfit residuals are close to unity, indicating a reasonable fit, but this statistic may not be an accurate measure of the accuracy of the rate estimates.In the next section, we assess these standard errors and the accuracy of the rate estimates.
[36] The solutions for the three-dimensional crustal deformation velocities and their standard errors for the different analyses are presented in Table 2.The solution used data prior to 1 June 2000.In sections 3.2.1 -3.2.3 we describe the different solutions.

Standard solution
[37] In the ''standard solution'' the model for the position estimates includes a mean value, a constant rate, an admittance parameter for atmospheric loading [vanDam and Wahr, 1993], and periodic terms with frequencies of one, two, and three cycles per year.(The periodic terms are meant to approximately model the effects of snow accumulation and other climatic effects.)Offsets at each antenna change (Table 1) were included.No editing was performed, except that for all sites the February 2000 data were

Edited solution
[38] The parameterization is identical to that of the standard solution, with the following exceptions.No admittance parameters for atmospheric loading were estimated, and only the annual periodic term was included.Offsets for antenna changes were estimated for all components.An editing loop was included that deleted data whose postfit residual was greater than 3 times the weighted root-meansquare (WRMS) residual.This loop was repeated three times.Data from time periods 1995.000 -1995.104, 1995.370 -1995.520, and 1996.438-1996.616were automatically deleted.These periods represent time spans during which radome change/replacement was occurring across the network.

Empirical orthogonal function solution
[39] An empirical orthogonal function (EOF) analysis for the radial rates was performed on a subset of the data.This analysis is described fully in section 5.1.Table 2 includes the radial rates resulting from this solution.

Comparison of Solutions
[40] In general, the edited and EOF solutions compare quite well to the standard solution, with several exceptions that may indicate weaknesses in the solution.The rootmean-square (RMS) differences of the edited solution relative to the standard solution are 0.3 mm yr À1 (east), 0.2 mm yr À1 (north), and 0.8 mm yr À1 (up).The RMS differences of the EOF solution are 0.5 mm yr À1 (east), 0.4 mm yr À1 (north), and 1.4 mm yr À1 (up).The RMS differences are dominated by a few large differences, the largest of which are for the Leksand site and several of the Finnish sites.The rate solutions for Leksand are significantly weaker than those for the other Swedish sites due to the large number of breaks (see section 5.3 and Table 1).There are several reasons why the Finnish sites may have weaker rate solutions.Most of these sites have been observing for a significantly shorter time (Table 1).
[41] To separate these causes, we have calculated RMS differences for the edited solution by dividing the complete data set into two data sets.(To remove the effects of the Leksand breaks, we have not included Leksand in these analyses.)First, we calculated RMS differences for the Swedish and Finnish subsets separately.The results (Table 3) indicate differences mostly for the up component, with the Swedish subset significantly smaller.We then calculated RMS differences separately for sites observing prior to June 1996 and sites first observing after this date.Several Finnish sites fall into the former category, whereas Swedish site Bora ˚s falls into the latter.This division yields a much more dramatic difference for the up component, with the RMS difference for the earlier sites being nearly a factor of 3 smaller.It seems clear, therefore, that apparent differences between the Swedish and Finnish networks are due mainly to the difference in operating time.
[42] For the EOF solution, comparison between the different subsets yields a different result, especially for the horizontal components, which can yield relatively large differences (up to 1 mm yr À1 ).The most striking result seems to be the very small RMS difference for the ''early'' subset, which seems conclusive until it is realized that this subset produces the largest RMS differences for the horizontal components.We do not presently understand these large differences for the solution and leave it to a future study to investigate this issue more thoroughly.

Interpretation of Observed Deformation Rates
[43] As stated in section 1, a primary goal of the BIFROST Project is to provide a new and useful GIA observable with which to constrain models of the GIA process in Fennoscandia.In order to achieve this goal the observations must exhibit a coherent signal that is clearly related to the regional GIA process.In this section we test this requirement by comparing the observed three-dimensional deformation rate signal to numerical predictions of this field and to the apparent sea level signal that has long been associated with the GIA process.We also consider several other geophysical effects that may produce temporal variations in site position.

Glacial Isostatic Adjustment
[44] A number of publications, some dating back to the 1930s [e.g., Haskell, 1935;Vening Meinesz, 1937], have employed sea level observations to infer Earth viscosity and ice sheet parameters in the Fennoscandian region [e.g., Fjeldskaar, 1994;Mitrovica, 1996;Lambeck et al., 1998a;Davis et al., 1999].In the recent study of Lambeck et al. [1998a] a three-layer Earth viscosity model and a regional ice model were proposed that provide a good fit to a carefully compiled and extensive data set based on geological sea level markers.The preferred Earth models are defined by a lithospheric thickness of 65-85 km, an upper mantle viscosity of 3 -4 Â 10 20 Pa s and a lower mantle viscosity that is a factor of 10 or more greater than the upper mantle value.This range of three-layer Earth models and an ice model were also found to produce a good fit to recent instrumented sea level and lake level records [Lambeck et al., 1998b].(These more recent, shorter time scale data apparently did not allow a robust inferenceof lower mantle viscosity.) [45] A number of previous inferences that appear to disagree with the above described viscosity profile [e.g., Wolf, 1987;Fjeldskaar, 1994] are, in fact, found to be compatible when the resolving depth of the various data sets is considered [Mitrovica, 1996].The inference of Mitrovica and Peltier [1993], which is based on the so-called Fennoscandian relaxation spectrum [McConnell, 1968], is not consistent with the Lambeck et al. [1998a] result.However, recent studies show that the paleoshoreline data upon which this spectrum is based require some revision [Wolf, 1996].Indeed, a recent reanalysis of the spectrum was found to eliminate the inconsistency between these two inferences [Wieczerkowski et al., 1999].We have limited the above discussion to recent GIA analyses that considered data from northwestern Europe.We have not considered recent analyses based on global sea level data sets (which may include data from northwestern Europe) in order to avoid the potential bias introduced to these inferences from lateral variations in viscosity structure.
[46] The following predictions are based on a spherically symmetric, compressible, Maxwell viscoelastic Earth model.We choose a three-layer viscosity model defined by a lithospheric thickness of 70 km, an upper mantle of 4 Â 10 20 Pa s, and a lower mantle viscosity of 5 Â 10 21 Pa s. (These values are consistent with the sea level constraints discussed above.)The elastic structure of our Earth model is taken directly from the seismically constrained Preliminary Reference Earth Model (PREM) [Dziewonski and Anderson, 1981].The ocean component of the surface load is computed via a revised sea level algorithm that solves the sea level equation [Farrell and Clark, 1976] in a gravitationally self-consistent manner while incorporating the effects of GIA-induced perturbations to the Earth's rotation vector [e.g., Milne, 1998] and the postglacial influx of ocean water/meltwater to onceice-covered regions [Milne, 1998].Predictions of the loadinduced three-dimensional deformation rate signal are calculated via the theory of Mitrovica et al. [1994a].This theory has also been extended to incorporate the influence of GIA-induced perturbations in the Earth's rotation vector [Mitrovica et al., 2001].The relative importance of the different components of the model will be described in a future publication.
[47] We require an ice model that provides a good fit to the sea level observations for our choice of Earth model.This criterion is met by the model proposed by Lambeck et al. [1998a].However, in order to accurately solve the sea level equation and to realistically compute GIA-induced perturbations to the Earth's rotation vector, we require a global ice model.To meet both of these requirements we remove the Fennoscandian and Barents Sea components of the lower resolution, global ICE-3G [Tushingham and Peltier, 1991] model and replace these by the high-resolution, regional model proposed by Lambeck et al. [1998a].The contours of uplift that we calculate using this ice model and the Earth model described above are shown in Figure 3a.
[48] To illustrate the pattern of uplift that we observe from the BIFROST network and to make a qualitative comparison with model predictions, we have fit a simple surface to the vertical rates from the standard solution.The model that we have chosen for the radial rate _ u at longitude l and f latitude is a two-dimensional Gaussian model: This mathematical representation was chosen because it has a number of parameters that are useful in describing the general amplitude, shape, and location of the observed uplift.The Gaussian is centered at (l o , f o ) and has maximal value uo .The parameters w 1 and w 3 control the widths of the Gaussian and w 2 controls the azimuth of its primary axis with respect to the north direction.In fitting for the seven parameters of the model, we have used for the observation uncertainties , where s _ u is the standard error of the uplift (see above) and s o = 0.8 mm yr À1 , chosen to achieve a c 2 residual rate per degree of freedom of unity.
[49] The resulting model is shown in Figure 3b, and it can be seen that GPS-derived model (which we will herein after refer to as the ''observed rates'') shares much in common with the uplift calculated from the Earth/ice model combination (the ''model rates'') described above.The estimated center of the uplift for the observed rates (l o = 22.5°, f o = 64.6°) is several degrees east of the center of uplift for the model rates.(The observed location of maximum of uplift is closer to that based on apparent sea level rates [Ekman, 1993].)The values of the observed and model maximum uplift rates are nearly the same, however.The areas undergoing subsidence for the observed field differ slightly from those for the model rates, but these areas are outside of the network and good agreement is not to be expected.Finally, the orientation and the amount of elongation agree quite well, although the model deformation is not so symmetric as the two-dimensional Gaussian used to represent the observed rates.This comparison is strong evidence that the vertical crustal motions observed with GPS are associated with the GIA process but nevertheless indicates that there exists a significant disagreement between the GIA model and the observed rates.
[50] A direct comparison between observed rates and those predicted for the Earth/ice model combination described above is shown in Figure 4. Excellent correlations for each component are clearly evident, but systematic differences can be observed.Assuming a simple linear model for the error in the predicted rates, we can determine the best fit lines shown in Figure 4, indicated by solid lines.The RMS residuals to this line are then a measure of the accuracy of our rate determinations.We find a value of 0.4 mm yr À1 for east, 0.3 mm yr À1 for north, and 1.3 mm yr À1 for up.
[51] As a further check on our vertical rates, we compare observed sea level rates from Baltic tide gauge data to rates calculated from the Gaussian model.The tide gauge data consisted of annual averages obtained from the Permanent Service for Mean Sea Level (PSMSL) [Pugh et al., 1987] for tide gauges with time spans of 40 years or longer after 1930.The sea level rate _ s at a tide gauge located within the Baltic at (l, f) is related to the land uplift _ u by where _ g is the rate of change of geoid and Á m is the eustatic sea level rate.The value of _ g is small, <0.5 mm yr À1 for these data, and so may be calculated from a model.In Figure 5 we plot _ s -_ g versus _ u, calculated by interpolating to the tide gauge position using our Gaussian fit to the GPS vertical crustal rates.(An uncertainty of 0.8 mm yr À1 for this interpolated rate yields a c 2 residual of unity.)A strong correlation is evident in Figure 5.This correlation is clear evidence of the relationship between the large apparent sea level rates observed for several centuries in the Baltic and the observed ongoing vertical crustal motions determined from the GPS data.These apparent sea level rates have long been interpreted as indications of GIA and have even been used to refine Earth and ice models [e.g., Davis et al., 1999;Lambeck et al., 1998b].The dotted line in Figure 5 is the line for _ u = 0, while the solid line is for the value of _ u that achieves the best fit, Á m = 1.9 ± 0.2 mm yr À1 , very similar to Model Rate (mm/yr)  inferences obtained with a global network of tide gauges [e.g., Douglas, 1997].
[52] We conclude that the secular vertical crustal rates that we observe using the BIFROST GPS data are mainly associated with the ongoing GIA process but that there is a significant difference between the observed rates and those predicted by our GIA model.In sections 4.2 and 4.3, we consider some other processes that may also contribute to the observed rates.In a future paper, we will present an analysis of the BIFROST observations in which we use the observations to improve the GIA model.

Ocean Tide Loading
[53] The effects of global ocean tide loading as well as the solid Earth tides are treated at the stage of GPS carrier phase data analysis.These motions are predominantly diurnal and semidiurnal; aiming for one site position estimate per day, it appears more advantageous to account for rapid station position variations in the early processing stages rather than having to remove the effects a posteriori.
[54] The ocean loading coefficients were computed with the same method as the one used for the International Earth Rotation Service (IERS) Conventions 1996 [McCarthy, 1996;Scherneck, 1991].The ocean tide model adopted for the processing is taken from Le Provost et al. [1994].This model does not contain the Baltic Sea, a sea area that is central to our region and the possible loading effects of which need a careful account and discussion.It is well established that the diurnal and semidiurnal tides in the Baltic Sea are <20 mm almost everywhere, and therefore their loading effects are to be expected at only submillimeter amplitudes; they can be neglected.The largest seiche mode of the Baltic Sea, an east-west oscillation involving the Bay of Kiel, the Baltic Proper, and the Gulf of Finland occurs at a 36-hour period [Wu ¨bber and Krauss, 1979].They can be excited by fast passing low-pressure areas.Significant amplitudes are found only in the bays at either end, lasting a couple of days.The existence of these modes argues for including either time series of water level at diurnal if not more rapid sampling rates of near-by tide gauges or predicted loading effects based on such observations by means of a hydrodynamic model interfaced with an elastic deformation model [Scherneck, 1991].
[55] On the timescale of days to years the situation for Baltic Sea is radically different.The geometry of the Baltic Sea basin, being well enclosed and connected to the open ocean only through narrows in Denmark and between Denmark and Sweden, causes the mass exchange with the world ocean to be retarded.Seasonal variations of the water level can reach ±0.5 m as a combined response of the air pressure and wind variability and due to the role of the Baltic Sea as a large catchment area where precipitation and evaporation are highly variable on seasonal to interannual timescales.
[56] In order to obtain rough estimates of the impact of variations in the hydrology of the Baltic Sea area on ground deformation we have conducted a simple pilot study.First, we assume that GPS monuments, the locations at which we aim to predict vertical crustal motion, are exactly following with the movement of a solid, homogeneous, elastic crust; that is, we neglect porosity-related deformations of soils and surface layers.We then devise a grid of 5-km mesh width covering the area of Sweden, Finland, and most of Norway, on which we distinguish the type of land and water coverage: (1) land, (2) open ocean, (3) Baltic Sea, (4) great lakes, and ( 5) large hydropower reservoirs.The deformation is modeled using integrated point load Green's functions in the usual way [e.g., Scherneck, 1991].Going through all cases separately, assuming a unit height slab of water in each of the land types, we thus can model admittance coefficients for the impact of loading due to accumulated snow, rain, and soil moisture for type 1 coverage and water level variations in each of the bodies for type 2 -5 coverages.We leave a detailed account of these studies for future publications.Here we use typical maximum values for the amplitudes of the loading processes in order to get an idea of the importance of the effects.The results are summarized in Table 4.
[57] We have also investigated the potential effects of a sea level tilt rate, observed from an analysis of data from the TOPEX satellite.The maximum vertical deformation occurs near the south of the Baltic, at latitude 56°, where we calculated a subsidence rate of $0.25 mm yr À1 .This effect is nearly negligible for our purposes.Furthermore, one would expect that such a tilt, if caused perhaps by wind stress, is quite variable over longer times, so that the effective secular elastic loading associated with this effect will on average decrease with time.

Atmospheric Loading
[58] In the EOF mixed regression we model a time series of atmospheric loading for every station.Previous work on this problem [vanDam and Wahr, 1993] showed generally low air pressure admittance at Onsala and small reductions of postfit c 2 .The timescale of the variations in the pressure field is, unlike in the case of ocean tide loading, predominantly in the range of more than 1 day.Also, the presently available processing software is not prepared for the input of three-dimensional time series of a priori site displacement information.
[59] We computed the atmospheric loading effect analogously to ocean tide loading, with the major difference being that we assume the loading effect is zero at the bottom of the open ocean assuming an inverse barometric response.Global air pressure fields at 1°Â 1°spatial and 6-hour temporal resolution are obtained from the European Centre for Medium-range Weather Forecasts (ECMWF) and convolved with elastic loading Green's functions for vertical and horizontal displacement.Although it can be shown that loading beyond 2000 km distance contributes little to station displacement, we use the entire global field.In this case the global mass balance is easy to maintain, and annual oscillations between the hemispheres do not offset the displacements.We compute an average pressure field for the entire time span in order to subtract the displacement due to the average atmosphere.Thus, for most of the stations a near zero mean for the computed pressure loading time series is obtained.
[60] In the eigenvector analysis the air pressure loading information that is orthogonal to the station residuals is retained, while the common mode will preserve coherent (correlated) residual signal power from this source.Such information, however, is expected to be greatly suppressed since the station residuals result from a regression that already includes air pressure loading.
[61] The reason why we estimate an admittance parameter of the predicted loading effect rather than applying the effect as a correction is as follows.The admittance coefficients that we obtain are systematically and significantly lower than unity.We suspect that the GPS orbits induce regional perturbations since atmospheric loading is not applied at the stage of orbit computation.Since only a small number of stations in northern Europe are used by the orbit centers, only a certain fraction of displacement is conveyed into the orbit.
[62] A second issue is the possible mapping of air pressure related information into other parameters of the GPS analysis.Here, the most probable candidate is the atmospheric delay parameter, and a particular reason to suspect it is the way the hydrostatic and the water vapor related delays are parameterized [Segall and Davis, 1997].However, we can show that only submillimeter vertical site offsets can be expected when atmospheric pressure varies as much as ±30 hPa.
[63] Bottom pressure equilibrium requires days to weeks to establish in shallow waters.Since the Baltic Sea is a nearly enclosed basin, the inverse barometer response is delayed.Thus it appears more promising to neglect the water response in the atmospheric loading model and add another signal channel in the linear regression representing the water level of the Baltic Sea at a tide gauge station nearby.

Discussion of Errors
[64] In this section we carefully examine the possible influence of a number of errors on our main geophysical observable, the site velocities.We include this study for several reasons.The expected magnitude of GIA contribution to the velocities is fairly small, typically <3 mm yr À1 for the horizontal component and <10 mm yr À1 for the vertical.The technique of continuous GPS is rather new, and no careful analysis of errors has yet been performed for determinations of velocity obtained from these data.We begin by assessing the spatial dependence of variations observed in our position determinations.In a future paper, we will examine spectral analyses of our time series.

Correlation Analysis
[65] The time series of Figure 2 display a great deal of variability that is clearly not white noise.Aside from the seasonal variations that are quite large at several sites, lowfrequency noise is also evident.We anticipate that there may be a number of noise sources affecting the stations in a region in a similar way, such as satellite orbit errors, reference frame errors, errors at one or possibly more sites that are propagated through the network due to the type of network solution we perform, environmental conditions that change over a region in a coherent way (for example, soil surface reflectivity affected by climatic factors, snow covering antenna and radomes), and short-lived nonsecular crustal deformations due to predominantly atmosphere and hydrosphere loading.
[66] Since we have time series of perturbed site positions, we seek a method that takes advantage of the statistics inherent in the large amount of information and relate to the separation of local and regional signals and noises.We have chosen to represent the degree of correlation of perturbations to site position (eventually including coherent, transient motion) between stations using an empirical orthogonal function (EOF) type of analysis.(Davis and Elgered [1998] used an EOF method with estimates of water vapor determined from BIFROST data.)It is desirable to utilize this information in the adjustment when we solve for rates, offsets, and other locally relevant parameters, thus attempting to discriminate between local deterministic processes and correlated transient signals.
[67] Secular GIA motion is, of course, correlated between the sites.Thus we need a first stage where parameters for this process (and others we are aware of) are estimated and residuals formed, followed by a second stage where we then use these residuals to obtain improved values for the rate and other parameters.
[68] In the locally relevant parameter set we include atmospheric loading since this perturbation can be predicted on a per-site basis.However, a certain fraction of atmospheric loading perturbations might actually be transferred into satellite orbit perturbations, since corrections of site positions due to this effect are not routinely estimated in the precise orbit generation phase.Thus, as we will see, the atmospheric loading signal is attenuated in the GPS time series, although some signal apparently leaks in during the second (EOF) stage due to correlation between the common mode and the atmospheric loading time series.5.1.1. Procedure [69] This section describes the two-stage least squares procedure we have developed that uses the GPS time series of station positions to determine a combined set of parameters that include (1) the usual set of rate and offset parameters; (2) admittance parameters for geophysical signals such as pressure and precipitation that we might expect to be correlated with the estimated GPS time series; and (3) between-site correlation parameters.We begin in the first step by looking for the solution to the linearized system The unknown parameter vector p contains both the usual set of rate and offset parameters and the admittance parameters for geophysical signals for one station.The vector d is a tobe-modeled time series (for example, north, east, or up components of position for one station) and e its vector of errors.The matrix G is the design matrix defined in the usual way.If we are including surface pressure in the model, then one of the columns of G will be the surface pressure at the site corresponding to the epochs for which the determinations of d were obtained.
[70] As usual, the solution to (3) is obtained by minimizing the mean square residual and can usually be determined using the standard least squares formulation.However, since this system can in general be underdetermined or otherwise singular, and for computational compatibility with the second step of the process, we solve (3) using the generalized inverse method of Lancxos [Aki and Richards, 1980].In this method we first normalize the system of (3) such that the expected variance of the e is the unit matrix.(We assume for this analysis that the observations are uncorrelated.)We denote the normalized quantities with a tilde ( G, d, and so on).The generalized inverse of G is expressed as a singular value decomposition (SVD), wherein we seek the solution to the eigenvalue equations where u is a matrix of eigenvectors that span the parameter space, v is a matrix of eigenvectors that span the data space, and l represents the set of eigenvalues.Aki and Richards [1980] demonstrate that u and v share the first m eigenvalues, where m is the minimum of the dimensions of u and v; the remaining eigenvalues are zero.
[71] In our analysis for this study we modeled the time series of vertical components of station position for the BIFROST sites.In addition to rates and biases for each site we estimated site-dependent sinusoidal amplitudes (in and out of phase) with annual, semiannual, terannual, and quarterannual periods and (site-dependent) admittance parameters for surface pressure.
[72] In the second stage of the analysis the solution was repeated with the modification that for each site we also estimated parameters that represent admittance for the residuals calculated in the first stage for all the BIFROST sites.For example, if using the data from the Umea ˚site yields a nonzero admittance for the stage-one residuals from Kiruna site, that result indicates a correlation between the unmodeled signal from the first stage for Kiruna and the observed signal for Umea ˚.In this sense, the second stage is equivalent to a traditional (unweighted) EOF analysis.Moreover, in a manner analogous to a traditional EOF analysis we keep only those parameters whose eigenvalues signify that they convey a significant amount of information.
[73] In the second (EOF) stage we define two subsets of the parameter space.The subset of the parameters that were estimated in the first stage we name C; these parameters will be estimated regardless of their eigenvalue.The subset consisting of the new (residual admittance) parameters of the EOF construct we name E. Sorting the eigenvalues of E by decreasing magnitude yields one which is the largest, l c .Associated with it is eigenvector v c , which is to be retained; the remaining subspace of E is to be ignored.In a traditional space-time EOF analysis, v c is known as the first temporal eigenvector.
[74] The eigenvector v c retained from subspace E extends the data space to comprise a common mode.The common mode time series is obtained from the parameter solution by simply projecting it on the subspace E and inserting it into the model (3).The common mode construction is a form of a spatial filter performing a weighted mean, where the weights are devised by the eigenvalue process and hence by the data themselves.The procedure that we devised may be contrasted to the filter of Wdowinski et al. [1997], which was formed by taking a mean of the time series for the stations in the network.Each observation is equally weighted regardless of uncertainty or station.In our procedure we use the uncertainties, and the station weights are determined by the station admittance parameters, calculated to minimize the mean-square residual.Assuming that the first eigenmode admits the station residuals with equal weight, our procedure will, in fact, perform exactly the same as the Wdowinski method.
[75] As in a traditional EOF analysis, the relative signal power propagated through each eigenvector in the inverse solution is proportional to the squared magnitude of the eigenvalue.In our solution using the vertical component of site position, the ratio l 1 2 /AE j l j 2 is 10%.Thus a filter based on the Wdowinski et al. [1997] approach would decrease the variance of the time series by this amount.In Figure 6 we show an example residual time series, for site Umea ˚, for the first-and second-stage analysis, as well as for the common mode.In this example the RMS of the ''corrected'' residuals are $50% of the ''uncorrected'' residuals.(The variance reduction is 66%.)

Spatial correlations
[76] The eigenvectors and eigenvalues carry information about the correlation of the time series among the sites, accounting for the model of the first adjustment stage.This information stems from the similarity of the scalar product used in the eigenvalue solution process and the manner in which correlations are computed.Figure 7, for example, contains a plot of the parameter eigenvectors for an analysis of the vertical components for a subset of 11 SWEPOS sites: Metsa ¨hovi, Skelleftea ˚, Umea ˚, Sundsvall, Ma ˚rtsbo, Leksand, Sveg, O ¨stersund, Vilhelmina, Arjeplog, and Kiruna.The 11 admittance parameters representing E are in positions 8 -18.It is clear from the primary eigenvector in Figure 7 that the contribution from the different sites is fairly large and constant.The exception might be the first site, Metsa ¨hovi.This site is used to constrain the reference frame (see section 3).The correlation between regional stations and a station which has been used for mapping into the reference frame is generally low because the average rotations and translations of the the small number of constraining stations that have been used are only slightly overdetermined; thus they are propagated throughout the network.
[77] From Figure 7 it can be seen that the first eigenvector is associated primarily with the admittances.The second eigenvector is associated mostly with the offset and rate terms, the third and fourth are associated with the seasonal parameters, and the fifth is associated with pressure loading.If we constructed a filter based on the first five eigenmodes, we would reduced the variance of the residuals by 34%, compared to 10% for the first eigenmode only.In fact, the EOF method enables one to construct filters that reduce the variance by any amount that is useful; if all eigenmodes are used, the variance is reduced 100%.The criterion that we use for acceptance of a mode is that the time series predicted by the mode should correlate with one of the residual time series by >50%.Using this criterion, we accept the first 17 modes, which produces a variance reduction over all the sites of 73%.
[78] It is straightforward to show that the correlation coefficient g jk between time series from the jth and kth stations can be computed using the first spatial eigenvector of E, denoted u c : where N is the number of data.The expected low correlations with constraining stations, in this case Metsa ¨hovi   (parameter 8), is clearly seen in Figure 8, which shows the correlation as a function of intersite distance.The remaining rate parameters exhibit great coherence (g ' 0.5) on a regional scale.Figure 8 does seem to indicate, however, that there is a clear but weak dependence of correlation on intersite separation.(The best fit line to the non-Metsa ¨hovi correlations of Figure 8 yields a slope of À0.2 per 1000 km.)This result implies that the cause of the correlation is network-wide, indicating perhaps a reference frame or orbital-type effect.

Effect of reference frame errors
[79] There are two types of reference frame errors to consider: errors formally accounted for by the error covariance matrix and biases.In the ideal case of independent Gaussian measurement errors with perfect geodetic models, linear propagation of the measurement error covariance would yield an accurate characterization of the errors in geodetic model parameter estimates.Uncertainties in velocity estimates derived from such position determinations would similarly be accurate and well understood.In the real world, however, geodetic models are not perfect, measurement noise processes are not normally (i.e., Gaussian) distributed, and the measurement errors are not generally independent.Each of these factors complicates the parameter estimation problem.The resulting parameter estimates are likely to be contaminated by systematic errors which are difficult, perhaps impossible, to fully assess.
[80] An exhaustive list of the mechanisms through which systematic errors could manifest themselves as reference frame biases and their respective importance will require continued research.However, it is not difficult to list some important potential mechanisms.Deficiencies in the geodetic models describing satellite orbital dynamics or the dynamics of Earth's orientation, for example, have obvious implications for an accurate reference frame realization.For local networks, errors in the atmospheric modeling or other spatially correlated errors, such as those due to similar scattering at similar monuments, could result in local reference frame biases.The same is also true of overconstrained a priori parameter estimates, such as for satellite ephemerides, satellite clock variations, or site positions, particularly if these a priori estimates were themselves correlated.These ''fiducial errors'' have been appreciated for some time [e.g., Larson and Agnew, 1991].Site-specific errors, associated with such phenomena as multipath or antenna phase center variations, may manifest themselves indirectly as reference frame biases as we now discuss.
[81] Trade-offs between parameters of the geodetic model, such as between the implicit specifications of the satellite and GPS network orientations, render the site position determination problem ill-posed.That is, in the absence of additional constraints the matrix of partial derivatives relating the basic prefit GPS carrier phase and pseudo range residuals (i.e., the data) to the first-order corrections to the a priori parameter estimates (including site positions, satellite parameters, and Earth orientation parameters) is singular.The degree of singularity depends on the geometry of the tracking network; regional scale networks are ''less singular'' than global-scale networks.A priori information, such as knowledge of the positions of certain sites in the network, is often effectively employed to regularize this singularity.The resulting position estimates implicitly define a reference frame.Depending on the level of constraints imposed and assuming, for now, that this level is consistent with the true errors in the a priori estimates, some linear combinations of the model parameter estimates, such as the differences in positions of the sites, may be significantly better determined than the This site was used to establish the reference frame, so the correlated variations in position might have been removed to a greater extent.The best-fit line is g = 0.6 À 0.2 B, where B is the line length in 1000 km.
''absolute'' model parameters themselves.However, sitespecific errors not formally accounted for by the covariance matrix will affect all components of the vector of linear combinations involving the contaminated site.The importance of any reference frame errors resulting from this transformation will thus depend on the transformation itself (the particular linear combination), the number of sites suffering from systematic errors, and the size of these systematic errors.
5.2.Atmosphere [82] As discussed in section 3, the GIPSY analysis software uses a stochastic filter to model the temporal variability of the zenith atmospheric propagation delay.The a priori model for the propagation delay for each site consists of a constant ellipsoidal height dependent term representing the hydrostatic zenith delay [Davis et al., 1985] that is ''mapped'' to the elevation angle of the GPS signal using the Lanyi dry mapping function [Lanyi, 1984].The a priori wet zenith delay is taken to be 100 mm.Corrections to this a priori total zenith delay are estimated by adopting the Lanyi wet mapping function to describe how the correction maps with elevation angle.In effect, this procedure maps the combined corrections to both the hydrostatic and the wet zenith delays using the wet mapping function.Using a simulation procedure similar to that of Elo ´segui et al.
[1995], we find that the error in the vertical component of site position (the primary geodetic parameter affected) caused by mapping the correction to the hydrostatic zenith delay using a wet mapping function is less than ±2 mm, assuming an annual pressure variation of ±30 mbar.
[83] The mapping functions used assume that the atmosphere is azimuthally symmetric as viewed from the station.Horizontal structure in the atmosphere nevertheless generally leads to an effective azimuth-dependent mapping function.In particular, horizontal gradients in the refractive index lead to sinusoidal variations in azimuth [e.g., Davis et al., 1993].We have not incorporated this model into our phase model or estimated any parameters associated with horizontal structure, leaving this to future studies.Although such errors can be quite large at low-elevation angles [Davis et al., 1993], we do not expect this error to systematically bias the estimate of the site rates.Of much greater concern are the azimuthaldependent effects of snow accumulation on the radomes.

Signal Effects
[84] In the course of our experiment we have investigated several errors involving the propagation of the GPS signals in the nearby environment of the GPS antenna.Near-field scattering [Elo ´segui et al., 1995] involves the reflection of the GPS signal off surfaces close to the GPS antenna.This effect is different from multipath in that the reflecting surfaces are well within the near field of the antenna.These surfaces thus become electromagnetically coupled to the antenna, effectively creating a new antenna.One such reflecting surface is the pillar directly beneath the GPS antenna, which often (and in the case of our Swedish sites) has a metal plate embedded onto its top surface.Jaldehag et al. [1996b] investigated this effect for sites of our Swedish GPS network and found that for the most part the scattering effects are equal for all the sites (and therefore cancel in the baselines), since with the exception of the Onsala GPS site all the pillars are identical and the local vertical direction at each of the sites is nearly parallel.It is possible that this effect is in part responsible for the common mode site variations that we find (section 5.1).In the principal component analysis that we performed, we found that the component of the spatial eigenvector for the Onsala site had the smallest value of all the Swedish sites.
[85] The error associated with the signal scattering effect is quite large, even when the size of the phase error is small.In the least squares estimation procedure the errors in each of these parameters can be ''magnified'' and yet the summed contribution to the phase can be fairly small.The effect on the horizontal components is smaller than that for the vertical but not zero.Elo ´segui et al. [1995] found that the error in the estimated vertical component was therefore dependent on the geometrical distribution of observations on the sky.For a given fixed constellation and sequence of observations the error in the estimate of site position is dependent on the minimum elevation angle accepted for use in the data analysis.In practice, however, it is not possible to maintain a fixed geometrical distribution of observations on the sky because of data dropouts.
[86] We also found that snow and ice, which adheres to the radome and accumulates, is a significant source of error simply by causing an additional propagation delay.This effect was first noted by Webb et al. [1995] when several meters of snow buried a GPS antenna near Long Valley.The amount of snow and ice which accumulates on our antennas is much smaller, perhaps several decimeters maximum, and the problem is therefore more insidious.(The lower amount of snow also tends to be distributed unevenly on the radome.)Jaldehag et al. [1996a] found that they were able to approximate the results from elevation angle cutoff tests using a conically symmetric distribution of snow, but visual inspection of our sites indicates that real conditions are not so ideal.
[87] An effect that we have already mentioned is the propagation delay due to the radome.In theory, an expression can be developed for this contribution, and in the future we plan to test such a model.For this paper, however, we have simply introduced an offset to the time series on the epochs which we have changes radomes.In the same way we have treated changes to microwave absorbing material placed around the antenna.The dates of these changes are also noted in Table 1.

Monument Stability
[88] As described in section 3, we monitor the relative position of the pillar reference point within a local (10 -15 m area) network of reference pins.The determination of the pillar position is given in Table 5.The expected uncertainty of these determinations, based on propagation of the instrumental errors through the resection technique, is $0.1 mm.On the basis of the RMS analysis of all the data the standard deviation of a single measurement is 0.2 mm for the north, 0.3 mm for the east, and 0.4 mm for the vertical.However, there are several occasions when repeated measurements are obtained in a single day.The five repeated measurements from the Leksand site, for instance, yield RMS differences of 0.01 mm for north, 0.04 for east, and 0.06 mm for vertical.The single repeated measurement for Kiruna, on the other hand, yields an RMS difference (based on two measurements) of 0.2 (north), 0.3 (east), and 0.1 mm (up), roughly consistent with the same statistics for the single repeated measurement for the O ¨verkalix site, 0.03 (north), 0.1 (east), and 0.1 mm (up); and for the Norrko ¨ping site, 0.1 (north), 0.1 (east), and 0.1 mm (up).If we combine the Kiruna, O ¨verkalix, and Norrko ¨ping repeat measurements, we find an RMS error of 0.10 (north), 0.18 (east), and 0.10 (north).We will therefore adopt the average 0.13 mm as our uncertainty in each of the components.We have no explanation as to why the differences in the repeated Leksand measurements are so small relative to these values.Eliminating the multiple measurements by replacing them with their average, the overall RMS variations for all sites are 0.20 (north), 0.31 (east), and 0.38 mm (up).Thus, overall, the average RMS variations are consistent with the adopted measurement error at the level of 3s.
[89] Nevertheless, for individual pillars, most notably Skelleftea ˚S (east component) and Umea ˚(up component), the variations are greater than 1 mm.In the case of Skelleftea ˚S, the variations for the north and vertical component are small, so unless the pillar happened to have moved on a nearly north-south axis this result is spurious.The Umea ˚variation is in the vertical component, which, one can surmise, could experience motions independent from the horizontal components.The most worrisome effect would be settling of the pillar.However, the measurements indicate that the pillar was 2.3 mm higher relative to the nearby network in June 1995 than in June 1993.In order to resolve this and other issues, further pillar measurements are planned for the near future.However, these future measurements will be obtained in a slightly different manner.Reflectors will be fixed to the pillars, and the antenna will not have to be removed.
[90] More frequent measurements were obtained on the Leksand site than for others for test purposes.The time series of measurements is shown in Figure 9. Langbein and Johnson [1997] have investigated temporal variations in line lengths determined from laser distance measurements; their analysis indicates that the power spectra of the variations follows a f À2 behavior, indicative of a random walk process.They argue that these variations are associated with pillar motions.Although our Leksand data are sparse, we can perform a maximum-likelihood analysis on the assumption that the temporal variations of the pillar position can be described by a random walk and that the measurement uncertainties are 0.13 mm, as described above.For the three components the maximum-likelihood estimate of the random walk variance was 0.3 (north), 0.2 (east), and 0.1 mm 2 yr À1 (up), values much less than the average value of 1.7 mm 2 yr À1 found by Langbein and Johnson [1997].The monuments they investigated are very similar to those we used for the Swedish GPS network: for instance, the sites of the Parkfield network used galvanized steel pipe of diameter 20 mm placed into augered holes 2 m deep.Near two of these Parkfield sites, special monuments were installed to a depth of 10 m; the average random walk variance for these two deep-anchored monuments was $2 mm 2 yr À1 , whereas that for the shallow monuments was $15 mm 2 yr À1 [Langbein and Johnson, 1997].The differences between the SWEPOS and the Parkfield monuments may be in the geology; all of the SWEPOS monuments are located on exposed bedrock (see section 2.1).
[91] The average random walk variances, based on the difference between end-point measurements for the 20 sites with pillar measurements spanning longer than 6 months, is 0.08 (north), 0.21 (east), and 0.26 mm 2 yr À1 (up).

Discussion and Future Work
[92] We have used a regional continuous GPS network to measure three-dimensional crustal deformation rates in Fennoscandia.The observed rates correlate highly to rates predicted for a viscoelastic Earth undergoing present-day glacial isostatic adjustment (GIA) due to the rapid melting of late Pleistocene glaciers, but significant systematic differences from the predictions can be observed.The observed rates correlate well also with observed sea level rates obtained from tide gauge data over the last 30 years.
[93] We see no evidence for the regional shear inferred by Pan and Sjo ¨berg [1999] using a much smaller GPS network and only two campaign measurements.Pan and Sjo ¨berg [1999] found relative horizontal motions across the north-south extension of the Baltic Sea to be 2 -3 mm yr À1 , whereas using the Table 2 rates, we find the motions of these sites to be generally less than 1 mm yr À1 , probably within the uncertainties.Our velocities are roughly consistent with the expected velocities near the center of uplift.[94] We get much better agreement with theoretical predictions of crustal velocity based on GIA, relative to the amplitude of the predictions, for the radial (vertical) component than for the horizontal components of station velocity.This result may seem counterintuitive since radial velocity estimates obtained from GPS data are well known to be less accurate that horizontal velocity estimates.Errors in the ice model, however, would tend to have a proportionally greater effect on the horizontal components of velocity.For example, a slightly misplaced maximum ice thickness could easily change the predicted horizontal rates by 100% near the maximum, while changing the radial rate by only $10%.
[95] Although we have much evidence for time-correlated behavior in the site position time series both from power spectrum calculations and from visual inspection of the time series themselves, our measurements for pillar motion indicate that the top of the pillar moves less than 1 mm relative to a very local network of pins.This result may indicate that ''monument instability,'' if it exists for our sites, may not be the result of local deformation of the crust.If so, then the ''footprint'' would have to be extended in order to measure this phenomenon.
[96] The empirical orthogonal function (EOF) analysis that we performed revealed a correlated noise in radial position time series across a large region ($1000 km).This analysis also indicated that these correlations decrease slightly with distance.Typical RMS residuals for our radial time series are 7 mm.The EOF analysis indicated that over half of the noise contributing to the RMS residual is from this correlated noise.The source of the correlated noise, unfortunately, cannot be revealed by an EOF analysis.We and others have speculated that the correlations are due to systematic errors in the satellite orbits, terrestrial reference frame, or both.If the source of this error can be identified and eliminated, then there is the possibility of reducing velocity uncertainties from GPS by $50%.
[97] We have reported the ''standard errors'' so that users of the rates may have the straightforward results of the least squares analysis.We have preliminarily attempted to determine our ''true'' uncertainties in several ways.The most conservative method would be our comparison to a model (section 4.1).This comparison would indicate uncertainties of 0.4 mm yr À1 for east, 0.3 mm yr À1 for north, and 1.3 mm yr À1 for up.Our comparison of the different solutions (section 3.3), however, would indicate that the sites that have been observing the longest may have uncertainties $50% of these values.
[98] The results presented here provide the first dense regional geodetic determinations of ongoing crustal deformation associated with GIA.The rates estimated herein are intended to be used as a new observable in the study of this phenomenon.In a future work, we will use the differences between the GIA predictions and the observed rates to improve the ice and Earth models used to calculate the predictions.We will also use the radial rates to correct tide gauge data for vertical crustal motion to obtain absolute sea level rates.
[99] Acknowledgments.This work was supported by NASA grants NAG5-1930 and NAG5-6068, NSF grant EAR-9526885, the Natural Sciences and Engineering Research Council of Canada, the Swedish Natural Science Research Council, the Swedish Council for Planning and Coordination of Research (FRN), the Knuth and Alice Wallenberg Foundation, Posiva Oy of Finland, and the Smithsonian Institution.GPS data in noncontinuous mode were obtained with the support of the University Navstar Consortium (UNAVCO) Facility.We thank S. Nerem for providing Baltic sea level rates estimated from his analysis of TOPEX data.We greatly appreciate the comments and suggestions of M. Ekman.B. Bills, T. James, and an anonymous referee provided comments that substantially improved the manuscript.In the analysis, data from sites of the IGS were used.These data were obtained from the NASA Crustal Dynamics Data Information System (CDDIS) and University of California, San Diego, SOPAC archives.Some figures were generated using Generic Mapping Tools (GMT) version 3 [Wessel and Smith, 1995].

Figure 1 .
Figure 1.Map showing GPS sites used in the study.Subnetworks are indicated by the symbols used: SWEPOS (triangles) and FinnRef (circles).IGS sites are indicated by diamonds.The Kiruna IGS site, located close to the Kiruna SWEPOS site, is not shown.The Onsala and Metsa ¨hovi sites are also IGS sites.

Figure 2 .
Figure 2. Time series of site positions for the BIFROST sites and Tromsø from the standard solution.Differences (in mm) from a series mean for the (left) east, (middle) north, and (right) radial components of site position.The light line shows the model fit to the time series, and the vertical lines indicate epochs at which an offset was introduced.

Figure 3 .
Figure 3. (a) Vertical crustal rates from GIA predictions calculated using a method outlined in the text.(b) Simple Gaussian model produced by a least squares fit to the observed GPS vertical crustal rates.See color version of this figure at back of this issue.

Figure 4 .
Figure 4. Comparison between the rates observed from the GPS data and those predicted by the Earth/ice model combination discussed in the text for the (a) east, (b) north, and (c) radial components of velocity.

Figure 5 .
Figure 5. Observed sea level rates from Baltic tide gauges versus negative crustal vertical rates obtained from the Gaussian model of Figure3b.The mathematical relationship between these quantities is given in equation (2).

Figure 6 .
Figure 6.Residuals for radial component for Umea ˚from EOF analysis (a) first and (c) second stages, after the recycling of the (b) ''common mode'' based on the residuals for the network of sites.For clarity the error bars are not shown.The WRMS residuals are 8 mm (Figure6a); 6 mm (Figure6b), and 4 mm (Figure6c).

Figure 7 .
Figure7.Eigenvectors in model space.The signal channels through which residuals at neighbor stations have been recycled are numbered 8 -18.The associated GPS stations are, in order of their parameter number, Metsa ¨hovi, Skelleftea ˚, Umea ˚, Sundsvall, Ma ˚rtsbo, Leksand, Sveg, O ¨stersund, Vilhelmina, Arjeplog, and Kiruna.The residuals contribute to the set of signals from which the (orthonormal) eigenvectors are constructed.Eigenvector 1 conveys a common mode of the recycled noise.Deleting from the model space all eigenvectors that convey the residual information except the one that is associated with the largest eigenvalue constitutes the extended EOF method used in this paper.

Figure 8 .
Figure8.Correlation between vertical time series from pairs of GPS sites as a function of intersite line length.The smaller ''branch'' at a correlation of $0.2 contains correlations with site Metsa ¨hovi.This site was used to establish the reference frame, so the correlated variations in position might have been removed to a greater extent.The best-fit line is g = 0.6 À 0.2 B, where B is the line length in 1000 km.

Figure 9 .
Figure 9.Time series of the pillar position measurements for the Leksand site relative to a local network of monuments.The error bars are equal to the RMS based on the repeated measurements from other sites (see text).The symbols represent the different components of position: square with solid line, north; circle with dashed line, east; and triangle with short dashed line, up.Each of the time series have been de-meaned, and slightly offset in the abscissa for clarity.

Figure 3 .
Figure 3. (a) Vertical crustal rates from GIA predictions calculated using a method outlined in the text.(b) Simple Gaussian model produced by a least squares fit to the observed GPS vertical crustal rates.

Table 1 .
(continued) the dates on which significant modifications were made to the site hardware.(Minorchanges,suchasthose affecting communication only, are not indicated in Table1.)Thehardwaredescribed in Table1includes the GPS receiver type, GPS antenna type, radome type, monument type, and approximate positions.IERS Domes numbers for those sites, which are IGS sites, are given as well.

Table 1 .
(continued)The Domes number is a unique station identifier issued by the International Earth Rotation Service.Only some of the SWEPOS sites have been so registered.DM Ash.''ARR indicates that the antenna was removed and replaced.The ''a'' indicates a change to microwave absorbing material placed near the antenna.

Table 2 .
Site Velocities in the Eurasia-Fixed System a See text.The formal uncertainties are not shown as they are nearly all equal for a given component: east and north, 0.1 mm yr À1 , and up, 0.2 mm yr À1 .The exceptions are Leksand east, 0.2 mm yr À1 , and up, 0.4 mm yr À1 ; Onsala up, 0.4 mm yr À1 ; and Romuvaara up, 0.3 mm yr À1 .Tromsø was not included in the EOF analysis.

Table 3 .
Comparison of Edited and EOF Solutions With the Standard Solution a Leksand included in the ''all sites'' data set only (see text).

Table 4 .
Tidal Loading Effects b Also Ha ¨ssleholm.