A dynamic model for the Iceland Plume and the North Atlantic based on tomography and gravity data

(1)

GJI Geo desy , p otential ﬁeld and applied geophysics

A dynamic model for the Iceland Plume and the North Atlantic based on tomography and gravity data

Gabriele Marquart

¹

and Harro Schmeling

²

1SRON and Institute of Earth Science, University of Utrecht, the Netherlands. E-mail: marquart@geo.uu.nl

2Institute of Meteorology and Geophysics, J. W. Goethe University Frankfurt, Germany

Accepted 2004 June 18. Received 2004 April 20; in original form 2003 May 10

S U M M A R Y

The North Atlantic around Iceland is characterized by a large geoid high and an anomalous shallow ocean floor; both observations are presumably due to upwelling of hot material in the mantle. We extracted this region from a global tomography data set to study the structure of the mantle in more detail. The tomography data reveal a confined low-velocity structure from the core–mantle boundary (CMB) to the upper mantle, stretching from a location at the CMB below southwest Greenland towards the upper mantle in a strongly eastward inclination. In the present study we compare the observed gravity potential field in the North Atlantic with a modelled field based on mantle temperature variations estimated from tomography. Seismic traveltime residuals are converted to temperature variations assuming a linear relation between seismic velocity and density and a pressure-dependent thermal expansivity. We found a maximum excess temperature in the plume conduit at the CMB of∼250^◦C, weakly decreasing towards the phase transition zone (PTZ). In the PTZ a temperature rise of∼50–70^◦C is found, which is in agreement with the latent heat release by the olivine phase transitions. The 3-D temperature field is then used as the driving force for viscous flow in a Cartesian convection model. The model dimensions are chosen four times as large as the tomographic section to allow resolution for long wavelengths and convective return flow. For a number of constant viscosity and temperature- and depth-dependent viscosity cases the dynamic topography, gravity anomaly and geoid undulations are calculated and compared with the EGM96 potential field coefficients in the wavelength range of 400 to 4000 km. The observation data were also corrected for ocean lithosphere cooling and isostatic compensation of continental crust. The best agreement between observation and modelled data (78 per cent fit for geoid and 47 per cent for gravity) is obtained for a temperature-dependent viscosity of about one order of magnitude for 500

◦C temperature variation and an increase of viscosity with depth by no more than a factor of 50 from the upper to the lower mantle. The generally good spatial agreement supports the tomographic model, at least for the upper mantle, and indicates that the East Greenland margin as well as the outer Faroer Ridge are dynamically supported. Low-density material west of the Kolbeinsey Ridge might be linked to low-density anomalies below the Greenland Shield. The presence of lower mantle anomalies causes a large-scale geoid high of∼3 to 5 m in agreement with observations, but our approach cannot further constrain the spatial distribution of anomalies in the lower mantle.

Key words:geoid anomalies, gravity anomalies, mantle viscosity, Iceland, North Atlantic, tomography.

1 I N T R O D U C T I O N

The opening of the North Atlantic started at about 60 Ma with the extrusion of huge amounts of flood basalts. Anomalous volcanic activity has continued since then. The entire region is characterized by a distinct geoid anomaly, with an equipotential surface elevation of∼40 m above the hydrostatic reference figure of the Earth, and

an anomalously shallow ocean depth compared with other ocean floors of the same age. No other ridge-centred or intraplate volcanic province has similar large-scale characteristics. To explain these general observations, it has been assumed that the mantle beneath the North Atlantic and the adjacent continental regions is anomalously warm, giving rise to large-scale dynamic topography and positive gravity potential field anomalies. To test this assumption, we studied

(2)

tomography and gravity data for the North Atlantic in a geodynamic context.

Tomography data show the variations of seismic wave velocity in the Earth’s interior. Under the assumption of a chemically homogeneous mantle, tomography data reflect the temperature variations in the Earth’s mantle and are directly related to density variations.

The seismic velocity structure in the North Atlantic has been investigated by different groups. Regional models, mainly based on the data set of the ICEMELT campaign (Wolfeet al.1997; Allenet al.

1999) show a low-velocity anomaly beneath Iceland extending from at least 400 km depth up to about 150 km with a rather small maximum decrease inP-wave velocity of about 2 per cent. The upper 150 km seem to have an ambiguous structure with regions of high velocity over others with lower velocities (Allenet al.1999). However, surface wave dispersion data along north–south profiles across Ice- land indicate very lowS-wave velocities at around 80 km depth for central Iceland (Bjarnason, personal communication; Allen et al.2002). The large-scale tomographic model used in this study (Bijwaard & Spakman 1999) indicates aP-wave velocity reduction of around 7 per cent in the upper 100 km below Iceland. The ap- parent discrepancy is at least partly related to the seismic reference mantle; for the large-scale model it was a global average mantle model, while for the regional studies an already hotter reference state was assumed, as it might be meaningful for the North Atlantic.

Density variations at depth produce gravity and geoid anomalies at the surface depending on the magnitude and the depth of the anomalies. Negative density anomalies, indicating a hotter mantle, cause negative gravity potential field anomalies. In addition, the driving buoyancy of the density anomalies forces solid state flow, which leads to deflection of the surface (and internal interfaces) and produces an additional gravity effect. The observed gravity potential field anomaly is the result of the superposition of both effects.

Even if the internal density variations are known, the magnitude and geometry of the resulting flow, and thus the deflections, are not, since the flow greatly depends on the viscosity structure of the Earth’s mantle. Here we start from a tomographic data set and infer internal densities and mantle temperatures. Based on these data, we determine by geodynamic modelling the gravity, geoid and dynamic topography in the North Atlantic. By comparing these ‘model predictions’ with the real surface observables we can first determine the laterally averaged vertical viscosity structure of the mantle. By allowing additional temperature-dependent viscosity variations and estimating the overall spatial fit, we can find the ‘best-fitting’ model, allowing further implications for mantle flow in the North Atlantic region.

While most dynamic models for the Iceland Plume (as for plumes in general) are fully dynamic in the sense that the Navier–Stokes equation and the heat transport equation are solved simultaneously and an idealized plume rise is studied in its temporal evolution (e.g.

Ribeet al.1995; Itoet al.1999; Ruedaset al.2004), a different approach is used here. We start with a given buoyancy field determined from a tomography model and calculate the dynamic response in a limited area. Those approaches have been widely used for global studies (e.g. Richards & Hager 1984) but hardly for regional studies and temperature-dependent viscosity. The disadvantage of the approach is that we cannot determine any temporal evolution; the advantage is that the model is closer to observations.

We first describe the physics and layout for the regional geodynamic model of the North Atlantic, then how we prepared the tomography and gravity data to be used for this model, followed by an estimate of the mantle viscosity structure for this area, to give a general fit between modelled and observed gravity potential

field anomalies. Finally we discuss the spatial agreement and some implications.

2 P H Y S I C S A N D L AY O U T O F T H E G E O D Y N A M I C M O D E L

The geodynamic model is based on the assumption of a viscous Earth where solid state creep is driven by thermal buoyancy. The physi- cal formulation includes the Navier–Stokes and the mass conserva- tion equations for infinite Prandtl number and extended Boussinesq approximation. The equations are non-dimensionalized and solved numerically in 3-D Cartesian coordinates using a hybrid Fourier/FD approach (the method is described in more detail in Ruedaset al.

2004). The formulation includes the buoyant effects of the phase boundaries at 410 km and 660 km depth. Denoting non-dimensional variables by primes, the governing equations are

0= −∇p+∂(2µe˙_{i j})

∂x_j −

RaT−RacρO/SfO/S−RacρS/PfS/P

ez (1)

0= ∇ ·u (2)

The subscript O/S and S/P denote the olivine–spinel and spinel–

perovskite phase change andf is the percentage of transformed material. The scaling numbers and all variables used in this paper, as well as their numerical values, are defined in Table 1. The thermal expansivity was assumed as pressure (depth) dependent in a way described below. Thus, alsoRaandRacare functions of depth. The depth- and temperature-dependent viscosity is modelled by

µ(T,z)=µ0w(z)[exp(−c1T)] (3)

wherew(z) can be a function of the type exp(c2 z) describing a 10th to a 100-fold increase of viscosity with depth or a function describing a sudden jump in viscosity in the mantle transition zone.

c1 is chosen in the way that one order of magnitude in viscosity change can be related to 100 to 800 K temperature change.

For the model we assumed dimensions of 1.6×10⁴km in both horizontal directions and 2880 km depth and used 256×256× 120 gridpoints. This implies a resolution of∼63 km horizontally and 25 km in depth. The phase transition of olivine to spinel at 410 km was assumed to occur linearly over a depth interval of 100 km (equivalent to four gridpoints) and the spinel–perovskite transition at 660 km over 50 km (i.e. two gridpoints). The Clapeyron slopes for olivine were chosen from two papers by Akaogiet al.

(1989, 1998) and treated as discussed explicitly in Marquartet al.

(2000). The anomalous temperature distribution was derived from tomography and will be described further below, and is imposed in a central area of 4×10³km by 4×10³km. For the undisturbed temperature profile we assumed a conducting half-space model after 30 Myr of cooling, and an adiabatic temperature profile throughout the mantle. 30 Myr of cooling for the lithosphere was assumed as a rough lateral average for this part of the North Atlantic. We used this approximation for the lithosphere instead of the real ridge geometry and age distribution to avoid ridge-related anomalies in the upper 80 km of our model which are not resolved by the tomography data. To be consistent we also removed lithosphere age-related signals from the geoid and gravity observations (explained in more detail in Section 4). At the lower boundary we assumed a constant heat flux of 20 mW m⁻². This is motivated by the fact that the thermal conditions at the CMB are not well constrained on a global

(3)

Table 1. Definition of variables used in this study.

Symbols Parameter Values and units

α Thermal expansivity K⁻¹

γ Gr¨ueneisen parameter 1.2

κ Thermal diffusivity 10⁻⁶m²s⁻¹

µ Viscosity Pa s

µ0 Scaling viscosity 10²¹Pa s

ρ⁰ Mean mantle density 4250 kg m⁻³

ρc Mean crustal density 2780 kg m⁻³

ρw Water density 1000 kg m⁻³

ρO/S Density jump at the olivine–spinel transition 196 kg m⁻³ ρS/P Density jump at the spinel–perovskite transition 253 kg m⁻³

ρ1 Density jump at the CMB 5000 kg m⁻³

cp Specific heat at constant pressure 1.3×10⁵J kg⁻¹K

cv Specific heat at constant volume 1.3×10⁵J kg⁻¹K

˙

eij Strain rate tensor s⁻¹

fO/S Material transformed from olivine to spinel per cent

fS/P Material transformed from spinel to perovskite per cent

g Gravity acceleration 9.87 m s⁻²

G Gravitational constant 6.6732×10⁻¹¹N m²kg⁻²

H Height of model 2880 km

h Topography m

Hh Mean isostatic compensation depth 304 or 40 km

KT Bulk modulus at constant temperature above 410 km 135×10⁹Pa Bulk modulus at constant temperature below 410 km 200×10⁹Pa

p Dynamic pressure Pa

P Hydrostatic pressure Pa

Ra Rayleigh number (ρgα TH³)/(κ µ0)

Rac Chemical Rayleigh number (ρgH³)/(κ µ0)

T Temperature ^◦C or K

T Scaling temperature difference 1200^◦C

u Flow velocity m yr⁻¹

VP SeismicP-wave velocity anomaly km s⁻¹

adiabatic gradient

no slip, prescribed heat flux 2880 km

1.6 e3 km 4000 km

4000 km

1.6 e3 km free slip, prescribed temperature

gradient plus

"tomographic"

adiabatic

temperature

periodic

Figure 1. General layout of the numerical model.

average and the tomography model does not have sufficient resolution at the CMB. The background temperature at the lithosphere base was thus 1200^◦C and 1550^◦C just above the thermal boundary layer at the CMB at a depth of 2880 km. The general model layout is shown in Fig. 1. This set-up with a large numerical box and a smaller internal area with prescribed temperature variations was used to allow long-wavelength flow to occur in the model box.

Since we used a Cartesian model but intended to compare model quantities such as geoid and gravity with observations, one has to consider this long-wavelength flow. If a numerical box as small as the central data area were to be used, return flow would be forced to occur within the box and the flow field would strongly be biased.

Alternatively one may assume an open side wall with in and out flux, but this would only partly solve the problem, since the lateral flow field has to be prescribed in some artificial way. Therefore we used a large box, but, since the long-wavelength components of geoid

and gravity anomalies are not comparable to observations on the (spherical) Earth, the resulting numerical data sets were filtered by a 2-D Fourier method for wavelengths of less than 4000 km (an ade- quate procedure was applied to the observation data, as explained in Section 4). We checked that the amplitude of long-wavelength flow was small in our model (nearly two orders of magnitude smaller than the velocity in the central plume) and thus the flow dynamics was virtually unaffected by the choice of the model box. However, the effect, especially on the geoid, is large and therefore this filter procedure is necessary.

3 T H E T O M O G R A P H Y M O D E L A N D R E L AT E D M A N T L E T E M P E R AT U R E S For this study we use the high-resolution global tomography data set of Bijwaard & Spakman (1999). This data set consists of

(4)

P-wave traveltime residuals compared with an average preliminary reference earth model (PREM) mantle. From this data set a subset for the North Atlantic is extracted and represented in blocks of 0.6^◦ length in the region between 49.8^◦N, 50^◦W to 85.2^◦N, 15.4^◦E and in 26 non-equidistant depth layers. Bijwaard & Spakman (1999) have already reported on a clear indication of a plume conduit related to the Iceland hotspot in this segment of the Earth.

Seismic traveltime residuals are often stronger in the upper mantle than in the lower mantle. For the data set used here the maximum residual for the upper 400 km is 7.0 per cent, for the depth interval 400 to 1000 km it is 1.5 per cent and for the lower mantle to the CMB it is 0.74 per cent. Similar variations have been observed in other data sets as well and are explained by the increase in the elastic moduli with depth.

3.1 Preparing the tomography data for the geodynamic model

To prepare the input for the Cartesian convection code the tomography data set was regridded on 64×64×120 equidistant values, using a Mercator projection. First, to investigate the seismic structures of the mantle beneath the North Atlantic in more detail we scaled the residuals by their maximum values in the three depth intervals introduced above to allow for comparison of structures in the upper and lower mantle. Fig. 2 shows a 3-D isosurface representation of the data set under four different viewing angles to give a 3-D impression of the structures. The figure shows the contour plane where 75 per cent of the maximum negative velocity anomaly scaled in the way explained above was reached. The contour plane

Figure 2. Regridded tomographic model after Bijwaard & Spakman (1999). Notice the clear presence of an uprising low-seismic-velocity anomaly, originating at the CMB at a position beneath the southern tip of Greenland. This anomaly rises, inclined strongly northeastwards, through the lower mantle and interacts with the mantle transition zone beneath the western European margin at the latitude of Great Britain. Through the transition zone itself, no clear continuation could be identified, but a strong upper mantle anomaly is present below Iceland.

representation also reveals the original block structure of the tomography inversion grid. A plume structure is clearly visible, stretching from the CMB to the northeast through the entire lower mantle.

The structure seems to split up when reaching the upper mantle with one branch continuing below western Europe (outside the region under investigation) and a second branch stretching towards the North Atlantic Ridge and Iceland. A second strong anomaly is located below the Greenland Shield and, by following the base of the continental lithosphere of Greenland below the northwestern part of the North Atlantic, reaches the Kolbeinsey Ridge north of Iceland.

This strong tilting and branching of low-seismic-velocity material in the lower mantle has not been expected from numerical modelling of mantle plumes and has initiated a controversial discussion about whether this structure should be identified as a lower mantle plume (Bijwaard & Spakman 1999) or not (Foulgeret al.2001). As a kind of compromise, the recent finite-frequency tomography model of Montelliet al.(2004) reveals a−0.3 per centVP anomaly in the lower mantle similar to that of Bijwaard & Spakman (1999) which can be traced down to a depth of 1900 km, but in contrast to that study this anomaly is lost at depths greater than 2350 km.

3.2 Mantle temperatures from the tomography model For the dynamic model we need to relate theP-wave variations to temperature, as temperature controls the buoyancy forces and the viscosity variations in our model. Furthermore, we are also interested in the excess temperature, predicted from the tomographic model for the North Atlantic mantle. Under the assumption of a chemically homogeneous mantle, the conversion from seismic

(5)

velocity anomalies to density can simply be done by a constant factor (Karato 1993), here we used a factor of 0.3 fordln ρ/dln VP. For the conversion to temperature we used the linearized relation between temperature and density, with the thermal expansivityα given by

α= γc_vρ KT

(4) withγ as the Gr¨uneisen parameter,cvas specific heat at constant volume andKT as the bulk modulus at constant temperature. All parameters on the right-hand side depend on the pressureP, but by far the strongest dependence is found forKT. Therefore onlyKT(P) is considered in the form

KT(P)=KT0+ d K

d pP(z). (5)

All other parameters in eq. (4) vary only slightly with depth (e.g.

Walzeret al.2003) and were treated as constant. For a depth above 410 kmKT0was assumed as 135×10⁹Pa anddK/dpas 6, and below KT0was 200×10⁹Pa anddK/dpequal to 3.5. These values are in accordance with a PREM mantle (Stacey 1992). This parametriza- tion results in a dependence ofαon pressure (depth) in the way that αhas a value of 4.3×10⁻⁵K⁻¹at the top, decreases to about half this value in the mantle transition zone and reaches a value of 1.1× 10⁻⁵K⁻¹close to the CMB. Thus, this treatment ofαenhances the effect on temperature in the deep mantle and damps it in the upper mantle.

The relation between seismic velocity and temperature was thus obtained by

∂lnT

∂lnV_P = ∂lnT

∂lnρ

∂lnV_P = −0.3 1

αT. (6)

The resulting maximum excess temperature with depth (temperature above the adiabatic gradient) is shown in Fig. 3. The excess temperature in the lower part is about 200 to 250^◦C; about 100

◦C is lost on rising through the lower mantle. This value is in good agreement to those obtained in other studies. While the excess temperatures in the lower mantle are lower than estimated by numerical plume models (e.g. 400 K or more by Farnetani 1997), the relative

Figure 3. Maximum excess temperature profile (reduced by the adiabatic gradient) in the North Atlantic based on tomography data and eq. (6).

loss of about 40 per cent between the bottom and top of the lower mantle matches that of a plume with moderate buoyancy flux (Al- bers & Christensen 1996). In the transition zone the temperature rises again by 50 to 80 K. Note that this increase is relatively small compared with the strong increase of the seismic anomaly due to the depth dependence ofα. We attribute the moderate temperature increase to the effect of latent heat of the perovskite–spinel–olivine transitions. It should be noted that this interpretation implies that in the surrounding mantle vertical velocities are considerably smaller, resulting in less production and consumption of latent heat there.

In the uppermost mantle, above 300 km, the temperature appears to increase again; however, here we expect that effects such as mantle depletion and melting occur and bias the relation. From geochem- ical studies on Icelandic lavas (Schilling 1991; Nicholson & Latin 1992) a plume excess temperature around 250^◦C was estimated; a similar or slightly lower excess temperature had been estimated from seismic velocity variations in the uppermost mantle below Iceland (Allenet al.1999; Wolfeet al.1997; Foulgeret al.2001). Fully dynamic models of plume rise including melting in the plume centre (Ruedaset al.2004; Itoet al.1999; Keen & Boutilier 2000) gave somewhat lower excess temperatures of 100 to 180^◦C.

Since only the maximum temperature is shown in Fig. 3, we are not sure if we are always inside the ‘plume conduit’. It might therefore be of interest to discuss the 3-D temperature field for the area of investigation with respect to Fig. 2. We found that temperatures above 200^◦C are present in the ‘deep mantle plume conduit’, in the deeper upper mantle below Greenland, and along the entire ridge axis, but more strongly in the northern branch. In the very uppermost mantle, above 200 km, the maximum excess temperature focuses again at Iceland. With this temperature field, obtained in the way described above, we solved the equations for convective flow in the entire mantle below the North Atlantic region.

4 S U R FA C E O B S E RVA B L E S : G R AV I T Y, G E O I D A N D T O P O G R A P H Y

The Earth gravity field is a sensor for density distribution at depth.

To further test and constrain the tomographic model and the derived temperature field, we model the geoid and gravity field and compare them with observational data. We used the spherical harmonic expansion EGM96 (Lemoineet al.1998) of the potential field.

Since the maximum wavelength which we can resolve in the numerical model is constrained by the size of the numerical box (see Fig. 1) we need to filter the surface observational data for the same wavelengths. Therefore we used degree 6 as the cut-off degree for long wavelengths: this means that at a latitude of∼65^◦N all wavelength longer than∼4000 km are filtered out. In relation to geodynamic processes this is equivalent to ignoring large-scale, global- wide mantle circulation. Furthermore, since the geodynamic model allows only viscous rheology, we cannot model elastic support of any load, neither is the contribution of continental crust modelled.

Both effects are mainly reflected in the short-wavelength anomalies. To eliminate these effects, we used a short-wavelength cut-off degree of 55, which is about 400 km at∼65^◦N. Thus, the gravity potential and the topography data which we used in the further investigation contain wavelengths between∼400 and 4000 km. To avoid artificial effects of tapering we used a cosine slope at both edges of the filter band. We tested the sensitivity of the geoid, gravity and topography data on the tapering by using different slopes and edge points for the filter. For the long-wavelength range we shifted the

(6)

EGM96 Gravity L=6-55

310˚ 320˚ 330˚ 340˚ 350˚ 0˚ 10˚

60˚

70˚

80˚

-40 -30 -20 -10 0 10 20 30 40

mgal

EGM96 Geoid L=6-55

310˚ 320˚ 330˚ 340˚ 350˚ 0˚ 10˚

60˚

70˚

80˚

-30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30

m

Residual Gravity

310˚ 320˚ 330˚ 340˚ 350˚ 0˚ 10˚

60˚

70˚

80˚

-40 -30 -20 -10 0 10 20 30 40

mgal Residual Geoid

310˚ 320˚ 330˚ 340˚ 350˚ 0˚ 10˚

60˚

70˚

80˚

-30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30

m

(a) (b)

Figure 4.(a) Filtered fields of surface observables. Upper figure: geoid height anomaly; lower figure, gravity anomaly. All data sets are filtered for spherical harmonics between degree 6 and 55, the EGM96 model was used for the potential field coefficients and a spherical harmonic representation of ETOPO5 for topography (National Geophysical Data Center 1988). (b) Fields of residual of surface observables after subtracting the effect of the cooling oceanic lithosphere and the isostatic effect of continental topography. Upper figure: geoid height anomaly; lower figure, gravity anomaly.

edge point from spherical harmonic degree 5 to 6 and varied the slope of the cosine over 0 to 3 spherical harmonic degrees. At the short-wavelength side we shifted the edge point between degree 40 and 60 and studied slopes over 2 to 10 degrees. The effect of the different filter bands resulted in variations for the extrema of the gravity field of about 6 mGal, and of about 5 m for the geoid undu- lation. Although, these values are about 20 per cent of the maximum anomalies, and therefore not small, they are, however, variations of peak value amplitudes; the regional pattern of the fields does not change considerably with the filter band. The fields which we used for further investigation are shown in Fig. 4.

The geoid in this region is simply characterized by a broad high. The long-wavelength gravity field, however, shows much more structure. Beside the ocean–continent distribution, the Faroer–

Iceland–Greenland Ridge, the Reykjanes Ridge to the south and the elevated topography of the Rockall Bank are clearly exposed, and additionally strong positive anomalies over eastern and southern Greenland and along a band north of Iceland, trending from southwest to northeast are visible. This positive anomaly is clearly related to the Kolbeinsey spreading axis which, however, has only a weak expression in the long-wavelength topography.

Since our numerical modelling does not include a spreading axis, we reduce the effect of the cooling oceanic lithosphere from the data.

To do this, we first used the age data for the oceanic crust (Mueller et al.1996) and filtered the age data for the same wavelength range as the gravity potential field, setting the continental areas to an artificial age of 60 Ma. Then we calculated the effect on geoid and bathymetry by using the formulae for a cooling lithosphere as given by Haxby & Turcotte (1978). Finally we reduced the effect on the gravity anomalies by applying the spectral relation between geoid and gravity. This was done by a 2-D Fourier transformation for the gridded geoid data of our area under investigation. IfNâis the geoid effect due to age andNâ_lmare the spectral modes, then the gravity anomaly due to age was determined bygâ_lm=g0(k²_l +k²_m)^1/2Nâ_lm, withklandkm as wavenumbers. This formula is only an approximation, but the error mainly effects the longest wavelengths which get strongly damped in gravity anomalies. The same procedure was used by Heller & Marquart (2002).

The influence of continental topography and isostatic crustal thickness variations is not of major concern for this study since we are mainly interested in the oceanic area. However, the Green- land Shield is a large part of the area of investigation and here we also find strong contributions to geoid and gravity anomalies which are partly due to crustal isostasy. To reduce the isostatic effect we simply assumed the topography to be balanced by Moho depth variations in relation to a reference depthHhwhich was set to 40 km for

(7)

the Greenland Shield and 30 km for the European part of the area.

For that case the isostatic geoid anomaly is given by (e.g. Marquart 1989)

NT(h)= πG

g ρcHhh(λ, φ). (7)

The long-wavelength topography as a function of longitude and latitude,h(λ,φ), was derived from a spherical harmonics expansion of ETOPO5 (Heller & Marquart 2002). Iceland itself, with its thick crust (Darbyshireet al. 2000) of up to∼40 km, was treated in the same way as Europe. This is also in agreement with results from admittance studies for Iceland and the North Atlantic which show that a considerable part of the positive gravity potential field anomaly across Iceland is related to shallow compensated masses (Marquart 1991; Heller & Marquart 2002). The effect of isostatic topography on the potential field is strongest in eastern Greenland where topography is high and reaches a maximum value of 10 m for the geoid and 25 mGal for gravity. The residual geoid and gravity anomalies after reducing the effects of a cooling oceanic lithosphere and of continental topography, assuming crustal isostasy for the area under investigation, are given in Fig. 4(b). A comparison of Figs 4(a) and (b) clearly shows that strong positive anomalies are still present, mostly located to the west of the spreading ridge. The anomaly has a few branches, one running along the ridge with an offset to the west and strongest in the area north of Iceland, another one along eastern Greenland, and a smaller one related to the outer Iceland–Faroer Ridge.

5 A D J U S T I N G T H E AV E R A G E V I S C O S I T Y P R O F I L E

For the dynamic model, the viscosity of the mantle material and its dependence on pressure and temperature are essential. The combined effect of increasing hydrostatic pressure and temperature with depth leads to an increase in viscosity. During the last 20 years many studies have been performed trying to relate global-scale geoid anomalies and mantle density structures (e.g. Richards & Hager 1984; Ricardet al.1993; Bunge & Richards 1996; Forte & Mitro- vica 2001). These models have been successively refined by adding surface plates (Ricard & Wuming 1991), compressibility (Forte &

Peltier 1991) or phase boundaries (Cadek & Fleitout 1999). Some studies try to use as few viscosity layers as possible (e.g. Ricardet al.

1993), others involve as many as 20 layers (e.g. Steinbergeret al.

2001). As a measure for the spatial agreement between modelled and observed data the fitfis suitable. Ifaiare the observations atn data points withi=1,. . .,nandbiare the modelled values at the same locations, and ¯ais the arithmetic mean of allai, the varianceδ of the observations is given by

δ(an)= 1 n

n i=1

(ai−a)² (8)

and

f(an,bn)=1− δ(an−bn)

δ(an) . (9)

This describes the fit between the two data sets.

The fit is generally found to be as good as 0.8 for very low degrees of harmonics but this decreases quickly with higher harmonic degrees. From the above studies we know very well that the viscosity increases with depth by about a factor 30 to 100 from the surface to the CMB. Changes of viscosity due to variations in dynamic pressure are negligible, but the temperature dependence might be

Figure 5. Two test cases of viscosity–depth distribution, one with a gradual increase throughout the mantle and one with a sudden increase in the mantle transition zone.

important, as is well documented from many small-scale convection experiments (e.g. Larsen & Yuen 1997). Thus, as a first step, we determined the dependence of viscosity on depth for the North Atlantic region and, in a second step, allowed for the effect of temperature, which is important for our medium-scale analysis. We tried two assumptions: in the first the viscosity increases gradually with depth and in the second the increase is restricted to the mantle transition zone (Fig. 5). In addition to cases with a high viscous upper lid, as described in Section 2, we also modelled cases without a lid where the upper mantle viscosity of 10²¹Pa s stretched up to the surface.

For the temperature field determined from tomography and for different viscosity depth functions we calculated the geoid and gravity at the surface from the dynamic topography and the density variations at depth due to temperature anomalies and to phase transitions outside the average transition regions, taking into account 60 per cent of olivine in the mantle rocks. From the filtered observational data (Fig. 4b) we estimate the maximum geoid high in the central part of the region under investigation to be 25 m, and the gravity anomaly to be 25 mGal; we compared these values with the peak geoid and gravity anomalies in the central part of the models. Fur- thermore we determined the fit between observation and models by regridding the Mercator-projected filtered observation data (Fig. 4) onto a 64×64 data point grid to have the same resolution as the output of the dynamic modelling. Fig. 6 shows the maximum amplitudes (upper figures) and the fits (lower figures) for geoid (left) and gravity (right) depending on the viscosity ratio between the upper and lower mantle for the two different cases.

6 R E S U L T S

Inspecting the curves in Fig. 6 we can conclude the following. In- creasing viscosity with depth reduces the amplitude of geoid and gravity anomalies, the effect is more pronounced for a viscosity

(8)

Figure 6. Peak geoid (left) and gravity (right) anomaly in the central part of the area under investigation depending on the viscosity ratio between the upper and lower mantle. The triangles indicate cases with a viscosity jump in the transition zone and the squares those with a gradual increase. Full lines indicate cases without a high viscous lid and dashed lines those with a lid on top. The thick dashed-dotted line in the upper figures give the maximum amplitudes from the observations.

jump in the transition zone (triangles in Fig. 6) due to stronger decoupling of upper and lower mantle flow. The presence of a high viscous lid (Fig. 6, dashed lines) reduces the amplitude for an other- wise constant-viscosity mantle, since it damps dynamic topography.

If viscosity increases with depth the high viscous lid mainly damps the dynamic topography due to upper mantle flow, which for the mantle density distribution considered here leads to higher amplitudes of geoid and gravity if the increase in viscosity is above a certain value.

The fits and maximum amplitudes are considerably better for a sudden increase in viscosity at the mantle transition zone (Fig. 6, triangles), as preferred in global modelling of geoid anomalies (e.g.

Cadek & Fleitout 1999). The best fit (of 0.66 to 0.78 for geoid and 0.47 and 0.45 for gravity) is found for a viscosity increase of between a factor of 20 and 50 from the upper to the lower mantle for cases without a lid. That fact, that the gravity fit is generally worse than the geoid fit is explained by the higher sensitivity of gravity to short wavelengths which are less well constrained in the tomography model. Obviously, for the wavelength range considered the lithosphere of the North Atlantic does not behave as a homogeneous viscous lid resting on an constant-viscosity upper mantle and accumulating viscous stresses. Possible explanations for this find-

ing may include the presence of a spreading zone or of a shallow low-viscosity asthenosphere as decoupling processes reducing such lithospheric stresses.

For the successful case with a viscosity jump at the transition zone and no lid we further investigated the temperature dependence of viscosity following eq. (3), for a factor of 10 to 30 increase in viscosity. The effect of temperature-dependent viscosity is much stronger in the upper than in the lower mantle, since here the temperature variations have stronger amplitudes and the background viscosity is already lower. By inspecting the flow pattern in the model, we found that in the case where the viscosity drops by one order of magnitude for a temperature rise of 300 K or less, small-scale flow starts in the upper mantle. This small-scale flow affects gravity more strongly than the geoid; as a result (Fig. 7) the geoid fit (left figure) is not much affected. In the case of a stronger increase in viscosity (full line) only a weak dependence on temperature (change of one order in viscosity for a∼500 to 600 K temperature change) leads to a slightly better fit. The gravity fit, however, is strongly reduced in the case of a stronger temperature-dependent viscosity and the presence of small-scale flow. A slight temperature dependence leads to some focusing in the lateral distributions of anomalies which are already present and weakly increases the fit.

(9)

Figure 7. Geoid and gravity fit for models with temperature-dependent viscosity. The temperature scale at the horizontal axis indicates the temperature change necessary for a viscosity change of one order of magnitude. The symbols at the right side of the figures give the fits for a purely depth-dependent viscosity variation. Full lines are for a viscosity jump of a factor 30, dashed lines for a factor 10, respectively.

Figure 8. Comparison between modelled and observed geoid (left) and gravity (right) anomalies.

To allow a visual comparison of the modelled and observed gravity potential field, Fig. 8 shows the observation data (right-hand side) and the modelled fields for the case of a 30-fold increase in viscosity from the upper to the lower mantle with a temperature dependence of one order of magnitude change in viscosity over a 500 K temperature change.

The strongest differences between the observed and modelled gravity fields are found in the continental regions. In the upper 300 km below Greenland the tomography model indicates high seismic velocities in accordance with an old shield region. These high velocities are transferred to high densities with the tendency to drive flow downward. This explains the slightly negative geoid values in the model. In nature, high seismic velocities might, at least partly, be caused by chemical depletion of mantle material. The observed strongly negative geoid values in northwest Europe do not show up in the model and we do not see a simple way to explain them. The

gravity field, showing the small-scale agreement between observation and model, gives more interesting findings. The elongated positive anomaly running from southwest to northeast through the entire region of investigation can clearly be linked to temperature anomalies in the upper mantle unrelated to average lithospheric cooling.

These anomalies are stronger and situated more deeply in the area north of Iceland. The maximum of the gravity anomaly is neither in the observed nor in the modelled gravity field centred on Iceland, but shifted to the west, close to the eastern Greenland margin. This proves a dynamic support for the eastern Greenland margin, which is in agreement with observations in drilling cores from the east Greenland shelf area indicating a delayed subsidence of the seafloor (Cliftet al.1995). Furthermore, the western part of the Faroer Ridge seems to be supported by lighter (hotter) mantle material, which is confirmed by an admittance study of North Atlantic plateaus (Heller

& Marquart 2002) where high admittance values of 2.44 m km⁻¹

(10)

Figure 9. Dynamic topography for the model with a factor of 30 increase in viscosity within the mantle transition zone and a temperature dependence of viscosity of one order of magnitude over 500 K.

(long-wavelength geoid to topography) were found for the western Faroer Block.

Fig. 9 shows the dynamic topography of our best model. A maximum of almost 1 km is predicted close to Iceland, at the same position as the maximum of the gravity and geoid anomalies. This value is larger than the dynamic topography of about 0.6 km predicted by the dynamic flow models of Conradet al.

(2004). In contrast to our models they did not take into account mantle density anomalies at depths of less than 325 km. On the other hand, the residual seafloor topography (topography corrected for age and sediment cover) in a region with a diameter of∼1000 km around Iceland shows amplitudes of between 1.5 and 3 km, reaching locally 4.5 km on Iceland and the Greenland–Faeror–Iceland Ridge (White 1997; Conradet al.2004). Since crustal thickening may only partly be responsible for this seafloor height anomaly, as much as half of this topography may be attributed to dynamic support (White 1997). This is in good agreement with our predictions.

The flow characteristics of this model for the whole mantle are shown in Fig. 10. The vertical flow (lower left) indicates that a weak upwards flow of 2 to 3 mm yr⁻¹ occurs in the entire region. The rising velocity related to the Iceland Plume (dashed line) is about 3.2 cm yr⁻¹in the upper mantle and∼2 cm yr⁻¹in the lower mantle. Note that our vertical profile does not precisely catch the ‘plume conduit’ in the lower mantle, due to its strongly inclined geometry, and the velocity might be a bit higher in the deeper parts of the lower mantle. For the overall model, horizontal velocities are negligible in the lower mantle but have values of 3 to 4 mm yr⁻¹in the upper mantle, indicating a drift of the upper mantle in a roughly westward direction. This is well in accordance with the predicted velocities for the North American–Eurasian plate system from the NUVEL hotspot reference frame (Gripp & Gordon 1990). At the Iceland location (dashed lines) pronounced horizontal velocities are confined to the upper 300 km where the spreading of the plume occurs, which gives rise to horizontal velocities in thex-direction (roughly east–

west) of up to 6 cm yr⁻¹, but considerably smaller velocities in the y-direction (roughly north–south) of 2.5 cm yr⁻¹. This difference in horizontal velocities corresponds well to that which is expected for a plume head below a spreading ridge, where the ridge perpendicular velocities are expected to be greater than the along-ridge velocities (Feighneret al.1995). Note that this velocity field is the result of the tomographic buoyancy field without prescribed kinematic velocities at the upper boundary.

7 C O N C L U S I O N

Based on the tomography model of Bijwaard & Spakman (1999) we determined for the North Atlantic region mantle temperatures related to the Iceland Plume. Relating this temperature field to a buoyancy field we estimated the gravity potential field at the surface. The spatial agreement between the observed and modelled medium-wavelength-range geoid and gravity is quite good, taking into account all the filtering procedures and assumptions entering this study. The presence of a low-seismic-velocity (low-density) anomaly in the lower mantle makes a contribution to the geoid of a few metres in the wavelength range considered. This supports the idea of a high-temperature anomaly in the lower mantle in this region, but our study does not allow us to confirm its spatial distribution in any more detail.

Our investigation suggests a rather weak viscosity increase throughout the mantle, by a factor of 10 to 50 between the upper and lower mantle and a preference for a relatively sudden increase at a depth of 500 to 1000 km. Also the dependence on temperature turned out to be weak in our modelling, however, the numerical grid has a resolution of only 50 km and cannot resolve small-scale variations in flow field and viscosity. The decrease in fit for stronger temperature- dependent viscosity may also be due to enhancing small-scale tomographic errors in the upper mantle. It still remains an open question whether uncertainties in seismic tomography models are enhanced by variable viscosity flow models in a way as to artificially reduce the fit. Higher-resolution seismic tomography and flow modelling would both help to resolve this question.

Furthermore, our study points out that in the North Atlantic region all the mantle material has a slow upward movement, indicating a slow exchange between lower and upper mantle material.

Whether this broad upwelling may be called a lower mantle plume or not may be a question of definition. We also estimated the volume flux in the plume-like structure in the lower mantle and obtained a value of 5 km³ yr⁻¹ for the mid-lower mantle and 3.8 km³ yr⁻¹ for the upwelling region in the upper mantle below Iceland. The value for the Iceland plume is in good agreement with other esti- mates based on chemical observation (a value of 1.43 km³ yr⁻¹; Schilling 1991), seismological observations (a value of 2.5 km³ yr⁻¹; Allenet al.2002) and numerical modelling (Ribeet al.1995;

Feighneret al.1995; Itoet al.1996, 1999) where values between about 2 and 6 km³ yr⁻¹were reported. However, one should keep in mind that due to the relatively high viscosity in the lower mantle, the vertical upward flow is considerably more widespread than the temperature anomaly. For our estimated flux values, we inte- grated the flux over the area with an excess temperature of more than 50^◦C, but the entire vertical flux through the lower mantle is considerably higher. In cases with a confined thermal anomaly but widespread flow the term plume flux might not be appropriate at all.

We also found a relative horizontal movement between the upper and lower mantle. Inspection of the flow pattern showed this to be a westward migration, which is well in accordance with the north- westward movement of the North American–Eurasian plate system in a hotspot reference frame.

It should be emphasized that the problem of fitting geoid and gravity by a dynamic model based on tomography with variable viscosity is a highly non-linear problem. While we focused par- ticularly on the temperature and depth dependence of viscosity, other combinations of model parameters such as different proper- ties of the lid or other parameters controlling the Rayleigh number might further improve the fit to the observed fields. However, such a

(11)

Figure 10. Depth profiles for the model with a factor of 30 increase in viscosity within the mantle transition zone and a temperature dependence of viscosity of one order of magnitude over 500 K. The full line gives the average of the part of the model including the tomography model and the dashed line is a vertical profile at the location of Iceland.

multidimensional search of the parameter space is beyond the scope of this paper.

Furthermore, it cannot be excluded that chemical heterogeneities might be responsible for some of the velocity anomalies. Our interpretation is based on the assumption that theVP-anomalies in the lower mantle are thermal in nature and are related to density anomalies and temperature in the way described in Section 3. To resolve chemical heterogeneitiesVS-anomalies and seismic attenu- ation are needed as independent observations. Some recent studies indicate that chemical heterogeneities are important below continental shields (Deschampset al.2002), in the deep mantle (Ishii &

Tromp 1999) and to a lesser degree in the mantle transition zone (Marquart 2004). Reduced buoyancy in the deepest mantle due to chemical heterogeneities would not have a major influence on the results shown here, since the effects due to the longest-wavelength flow components have been removed. However, it is very likely that the highVPvelocities in the upper 250 km below Greenland are to some extent caused by depletion in iron and not only by increased density due to cooling. Under these circumstances the (modelled) geoid and gravity anomalies for the area of the Greenland Shield, as shown on

Fig. 8, can be expected to be weaker. This would allow a better fit to the observations but would only have a minor effect on the viscosity structure. A stronger impact on the viscosity structure can be expected if chemical separation and water release occur in the mantle transition zone as proposed by Bercovici & Karato (2003); however, the consequences are hard to predict. Our modelling showed that the tomography model of Bijwaard & Spakman (1999) is able to explain important features of the medium-wavelength geoid and gravity field. It might be interesting to test whether this fit may also be obtained with alternative tomography models, e.g. with the finite- frequency tomography model of Montelliet al.(2004), in which the low-velocity anomaly in the lowermost mantle has not been verified.

A C K N O W L E D G M E N T S

Our colleagues in the ‘Frankfurt-Mainz Iceland Plume Working Group’, especially Bernhard Steinberger, made an essential contribution to this work with their interest and discussions. Thomas

(12)

Ruedas gave us access to his comprehensive literature overview of Iceland plume research. Wim Spakman was so kind as to provide us with the North Atlantic part of his tomography model. The com- ments and suggestions of two anonymous reviewers were very help- ful in improving the original manuscript. The study was supported by a research grant of the Deutsche Forschungsgemeinschaft; the Space Research Organization of the Netherlands (SRON) provided the possibility to finish the manuscript. Figs 4, 8 and 9 are produced withGMT(Wessell & Smith 1995).

R E F E R E N C E S

Akaogi, M., Ito, E. & Navrotsky, A., 1989. Olivine-modified spinel–spinel transition in the system Mg2SiO4–Fe2SiO4—calorimetric measurements, thermochemical calculations, and geophysical application,J. geophys.

Res.,94,15 671–15 685.

Akaogi, M., Kojitani, H., Matsuzaka, K., Suzuki, T. & Ito, E., 1998. Post- spinel transformations in the system Mg2SiO4–Fe2SiO4—element par- titioning, calorimetry, and thermodynamics calculations, inProperties of the Earth and Planetary Materials at High Pressure and Temperature, AGU Geophysical Monograph 101, pp. 373–384, eds Manghnani. H.H.

& Yagi, T., Washington, DC, American Geophysical Union.

Albers, M. & Christensen, U.R., 1996. The excess temperature of plumes rising from the core–mantle boundary,Geophys. Res. Lett.,23,3567–

3570.

Allen, R.M. et al, 1999. The thin hot plume beneath Iceland,Geophys. J.

Int.,137,51–63.

Allen, R.M. et al, 2002. Imaging the mantle beneath Iceland using in- tegrated seismological techniques,J. geophys. Res., 107(B12), 2325.

doi:10.1029/2001JB000595.

Bercovici, D. & Karato, S., 2003. Whole-mantle convection and the transition-zone water filter,Nature,425,39–44.

Bijwaard, H. & Spakman, W., 1999. Tomographic evidence for a narrow mantle plume below Iceland,Earth planet. Sci. Lett.,166,121–126.

Bunge, H.-P. & Richards, M.A., 1996. The origin of long-wavelength structure in mantle convection: effects of plate motions and viscosity stratification,Geophys. Res., Lett.,23,2987–2990.

Cadek, O. & Fleitout, L., 1999. A geoid model with imposed plate velocities and partial layering,J. geophys. Res.,104,29 055–29 075.

Clift, P.D., Turner, J., Ocean Drilling Program Leg 152 Scientific Party, 1995. Dynamic support by the Icelandic Plume and vertical tectonics of the northeast Atlantic continental margin,J. geophys. Res.,100,24 473–

24 486.

Conrad, C.P., Lithgow-Bertelloni, C. & Louden, K.E., 2004. Iceland, the Farallon slab, and dynamic topography of the North Atlantic,Geology, 32,177–180.

Darbyshire, F.A., White, R.S. & Priestley, K.F., 2000. Structure of the crust and uppermost mantle of Iceland from a combined seismic and gravity study,Earth planet. Sci. Lett.,181,409–428.

Deschamps, F., Trampert, J. & Snieder, R., 2002. Anomalies of temperature and iron in the uppermost mantle inferred from gravity data and tomographic models,Phys. Earth planet. Inter.,84,245–264.

Farnetani, C.G., 1997. Excess temperature of mantle plumes: the role of chemical stratification across D,Geophys. Res. Lett.,1583–1586.

Feighner, M.A., Kellogg, L.H. & Travis, B.J., 1995. Numerical modeling of chemically buoyant mantle plumes at spreading ridges,Geophys. Res.

Lett.,22,715–718.

Forte, A.M. & Mitrovica, J.X., 2001. Deep-mantle high viscosity flow and thermomechanical structure inferred from seismic and geodynamic data, Nature,410,1049–1056.

Forte, A.M. & Peltier, R., 1991. Viscous flow models of global geophysical observables; 1, Forward problems,J. geophys. Res., 96,20 131–

20 159.

Foulger, G.R. et al, 2001. Seismic tomography shows that upwelling beneath Iceland is confined to the upper mantle,Geophys. J. Int.,146,504–

530.

Gripp, A.E. & Gordon, R.G., 1990. Current plate velocities relative to the hotspots incorporating the NUVEL-1 global plate motion model,Geo- phys. Res. Lett.,17,1109–1112.

Haxby, W.F. & Turcotte, D.L., 1978. On isostatic geoid anomalies, J. geo- phys. Res.,83,5473–5478.

Heller, D.-A. & Marquart, G., 2002. An admittance study of the Reykjanes Ridge and elevated plateaux between the Charlie-Gibbs and Senja fracture zones,Geophys. J. Int.,148,65–76.

Ishii, M. & Tromp, J., 1999. Normal mode and free air gravity constraints on lateral variations in velocity and density in the Earth’s mantle,Science, 285,1231–1236.

Ito, G., Lin, J., & Gable, C.W., 1996. Dynamics of mantle flow and melting at a ridge-centered hotspot: Iceland and the mid-Atlantic ridge,Earth planet.

Sci. Lett.,144,53–74.

Ito, G., Shen, Y., Hirth, G. & Wolfe, C.J., 1999. Mantle flow, melting and dehydration of the Iceland mantle plume,Earth planet. Sci. Lett.,165, 81–96.

Karato, S., 1993. Importance of anelasticity in the interpretation of seismic tomography,Geophys. Res. Lett.,20,1623–1626.

Keen, C.E. & Boutilier, R.R., 2000. Interaction of rifting and horizontal plume sheets at volcanic margins, J. geophys. Res., 105,13 375–

13 387.

Larsen, T.B. & Yuen, D.A., 1997. Fast plumeheads—temperature dependent versus non-Newtonian rheology,Geophys. Res. Lett.,24,1995–1998.

Lemoine, F.G. et al, 1998. The Development of the Joint NASA GSFC and NIMA Geopotential Model EGM96,NASA/TP-1998-206861, NASA Goddard Space Flight Center, Greenbelt, MD.

Marquart, G., 1989. Isostatic topography and crustal depth corrections for the Fennoscandian geoid,Tectonophysics,169,67–77.

Marquart, G., 1991. Interpretations of geoid anomalies around the Iceland hotspot,Geophys. J. Int.,106,149–160.

Marquart, G., 2004. Inferring mantle viscosity andS-wave-density conversion factor from new seismic tomography and geoid data,Geophys. J. Int., submitted.

Marquart, G., Schmeling, H., Ito, G. & Schott, B., 2000. Conditions for plumes to penetrate the mantle phase boundaries,J. geophys. Res.,105, 5679–5694.

Montelli, R., Nolet, G., Dahlen, F.A., Masters, G., Engdahl, E.R. & Hung, S.-H., 2004. Finite-frequency tomography reveals a variety of plumes in the mantle,Science,303,338–343.

Mueller, R.D., Roest, W.R., Royer, J.-Y., Gahagan, L.M. & Sclater, J.G., 1996.A Digital Age Map of the Ocean Floor,SIO References Series 93-30, Scripps Institution of Oceanography, La Jolla, CA.

National Geophysical Data Center, 1988. Digital Relief of the Sur- face of the Earth,Data Announcement 88-MGG-02, NOAA, Boulder, CO.

Nicholson, H. & Latin, D., 1992. Olivine tholeiites from Krafla, Iceland:

evidence for variations in melt fractions within a plume,J. Petrol.,33, 1105–1124.

Ribe, N.M., Christensen, U.R. & Theissing, J., 1995. The dynamics of plume- ridge interaction, 1: Ridge-centered plumes,Earth planet. Sci. Lett.,134, 155–168.

Ricard, Y. & Wuming, B., 1991. Inferring the mantle viscosity and its three- dimensional structure from geoid, topography and plate velocities,Geo- phys. J. Int.,105,561–571.

Ricard, Y., Richards, M.A., Lithgow-Bertelloni, D. & LeStnuff, Y., 1993. A geodynamic model of mantle density heterogeneity,J. geophys. Res.,98, 21 895–21 909.

Richards, M.A. & Hager, B.H., 1984. Geoid anomalies in a dynamic earth, J. geophys. Res.,89,5987–6002.

Ruedas, T., Schmeling, H., Marquart, G., Kreutzmann, A. & Junge, A., 2004. Temperature and melting of a ridge-centred plume with application to Iceland. Part I: Dynamics and crust production,Geophys. J. Int.,158, 729–743.

Schilling, J.-G., 1991. Fluxes and excess temperature of mantle plumes inferred from their interaction with migrating mid-ocean ridges,Nature, 313,187–191.

Stacey, F.D., 1992.Physics of the Earth,3rd edn, Brookfield Press, Brisbane.

(13)

Steinberger, B., Schmeling, H. & Marquart, G., 2001. Large-scale lithospheric stress field induced by global mantle circulation,Earth planet.

Sci. Lett.,186,75–91.

Walzer, U., Hendel, R. & Baumgardner, J., 2003. Variation of non- dimensional numbers and a thermal evolution model of the Earth’s mantle, inHigh Performance Computing in Science and Engineering ’02,eds Krause, E. & Jaeger, W., Springer, Berlin.

Wessel, P. & Smith, W.H.F., 1995. New version of the Generic Mapping Tools released,EOS, Trans. Am. geophys. Un.,76,329.

White, R.S., 1997. Rift-plume interaction in the North Atlantic.,Phil. Trans.

R. Soc. Lond.,A,355,319–339.

Wolfe, C.J., Bjarnarson, I.T., VanDecar, J.C. & Solomon, S.C., 1997.

Seismic structure of the Iceland mantle plume, Nature, 385, 245–

247.