Sensitivity of glacial isostatic adjustment models with shallow low-viscosity earth layers to the ice-load history in relation to the performance of GOCE and GRACE

Earth and Planetary Science Letters 236 (2005) 828 – 844 www.elsevier.com/locate/epsl

Sensitivity of glacial isostatic adjustment models with shallow low-viscosity earth layers to the ice-load history in relation to the performance of GOCE and GRACE H.H.A. Schotmana,b,*, L.L.A. Vermeersena a

DEOS, Astrodynamics and Satellite Systems, Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands b SRON, National Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands Received 8 March 2005; received in revised form 23 March 2005; accepted 8 April 2005 Available online 23 May 2005 Editor: S. King

Abstract The GOCE satellite mission, which is planned by ESA for launch in August 2006, is designed to map the static global gravity field with centimeter accuracy in geoid height at 100 km or better resolution. Such a global high resolution gravity field might be able to constrain properties of shallow low-viscosity zones (LVZs) using glacial isostatic adjustment (GIA) models. In (L.L.A. Vermeersen, The potential of GOCE in constaining the structure of the crust and lithosphere from post-gracial rebound, Space Sci. Rev. 108 (2003) 105–113.) and (W. van der Wal, H.H.A. Schotman, L.L.A. Vermeersen, Geoid heights due to a crustal low viscosity zone in glacial isostatic adjustment modeling; a sensitivity analysis for GOCE, Geophys. Res. Lett. 31 (2004) 10.1029/2003GL019139.) it is shown that a crustal low-viscosity zone (CLVZ) introduces variations in geoid height up to several decimeters with spatial scales down to hundred kilometers underneath and just outside formerly glaciated areas. In (W. van der Wal, H.H.A. Schotman, L.L.A. Vermeersen, Geoid heights due to a crustal low viscosity zone in glacial isostatic adjustment modeling; a sensitivity analysis for GOCE, Geophys. Res. Lett. 31 (2004) 10.1029/2003GL019139.) it is shown that the response is sensitive to both changes in the properties of the CLVZ and the Late Pleistocene ice-load history. In this study we quantify the sensitivity to ice-load history, and investigate the effect of an asthenospheric low-viscosity zone (ALVZ) just below the lithosphere. We show, using spherical harmonic degree amplitudes, that GOCE is predicted to be sensitive to differences in the load history up to degree 130 for a CLVZ and degree 70 for an ALVZ. The sensitivity of GRACE, using the realized performance over a 111-day period (GGM01S, (B.D. Tapley, S. Bettadpur, M. Watkins, C. Reigber, The gravity recovery and climate experiment: Mission overview and early results, Geophys. Res Lett. 31 (2004) 10.1029/2004GL019920.)) is limited to lower degrees. This means that GOCE, and to a lesser degree GRACE, could provide information on the ice-load history in the presence of an LVZ, provided that our knowledge of the earth mode is sufficient. To extract information about the LVZ from the degree amplitudes, we show that it is possible to largely remove the influence of the ice-load history and thus

* Corresponding author. DEOS, Astrodynamics and Satellite Systems, Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands. Tel.: +31 15 278 2086; fax: +31 15 278 5322. E-mail address: [email protected] (H.H.A. Schotman). 0012-821X/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.epsl.2005.04.008

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

829

compute spectral signatures for different properties of the LVZ. We focus on the continental shelf-areas of western Europe to show the sensitivity of estimates of present-day sea level change to the inclusion of an LVZ. D 2005 Elsevier B.V. All rights reserved. Keywords: glacial rebound; low-viscosity layers; satellite gravity missions

1. Introduction The European Space Agency’s upcoming Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite mission, which is planned for launch in August 2006, is designed to map the static global gravity held with centimeter accuracy in geoid height at a resolution of 100 km or better [1]. Features in such a high-resolution gravity field can be associated with geophysical processes that act on a regional scale. An example of such a process is the response of shallow low-viscosity layers to loading and unloading of the crust, as for example in glacial isostatic adjustment (GIA). In GIA studies it is generally assumed that the lithosphere has very high viscosity and is therefore effectively elastic. This is a good approximation for oceanic lithosphere. In continental lithosphere however, intraplate earthquakes are mainly confined to the upper crust (0–20 km, [2,3]) and the subcrustal mantle part of the lithosphere [3], which are thought to be brittle. The lack of seismicity in the lower crust is associated with ductile flow, with a viscosity that depends on composition, fluid content and geotherm [2,3]. From reflection seismology it is found that continental regions with relatively high (N 70 mW m
1) heat flow show lower crustal reflectivity [4], which indicates lamination supported by ductile flow. From mine-induced deformation studies [5] and postseismic deformation studies (e.g. [6]) viscosities as low as 1017 Pa s are deduced. This indicates that crustal low-viscosity zones (CLVZs) can be expected in certain regions with relatively high heat flow, which excludes for example the old cratons of Canada and Scandinavia. Some post-seismic deformation studies prefer relaxation in the asthenosphere instead of the lower crust (e.g. [7]) with an effective viscosity of 1018–1019 Pa s. The possibility of ductile flow in the upper mantle is confirmed by seismic studies, which show a low-velocity zone just below the lithosphere

(see e.g. [8], p. 170). These asthenospheric low-viscosity zones (ALVZs) can be expected on a more global scale than CLVZs. Some GIA studies, focussing on the vertical displacement [9], present-day [10] and late Holocene [11] sea level change, and the global gravity field as expected from GOCE [12], have included a CLVZ at a depth of 20–35 km, with a thickness of 10–15 km and a viscosity of 1017–1019 Pa s. In [13] we have shown that a CLVZ introduces variations in geoid heights typically up to a few decimeters with spatial scales down to hundred kilometers underneath and just outside former glaciated areas. The response to changes in the properties of the CLVZ is shown to be wavelength-dependent. We have also indicated that the response might be very sensitive to the ice-load history. In this paper we investigate the sensitivity of the response to a CLVZ and an ALVZ by using three different ice-load histories: a modified version of ICE3G [14], this model with the Fennoscandian ice sheet replaced by a more recent model [15] and an ice sheet-dynamical model [16]. We will also test the sensitivity to the ocean-load history by implementing a time-dependent ocean function, and modelling meltwater in flux (MWIF) in formerly glaciated areas that are now below sea level. The effect of a time-dependent ocean function is especially important in shallow seas such as the North Sea, which fell dry during glaciation, and thus experienced a smaller load than when using present-day fixed coastlines. The effect of meltwater influx is important in areas that were once glaciated and are now below sea level, as the Gulf of Bothnia. For these areas, the total ocean load at a certain time is required, which is the negative of the topography at a certain time, and not the ocean load perturbation due to GIA only. We focus on the geoid signal in the continental shelf-areas of western Europe because these areas have relatively high heat flow [17,18], making the

830

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

presence of a CLVZ possible. Moreover, this region is just outside the margins of the maximum ice extent and is often used in studies of present-day sea level change (e.g. [19,20]), and corrections for GIA are sensitive to the inclusion of a CLVZ [10]. We show the effect on the corrections at the end of this paper. We start with an introduction of the normal-mode theory that is used to compute viscoelastic Love numbers, and of the pseudo-spectral approach to solve the sea-level equation. This is followed by a description of the used earth stratifications, and ice- and ocean-load histories. We then show that, whereas the relaxation due to a CLVZ is dominated by an extra buoyancy mode, the relaxation due to an ALVZ is guided by changes in the M0 and L0 modes. After this we will look at features of the geoid height perturbation field due to an LVZ in western Europe. We show in the spectral domain, using spherical harmonic degree amplitudes, that for a CLVZ the differences are above the expected GOCE performance up to spherical harmonic degree 130 and for an ALVZ up to degree 70. The differences are also above the realized GRACE performance over a 111-day period (GGM01S, [21]), up to degree 70 for a CLVZ and degree 50 for an ALVZ. The results for an ALVZ are comparable to GOCE if we use the more recent GGM02S (363-day period). This means that GRACE, but especially GOCE, could provide information on the ice-load history in the presence of an LVZ. In practice things are more complicated, as the measured geoid is the sum of all mass anomalies in the earth. To extract information on the ice-load history or on LVZs, all other processes that give rise to mass anomalies have to be modelled and removed. This is however not feasible at the moment, as most processes are not sufficiently well modelled (see e.g. [22]). Even in the presence of unmodelled errors, as uncertainties in the ice-load history, it is possible to extract information from the static field by some sort of spatio–spectral filtering (e.g. [23]). In this paper, we will look at a part of the spectral content of the gravity field signal induced by a CLVZ using degree amplitudes. We show that it is possible to largely remove the influence of (uncertainties in) the iceload history and thus compute spectral signatures for different properties of the CLVZ. In future studies, we will concentrate on extracting information on CLVZs in the presence of other error sources.

2. Theory 2.1. Earth relaxation model We use a semi-analytical viscoelastic relaxation model as described in [24] based on the normal mode formalism [25,26] to compute the response of a model earth to a surface load. This response is expressed in viscoelastic load Love numbers defined as [25]: X   kl ðt Þ ¼ klE dðt Þ þ ril exp sli t ð1Þ i

kEl

where is the gravitational potential perturbation elastic Love number for degree l, the summation is over all modes i of the particular earth stratification, and rli is the potential perturbation residual for the inverse relaxation time sli. For a Heaviside kind of loading, the viscoelastic Love numbers become [26]: klH ðt Þ ¼ klE þ

X rl    i 1
exp sli t l i si

ð2Þ

where the ratio rli / sli is called the modal strength. 2.2. Sea level equation We use the pseudo-spectral sea level theory of [27] to solve the sea level equation [28]. To include a timedependent ocean function and meltwater influx (MWIF) into formerly glaciated areas that are now below sea level, we use the theory as described by [29] and recently by [30]. From [30] we use Eq. (39), which is in somewhat different notation (and where the geographical coordinates are omitted):     dSj ¼ DSLj
DSLj
1 Cj Fj þ Tj
1 Cj
1 Fj
1
Cj Fj ð3Þ with dS j the incremental change in ocean height at timestep j, DSL j the change in sea level (defined globally, so also on land) from the onset of loading to time t j , T j the topography at time t j , and where C j describes the time-dependent ocean function and F j the grounding of ice. In Table 1 we have shown possible combinations of the latter for subsequent times t j
1, t j . To compute the time-dependent ocean function C j we need the paleotopography at time t j . This is estimated from the present-day topography and the sea level change from that time to the present. As

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844 Table 1 Operator F j that describes the grounding and floating of ice Fj = 0 Fj = 1

Fj
1 = 0

Fj
1 = 1

Ice growth/melt Influx of water

Grounding of ice No ice growth/melt

the latter is not known, but is computed in the process, the computation is iterated over a full glacial cycle. The operator F j is determined by comparing the mass per unit area of the ocean at a certain grid point with that of the prescribed ice-load (as taken from a certain ice-load history) at that point.

3. Model parameters 3.1. Earth stratification The continental shelf areas in northwestern Europe have a crustal thickness of 30–35 km [31]. The North Sea area has a relatively high heat flow of 60–80 mW/ m2 [17,18], which makes a crustal low-viscosity zone (CLVZ) of 1018 Pa s possible. We choose a model with a lower crust starting at a depth of 20 km, with a thickness of 12 km and a viscosity of 1018 Pa s (model Md20t12v18). The total lithospheric thickness is 80 km. This model seems to be more appropriate for this coastal shelf area as the models used by [10–12], which have a deeper and thicker CLVZ. We have to remark that this CLVZ is an extreme end-member in other parts of this area, as earthquakes as deep as 25–30 km are found in the continental shelf area off the coast

831

of Norway [32]. Moreover, a CLVZ in the neighbouring Fennoscandian shield, with its cold and thick crust (see e.g. [18]) is highly unlikely. In oceanic lithosphere, earthquakes occur from a depth of 5–40 km [2,33], which indicates brittle behaviour of both the crustal and mantle part of the lithosphere. Below the mantle part of the lithosphere, i.e. the asthenosphere, seismic data show a low-velocity zone (see e.g. [8], p. 170). This low-velocity zone is associated with ductile flow, as deformation laws, the convergence of the geotherm and the solidus, and the strong vertical advection of deep heat from mantle convection predict low-viscosity material in this area (e.g. [8], p. 204). Post-seismic deformation, levelling and shoreline elevation studies in the western US, which has a slightly higher heat flow than western Europe, give upper mantle viscosities in the range of 1017–1019 Pa s [34]. We have modelled an asthenospheric low-viscosity zone (ALVZ) below a fully elastic lithosphere of 80 km, with a thickness of 35 km and a viscosity of 1018 Pa s (model Ml80t35v18). As we are only interested in geoid height perturbations due to an LVZ, we subtract the results of a reference model that has the same number of layers as the model with an LVZ, but with a viscosity value that is either very large (the 7-layer model Md20t12ela, in case of a CLVZ) or equal to the upper mantle viscosity (the 6-layer model Ml80t35vum, for an ALVZ). Both these models are comparable to the 5layer volume-averaged models as used by [24] and [35]. The viscosity stratification of our models is given in Table 2. We take volume-averaged elastic parameters from PREM [36].

Table 2 Viscosity stratification of our earth models Layer

Depth (km)

Viscosity (Pa s) Md20t12v18

Lithosphere

Upper mantle

Lower mantle Core

0–20 20–32 32–80 80–115 115–400 400–670 670–2891 2891–6371

1d 1d 1d 5d

50

10 1018 1050 1020

5 d 1020 5 d 1021 0

Md20t12ela 1d 1d 1d 5d

50

10 1050 1050 1020

5 d 1020 5 d 1021 0

Ml80t35v18 50

1 d 10

1d 5d 5d 5d 0

1018 1020 1020 1021

Ml80t35vum 1 d 1050

5d 5d 5d 5d 0

1020 1020 1020 1021

832

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

3.2. Loading history We use three different ice-load histories: a modified version of ICE3G of Tushingham and Peltier [14], this model with the Fennoscandian and Barents Sea part replaced by a model developed by Lambeck et al. from the Australian National University (ANU) [15] and an ice sheet-dynamical (ISD) model of the Institute for Marine and Atmospheric Research Utrecht (IMAU, Utrecht University), comparable to the model described in [16]. Both ICE3G and the ANU model are constrained using geomorphological data, sea level curves and assumptions on the viscosity structure of the earth [14,15]. We have smoothed the ICE3G model to remove holes that arise because of the finite disc definition of this model. The smoothing is accomplished by replacing zero ice heights with the average of the nearest neighbours if at least three of the neighbouring points are non-zero. In this way, the extent of the ice sheet is preserved, which is not the case when a filter is applied in the spectral domain (as e.g. in [13]). A disadvantage is the increase in volume, which is however of minor importance, as we have scaled the model, except for the Fennoscandian and Barents Sea parts, to include a total ice-equivalent sea level of 120 m at the last glacial maximum (LGM). We will call this modified ICE3G model dI3GT to distinguish it from the original, and have plotted the ice heights at LGM of this modified ICE3G model for the northern hemisphere in Fig. 1a. We have complemented the ANU model for Fennoscandia and the Barents Sea with scaled British Isles, Kara Sea, East-Siberia, North-America, Greenland and Antarctica parts of ICE3G to have a total volume at LGM of 120 m (Fig. 1b). We have to remark that for ease of implementation and comparison, we have made some assumptions which alter the ANU model substantially. Firstly, we have assumed that ICE3G and ANU are given as uncalibrated C14 years before present (kyr bp). This is however not the case, as the ANU model is constructed using a calibrated C14 time scale, which leads to an error in timing. In this paper, LGM for both models is assumed to be at 18 kyr bp. Moreover, we have implemented the ANU model as a series of Heaviside step functions and not as stepwise linear, as originally developed. Finally, in the complete (i.e. global) model produced by Lambeck et al., the Kara

Sea area is predicted to be largely ice-free and the total ice-equivalent sea level is estimated to be 140 m (Tony Purcell, personal communication). The ISD model is constructed using paleo-temperatures for the atmosphere, a global circulation model, eustatic sea level data and assumptions on the flow of ice [16]. We have used only the North-American and Eurasian part of this model, complemented with scaled Greenland and Antarctica parts of ICE3G to have a total volume of 120 m. The Eurasian part includes the British Isles, Fennoscandia, and the Barents and Kara Sea area (Fig. 1c). In Fig. 1d we show the time development of the Fennoscandian and Barents Sea ice sheet (FEN) for ANU and I3G (modified ICE3G), and the total of the Eurasian ice sheets (EAS, includes FEN) for all models. The large volume of the ISD model is mainly confined to the Barents Sea and Kara Sea area, where the ice-dynamical model predicts a huge ice sheet, see Fig. 1c and [16]. The large ice sheet over the Kara Sea is probably not very realistic, as this area seems to have been largely ice-free at LGM [37]. The ANU and I3G model are very similar, with however a smaller mass for Fennoscandia in the ANU model, and an earlier ending of deglaciation (note that both might be due to the assumption that ANU is given in an uncalibrated C14 time scale). As the ICE3G and ANU ice-load histories are computed using the same type of information, we regard differences in the results between the modified ICE3G (I3G) and ANU models as an indicator of the upper limit of the sensitivity of the GIA model, or a lower limit of the error introduced by the use of a certain ice-load history. The difference between I3G and the ISD model can then be regarded as an indicator of the lower limit of the sensitivity of our GIA model, or an upper limit of the error introduced by the use of a certain ice-load history, because they are based on mutually independent information and methods. In this study, all ice-load histories are defined from 18 kyr bp with a time step of 1 kyr. We interpolate to a time step of 0.5 kyr in the deglaciation phase, and precede the given load history with a linear glaciation phase of 90 kyr with a time step of 5 kyr, except for 10 steps of 0.5 kyr just before LGM. We solve our models on a grid of 384 lines in latitude and 512 lines in longitude up to degree 256, which is sufficient for

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

(b) ANU + I3G

(a) I3G (modified ICE3G)

0

500

833

1000 1500 2000 2500 3000

0

500

ice height at LGM [m]

1000 1500 2000 2500 3000

ice height at LGM [m]

(c) ISD + I3G

(d) deglaciation history

ice–equivalent sea level (m)

60

40 30 20 10 0 18

0

500

1000 1500 2000 2500 3000

I3G, FEN I3G, EAS ANU, FEN ANU, EAS ISD, EAS

50

16

14

12

10

8

6

4

2

time (kyrs bp)

ice height at LGM [m] Fig. 1. Ice heights at LGM in northern Europe of the modified ICE3G model (I3G, a), the modified ICE3G model with the Fennoscandian part replaced by the ANU ice-load model (ANU + I3G, b) and the ice sheet-dynamical model with modified ICE3G outside NAM and EAS (ISD + I3G, c), and the deglaciation history for the three ice models (d).

our purposes. The land–sea mask and the topography at present, as needed for the computation of MWIF, are derived from the ETOPO5 data set.

4. Effect of an LVZ on relaxation spectra and geoid heights 4.1. Changes in relaxation spectra The modal relaxation times and strengths for model Md20t12v18 (dots) and Md20t12ela (squares)

are plotted in Fig. 2a. The introduction of a CLVZ induces 1 extra buoyancy mode and 3 extra viscoelastic modes, as already demonstrated by [9] for their halfspace model. The extra buoyancy mode, labelled MC, has significant strength and is even stronger than M0 above degree 50, where the associated relaxation time is tens of kiloyears. The typical behaviour of the MC mode, with large relaxation times for small harmonic degree, is probably because adjustment for long-wavelength loads to isostasy is slower for channel flow than for mantle flow (see e.g. [38], p. 158). Note also the increase in

834

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

(a) Modal relaxation times, CLVZ

(b) Modal relaxation times, ALVZ

6

6

10

dots: Md20t12v18 solid line: total Md20t12v18 squares: Md20t12ela dashed line: total Md20t12ela

MC

4

10

10

M3

4

10

M2

M2 2

10

relaxation time [kyr]

relaxation time [kyr]

M1

C0 M0

0

MC

L0

10

dots: Ml80t35v18, solid line: total Ml80t35v18 squares: Ml80t35vum, dashed line: total Ml80t35vum

2

C0

10

M1 M0 (vum) 0

10

T1(2),T2(2) M0 (v18)

T1(2),T2(2) LC

-2

10

-2

L0

10

TC(2) -4

10

T3(2)

-4

0

50

100 150 harmonic degree

200

10

250

0

(c) Modal strengths, CLVZ

50

100 150 harmonic degree

200

250

(d) Modal strengths, ALVZ

0

0

10

10 dots: Md20t12v18 squares: Md20t12ela

MC

M0

dots: Ml80t35v18 squares: Ml80t35vum

L0 -1

10

modal strength

modal strength

-1

M0 (v18)

M0 (ela) -2

10

50

L0 (v18) M0 (vum)

-2

T3

solid line: difference in total strength between Md20t12v18 and Md20t12ela 0

L0 (vum)

10

LC

-3

10

10

100 150 harmonic degree

200

M0 (v18)

-3

250

10

0

50

100 150 harmonic degree

200

250

Fig. 2. Modal relaxation times and modal strengths for a CLVZ (a, c) and an ALVZ (b, d).

the relaxation time of the M0 mode. The viscoelastic LC mode, due to the difference in Maxwell time between the upper and lower crust, and the two transient TC modes, due to the difference in Maxwell time of the lower crust and the lithosphere, have strengths up to an order of magnitude smaller than the dominant modes and might be interesting for fast, high degree relaxation as in post-seismic deformation. We have computed the total relaxation time as a function of harmonic degree by summing the modal relaxation times weighted with their relative modal strengths. The total relaxation time for Md20t12v18 (solid line) and Md20t12ela (dashed line) is also

plotted in Fig. 2a and shows that the relaxation time is significantly increased for all harmonic degrees. In Fig. 2c we have plotted the absolute value of the difference in total strength between Md20t12v18 and Md20t12ela, which shows that the difference in response is mainly confined to degree 50 to 150. It also shows that our normal-mode method has very high accuracy, as the curve is smooth. This is not trivial because for large relaxation times the system of differential equations is unstable, as the sixth-order system degenerates to a second order system in the fluid limit [24,26]. This fluid limit is equal to the plotted total strength plus the elastic limit; compare Eq. (2).

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

10˚

0˚

–10

30˚

20˚

–2 0˚

˚

(a) Total GIA induced 65˚

60

˚

60˚

–1

1

–0

55˚

.5

0

55˚

50˚ 0.5

50˚

45˚ 45˚ 1

The modal relaxation times and strengths for Ml80t35v18 (dots) and Ml80t35vum (squares) are given in Fig. 2b and d. In contrast to the model with a CLVZ, the extra buoyancy mode (M3) has no significant strength. The response due to an ALVZ seems to be guided by changes in both the relaxation time and strength of the M0 and L0 modes, leading to a smaller total relaxation time (Fig. 2b). The difference in total strength between Ml80t35v18 and Ml80t35vum is zero as the models show identical behaviour in the fluid limit, and is therefore not shown in Fig. 2d.

835

40˚

4.2. Perturbations in geoid height -6

-5

-4

-3

-2

-1

0

1

2

geoid height [m]

10˚

65˚

–1 0

0˚

–10

–2

30˚

60

20˚

0˚

˚

(b) CLVZ induced perturbation ˚

0

0

10

60˚

10

55˚

10

10

55˚

0

–10

50˚

0 –1 0

10

0

10 0

0

–

50˚

45˚ 45˚ –10

40˚

40˚

-60

-50

-40

-30

-20 -10

0

10

20

30

geoid height perturbation [cm]

10˚

0˚

–10

30˚

65˚

˚ 0.5

60

20˚

–2 0˚

˚

(c) ALVZ induced perturbation

–1

.75

–0

60˚

–0.5

55˚

–0.25

55˚

0.25 0

The results of our GIA model show that the total GIA induced geoid height (using I3G) is dominated by a deep low at the center of rebound, see Fig. 3a. The white line in this figure indicates the continental margin. The North Sea continental shelf area is located at the collapsing forebulge, which has roughly its peak where the geoid height changes its sign. The geoid height perturbations induced by a CLVZ (model Md20t12v18, with the results of model Md20t12ela subtracted) are shown in Fig. 3b. The induced signal shows amplitudes of tens of centimeters with spatial scales down to 100 km. The pattern is mainly explained by the extra flow away from the formerly glaciated areas through the low-viscosity channel during the glaciation phase, leading to differential forebulges around the ice sheets. Due to the relatively long relaxation time of the MC mode (Fig. 2a) and the short deglaciation period, these areas have only slightly adjusted to isostasy and thus show still a positive perturbation. The position of the ice sheets can clearly be distinguished in these figures, which already indicates the sensitivity to ice-sheet geometry. An important farfield phenomenon is the extra tilting of the coastlines due to enhanced material flow forced by the ocean load [10,11]. This leads to negative perturbations on the land-side of the coastline, due to the negative sea-load in the glaciation period and the much shorter deglaciation period.

40˚

50˚

50˚

–0.25

45˚ 0

45˚

40˚ 40˚

Fig. 3. Total GIA induced geoid heights (a) as predicted by our GIA model, and CLVZ (b) and ALVZ (c) induced geoid height perturbations. The white line indicates the continental margin.

-2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

geoid height perturbation [m]

1.0

1.5

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

10˚

0˚

–10

0˚

0˚

10˚

˚ –10

–2

(b) ALVZ, difference between ANU and I3G 30˚

0

60 ˚

65˚

0

–1 0

55˚

0

60˚

10

–1

20

10

60˚

0

10 0

55˚

We test the sensitivity of our GIA model with a CLVZ to ice-load scenarios by subtracting the geoid height perturbations induced by using ice-load history I3G from the geoid height perturbations using ice-load history ANU or ISD. Results are then geoid height perturbation differences as plotted in Fig. 4a for the difference between ice-load histories ANU and I3G,

–10

65˚

–10

60 ˚

5.1. Ice-load sensitivity

20˚

30˚

20˚

0˚

(a) CLVZ, difference between ANU and I3G

5. Sensitivity to the load history

˚

The effect of an ALVZ, using Ml80t35v18 with the results of Ml80t35vum subtracted, is shown in Fig. 3c. Due to the larger depth of the LVZ, the pattern is much less dependent on the fine-scale structure of the ice load. Moreover, instead of extra material flow away from the formerly glaciated areas, the differential geoid signal is dominated by extra material flow towards the center of rebound, thus accelerating the process of isostatic adjustment. This can be explained by the shorter relaxation times of the M0 and L0 mode for model Ml80t35v18 compared with Ml80t35vum, see Fig. 2b.

–2

836

20

0

50˚

0

55˚

10

50˚

20

55˚

0

–10

50˚

0

45˚

50˚

45˚

20

45˚

40˚

45˚

10

40˚

40˚

-10

0

10

20

-100

–2 0˚

0˚

10˚

˚ –10

0

–10

0 –1

0

60˚ –10

55˚

55˚

0 –2

50˚

65˚

0

10

20

10

55˚

0

0 –1 0

0 –1

20

10

0

10

0

60 ˚

60˚

–10

-20

(d) ALVZ, difference between ISD and I3G

0 –1

55˚

-40

30˚

65˚

60 ˚

-60

20˚

30˚

20˚

–2 0

˚

(c) CLVZ, difference between ISD and I3G

-80

geoid height perturbation difference [cm]

10˚

-20

0˚

-30

˚

-40

geoid height perturbation difference [cm]

–10

-50

40˚

50˚ 50˚

45˚

50˚

45˚ 45˚

0

40˚

45˚

40˚

40˚

-50

-40

-30

-20

-10

0

10

20

geoid height perturbation difference [cm]

40˚

-100

-80

-60

-40

-20

0

20

geoid height perturbation difference [cm]

Fig. 4. CLVZ induced geoid height perturbation differences between ANU and I3G (a), and ISD and I3G (c), and ALVZ induced differences between ANU and I3G (b), and ISD and I3G (d).

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

and in Fig. 4c for the difference between ISD and I3G. A common feature in both plots is the high in areas where the ANU and ISD model define less ice than ICE3G (e.g. southern tip of Norway and Sweden, compare Fig. 1) and just outside the maximum extent of ANU and ISD (e.g. Denmark), where the differential forebulge is pushed southward. Another common feature is the large slightly negative area around the former ice sheets, which can be explained by the smaller mass of the Fennoscandian ice sheet for ANU and ISD, see Fig. 1d, and the slower collapse of the forebulge as explained in Section 4.2. The sensitivity for a model with an ALVZ is shown in Fig. 4b and d. The difference in perturbations between ANU and I3G is showing that the lows are less deep and the highs are less high for the ANU model, mainly due to its smaller mass for Fennoscandia; compare Fig. 1. Outside the once-glaciated areas we still find amplitudes larger than 20 cm. For the difference between the ISD and I3G model the picture is completely different because of the large ice mass in the Barents and Kara Sea areas in the ISD model. This gives a differential high in those areas, and a differential low over large parts of Europe.

837

and became gradually flooded during the deglaciation period. This effect is described by the time-dependent ocean function in the implementation of MWIF and results in negative differences in CLVZ induced geoid height perturbations between a model that includes MWIF and a model that does not. This can be seen in Fig. 5a, which is produced using I3G. Due to the prevailing negative ocean-load, the area will show a negative geoid height perturbation when defined as land and a positive perturbation when defined as ocean; compare with the tilting effect as described in Section 4.2. The positive patch in the Gulf of Bothnia is due to the fact that in the model without MWIF, it is possible that areas are both ice and ocean loaded, which has a large effect in the glaciation phase. As a result, the model with MWIF gives smaller negative perturbations at LGM than the model without MWIF. The picture for an ALVZ (Fig. 5b) is smoother due to the greater depth of the LVZ, and more or less of opposite sign, as explained in Section 4.2.

6. Results in the spectral domain

5.2. Ocean-load sensitivity

6.1. Comparison with the GOCE and GRACE performance

During the last glacial cycle, shallow seas such as the North Sea fell dry at a certain moment in the past,

The results we have shown so far are regional results from global computations. To compare our

0 0

10˚

0˚

0

0

50˚

50˚

5

45˚ 0

0

45˚

45˚

0

40˚

40˚ 0

-20 -15 -10

55˚

50˚

0 0

˚

55˚

10

0

45˚

–10

0˚

5 –5

0

0

60˚ 55˚

–5

–5

50˚

65˚

60 ˚

60˚

0

–5

55˚

–2

10˚

0˚

–10 0

30˚

65˚

60 ˚

20˚

–2

(b) ALVZ, effect of MWI

30˚

20˚

0˚

˚

(a) CLVZ, effect of MWI

-5

0

40˚

40˚

5

10

15

20

geoid height perturbation difference [cm]

-20 -15 -10

-5

0

5

10

15

20

geoid height perturbation difference [cm]

Fig. 5. CLVZ (a) and ALVZ (b) induced geoid height perturbation differences between a model that includes MWIF and a model that does not.

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

results with the performances of the GOCE and GRACE satellite mission we need to use global results, as the performances are defined in spherical harmonics, which give only global information. We therefore transform global geoid height perturbations and perturbation differences to the spectral domain and compute degree amplitudes, defined as the square root of the degree variances: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u l uX  2 þ S2 cl ¼ t Clm lm

ð4Þ

(a) perturbations compared with gravity models 2

10

geoid height degree amplitude (cm)

838

EGM96 GGM01S GOCE CLVZ AVLZ

1

10

0

10

-1

10

-2

10

m¼0 -3

10

0

50

100 150 harmonic degree

200

250

(b) CLVZ induced perturbation differences 2

geoid height degree amplitude (cm)

10

ANU – I3G ISD – I3G MWIF – no MWIF EGM96 GGM01S GOCE

1

10

0

10

-1

10

-2

10

-3

10

0

50

100 150 harmonic degree

200

250

(c) ALVZ induced perturbation differences geoid height degree amplitude (cm)

with (C lm , S lm ) a set of spherical harmonic (Stokes) coefficients of the geoid height perturbations or perturbation differences. Note that we are using only part of the information that is available in our results, because for each degree we sum over all orders. This simplifies comparisons but neglects a large amount of longitudinal information. The performances of GOCE and GRACE as given below are for the same reason too pessimistic. In Fig. 6a we have plotted the CLVZ (circles) and ALVZ (pluses) induced geoid height perturbation degree amplitudes together with the error characteristics of EGM96 ([40], dotted), the realized performance of GRACE for a 111-day period (GGM01S, [21], dashed) and the expected performance of GOCE based on formal errors (Pieter Visser, personal communication, solid). The degree amplitudes of the CLVZ induced perturbations are above the performance of GOCE up to degree 130 (above GGM01S up to degree 80) and of the ALVZ induced perturbations up to degree 60 (GGM01S: 50). This means that a large part of the information about LVZs is contained in the GOCE signal and, to a lesser degree, in the GRACE signal. The more recent GRACE model GGM02S improves the results for an ALVZ to the level of GOCE and slightly improves the results for a CLVZ. We are interested if GOCE and GRACE are also sensitive to changes in the used ice- and ocean-load history in the presence of a CLVZ. In Fig. 6b we have plotted the degree amplitudes of CLVZ induced geoid height perturbation differences between ANU and I3G (circles), ISD and I3G (pluses), and a model with and without MWIF (crosses). We see that GOCE is sensitive to the difference between ISD and I3G up to degree 130 (GGM01S: 80) and for the difference between

ANU – I3G ISD – I3G MWIF – no MWIF EGM96 GGM01S GOCE

1

10

0

10

1

10

-2

10

0

50

100 150 harmonic degree

200

250

Fig. 6. Geoid height perturbation degree amplitudes for a CLVZ (Md20t12v18–Md20t12ela) and an ALVZ (Ml80t35v18–Ml80t35vum) compared with existing and future gravity models (a), and CLVZ (b) and ALVZ (c) induced perturbation differences for different load histories.

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

(a) degree correlation coefficients 1

clvz(I3G), I3G clvz(ANU), I3G clvz(ISD), I3G

0.9 degree correlation coefficient

ANU and I3G up to degree 100 (GGM01S: 60). The lower sensitivity to the difference between ANU and I3G is due to the fact that degree amplitudes give global spectral information and ANU and I3G differ only significantly for northern Europe. The sensitivity to the inclusion of MWIF is comparably low due to the mainly local effect near former ice sheets. GOCE and GRACE are however sensitive to the inclusion of MWIF, which means that this effect has to be taken into account in the computations. In the presence of an ALVZ (Fig. 6c) the sensitivity is confined to lower degrees, where GOCE is sensitive up to degree 70 and GRACE up to degree 50.

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

6.2. Spectral signatures

100 150 harmonic degree

200

250

normalized geoid height degree amplitude

-1

10

-2

10

-3

10

-4

I3G ANU ISD best estimate

alm blm ð5Þ

m¼0

with a lm a set (C lmA , S lmA ) of spherical harmonic coefficients and b lm a set (C lmB , S lmB ). The absolute value is equal to 1 if set A and set B are linearly dependent for a certain degree, and 0 if there is no correlation between the two. We can now use the degree correlation coefficient to determine the most likely actual ice-load history. The degree correlation coefficient between the Stokes coefficients of CLVZ induced geoid height perturbations (using I3G) and the different ice-load histories (I3G, ANU, ISD) at LGM is shown in Fig. 7a. It is

0

50

100 150 harmonic degree

200

250

(c) spectral signatures, different earth models normalized geoid height degree amplitude

m¼0

10

10

m¼0 ql ¼ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u l l X uX t alm aTlm blm bTlm

50

(b) spectral signatures, different ice loading histories

We have shown in Fig. 6b that if we know the properties of the CLVZ, we should be able to discriminate between different ice-load histories using GOCE, and to a lesser extent, GRACE. However, we do not know the properties of the CLVZ and moreover, from a geophysical point of view, we are more interested in constraining the properties of the CLVZ using GOCE. The differences between ice-load histories in CLVZ induced geoid height perturbations can then be regarded as an error source. We show that it is possible to extract information from GOCE measurements about a possible CLVZ without exact information on the ice-load history. For this purpose we will use the degree correlation coefficient (see e.g. [39]) defined as: l X

839

-1

Md20t12v18 Md20t20v18 Md20t12v19 Md32t08v18

10

10

10

-2

-3

-4

0

50

100 150 harmonic degree

200

250

Fig. 7. Degree correlation coefficients (a) between CLVZ induced geoid height perturbations (using I3G) and different ice-load histories (I3G: solid line, ANU: circles, ISD: pluses), spectral signatures of the CLVZ for different ice-load histories (b), and best estimates of the spectra signature for different properties of the CLVZ (c).

840

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

clear that the degree correlation coefficients always prefer the ice-load history that is used in the computation of the CLVZ perturbations, even if the spectra are close together as for I3G and ANU. In reality, we do not know the actual ice-load history, but we can state that the history which has the highest correlation with the measured geoid heights is the most likely actual history. We can then normalize the geoid height perturbations with the degree amplitudes of the most likely actual history to obtain spectral signatures. In Fig. 7b we show that the spectral signatures of a CLVZ for different ice-load histories are very close together, i.e. we have largely removed the influence of (uncertainties in) the ice-load history. Moreover, if we use the degree correlation coefficients as a measure for the error in the computed signature, we can compute a least-squares fit to the curves, weighted with the degree correlation coefficients. The errors in the signatures are defined in such a way that the error is 0 if the correlation is 1, equal to the standard deviation of the signatures if the correlation is 0.5, and infinitely large if the correlation is 0. We have also plotted this best estimate of the spectral signature in Fig. 7b for values that are above the propagated error. This propagated error is computed from the error in the spectral signatures as defined above, using standard leastsquares techniques. Finally, we have plotted estimates of spectral signatures for different properties of the CLVZ in Fig. 7c, where Md20t20v18 has a thicker CLVZ (20 km), Md20t12v19 has a higher viscosity CLVZ (1019 Pa s) and Md32t08v18 both a deeper (at 32 km) and thinner (8 km) CLVZ. As the signatures differ significantly, it should be possible to extract information about the properties of a CLVZ from GOCE measurements by following the above-mentioned procedure. That is, estimate a most likely actual ice-load history from degree correlations with the measured geoid heights, use this information to largely remove the influence of the load history to obtain a spectral signature, and compare this measured signature with a library of signatures for different properties of the CLVZ. Note that the spectral signature for our reference model (Md20t12v18) is roughly the same as the total difference in modal strength, compare Fig. 2c. Actually, if we put a load on the earth for an infinitely long time, the spectral signature will be equal to the fluid Love number. Equally, for transient time scales the

spectral signature will be equal to the Love number on transient time scales. This can be understood if we view our GIA model as a linear shift-invariant, discrete-space system (compare [41], p. 11): yðnÞ ¼ hðnÞTxðnÞ

ð6Þ

where * denotes convolution of the unit-sample response h(n) (i.e. the Green’s function) with the input x(n) (i.e. the load history) to give the output y(n) (e.g. geoid heights), with n a set of geographical coordinates. The spectral form of this equation is ([41], p. 68): Y ðmÞ ¼ H ðmÞX ðmÞ

ð7Þ

where m is a frequency associated with a certain base function (in this case degree and order of a spherical harmonic expansion) and H(m) is called the system function. As our system is both a function of geographical coordinates and of time, we can regard the Love number as the system function of our GIA model for a certain time scale. For a load that is constant in time, the spectral signature is also equal to the system function for a certain time scale as it is the square-root of the response power density (i.e. degree variance of the geoid height perturbations) divided by the forcing power density (i.e. degree variance of the load); compare [41] p. 393: jH ðmÞj2 ¼ /yy ðmÞ=/xx ðmÞ

ð8Þ

As the load that we have used in this study is phased in time, we can regard the spectral signature as an approximation of the system function of our GIA model for a certain transient time scale. For an example how transient transitions take place in a layered earth we refer to [44].

7. Effect of an LVZ on present-day sea level change In global warming studies, estimates of presentday sea level change (SLC) from tide gauges are corrected for ongoing GIA, see e.g. [19,20,42,43]. If a certain GIA model decreases the standard deviation globally or in a certain area, then the corrected estimate is regarded as a better one than the original. In this way, the estimates can be improved without knowledge about other processes that influence esti-

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

10˚

0˚

–2 0˚

65˚

60 ˚

60˚ 55˚ 55˚

50˚ Ne

50˚

Br

45˚ Tr Ma

40˚

45˚

Ge

Ca

-0.75

40˚

-0.50

-0.25

0.00

0.25

0.50

0.75

rsl rate [mm/yr]

10˚

0˚

˚ 20

30˚

20˚

1 0˚

(b) CLVZ 65˚

60 ˚

60˚ 55˚ 55˚

50˚ Ne

50˚

Br

45˚ Tr Ma

40˚

45˚

Ge

Ca

-0.75

40˚

-0.50

-0.25

0.00

0.25

0.50

0.75

rsl rate [mm/yr] 10˚

0˚

0˚ –2

30˚

60

20˚

–10 ˚

(c) ALVZ 65˚

˚

60˚ 55˚ 55˚

50˚ Ne

50˚

Br

45˚ Tr Ma

40˚

Fig. 8. GIA induced rate of relative change level rise for an earth model without an LVZ (Md20t12ela, a), a CLVZ (Md20t12v18, b) and an ALVZ (Ml80t35v18, c). The white line indicates zero rsl rate.

30˚

–10 ˚

(a) no LVZ 20˚

mates of SLC as tectonics. The correction is very sensitive to the ice-load history and the earth model. Here we will not argue for a specific earth and ice model combination, but merely show the effect of an LVZ on the predicted rate of present-day SLC in western Europe. We show that for this area the standard deviation of the corrected tide gauge measurements can be smaller for a model with an LVZ than for a model without an LVZ. As we are using a relatively small set of tide gauge measurements, these results should be regarded as a possible way to explain short-scale differences in tide gauge measurements, not as better corrections than proposed by other authors (as e.g. [43]). In Fig. 8a we show the predicted rate of relative SLC due to GIA in western Europe, computed using earth model Md20t12ela (without an LVZ) and iceload history I3G. The main features are a large positive rate of relative SLC (N 1.5 mm/yr) due to collapse of the forebulge, a large negative rate of relative SLC (up to ~ 10 mm/yr in Scandinavia) due to rebound of the formerly glaciated areas, and a small positive or negative rate near the coasts depending on the location of a site relative to the long-wavelength trend of the coastline (see e.g. [11]). In Fig. 8 we have also indicated the 6 tide gauge stations that we use (Ne: Newlyn, Br: Brest, Ma: Marseille, Ge: Genoa, Tr: Trieste and Ca: Cascais). SLC values estimated from these tide gauges are taken from 2 studies of Douglas (1991: [19], 2001: [42]). As the standard deviation is a measure for the appropriateness of the GIA correction, we work with the mean (dMEANT) of the 6 estimates and the standard deviation (dSIGMAT), which are shown in the column dTrendT of Table 3. In the next column (dTrendGIAT) we show the estimates of SLC after correction for ongoing GIA using a certain earth and ice model (1991: ICE3G [14,19], 2001: ICE4G [43]). Note that for 2001 the standard deviation is increased, meaning that the correction is not appropriate for this region. In the third column (dNo LVZT) of Table 3 we use the values from Fig. 8a and see that our correction gives a smaller mean with a comparable (1991) or smaller (2001) standard deviation.

841

45˚

Ge

Ca

-0.75

40˚

-0.50

-0.25

0.00

0.25

rsl rate [mm/yr]

0.50

0.75

842

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

Table 3 Rates of present-day sea level change corrected for GIA (in mm/yr) in western Europe

MEAN SIGMA MEAN SIGMA

1991 1991 2001 2001

Trend

Trend-GIA

No LVZ

CLVZ

ALVZ

1.42 0.22 1.35 0.24

1.65 0.19 1.24 0.36

1.32 0.21 1.26 0.25

1.31 0.19 1.21 0.23

1.18 0.24 1.13 0.22

Fig. 8b shows the GIA induced relative SLC if a CLVZ is incorporated in our earth model (Md20t12v18), showing a pattern of relatively large gradients along coastlines due to extra tilting (see Section 4.2 and e.g. [11]). The effect on the rate of SLC in Europe is to lower both the mean and the standard deviation (dCLVZT, Table 3), the first with a maximum of 4% (from 1.26 to 1.21 mm/yr) and the latter with 10%. The effect of an ALVZ is smoother and tends to increase the GIA induced rate of relative SLC, due to the faster collapse of the bulge area (Fig. 8c). The effect of an ALVZ is larger than that of a CLVZ, with a reduction in the mean rate of about 10% (dALVZT, Table 3). For the tide gauge estimates of 2001 this leads to a mean SLC rate corrected for GIA of 1.13 mm/yr (1.26 mm/yr for the model without an LVZ) and a standard deviation of 0.22 mm/yr (0.25 mm/yr for the model without an LVZ), which is the best estimate for 2001 in terms of standard deviation.

8. Conclusions From seismology, laboratory experiments and postseismic deformation studies there is substantial evidence for the existence of regional crustal and asthenospheric low-viscosity zones (LVZs). The properties of these LVZs (depth, thickness, viscosity) and their lateral variation are uncertain. We have modeled a crustal low-viscosity zone (CLVZ) in western Europe, which might be reasonable for the continental shelf areas outside the Fennoscandian shield. In the oceanic lithosphere and the shield area, a CLVZ is unlikely, though the existence of an asthenospheric low viscosity zone (ALVZ) is possible. We have shown that the response of a model earth with a shallow LVZ to loading and unloading during the last glacial cycle is a complex function of sphe-

rical harmonic degree and time. For a CLVZ the relaxation time and strength are increased for all degrees due to the extra buoyancy mode MC. This leads to extra material flow away from the (formerly) loaded areas, especially for degrees 50–150. Coastlines experience additional tilting, thus lowering estimates of present-day sea level change. For an ALVZ the relaxation time is decreased for all degrees, leading to faster adjustment to isostasy and lower estimates of present-day sea level change. For both LVZs, the standard deviation of the estimates is decreased for a region like western Europe. The response is also very sensitive to the prescribed ice-load history and to a lesser extent the ocean-load history. We have tested this sensitivity by using three ice-load histories and the implementation of meltwater in flux (MWIF) in formerly glaciated areas. In the spatial domain we have shown that the sensitivity to the ice-load history is mainly due to ice sheet geometry and volume, and that the effect is confined to an area of a few hundred kilometers around the formerly glaciated areas. The implementation of MWIF gives mainly improvements for shallow seas such as the North Sea and sea areas that were once glaciated as the Gulf of Bothnia. We have compared our results with the expected performance of GOCE and the realized performance of GRACE (GGM01S). CLVZ induced geoid height perturbation degree amplitudes are above the former up to degree 130 and for the latter up to degree 80. For an ALVZ this is up to degree 60 (GOCE) and 50 (GGM01S), though the results are equal to GOCE for an ALVZ using the more recent GGM02S model. For a CLVZ, GOCE is sensitive to differences in the iceload history from degree 100–130 (GGM01S: 50–70), for an ALVZ from degree 50–70 (GGM01S: 40–50). This means that GOCE, and to a lesser degree GRACE, should be able to discriminate between different ice-load histories if we know the properties of the LVZ. We have also shown that both missions are sensitive to the inclusion of a more realistic oceanload history by implementing MWIF. As we are interested in constraining properties of shallow LVZs, the sensitivities become measures for errors. We have shown that it is possible to extract information on the CLVZ in the presence of these errors by using the degree correlation coefficient and

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

normalized degree amplitudes, or spectral signatures. These signatures show a distinct behaviour for different properties of the CLVZ and are largely independent of the ice-load history. They can be regarded as an estimate of the system function of our system for a certain transient time scale. We have shown that the modelled signal due to LVZs is above the performance of GOCE and GRACE. In reality however, we first have to extract the signal due to an LVZ from the measurements. This is not trivial, as most processes that give rise to mass anomalies are not very well modelled. In future work we will use simulated GOCE measurements to estimate how well the gravity signal due to a shallow LVZ can be extracted from these measurements in the presence of unmodelled other high resolution gravity signals, noise, and errors due to uncertainties in the ice-load history (see e.g. [22] for the GRACE mission). A second aspect will be the introduction of lateral heterogeneities in our earth model using the finite element method (see e.g. [45]). We then expect more realistic estimates of the effect of an LVZ on the global gravity field and on present-day sea level change.

Acknowledgements We thank Patrick Wu, Erik Ivins, Roberto Sabadini and Scott King for their constructive comments, Kurt Lambeck and Tony Purcell (ANU, Canberra) for their Fennoscandian ice sheet model, Richard Bintanja (IMAU, Utrecht) for his ice sheet-dynamical model, Pieter Visser (DEOS, Delft) for computing the expected GOCE performance, Mark-Willem Jansen (DEOS, Delft) for improvements to our normalmode code, and Gabriele Marquart and Radboud Koop (SRON, Utrecht) for discussions.

References [1] P.N.A.M. Visser, R. Rummel, G. Balmino, H. Su¨nkel, J. Johannessen, M. Aguirre, P.L. Woodworth, C. Le Provost, C.C. Tscherning, R. Sabadini, The European earth explorer mission GOCE: impact for the geosciences, in: J.X. Mitrovica, L.L.A. Vermeersen (Eds.), Ice Sheets, Sea Level and the Dynamic Earth, AGU Geodynamics Series 29, AGU, Washington DC, 2002, pp. 95 – 107.

843

[2] A.B. Watts, E.B. Burov, Lithospheric strength and its relationship to the elastic and seismogenic layer thickness, Earth Planet. Sci. Lett. 213 (2003) 113 – 131. [3] G. Ranalli, D. Murphy, Rheological stratification of the lithosphere, Tectonophysics 132 (1987) 281 – 295. [4] R. Meissner, N.J. Kusznir, Crustal viscosity and the reflectivity of the lower crust, Ann. Geophys. 5B (1987) 365 – 374. [5] A. Klein, W. Jacoby, P. Smilde, Mining-induced crustal deformation in northwest Germany: modelling the rheological structure of the lithosphere, Earth Planet. Sci. Lett. 147 (1997) 107 – 123. [6] E.H. Hearn, What can GPS data tell us about the dynamics of post-seismic deformation? Geophys. J. Int. 155 (2003) 753 – 777. [7] F.F. Pollitz, Transient rheology of the uppermost mantle beneath the Mojave Desert, California, Earth Planet. Sci. Lett. 215 (2003) 89 – 104. [8] S. Stein, M. Wyssesion, Introduction to Seismology, Earthquakes, and Earth Structure, Blackwell Publishing, Oxford, 2003, 540 pp. [9] V. Klemann, D. Wolf, Implications of a ductile crustal layer for the deformation caused by the Fennoscandian ice sheet, Geophys. J. Int. 139 (1999) 216 – 226. [10] G. Di Donato, J.X. Mitrovica, R. Sabadini, L.L.A. Vermeersen, The influence of a ductile crustal zone on glacial isostatic adjustment; geodetic observables along the U.S. East Coast, Geophys. Res. Lett. 27 (2000) 3017 – 3020. [11] R. Kendall, J.X. Mitrovica, R. Sabadini, Lithospheric thickness inferred from Australian post-glacial sea-level change: the influence of a ductile crustal zone, Geophys. Res. Lett. 30 (2003) 1461 – 1464. [12] L.L.A. Vermeersen, The potential of GOCE in constraining the structure of the crust and lithosphere from post-glacial rebound, Space Sci. Rev. 108 (2003) 105 – 113. [13] W. van der Wal, H.H.A. Schotman, L.L.A. Vermeersen, Geoid heights due to a crustal low viscosity zone in glacial isostatic adjustment modeling; a sensitivity analysis for GOCE, Geophys. Res. Lett. 31 (2004), doi:10.1029/2003GL019139. [14] A.M. Tushingham, W.R. Peltier, ICE3G: a new global model of late Pleistocene deglaciation based upon geophysical predications of postglacial relative sea level change, J. Geophys. Res. 96 (1991) 4497 – 4523. [15] K. Lambeck, C. Smither, P. Johnston, Sea-level change, glacial rebound and mantle viscosity of northern Europe, Geophys. J. Int. 134 (1998) 102 – 144. [16] R. Bintanja, R.S.W. van der Wal, J. Oerlemans, Global ice volume variations through the last glacial cycle simulated by a 3-D ice-dynamical model, Quatr. Int. 95–96 (2002) 11 – 23. [17] H.N. Pollack, S.J. Hurter, J.R. Johnson, Heat flow from the earth’s interior: analysis of the global data set, Rev. Geophys. 31 (1993) 267 – 280. [18] V. Pasquale, M. Verdoya, P. Chiozzi, Heat flux and seismicity in the Fennoscandian Shield, Earth Planet. Sci. Lett. 126 (2001) 147 – 162. [19] B.C. Douglas, Global sea level rise, J. Geophys. Res. 96 (1991) 6981 – 6992.

844

H.H.A. Schotman, L.L.A. Vermeersen / Earth and Planetary Science Letters 236 (2005) 828–844

[20] J.X. Mitrovica, J.L. Davis, Present-day post-glacial sea level change far from the Late Pleistocene ice sheets: implications for recent analyses of tide gauge records, Geophys. Res. Lett. 22 (1995) 2529 – 2532. [21] B.D. Tapley, S. Bettadpur, M. Watkins, C. Reigber, The gravity recovery and climate experiment: mission overview and early results, Geophys. Res. Lett. 31 (2004), doi:10.1029/ 2004GL019920. [22] I. Velicogna, J. Wahr, Postglacial rebound and Earth’s viscosity structure from GRACE, J. Geophys. Res. 107 (2002), doi:10.1029/2001JB001735. [23] M. Simons, B.H. Hager, Localization of the gravity field and the signature of glacial rebound, Nature 390 (1997) 500 – 504. [24] L.L.A. Vermeersen, R. Sabadini, A new class of stratified visco-elastic models by analytical techniques, Geophys. J. Int. 139 (1997) 530 – 571. [25] W.R. Peltier, The impulse response of a Maxwell earth, Rev. Geophys. Space Phys. 12 (1974) 649 – 669. [26] P. Wu, W.R. Peltier, Viscous gravitational relaxation, Geophys. J. R. Astron. Soc. 70 (1982) 435 – 485. [27] J.X. Mitrovica, W.R. Peltier, On postglacial sea level, J. Geophys. Res. 96 (1991) 20 053 – 20 071. [28] W.E. Farrell, J.A. Clark, On postglacial geoid subsidence over the equatorial oceans, Geophys. J. R. Astron. Soc. 46 (1976) 647 – 667. [29] G.A. Milne, J.X. Mitrovica, J.L. Davis, Near-field hydro-isostasy: the implementation of a revised sea-level equation, Geophys. J. Int. 139 (1999) 464 – 482. [30] J.X. Mitrovica, G.A. Milne, On post-glacial sea level: I. General theory, Geophys. J. Int. 154 (2003) 253 – 267. [31] W.D. Mooney, G. Laske, T.G. Masters, CRUST 5.1: A global crustal mode at 5050, J. Geophys. Res. 103 (1998) 727 – 747. [32] U. Byrkjeland, H. Bungum, O. Eldholm, Seismotectonics of the Norwegian continental margin, J. Geophys. Res. 105 (2000) 6221 – 6236. [33] D.A. Wiens, S. Stein, Age dependence of oceanic intraplate seismicity and implications for lithospheric evolution, J. Geophys. Res. 88 (1983) 6455 – 6468. [34] J.E. Dixon, T.H. Dixon, D.R. Bell, R. Malservisi, Lateral variation in the upper mantle viscosity: role of water, Earth Planet. Sci. Lett. 222 (2004) 451 – 467.

[35] G. Di Donato, L.L.A. Vermeersen, R. Sabadini, Sea-level changes, geoid and gravity anomalies due to Pleistocene deglaciation by means of multilayered, analytical earth models, Tectonophysics 320 (2000) 409 – 418. [36] A.M. Dziewonski, D.L. Anderson, Preliminary reference earth model, Phys. Earth Planet. Inter. 25 (1981) 297 – 356. [37] M.J. Siegert, J.A. Dowdeswell, Numerical reconstructions of the Eurasian Ice Sheet and climate during the Late Weichselian, Quat. Sci. Rev. 23 (2004) 1273 – 1283. [38] L.M. Cathles, The Viscosity of the Earth’s Mantle, Princeton Univ. Press, Princeton, 1975, 390 pp. [39] J.X. Mitrovica, W.R. Peltier, Pleistocene deglaciation and the global gravity field, J. Geophys. Res. 94 (1989) 13 651 – 13 671. [40] F.G. Lemoine, D.E. Smith, L. Kunz, R. Smith, E.C. Pavlis, N.K. Pavlis, S.M. Klosko, D.S. Chinn, M.H. Torrence, R.G. Williamson, C.M. Cox, K.E. Rachlin, Y.M. Wang, S.C. Kenyon, R. Salman, R. Trimmer, R.H. Rapp, R.S. Nerem, The development of the NASA GSFC and NIMA joint geopotential model, in: J. Segawa, H. Fujimoto, S. Okubo (Eds.), Gravity, Geoid and Marine Geodesy, International Association of Geodesy Symposia, vol. 117, 1997, pp. 461 – 469. [41] A.V. Oppenheimer, R.W. Schafer, Digital Signal Processing, Prentice-Hall, Englewood Cliffs, 1975, 784 pp. [42] B.C. Douglas, Sea level change in the era of the recording tide gauge, in: B.C. Douglas, M.S. Kearney, S.P. Leatherman (Eds.), Sea Level Rise: History and Consequences, International Geophysics Series 75, Academic Press, San Diego, 2001, pp. 37 – 64. [43] W.R. Peltier, Global glacial isostatic adjustment and modern instrumental records of relative sea level history, in: B.C. Douglas, M.S. Kearney, S.P. Leatherman (Eds.), Sea Level Rise: History and Consequences, International Geophysics Series 75, Academic Press, San Diego, 2001, pp. 65 – 95. [44] E.R. Ivins, C.G. Sammis, Transient creep of a composite lower crust: 1. Constitutive theory, J. Geophys. Res. 101 (1996) 27 981 – 28 004. [45] P. Wu, H. Wang, H. Schotman, Postglacial induced surface motions, sea-levels and geoid rates on a spherical, self-gravitating, laterally heterogeneous Earth, J. Geodyn. 39 (2005) 127 – 142.