## Abstract

Analysis of Global Positioning System (GPS) data demonstrates that ongoing three-dimensional crustal deformation in Fennoscandia is dominated by glacial isostatic adjustment. Our comparison of these GPS observations with numerical predictions yields an Earth model that satisfies independent geologic constraints and bounds both the average viscosity in the upper mantle (5 × 10^{20} to 1 × 10^{21} pascal seconds) and the elastic thickness of the lithosphere (90 to 170 kilometers). We combined GPS-derived radial motions with Fennoscandian tide gauge records to estimate a regional sea surface rise of 2.1 ± 0.3 mm/year. Furthermore, ongoing horizontal tectonic motions greater than ∼1 mm/year are ruled out on the basis of the GPS-derived three-dimensional crustal velocity field.

Previous maps of ongoing postglacial “rebound” in Fennoscandia have been based on tide gauge records of sea level change and on conventional geodetic leveling surveys from the past century (1). These maps show a vertical deformation field that is broadly correlated with maximum Late Pleistocene ice cover and a peak apparent uplift rate of 9 mm/year near the northern tip of the Gulf of Bothnia. Improving on the accuracy of these previous studies and extending them to incorporate three-dimensional (3D) crustal deformations are important in several respects. First, geophysical observables related to vertical glacial isostatic adjustment (GIA) have been used to constrain the viscosity structure of Earth's mantle and the space-time geometry of Late Pleistocene glaciation (2–10). Numerical models show that horizontal motions are sensitive to variations in ice geometry and mantle visocity (11–15), and therefore observations of these motions provide additional information to improve previous inferences (7). Second, efforts to determine present-day global or regional sea level trends by correcting tide gauge records for the GIA signal have eschewed the Fennoscandian tide gauge record because of the large GIA amplitude in this region (up to ∼10 mm/year) relative to the global trend (1 to 3 mm/year) (16). A map of radial crustal velocities permits a direct, rather than model-dependent, correction of the GIA contribution to tide gauge–derived secular sea level trends. Third, it has been suggested that Fennoscandia is subject to recent tectonic deformations oriented along localized zones of weakness thought to exist in the Baltic Shield (17). This possibility has implications for regional seismic hazards, and a map of 3D motions will clarify the relative importance of GIA and tectonics in the present-day Fennoscandian deformation field.

**BIFROST map of 3D crustal velocities**. The Baseline Inferences for Fennoscandian Rebound Observations Sea Level and Tectonics (BIFROST) project was initiated in 1993 to map 3D regional deformation associated with GIA (18, 19). The BIFROST GPS network consists of two subnetworks: the Swedish SWEPOS network (21 receivers operating continuously since 1993) and the Finnish FinnRef network (12 receivers operating continuously since 1995), which together cover the region that is thought to be presently uplifting (Fig. 1A). For each site, three orthogonal components of position are estimated from daily solutions of the GPS measurements (19, 20). Site-specific crustal velocities have been obtained for all 34 GPS sites, and these were used to construct maps of 3D regional deformation (Fig. 1, B and C).

The pattern and amplitude of 3D velocities are consistent with published theoretical predictions of postglacial crustal deformations (11–15). The GPS-derived radial velocity field shows a broad ellipsoidal uplift dome with a major axis oriented roughly southwest to northeast. The Fennoscandian region is in active uplift, with a maximum uplift rate of ∼11.2 ± 0.2 mm/year for the site Umeå. The uncertainties in the uplift rates are higher for inland sites relative to coastal sites. This may be due to electromagnetic propagation effects associated with snow accumulation on the GPS receivers, because coastal sites experience more moderate winter temperatures (21). The horizontal velocities are relatively low where the radial uplift rates are largest (such as the central Baltic Sea), and they are directed outward from this location on all sides. In further agreement with numerical predictions, these rates increase with distance away from the uplift center, and they reach ∼1 to 2 mm/year at sites marking the perimeter of the BIFROST network. The Fennoscandian region is thus subject to widespread present-day extension. The error bars for the horizontal rates at Finnish sites are higher than those at Swedish sites because of the shorter time span of data collection in Finland.

**Earth structure**. We performed numerical predictions based on the response of a spherically symmetric (Maxwell) viscoelastic Earth model (22) to a load composed of a model of Late Pleistocene ice cover and a gravitationally self-consistent ocean load (23). The ice model incorporated a recent high-resolution reconstruction of Fennoscandian ice cover (24). The elastic structure of the Earth model was derived from seismic data (25), and the viscous structure was represented by a simple three-layer model defined by an elastic lithosphere of thickness *L*_{T} and uniform upper and lower mantle viscosities (denoted by ν_{um} and ν_{lm}, respectively), where the boundary between the viscous layers coincided with the seismic discontinuity at a depth of 670 km. We predicted 3D crustal velocities in Fennoscandia for a suite of Earth models in which *L*_{T} was set to 120 km, and ν_{um} and ν_{lm} were varied from 10^{20} Pa·s to 5 × 10^{21} Pa·s, and 10^{21} Pa·s to 5 × 10^{22} Pa·s, respectively. A χ^{2} misfit (per degree of freedom) between the predictions for each model and the GPS-derived observations was then computed (Fig. 2).

The radial velocity estimates place tighter bounds on ν_{um} than on ν_{lm}, particularly as the viscosity contrast between these two regions is increased (Fig. 2A). The misfit between the observed and predicted radial velocities increases as the upper mantle viscosity decreases below ∼5 × 10^{20} Pa·s, and thus ν_{um} values less than ∼4 × 10^{20} Pa·s are ruled out by the GPS data. This constraint on ν_{um} results from a drop in the magnitude of the predicted radial velocity as the upper mantle is weakened. As an example, the peak predicted uplift in the region drops from ∼10 mm/year to ∼1 mm/year as ν_{um} is reduced from 5 × 10^{20} Pa·s to 10^{20} Pa·s.

Misfits for the horizontal rates reflect a distinct sensitivity to ν_{um} and ν_{lm} relative to the vertical rates. The combined 3D rates are best fit by a model with ν_{um} = 8 × 10^{20} Pa·s and ν_{lm} = 10^{22} Pa·s (Fig. 2C). The 95% confidence interval for these parameters, based on an *F*-test method, is mapped out in Fig. 2C by the ranges 5 × 10^{20} ≤ ν_{um} ≤ 10^{21}Pa·s and 5 × 10^{21} ≤ ν_{lm} ≤ 5 × 10^{22} Pa·s, where a tradeoff exists such that an increase in ν_{um} requires a decrease in ν_{lm}to achieve the same level of misfit. These ranges are consistent with viscosities estimated from analyses of geologic and tide gauge markers of relative sea level change in Fennoscandia (8,9).

We performed a resolving power analysis that indicated that the 3D rates are able to resolve upper mantle viscosity but that any estimates of lower mantle viscosity are correlated with values in the upper mantle. Furthermore, the lower mantle sensitivity of the 3D rates is greatest in the top ∼1000 km of the region.

We repeated the calculations in Fig. 2 for lithospheric thicknesses of 70, 95, and 145 km and found that the best fit to the 3D data set was obtained using a thickness of 120 km. We then considered the misfit between the predicted and observed radial, horizontal, and 3D rates for a suite of Earth models in which the lithospheric thickness was varied while ν_{um} and ν_{lm} were fixed to values that best fit the 3D rates for*L*
_{T} = 120 km. The χ^{2} values for the 3D rates indicate a 95% confidence interval for*L*
_{T} of 90 to 170 km (26).

We computed maps of crustal deformation using the numerical model that best fit the 3D GPS observations (Fig. 3, A and C). A comparison of these high-resolution numerical predictions with discrete GPS-derived 3D rate estimates can be misleading. Accordingly, we also provide maps (Fig. 3, B and D) generated from numerical predictions limited to the BIFROST sites to compare with the maps derived from our GPS observations (Fig. 1, B and C). The predicted amplitude and geometry of the postglacial uplift of Fennoscandia (Fig. 3, A and B) are consistent with the observed field (Fig. 1B). In both cases, a maximum uplift rate (of ∼11 mm/year) is obtained near Umeå. The 8 mm/year and 10 mm/year contours in Fig. 1B extend further toward the northeast than do the analogous predictions in Fig. 3B; however, the observational uncertainties are also largest in this region. The pattern of uplift contours at northern sites such as Tromsø and Kevo (Fig. 1B) is due to the poor sampling by the GPS network of the vertical deformation field in this region.

The amplitude and orientation of the observed horizontal velocity vectors in Fennoscandia are also consistent with the GIA model prediction. In both cases, the horizontal speeds are small in the vicinity of the Gulf of Bothnia, and they tend to radiate outward from this area. The observed and predicted horizontal rates (Figs. 1C and 3, C and D) have greater amplitude westward of the Gulf of Bothnia than eastward of this region. This asymmetry is a consequence of several factors. The surface mass (ice plus water) load is not symmetric about the Gulf of Bothnia. Furthermore, deglaciation of Late Pleistocene North American ice sheets and rotational effects produce long-wavelength present-day horizontal deformations in Fennoscandia directed, respectively, toward the northwest (12) and east.

The Fennoscandian relaxation spectrum (27,28) provides the relaxation time of the postglacial uplift in the region as a function of the spherical harmonic degree or spatial wavelength of the deformation. The spectrum has been estimated through spectral analysis of Fennoscandian strandline data of different ages, and it represents a constraint on Earth structure that is, at least in theory, independent of ice load geometry (27,28). The Earth model that best fits the GPS-derived 3D deformation rates has a relaxation spectrum (29) that skirts the upper bound of acceptable relaxation times for degrees less than ∼30 (Fig. 4). Many viscosity models that fall within the 95% confidence ranges for ν_{um} and ν_{lm} (such as ν_{um} = 6.5 × 10^{20} Pa·s and ν_{lm} = 8 × 10^{21} Pa·s) yield relaxation spectra that fall within the observational bounds (Fig. 4).

**Regional sea level rise**. Global sea level rise is an indicator of global warming. The estimated rate of present-day sea level rise is about 1 to 3 mm/year (30). One of the primary sources of uncertainty in these estimates is the radial motion of the tide gauges due to geophysical processes such as ongoing GIA. To account for GIA motions, researchers commonly “correct” the tide gauge rates using numerical predictions based on specific Earth models and glaciation histories (31–35). This method is sensitive to the Earth model, however, even in regions within the intermediate and far field of previously glaciated areas (35, 36). Moreover, the method cannot be applied to records from tide gauges in previously deglaciated areas because of the sensitivity to details in the glaciation history and Earth structure (16).

Geodesy with GPS affords us the possibility of correcting tide gauge records for contamination due to crustal deformation, using direct measurements of radial site motions. With this goal in mind, a number of the BIFROST GPS sites were positioned near Fennoscandian tide gauges. We can write the rate of change of sea level *s˙*(37) at a particular location (longitude λ, latitude φ) as(1)where *u˙* is the radial crustal velocity and the geoid (or sea surface) rate is separated into a geographically uniform signal (μ˙) and a geographically varying term (*g˙*). The terms *u˙* and*g˙* are to be associated with any geophysical processes including, but not exclusively limited to, GIA (38).

The term *g˙* is smaller than μ˙, and little error is introduced if we correct *s˙* using a numerical prediction for *g˙* due to GIA (39). Accordingly, if we rewrite Eq. 1 as (*s˙ −g˙*) = −*u˙ *+ μ˙, then on a plot of geoid (*g˙*)-corrected tide gauge rates versus radial crustal velocities determined from GPS, the *y*intercept is the sea surface rate μ˙. We paired GPS-determined vertical rates from 20 BIFROST sites with sea level rates from the closest tide gauge sites (Fig. 5). We used only tide gauge sites with records spanning 35 years or more and used only those data acquired during or after 1930 (40). The dotted line is the expected result for μ˙ = 0. The solid line is our best estimate, μ˙ = 2.1 ± 0.3 mm/year, where the 1σ uncertainty includes the effects of possible reference frame errors and the postfit scatter of the residuals to the model (41). This line, and in particular the location of the data relative to the null result (μ˙ = 0), reflect a regionally coherent residual sea surface trend. The crustal uplift rate correction amounts to a maximum of ∼10 mm/year, which is nearly a factor of 5 greater than the “signal” (the geographically uniform sea surface rate) we determined.

**Neotectonics**. In Fig. 6, we plot residual radial and horizontal rates determined by removing from the GPS observations (Fig. 1, B and C) the GIA prediction that best fits the full 3D deformation field (Fig. 3, B and D). Although limitations in the forward model will contribute to these residuals, the maps nevertheless provide a measure of the level at which other geophysical processes, in particular neotectonics, contribute to the observed 3D velocity field.

With the exception of a few sites in northeastern Finland (such as Kevo, Sodankylä, and Kuusamo), the uplift rate residuals over the GPS network (Fig. 6A) fall within the approximate range –1 to 1 mm/year. Furthermore, the residuals between observed and predicted horizontal rates are all less than 1 mm/year, with the exception of three sites (Kuusamo, Romuvaara, and Kevo) in eastern Finland. The rather systematic east-west geographic trend evident in Fig. 6A may be due to errors in the GIA model, such as, for example, an underestimation of Late Pleistocene ice cover toward the east. Alternatively, these trends may arise from a geophysical process producing long-wavelength tilting of the Fennoscandian platform. Either interpretation requires caution, because the largest residuals in Fig. 6 occur at sites in northeastern Finland that are characterized by the largest observational uncertainties in the 3D rates. In any event, we conclude that contemporary neotectonic deformations contribute less than ∼1 mm/year to the horizontal rates in Fennoscandia (42).

The residuals in Fig. 6B show no systematic evidence for zones of shear. A site of particular relevance in this regard is Kiruna, because it is located in the vicinity of two major shear zones (the Bothnian-Senja and the Bothnian-Kvænangen zones) that have had a history of reactivation (43). Any neotectonic activity at this site is limited to the sub–millimeter-per-year level. GPS campaign data have recently been used to argue for ∼2 mm/year of strike-slip motion along a north-south–oriented zone of megashear through the middle of the Baltic Sea (extending south from Umeå to a location east of the island Visby) (44). Our results (Fig. 6B) show no evidence for an active shear zone in this region.