Research Article

Space-Geodetic Constraints on Glacial Isostatic Adjustment in Fennoscandia

See allHide authors and affiliations

Science  23 Mar 2001:
Vol. 291, Issue 5512, pp. 2381-2385
DOI: 10.1126/science.1057022


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 × 1020 to 1 × 1021 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).

Figure 1

(A) Site map showing the locations of the 21 SWEPOS GPS receivers (triangles) and the 12 FinnRef GPS receivers (circles). The IGS (International GPS Service) site Tromsø in Norway is also indicated (diamond). The network was designed so that the GPS sites are relatively evenly distributed over the Fennoscandian region. An effort was also made to place the GPS receivers in proximity to tide gauge sites. (B) Map of present-day radial velocity in Fennoscandia, constructed by fitting a polynomial function to rates obtained by GPS measurements at the BIFROST sites. The color scale is defined at the base of the plot. The locations of the GPS sites are indicated by the solid dots [see (A)], and the associated 1σ error bars (19) represent the precision of the site-specific radial velocity estimate (following the scale at the bottom right of the panel). (C) Horizontal velocity vectors estimated at each of the BIFROST sites. The scale associated with each of these vectors, as well as with the associated 1σ error ellipses (19), is given at the base of the plot.

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 LT 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 LT was set to 120 km, and νum and νlm were varied from 1020 Pa·s to 5 × 1021 Pa·s, and 1021 Pa·s to 5 × 1022 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).

Figure 2

χ2 misfit per degree of freedom (d.o.f.) between GPS-derived crustal velocities (Fig. 1, B and C) and numerical GIA predictions based on a suite of Earth models. Misfit is shown as a function of νum (ordinate scale) and νlm (abscissa scale) for the (A) radial, (B) horizontal, and (C) full 3D velocity components, respectively (45). The lithospheric thickness of the Earth models was fixed to 120 km.

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 × 1020 Pa·s, and thus νum values less than ∼4 × 1020 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 × 1020 Pa·s to 1020 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 × 1020 Pa·s and νlm = 1022 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 × 1020 ≤ νum ≤ 1021Pa·s and 5 × 1021 ≤ νlm ≤ 5 × 1022 Pa·s, where a tradeoff exists such that an increase in νum requires a decrease in νlmto 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 forL T = 120 km. The χ2 values for the 3D rates indicate a 95% confidence interval forL 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.

Figure 3

(A) Map of numerically predicted present-day radial velocity in Fennoscandia due to GIA. The calculation is based on the Earth model that provides the best fit to the GPS-determined crustal velocities in Fig. 2C (namely, a lithospheric thickness of 120 km, νum = 8 × 1020 Pa·s, and νlm = 1022 Pa·s). (B) As in (A), except that the map was constructed, as in Fig. 1B, by fitting a polynomial function to rates predicted at only the 34 sites sampled by the BIFROST network. (C) As in (A), except for predicted horizontal crustal velocities. (D) As in (C), except that the predicted horizontal velocities are shown only for the 34 sites in the BIFROST network.

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 × 1020 Pa·s and νlm = 8 × 1021 Pa·s) yield relaxation spectra that fall within the observational bounds (Fig. 4).

Figure 4

Vertical dotted lines represent the Fennoscandian relaxation spectrum (28) over the spherical harmonic degree range from 11 to 73. The solid line shows the relaxation spectrum predicted (29) using the viscoelastic model that best fits the GPS-derived 3D velocities in the BIFROST network (a lithospheric thickness of 120 km, νum = 8 × 1020 Pa·s, and νlm = 1022Pa·s). The dashed line is identical to the solid line, with the exception that a model with νum = 6.5 × 1020 Pa·s and νlm = 8 × 1021 Pa·s was adopted in the predictions.

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 (37) at a particular location (longitude λ, latitude φ) asEmbedded Image(1)where is the radial crustal velocity and the geoid (or sea surface) rate is separated into a geographically uniform signal (μ˙) and a geographically varying term (). The terms and are to be associated with any geophysical processes including, but not exclusively limited to, GIA (38).

The term is smaller than μ˙, and little error is introduced if we correct using a numerical prediction for due to GIA (39). Accordingly, if we rewrite Eq. 1 as (s˙ −g˙) = −+ μ˙, then on a plot of geoid ()-corrected tide gauge rates versus radial crustal velocities determined from GPS, the yintercept 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.

Figure 5

Rates of sea level change determined from tide gauge records at 20 sites in Fennoscandia (40), corrected for regional geoid variations due to GIA (39) versus GPS-determined radial velocities at the same (or nearby) sites. The lower diagonal line is, following Eq. 1, the result for μ˙= 0. The upper diagonal line is the best estimate through the data, and it yields μ˙ = 2.1 ± 0.3 mm/year of regionally coherent sea surface rise.

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.

Figure 6

(A) Residual vertical rates in Fennoscandia determined by subtracting from the observations (Fig. 1B) the predictions obtained using our best-fit model (Fig. 3B). (B) As in (A), except for horizontal rates (Fig. 1C minusFig. 3D).

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.

  • * To whom correspondence should be addressed. E-mail: g.a.milne{at}

  • Present address: Institute of Geodesy and Cartography, Helsinki University of Technology, FI-02015, TKK, Finland.


View Abstract

Stay Connected to Science

Navigate This Article