# A first detailed look at the Greenland lithosphere and upper mantle, using Rayleigh wave tomography.

F.A. Darbyshire, T.B. Larsen, K. Mosegaard, T. Dahl Jensen, O.

Gudmundsson, T. Bach, S. Gregersen, Helle A. Pedersen, W. Hanka

► To cite this version:

F.A. Darbyshire, T.B. Larsen, K. Mosegaard, T. Dahl Jensen, O. Gudmundsson, et al. A first detailed look at the Greenland lithosphere and upper mantle, using Rayleigh wave tomography. Geophysical Journal International, 2004, 158 (1), pp.267-286. 10.1111/j.1365-246X.2004.02316.x. hal-00109931

HAL Id: hal-00109931

https://hal.science/hal-00109931v1

Submitted on 5 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A first detailed look at the Greenland lithosphere and upper mantle, using Rayleigh wave tomography

Fiona A. Darbyshire,<fnref n="1" />* Tine B. Larsen,<fnref n="1" /> Klaus Mosegaard,<fnref n="2" /> Trine Dahl-Jensen,<fnref n="3" /> Ólafur Gudmundsson,<fnref n="4" /> Torben Bach,<fnref n="5" /> Søren Gregersen,<fnref n="1" /> Helle A. Pedersen<fnref n="6" /> and Winfried Hanka<fnref n="7" />

Accepted 2004 April 20. Received 2004 April 18; in original form 2003 July 2

## SUMMARY

We present the results of a surface wave study carried out across Greenland as part of the ‘GLATIS’ (Greenland Lithosphere Analysed Teleseismically on the Ice Sheet) project. Rayleigh wave phase velocity dispersion curves were estimated for 45 two-station paths across Greenland, using data from large teleseismic earthquakes. The individual dispersion curves show characteristics broadly consistent with those of continental shields worldwide, but with significant differences across the Greenland landmass. Reliable phase velocity measurements were made over a period range of 25–160 s, providing constraint on mantle structure to a depth of ~300 km.

An isotropic tomographic inversion was used to combine the phase velocity information from the dispersion curves, in order to calculate phase velocity maps for Greenland at several different periods. The greatest lateral variation in phase velocity is observed at intermediate periods (~50–80 s), where a high-velocity anomaly is resolved beneath central-southwestern Greenland, and a low-velocity anomaly is resolved beneath southeastern Greenland.

The results of the phase velocity inversion were used to construct localized dispersion curves for node points along two parallel north–south profiles in southern Greenland. These curves were inverted to obtain models of shear wave velocity structure as a function of depth, again with the assumption of isotropic structure. A similar inversion was carried out for two two-station dispersion curves in northern Greenland, where the resolution of the phase velocity maps is relatively low.

The models show a high-velocity ‘lid’ structure overlying a zone of lower velocity, beneath which the velocity gradually increases with depth. The ‘lid’ structure is interpreted as the continental lithosphere. Within the lithosphere, the shear wave velocity is ~4–12 per cent above global reference models, with the highest velocities beneath central-southwestern Greenland. However, the assumption of isotropic structure means that the maximum velocity perturbation may be overestimated by a few per cent. The lithospheric thickness varies from ~100 km close to the southeast coast of Greenland to ~180 km beneath central-southern Greenland.

**Key words:** Greenland, lithospheric structure, Rayleigh waves.

## 1 INTRODUCTION

Continental Greenland has previously only been sparsely investigated, despite Greenland’s great size and the fact that it may hold the key to essential questions about the pre-history of the Iceland

<footnote n="1">¹National Survey and Cadastre, Rentemestervej 8, DK-2400 Copenhagen NV, Denmark</footnote>

<footnote n="2">²University of Copenhagen, Denmark</footnote>

<footnote n="3">³Geological Survey of Denmark and Greenland</footnote>

<footnote n="4">⁴Danish Lithosphere Centre, Denmark</footnote>

<footnote n="5">⁵University of Aarhus, Denmark</footnote>

<footnote n="6">⁶Université Joseph Fourier, Grenoble, France</footnote>

<footnote n="7">⁷Geoforschungszentrum Potsdam, Germany</footnote>

<footnote>*Now at: Geological Survey of Canada, 7 Observatory Crescent, Ottawa, ON, K1A 0Y3, Canada. E-mail: fiona@seismo.nrcan.gc.ca.</footnote>
plume and the opening of the North Atlantic. Most of Greenland is a Precambrian shield, including an Archaean craton in the south and several Proterozoic mobile belts. The continent has been affected by several orogenic and rifting episodes. The Greenland Shield is separated from the North American continent by the Labrador Sea, Davis Strait and Baffin Bay. The opening of the Labrador Sea at ~61 Ma (Chalmers & Pulvertaft 2001) created a complex passive margin on the southwest coast of Greenland. The east coast is also a passive margin, created by the opening of the North Atlantic at ~52 Ma, connected with break-up above the Iceland plume.

Roughly 80 per cent of Greenland is covered by a massive ice sheet, in places up to 3400 m thick, which poses a logistical challenge for seismological and other experiments. Seismic investigations have primarily focused on the continental margin and the continent-ocean transition (e.g. Chian & Louden 1992; Dahl-Jensen et al. 1998; Schlindwein & Jokat 1999), and only a few studies have provided information about the crust in the interior of Greenland (Gregersen 1970; Gregersen 1982; Gregersen et al. 1988; Dahl-Jensen et al. 2003). The Greenland lithosphere and upper mantle have not previously been the subject of a direct investigation, though some information can be discerned from global tomographic models (e.g. Su & Dziewonski 1994; Ritzwoller & Lovely 1995; Li & Romanowicz 1996; Masters et al. 1996; Ekström et al. 1997; Ritsema et al. 1999; Megnin & Romanowicz 2000; Shapiro & Ritzwoller 2002). In general, global models show higher than average shear wave velocities beneath Greenland, to a depth of ~200 km, with the highest velocities in the north of the continent. However, the modelled anomalies beneath Greenland may be influenced by the fast Canadian Shield to the west and the slow mid-Atlantic ridge/hotspot system to the east. Most global tomographic models

employ lateral smoothing constraints; in many cases, lateral scale resolution is of the order ~1000 km (e.g. Ritsema & Allen 2003). Lateral resolution is improved in the present regional study.

## 1.1 Project GLATIS

The GLATIS project (Greenland Lithosphere Analysed Teleseismically on the Ice Sheet) was initiated as a first detailed study of the Greenland Shield, contributing valuable information to the major tectonic questions concerning the region. The experimental set-up consisted of 16 temporary broad-band seismograph stations to supplement the four permanent broad-band stations already operating in Greenland (Fig. 1; Table 1). A description of data acquisition and instrumentation can be found in Dahl-Jensen et al. (2003). The station network layout was to a large degree dictated by logistics, and the remaining stations were placed in such a way that the network would be suitable for two-station surface wave analysis. Of particular interest to the study is the possible effect on lithospheric structure of the Iceland plume. Geological investigations place the region consisting of plume-like material between stations UPN and GDH on the west coast and stations SCO and SOE on the east coast, and link the appearance of flood basalts (the Lower Tertiary basalt sequence shown in Fig. 1) to the opening of the northern Atlantic (White & McKenzie 1989; Lawver & Müller 1994; Tegner et al. 1998; Larsen & Saunders 1998; Larsen et al. 1999). Knowledge of the structure of the Greenland lithosphere may provide information about the initiation and early history of the Iceland plume.

The objective of this paper is to describe and discuss the results of a fundamental-mode Rayleigh wave study using the GLATIS data set. The station configuration is too coarse to resolve a narrow

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Geological map of Greenland with seismograph stations used in the GLATIS project. Modified from Dahl-Jensen et al. (2003).</figcaption></figure>
<table><caption>**Table 1.** Seismograph stations used in the GLATIS project. Responsible institutions/projects are as follows: IRIS, Incorporated Research Institutions for Seismology; GLATIS, Greenland Lithosphere Analysed Teleseismically on the Icesheet; GEOFON, Geoforschungszentrum Potsdam, Germany; NEAT, northeast Atlantic Tomography (Danish Lithosphere Centre and University of Cambridge, UK); GSC, Geological Survey of Canada. The time periods shown in the final column of the table indicate the time periods from which teleseismic events have been extracted for this study, and do not necessarily reflect the full deployment times of the stations.</caption><thead><tr><th>Station name</th><th>Location</th><th>Project/ institution</th><th>Latitude</th><th>Longitude</th><th>Time period of data used</th></tr></thead><tbody><tr><td>ALE</td><td>Alert, Canada</td><td>IRIS</td><td>82.5033</td><td>-62.3500</td><td>1999-2002</td></tr><tr><td>ANG</td><td>Tasiilaq</td><td>GLATIS/NEAT</td><td>65.6160</td><td>-37.6370</td><td>2000-2002</td></tr><tr><td>DAG</td><td>Danmarkshavn</td><td>KMS/GEOFON</td><td>76.7714</td><td>-18.6547</td><td>1999-2002</td></tr><tr><td>DBG</td><td>Daneborg</td><td>GLATIS/NEAT</td><td>74.3080</td><td>-20.2139</td><td>2000-2002</td></tr><tr><td>IS1</td><td>Icesheet</td><td>GLATIS</td><td>71.9100</td><td>-31.0600</td><td>Not used here</td></tr><tr><td>IS2</td><td>Icesheet</td><td>GLATIS</td><td>69.1660</td><td>-44.7357</td><td>Summer 2000</td></tr><tr><td>IS3</td><td>Icesheet</td><td>GLATIS</td><td>68.9100</td><td>-31.5400</td><td>Not used here</td></tr><tr><td>DY2</td><td>Raven Camp</td><td>GLATIS</td><td>66.4740</td><td>-46.2639</td><td>Summer 2000</td></tr><tr><td>GDH</td><td>Qeqertarsuaq</td><td>GLATIS</td><td>69.2500</td><td>-53.5333</td><td>2000-2001</td></tr><tr><td>HJO</td><td>Hjørnedal</td><td>NEAT</td><td>70.3500</td><td>-28.1600</td><td>2000-2002</td></tr><tr><td>KAG</td><td>Kagssortoq</td><td>NEAT</td><td>63.2500</td><td>-42.0300</td><td>2001-2002</td></tr><tr><td>NGR</td><td>NorthGRIP</td><td>GLATIS</td><td>75.0010</td><td>-42.3148</td><td>Summer 2000</td></tr><tr><td>NOR</td><td>Station Nord</td><td>GLATIS</td><td>81.6047</td><td>-16.6602</td><td>2002</td></tr><tr><td>NRS</td><td>Narsarsuaq</td><td>KMS</td><td>61.1700</td><td>-45.4100</td><td>2000-2002</td></tr><tr><td>NUK</td><td>Nuuk</td><td>GLATIS/NEAT</td><td>64.1830</td><td>-51.7360</td><td>2000-2002</td></tr><tr><td>PAA</td><td>Paamiut</td><td>GLATIS</td><td>61.9914</td><td>-49.6613</td><td>2000-2001</td></tr><tr><td>SCO</td><td>Ittoqqortoormiit</td><td>KMS</td><td>70.4800</td><td>-21.9500</td><td>1999-2002</td></tr><tr><td>SFJ</td><td>Kangerlussuaq</td><td>KMS/GEOFON/IRIS</td><td>66.9967</td><td>-50.6156</td><td>1999-2002</td></tr><tr><td>SOE</td><td>Sødalen</td><td>NEAT</td><td>68.2000</td><td>-31.3800</td><td>2000-2002</td></tr><tr><td>SUM</td><td>Summit Camp</td><td>GLATIS</td><td>72.5763</td><td>-37.4538</td><td>2000-2002</td></tr><tr><td>TULE</td><td>Pituffik/Thule Air Base</td><td>GLATIS/GSC</td><td>76.4100</td><td>-68.5600</td><td>2000-2002</td></tr><tr><td>UPN</td><td>Upernavik</td><td>GLATIS/GSC</td><td>72.7848</td><td>-56.1408</td><td>1999-2000</td></tr></tbody></table>

plume; however, we have enough information to discern previously unresolved regional differences in the Greenland lithosphere and upper mantle. Surface waves are particularly well suited for a sparse network, and we have chosen to use teleseismic arrivals as the local seismicity is quite modest (Gregersen 1982). Our analysis is based on two-station phase velocity measurements of fundamental-mode Rayleigh waves, as the data did not allow for a reliable separation and analysis of higher modes. We study only the isotropic properties of the Greenland upper mantle; while this is obviously a severe assumption, work in Australia has shown that the main patterns of heterogeneity remain the same for both an isotropic and an anisotropic model based on surface wave tomography (Debayle & Kennett 2000).

The smaller and more difficult Love wave data set obtained from the GLATIS seismograph stations is currently under analysis, and the results will be presented in a future publication.

## 1.2 Surface wave studies of other continental shields

Large-scale broad-band seismological experiments around the world have provided valuable information about continental shield regions, e.g. Australia (van der Hilst et al. 1994; Debayle & Kennett 2000), Africa (James et al. 2001; Freybourger et al. 2001), North America (Grand 1994; van der Lee & Nolet 1997; Frederiksen et al. 2001), Finland (Bock et al. 2001; Bruneton et al. 2003) and Siberia (Priestley & Debayle 2003). The projects provide valuable information on the thickness and nature of the continental lithosphere. In Brazil a similar experiment revealed a structure interpreted as a fossil mantle plume beneath the shield (VanDecar et al. 1995).

In southern Africa, a model based on anisotropic surface wave inversion (Freybourger et al. 2001) gave shear wave velocities higher

than those of the preliminary reference earth model (PREM) to a depth of 220 km, though a corresponding study based on isotropic body wave tomography (James et al. 2001) suggested that the thickness of the lithosphere was 250–300 km. Beneath the Siberian Platform, Priestley & Debayle (2003) reported lithospheric thicknesses of 175–225 km, and maximum $S_v$ perturbations of ±7 per cent at 150 km depth. Debayle (1999) observed seismic heterogeneities of up to 13.5 per cent in Australia, though these amplitudes decreased when anisotropy was included in the analysis. The thickness of the Australian lithosphere was reported as ~200 km by Debayle & Kennett (2000) and 225 ± 50 km by Simons et al. (2002). A recent study of Canada by Frederiksen et al. (2001) showed a maximum range of $S_v$ heterogeneities of ±9 per cent, and a lithospheric thickness of 200–250 km depth. We compare our results with those from other shield regions, though we note that the differences in data types and modelling methods between the various studies will influence the values reported.

Studies of continental lithospheric thickness are influenced by the type of shear wave energy analysed. Recently, Gung et al. (2003) pointed out that models based purely on $S_v$ energy indicate continental roots extending to ~200–250 km depth, whereas $S_H$ and hybrid models generally yield continental roots to ~300–400 km depth. This is due to significant $S_v/S_H$ anisotropy in the 250–400 km depth range below most cratons. The discrepancies between $S_H$, hybrid and $S_v$ global models reviewed by Gung et al. (2003) indicate that it is important to take the type of shear wave energy into account when comparing different regional surface wave models. Our study is based on the analysis of $S_v$ energy alone.

In this study, we find lateral variations of up to ±4 per cent in Rayleigh wave phase velocity across Greenland. The amplitudes of the $S_v$ heterogeneities beneath the southern Greenland lithosphere are comparable with those reported in other shield regions, with fast anomalies of 4–12 per cent compared with global reference models.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Map of teleseismic events (stars) used in this study. Event magnitudes range from $M_s$ 5.5 to $M_s$ 8.0.</figcaption></figure>

However, the lithosphere appears to be a little thinner in Greenland than the average for shields worldwide, reaching a maximum thickness of <200 km across the region studied. The thickness and nature of the lithosphere appears to vary significantly across southern Greenland.

## 2 DATA COLLECTION AND PROCESSING

This study uses a data set made up of 121 large (magnitude 5.5 or greater) teleseismic events from the period 1999–2002 (Fig. 2). Initially, events were selected based on their location with respect to the Greenland seismograph stations. We use a two-station method for which the incoming surface wave energy should ideally lie exactly along the great circle path between a given pair of stations. In practice, events whose backazimuth was ≤5° from a given inter-station great circle path were considered to be appropriate for the analysis.

For each available station pair, the global earthquake database was searched for events within the correct backazimuth range, for the time period over which both stations were operating. Waveform data for the selected events were inspected visually for Rayleigh wave quality at a variety of different passbands. Events with clear Rayleigh waves visible on both stations were kept for further analysis. Fig. 3 shows a data example from an event in central America, recorded by one of the temporary stations in East Greenland.

Since the seismograph stations used in this study have a variety of digitizers and seismometers, it was important to correct the instrument responses for each pair of stations with different instrumentation. While the RefTek digitizers are zero-phase instruments, and simple to deal with in the data processing, care had to be taken when applying response corrections to the stations with Orion digitizers, as these instruments' DC filters caused significant phase distortions if not taken explicitly into account. At each station pair, the response of the broader-band station was removed from the seismograms by deconvolution and replaced with the response of the narrower-band station. The seismograms were then decimated, and filtered at a variety of passbands for phase dispersion analysis using a zero-phase Butterworth filter. Rayleigh wave data from the vertical-

component seismograms only are used for the phase velocity measurements.

The waveform data from the majority of GLATIS seismograph stations are stored in the archives of the GEOFON data centre, GFZ-Potsdam, in Germany. In addition, waveforms for stations SCO and NRS were obtained from the KMS archives in Copenhagen, and waveform data for stations ALE and SFJ were obtained from the IRIS Data Management Centre. The 2001–2002 data for the NEAT stations HJO, KAG and SOE were provided by the University of Cambridge. Data from ALE, NRS, SCO and SFJ are currently freely available to the scientific community, whereas data from the GLATIS and NEAT stations will be made available at a later date.

### 2.1 Quality control

In order to minimize errors in the dispersion calculations it was important to ensure that the parts of the waveform windowed for analysis were the result of simple, fundamental-mode Rayleigh wave energy arriving at the stations from the correct backazimuth. Initially, visual inspection of the waveforms was used to reject those which were not sampling the same structure between the source and the two stations making up the pair. For most of the paths crossing Greenland, the Rayleigh wave trains at the two stations were similar in appearance. However, for paths along the Greenland coasts, significant differences in the wave trains were visible. This implies that the Rayleigh wave energy was not travelling directly along the theoretical great circle path between the stations but was instead following different paths to each of the two stations. On the basis of these observations, it was decided to confine the analysis to paths crossing the Greenland mainland.

All the events used in this study originate at teleseismic distances from the Greenland stations, and the distance from the source to the stations is considerably greater than the interstation distance. The source-station group velocities for an event recorded at a given pair of stations should therefore be dominated by the source-to-Greenland path, with only minor contributions from the structure between the two stations. This assumption only holds true if the Rayleigh wave energy moves along the same great circle path to each of the stations. We calculated group velocity curves for the
<figure><img src="image_4.png" /></figure>

**Figure 3.** An example of the teleseismic data used in this study. The figure shows the three-component seismogram set (BHE, BHN, BHZ) for the 2000 July 21 (Ms 6.1) Costa Rica earthquake, recorded by stationDBG. The backazimuth of the earthquake from the station is 248° and the epicentral distance is 74°. Rayleigh waves on the vertical component seismogram clearly show their dispersive nature.

seismograms, using the frequency–time analysis method (FTAN) of Levshin *et al.* (1992), and compared the curves for the two stations. We found that those events for which the Rayleigh wave trains were visually similar also resulted in similar group velocity dispersion curves.

For events which were well-recorded on all three components of the seismometers, it was possible to assess the direction of the incoming Rayleigh wave energy and compare it with the theoretical backazimuth, calculated from the source and station coordinates. We used the polarization properties of Rayleigh waves, namely that the vertical component waveform is equal in phase to the inverse Hilbert transform of the radial component waveform. At incremental rotations of the horizontal components, the inverse Hilbert transform of the trial radial component was compared with the vertical component until the best waveform match was found. The analysis was carried out on filtered Rayleigh wave signals, bandpassed at a period of 30–60 s. This approach gave reliable results only for the clearest Rayleigh wave data. Records with reliable results showed small (generally <10°) deviations of the Rayleigh wave energy from the expected backazimuth. While it has been shown (e.g. Cotte *et al.* 2000) that deviations of Rayleigh wave energy from the expected great circle path can lead to errors in phase velocity dispersion curves, the small deviations measured here are unlikely to affect our results significantly.

When windowing the Rayleigh wave data for dispersion analysis, we were careful to include only that part of the seismogram for which simple dispersion was taking place. Multipathing of the Rayleigh wave energy results in a pattern of ‘beats’ in the later part of the

wave train. The analysis windows were chosen such that no energy arriving after the first ‘beat’ of the sequence was included.

# 3 ANALYSIS METHODS

## 3.1 Rayleigh wave dispersion

The dispersive properties of Rayleigh waves have been used in studies of Earth structure for several decades (e.g. Brune & Dorman 1963; Knopoff 1972; Gomberg *et al.* 1988; Debayle & Kennett 2000). Information may be gained from group velocity and phase velocity measurements, using single-station or two-station methods. In this study, we measure Rayleigh wave phase velocities between pairs of stations lying on a common great circle path from a given earthquake. Making measurements at two stations removes the amplitude and phase source terms, which are common to both seismograms, and thus eliminates the need for information about the source function of the earthquake. A two-station study also has the advantage that only the structure between the stations contributes to the dispersion curve, instead of the entire source–station path. This is a particularly important factor in this study, where we use data from events at teleseismic distances, several of which lie over 100° in epicentral distance from the Greenland seismograph network.

### 3.1.1 Theory

To obtain the phase velocity dispersion curves, we use a modified form of the transfer-function method of Gomberg *et al.* (1988). A
description of the method is given below (for further details see also Saunders *et al.* 1998).

If $S_1(t)$ and $S_2(t)$ represent the Rayleigh wave time-series at the stations nearest to and furthest away from the earthquake, then the ‘far’ seismogram can be represented as the ‘near’ seismogram convolved with an earth filter which describes the properties of the Rayleigh wave propagation between the two stations. In the frequency domain:

$$S_2(\omega) = F(\omega)S_1(\omega) \qquad (1)$$

where $F(\omega)$ is the earth filter:

$$F(\omega) = e^{-\gamma(\omega)\Delta r} e^{ik(\omega)\Delta r} = A(\omega)e^{i\phi(\omega)}. \qquad (2)$$

$A(\omega)$ is the amplitude term, comprising the attenuation parameter $\gamma(\omega)$ and the interstation distance $\Delta r$. $\phi(\omega)$ is the phase and $k(\omega)$ is the wavenumber, which is related to the phase velocity, $c(\omega)$, in the following way:

$$\phi(\omega) = k(\omega)\Delta r = \frac{\omega}{c(\omega)}\Delta r. \qquad (3)$$

A trial earth filter, $\tilde{F}(\omega)$, is constructed, based on an initial estimated dispersion relation. Use of the trial Earth filter gives a prediction for the form of the far seismogram, $\tilde{S}_2$, which is then compared with the observed data. The cross-spectrum $D(\omega)$ of the observed and predicted data is used to define and calculate a wavenumber correction vector, $\delta k(\omega)$, for the system, from which the phase term is obtained. The true earth filter may be written in the following form:

$$F(\omega) = D(\omega)\tilde{F}(\omega) \approx e^{-[\tilde{\gamma}(\omega) + \delta\gamma(\omega)]\Delta r} e^{[k(\omega) + \delta k(\omega)]\Delta r} \qquad (4)$$

where

$$D(\omega) = \frac{\tilde{S}_2^*(\omega)S_2(\omega)}{|\tilde{S}_2(\omega)|^2}. \qquad (5)$$

The tilde symbol denotes the trial values of the terms, and the asterisk in the cross-power spectrum denotes complex conjugation.

A waveform coherence estimate is made using the cross-spectrum $D(\omega)$, averaged over several frequency bins around the central frequency at which the measurement is made. The waveform coherence is used to weight the data equations at each frequency.

There is no linear approximation made to the $e^{i\delta k \Delta r}$ term in the modified version of the method of Gomberg *et al.* (1988) used here, unlike that used in the original method. The equation above can be solved for the wavenumber correction vector $\delta k(\omega)$ directly, which provides better stability and accuracy than the solutions calculated using the linear approximation, and enables the phase term to be calculated independently of the amplitude term.

The solution of the problem is subject to weighting of the data for their quality and coherence. Smoothing is imposed on the dispersion curve by requiring the group velocity associated with the results to lie close to the group velocity associated with the initial trial earth filter. The smoothing condition may be written as:

$$\delta k(\omega_{i+1}) - \delta k(\omega_i) = \pm \frac{\Delta \omega}{\tilde{U}(\omega_i)} \left( \frac{\delta U(\omega_i)}{\tilde{U}(\omega_i)} \right)_{\max} \qquad (6)$$

where $U$ is the group velocity and $(\delta U/\tilde{U})_{\max}$ is the maximum allowable deviation from the initial estimate of the group velocity, chosen by the user.

The wavenumber and hence the phase velocity for each frequency is found by solving the data equations for wavenumber perturbation, together with the smoothness constraints and weighting factors. The form of the equations is shown in detail in Gomberg *et al.* (1988).

The solution of the equations also produces a covariance matrix for the wavenumber, with diagonal elements $\sigma(k)_{ij}^2$ for the $i$th frequency and the $j$th seismogram pair.

For frequencies where the amplitude of the cross-power spectrum falls below a chosen fraction of the maximum amplitude there is no contribution made to the dispersion estimate. This reduces noise contamination and allows better control of the evaluation of the dispersion curve at any given bandwidth.

Several earthquakes lying within the correct range of great circle paths ($\le 5^\circ$ from the interstation backazimuth) for a given station pair may be analysed in this manner simultaneously, by minimizing the phase misfit of the waveforms of all events used.

Error estimates are also obtained for the dispersion curves. These are calculated from the covariance matrix of the wavenumber, and represent statistical errors which may be attributed to contamination of the Rayleigh waves by random noise. The standard deviation of the phase velocity at a frequency $i$ can be written as:

$$\sigma(c_i) = \frac{c(\omega_i)^2}{\omega_i} \sqrt{\sigma(k)_{ii}^2}. \qquad (7)$$

These statistical errors do not, however, take into account systematic errors such as those due to multipathing or off-azimuth arrivals, and should therefore be viewed as a minimum estimate of the errors in the dispersion curves.

### 3.1.2 Calculation of dispersion curves for the Greenland stations

Earthquake selection for the Greenland two-station paths resulted in coverage of 45 different paths from a total of 200 seismogram pairs. While some paths were constrained by data from just one event, others yielded high-quality seismograms from multiple events, with up to 16 seismogram pairs for a two-station path.

Each seismogram was filtered using a range of different passbands intended to enhance certain frequency bands. Typically, the following filters were used: 30–200 s, 50–200 s, 80–400 s, 100–500 s. Fig. 4(a) shows an example of the results of filtering the seismograms. The procedure was necessary because long-period signals, while significantly above the background noise levels, are much smaller in amplitude than the shorter-period (~20–50 s) signals. If all frequencies were analysed simultaneously, the relative contribution of the low-frequency signals would be insufficient to allow a meaningful calculation of the dispersion relation at long periods.

Consideration was given to the choice of starting model for the trial earth filter. Tests on the surface wave data using different starting models (typically with variation in the crustal thickness) were carried out; in general the resulting dispersion curves were most stable for a model based on the CANSD Canadian Shield model of Brune & Dorman (1963), with a 40 km thick crustal section. The results were generally insensitive to realistic changes to the mantle structure of the starting model.

For each passband, a dispersion curve was calculated for each seismogram pair on the given two-station path. Where multiple earthquakes were available, individual curves were calculated before simultaneous analysis, in order to check for, and remove, outliers and unstable results. The quality of the results was assessed by consideration of (1) the stability of the calculated dispersion curve, (2) the waveform fit between the seismogram at the near station, convolved with the best-fitting earth filter, to that of the seismogram at the far station (Fig. 4b), and (3) the value of the squared coherence. In addition, visual inspection of the dispersion curves was used to ensure
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Vertical component seismograms of the Costa Rica earthquake (2000 July 21, $M_s$ 6.1), recorded by station DAG and filtered at a range of different passbands. The seismograms show clear Rayleigh wave energy at periods greater than 100 s. (b) Illustration of waveform fitting using the transfer function method of Gomberg *et al.* (1988). The top two traces show the original seismograms and the bottom trace shows the waveform fit achieved when the near seismogram is filtered through the correct earth filter. The event used in this example is the same as that shown in (a) above, and a bandpass filter of 50–200 s has been applied to the data.</figcaption></figure>

that no $2\pi$ phase shifts ('cycle skipping') had occurred during the calculation. These shifts would result in unphysically high or low phase velocities in the dispersion curves.

Using the dispersion curve and the coherence values, the part of the dispersion curve that was well constrained by the procedure was extracted. This resulted in a set of dispersion curves for each two-station path, generally one or two per frequency band, which overlapped partially across the phase velocity/period space. Fig. 5 illustrates an example of this process.

The final dispersion curve for each two-station path was formed by taking the average of the dispersion curves from each frequency band. The error bars on the resulting curve take into account both the errors of the constituent dispersion curves as calculated from the covariance matrix and the spread of the sets of curves making up the final result. For each period, we calculated the average of the individual standard errors and the standard deviation of the phase velocity values, and selected the larger of these two values to represent the final error. Fig. 6 shows examples of the Greenland dispersion curves, plotted against the CANSD dispersion curve of Brune & Dorman (1963) for reference.

We were able to make further tests on the reliability of the dispersion curve results by calculating phase velocities for selected seismogram pairs using a different two-station method (Pedersen *et al.* 1994) in addition to the method described above. The method of Pedersen *et al.* (1994) applies a set of phase-matched filters to the seismogram pair, and uses the cross-spectrum of the two records to measure the interstation phase velocity. The Gomberg *et al.* (1988) method was used in this study primarily because it allows simultaneous analysis of several seismogram pairs. In all cases, the two resulting dispersion curves for each given seismogram pair were approximately the same, within error bars.

## 3.2 Mapping phase velocity variations across Greenland

### 3.2.1 Phase velocity inversion

We use an inversion routine to produce phase velocity maps from the information provided by the two-station dispersion curves at several different Rayleigh wave periods. An isotropic medium is assumed in the inversion. Phase velocity maps at eight different periods are shown in Fig. 7. At each period, the maps are coloured to indicate the deviation from the average phase velocity for that period, highlighting the lateral variations across Greenland.

The method used to calculate the phase velocity maps is a stochastic inversion with an assigned data covariance and model covariance (Franklin 1970). The inversion carries out the calculations subject to a reference velocity, $c_0$, supplied by the user which, in each case here, was chosen as the average phase velocity for each given period. A brief description of the method is given below.

Let $c_i$ represent the $i$th measurement of phase velocity along the $i$th path of length $L_i$. Assuming great circle path propagation, we can write

$$ \delta t_i \approx L_i(c_i^{-1} - c_0^{-1}) = \int_{\text{path}} \delta U(\mathbf{x}) dS \quad (8) $$

where $\delta U$ is a varying phase-slowness perturbation along the path, $S$ is the incremental path length and $\mathbf{x}$ is a position vector. The phase is described as a path integral of phase slowness (Woodhouse 1974). The local estimate of phase velocity is

$$ c^{-1}(\mathbf{x}) = c_0^{-1} + \delta U(\mathbf{x}). \quad (9) $$
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Example of dispersion curve computed for five earthquakes in the 50–200 s passband, for path ALE-ANG. Top: dispersion curve with error bars. The horizontal bar indicates where the results are well resolved and reliable. Bottom: waveform coherencies for the five earthquake pairs used in the analysis. (b) Composite dispersion curve for path ALE-ANG, composed of reliable results from several different passbands using several seismogram pairs in each case.</figcaption></figure>

If $c_0$ is the average phase velocity within the region of inversion then $\delta U$ has a zero expectation and the assumptions for the stochastic inversion hold. Discretization of the integral in eq. (8) gives

$$ \delta t_i = \int_{\text{path}} \delta U(\mathbf{x}) dS = G_{ij} \delta U_j(\mathbf{x}) \quad (10) $$

where $G_{ij}$ is the path length of the ith path in the jth model cell, $\delta U_j = \delta U(\mathbf{x}_j)$. The stochastic inverse estimate of the vectorized $\delta U$, $\mathbf{m}$, is

$$ \mathbf{m} = \mathbf{C}_m \mathbf{G}^T (\mathbf{G} \mathbf{C}_m \mathbf{G}^T + \mathbf{C}_d)^{-1} \delta \mathbf{t}. \quad (11) $$

Here, bold lower-case letters represent vectors and bold upper-case letters represent matrices. The matrices $\mathbf{C}_m$ and $\mathbf{C}_d$ are covariance matrices for the model, $\mathbf{m}$, and data, $\delta \mathbf{t}$, respectively. The matrix inversion in eq. (11) has the dimensions of data space, and one is free to choose a model-space discretization that is arbitrarily fine at a relatively low computational cost. In that sense, the model can be a continuum. Eq. (11) ignores the uncertainty of theory, i.e. the ray-theoretical integral in eq. (8).

We assume that data errors are independent and equal for all the data. Thus, the data covariance is a scaled identity matrix, $\mathbf{C}_d = \sigma_t^2 \mathbf{I}$, where $\sigma_t^2$ is the data variance. We also assume that the model covariance, $\mathbf{C}_m$, is stationary and Gaussian. Thus

$$ \mathbf{C}_m = \sigma_t^2 \mathbf{N}(\Delta, \alpha) \quad (12) $$

where $\sigma_t^2$ is the model variance, $\alpha$ is a correlation length for the model, $\Delta$ is the distance between two points in model space and $\mathbf{N}$ is a Gaussian function of $\Delta$ and $\alpha$ with maximum amplitude of 1 on the diagonal. Eq. (11) can then be rewritten as

$$ \mathbf{m} = \mathbf{N} \mathbf{G}^T (\mathbf{G} \mathbf{N} \mathbf{G}^T + \gamma^2 \mathbf{I})^{-1} \delta \mathbf{t} \quad (13) $$

where $\gamma = \sigma_t / \sigma_u$ can be regarded as a damping parameter whose units are in kilometres. Based on the relatively sparse path coverage throughout much of Greenland, we choose a value of $\alpha = 500$ km for the model correlation length.

The uncertainty of the individual phase velocity measurements is about 0.05 km s⁻¹ and the interstation distances of the station network in Greenland range from 350 to 2400 km. Taking a representative interstation distance of 1000 km, the uncertainty of the data, $\delta t$, is typically $\sigma_t = 3.25$ s for an average phase velocity of 4 km s⁻¹. If the variations in phase slowness in the model are of the order of 5 per cent, then we expect a standard deviation of the order of $\sigma_U = 0.0125$ s km⁻¹. We therefore expect a suitable choice of the damping parameter to be approximately 250 km. After experimenting with the trade-off between model uncertainty and resolution, we chose a constant value of $\gamma = 300$ km for the phase velocity inversion at all periods. The path coverage and the uncertainty of phase velocity measurements are similar throughout the period range of our study.

The phase velocity maps are parametrized as a set of cells, 1° in latitude by 2° in longitude, over a region covering the entire Greenland continent, from 60°N to 83°N and 70°W to 10°W. For each model cell, the inversion outputs a phase velocity and a measurement of resolution length associated with the cell. The resolution matrix, $\mathbf{R}$, for the inversion scheme may be defined as follows:

$$ \mathbf{R} = \mathbf{N} \mathbf{G}^T (\mathbf{G} \mathbf{N} \mathbf{G}^T + \gamma^2 \mathbf{I})^{-1} \mathbf{G}. \quad (14) $$

We collapse the complex information in the resolution matrix into simple measures of resolution length. This is defined as the half-width of a circularly symmetric Gaussian function with its maximum amplitude equal to the value of the diagonal of the resolution matrix, for which the integral over the model is unity. This is a simplistic measure of resolution because the resolution kernels may be centred away from the location in the model to which they correspond, and they may have a complex shape, e.g. be asymmetric about the model point. Nevertheless, this is a pragmatic way of representing resolution in a digestible way and will work well in regions of good resolution.
<figure><img src="image_7.png" /></figure>

**Figure 6.** Examples of the 45 phase velocity dispersion curves computed for Greenland. The label in the top left-hand corner of each plot refers to the two-station path (see Fig. 1 for station locations). In each case, the dispersion curve is plotted as a set of triangles, with error bars. The dashed line is the CANSD dispersion curve of Brune & Dorman (1963) and is included here as a reference to highlight the variations between the Greenland dispersion curves.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Phase velocity maps at eight different periods, created by a combination of phase information from the two-station dispersion curves. Model cells for which the resolution length exceeds 1000 km are masked in white. The colour scale shows deviations of up to ±4 per cent from the average phase velocity in each case.</figcaption></figure>
### 3.2.2 One-dimensional dispersion curves

The phase velocity maps (Fig. 7) provide us with important information about the lateral variations in structure across Greenland for a wide range of periods, but information about the shear wave velocity structure as a function of position and depth is more informative. To achieve this, phase velocity maps were calculated for a large number of periods (38 different periods between 20 and 170 s were evaluated). For each individual cell of the map, these data were used to create a new, 1-D phase velocity dispersion curve. For an ideal situation, where the two-station path coverage is both dense and homogeneous, it would then be possible to invert each of these new dispersion curves in turn and build up a 3-D model of the shear wave velocity structure. However, due to the differences in path coverage and resolution across Greenland, it was necessary to exercise caution in choosing which model cells were sufficiently well

resolved to provide reliable dispersion curves. Fourteen model cells in southern Greenland were chosen for further analysis after careful consideration of the resolution of the phase velocity maps. These are shown in Fig. 8, plotted on top of the 1-D resolution length for each cell, and the path coverage of the original two-station dispersion curves.

Two factors were considered when selecting individual model cells for further analysis:

(1) Model resolution length at different periods. Cells for which the model resolution length was large (>500 km) at several periods were immediately rejected. However, a good resolution does not necessarily imply a reliable model cell, due to the fact that the resolution length is a 1-D value describing a 2-D quantity. Several cells in northwestern Greenland between stations ALE and NGR, for example (Fig. 8), have resolution lengths of <300 km, but the phase

<figure><img src="image_9.png" /></figure>

**Figure 8.** Map showing dispersion curve path coverage (purple lines) and cell-by-cell model resolution for the phase velocity inversion procedure (squares, coloured according to resolution length). Seismograph stations are shown as red triangles. Each model cell is assigned an identification number. The numbered circles indicate the model cells which have been chosen for further analysis, based on the resolution length and the density of paths crossing within the cell. The numbers are coloured according to the cell resolution length. The heavy purple lines between stations ALE, DAG and NGR indicate two-station paths chosen for further analysis. No single-cell analysis was possible in northern Greenland due to the lack of crossing paths in the region.
<figure><img src="image_10.png" /><figcaption>Figure 9. Dispersion curves constructed for individual cells in the phase velocity maps of Fig. 7. The cells are identified by numbers which, for each of the two profiles, increase with increasing latitude. See Fig. 8 for exact cell locations.</figcaption></figure>

velocity values are in fact smeared out along the entire two-station path, as there are no crossing paths in this region.

(2) Density of crossing paths. Due to the smearing effect described above, it is important to take into consideration the directions of the dispersion curve paths passing through each model cell. Only cells which have a relatively short (<~400 km) resolution length and through which several dispersion curve paths cross can be considered to be reliable.

The new one-dimensional dispersion curves constructed for the chosen model cells are shown in Fig. 9.

### 3.3 Modelling shear wave velocity structure

We calculated models of shear wave velocity versus depth for the 1-D model cells described above, and also for two two-station paths in northern Greenland, where resolution was insufficient to justify the formation of cell-based dispersion curves. The models were obtained by the iterative linearized least-squares inversion method (surf96, Herrmann & Ammon 2002). The choice of the general form of the initial starting model was found to be important; in particular, a realistic increase in shear wave velocity with depth was necessary in the deeper sections of the model in order to prevent oscillations and edge effects. We chose to use the iasp91 model of Kennett & Engdahl (1991) as a basis for the inversion procedure, as this model has a smooth velocity-depth relationship throughout the depth range of interest. The iasp91 model is a global reference model, but because

of the predominance of continental seismic observatories over oceanic observatories it has a continental bias which makes it suitable for a continental region such as Greenland. The starting model was parametrized as a stack of 20 km thick layers down to a depth of 600 km, in which the shear wave velocity was allowed to vary. The layers were parametrized as spherical shells, using the Earth-flattening transformation of Schwab & Knopoff (1972). The only significant deviation from the iasp91 model was in the crust. The crustal section in the starting model consists of four layers, with a gradual increase in shear wave velocity from the surface to mantle-type velocities at 40 km depth. The 40 km crustal thickness represents an approximate average for Greenland (Dahl-Jensen et al. 2003). We avoided using a large velocity contrast at the Moho, as this tends to create artefacts in the sub-Moho layers of the model. Tests carried out indicated that the exact nature of the crustal section of the starting model had a negligible effect on the output velocity models below depths of 60–80 km. Further details are given in the Appendix.

In order to reduce the dependence of the inversion results on the initial model, we used the perturbation method of Ammon et al. (1990) to create 40 starting models, seeded from the iasp91 model. Each of these models was inverted in turn, and the convergence of the 40 results was inspected (Fig. 10). An average result was calculated, which was used as the basis for further analysis. The spread of the inversion results across the model space, compared with the spread of the starting models, does not supply formal statistics about model resolution and uncertainty but provides useful qualitative
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Example of inversion of phase velocity dispersion curves, using data from model cell 222 (see Fig. 8). Left: crust/mantle models, showing standard earth models (broken black lines), starting models formed from perturbation of the iasp91 model (Kennett & Engdahl 1991), the full set of solution models (solid grey lines) and the average solution (black line). The models are constrained by the phase velocity data in the depth range ~70–300 km. Right: corresponding dispersion curves.</figcaption></figure>

information about the depth range for which the models are well con-
strained. In addition, we take into account the approximate Rayleigh
wave sampling depths for the period range for which we have dis-
persion data. The models appear to be reliable in the depth range
~70–300 km, although tests (see Appendix) suggest that we have
some constraint to depths of approximately 400 km.

In most cases, the inversion results matched the dispersion data
well within the estimated data error bounds. However, there were
some model cells for which the inversion result did not match the
curvature of the data dispersion curves or where the synthetic curve
lay outside the data error bars. Using the averaged inversion results
as a basis, forward modelling was carried out for all the data sets.
The principal aim of the forward-modelling procedure was to find
the simplest velocity models which matched the dispersion data
within error bounds. The positions of dispersion curve maxima and
minima with respect to period in the synthetic curves were shifted,
where necessary, by altering the depths of shear wave velocity max-
ima and minima in the structural models, until an adequate match
to the curvature of the dispersion data was achieved. Forward mod-
elling was also used to test the reliability of the major features in the
models (see Appendix for details) and to demonstrate that these fea-
tures were necessary to provide an adequate match to the dispersion
data.

Rayleigh wave behaviour is influenced by both *P* and *S* waves,
but is primarily sensitive to the shear wave velocity structure of the

region over which the waves pass. The velocity models calculated by
surf96 are the result of inversion for shear wave velocity structure;
the $V_p/V_s$ ratio of the output models is constrained to be the same
as that of the input starting models. In this study, a simple ratio of
$V_p/V_s = \sqrt{3}$ is used for the starting models. surf96 then calculates
density from $V_p$ using the Nafe–Drake relationships.

4 RESULTS AND DISCUSSION

4.1 Two-station dispersion curves

In total, our phase velocity measurements resulted in 45 new two-
station dispersion curves across Greenland, each based on signals
from one to 16 earthquakes. The typical period range for the curves
is ~25–160 s, though some measurements for periods as short
as 20 s and as long as 200 s were possible in a few cases. A
subset of eight of the Greenland dispersion curves is shown in
Fig. 6. The initial observation is that all of the Greenland disper-
sion curves lie within approximately 0.1 km s⁻¹ of the CANSD
curve (Brune & Dorman 1963) for their entire period range. Com-
paring the form of the dispersion curves with those of ‘typical’
examples from a range of tectonic regimes (e.g. Knopoff 1972;
Priestley & Brune 1978) shows that the Greenland data set is best
described as ‘shield-like’ in its general characteristics, as expected
given the tectonic setting. Within this broad classification, however, there are distinct differences between the dispersion curves for different two-station paths. At long periods, the phase velocity differs by up to 0.2 km s⁻¹ from path to path (significantly above the estimated errors at these periods), and the period at which the curve flattens out after its steep rise at the shorter periods also differs significantly.

Due to the long paths sampled by most of the dispersion curves (interstation distance ranges from 360 to 2500 km), it is not possible to make detailed structural interpretations based on the variations in the characteristics of the curves. However, there are some indications in the data of large-scale structural variations. In particular, several of the paths sampling primarily the southwestern region of Greenland show phase velocities in excess of the CANSD reference at intermediate and long periods, whereas paths sampling between southwest Greenland and the central-eastern coastal region (stations SOE, HJO and SCO; see Fig. 1 for locations) tend to show lower velocities than CANSD for most of the period range of the dispersion curve.

## 4.2 Phase velocity maps

Maps of phase velocity variation are shown in Fig. 7. The patterns of isotropic lateral heterogeneity vary considerably with Rayleigh wave period, suggesting that the patterns of structural variation across Greenland vary with depth. The path coverage and resolution estimates associated with the maps (Fig. 8) indicate that the southern half of Greenland is well resolved in these maps, but the northern half of Greenland is poorly resolved.

At the shortest periods (<40 s) there is little lateral variation across the southern half of Greenland on a ±4 per cent scale. At a 40 s period variations are also small; the most prominent feature is a zone of relatively low phase velocity at the central east coast. Between periods of 50 and 80 s a strong pattern of phase velocity variation is observed in southern Greenland. A high phase velocity anomaly of amplitude 3–4 per cent, centred close to the west coast, dominates the pattern, with highest amplitudes in the 60–70 s period range. Low phase velocities are again observed in the eastern coastal region, and an east–west trending band of lower than average phase velocities can also be observed to the north of the high-velocity anomaly. In northern Greenland there is a general trend of higher than average phase velocities.

At periods >80 s the southern high-velocity anomaly is not observed. Lower than average phase velocities are present beneath the central-eastern coastal region and beneath the southern tip of Greenland, though the latter in particular lie on the edge of the model resolution. Some indications of high phase velocities in northwestern Greenland are visible, although model resolution is relatively poor in this region. At the longest periods (>140 s) phase velocities are relatively invariant across much of Greenland, the exception being a low phase velocity anomaly in the southeast.

The phase velocity maps do not show any strong northwest–southeast trending features in the southern half of Greenland, and hence no obvious correlation with the postulated plume track across the continent (e.g. Lawver & Müller 1994). Instead, there is a significant difference between the phase velocities beneath the regions of Tertiary flood basaltts on the west and east coasts of Greenland, despite the fact that comparable amounts of melt have been generated at approximately the same time on both sides of Greenland. This suggests that the flood basalt emplacement is not a dominant effect on the velocity structure of the upper mantle of southern Greenland.

The lack of any northwest–southeast track in the phase velocity maps might imply either that the impact of the Iceland plume had little effect on the Greenland lithosphere or that the lithosphere was affected by the plume but that the effect was spread out approximately evenly over the whole of southern Greenland.

## 4.3 One-dimensional dispersion curves in southern Greenland

We have constructed 1-D phase velocity dispersion curves for sets of model cells in southern Greenland, using results from the inversion of two-station phase velocity data at 38 different periods in the range 20–160 s. The aim is to use these localized curves to obtain velocity–depth information on the structures causing the lateral heterogeneities observed in the phase velocity maps. The model cells selected lie along two parallel profiles, one running south–north in southwestern Greenland and the other running south–north 8° in longitude to the east. Comparison of the dispersion curves for the selected model cells shows significant variation along each profile and between the two profiles (Fig. 9). Most of the along-profile variation in the west occurs in the period range 40–120 s, and can clearly be correlated with the high phase velocity anomaly observed in the phase velocity maps (Fig. 7). The maximum phase velocities in the 50–80 s period range are observed at cell 282, with values of ~4.23 km s⁻¹.

The eastern profile also shows some along-profile variation, though this is not as pronounced as that observed along the western profile. The phase velocities in the 50–80 s period range are somewhat lower than those observed on the western profile, with maximum values of ~4.15 km s⁻¹ at cell 286.

In the northern half of each profile we observe that the dispersion curves steepen with increasing latitude, with relatively low phase velocities in the 50–80 s period range and relatively high (up to ~4.4 km s⁻¹) phase velocities at periods >130 s. We observe a further steepening of the dispersion curves to the north of the cells displayed in Fig. 9, although the resolution of these cells is not sufficient to justify detailed analysis.

## 4.4 Shear wave velocity structure across Greenland

Inversion and forward modelling of the 1-D phase velocity dispersion curves in southern Greenland results in a set of models of shear wave velocity against depth for the region. All models assume isotropic structure. The models, plotted against the standard Earth models iasp91 (Kennett & Engdahl 1991) and PREM (Dziewonski & Anderson 1981), are shown in Fig. 11. In all cases, the models show a high-velocity ‘lid’ structure overlying a lower-velocity region, beneath which the shear wave velocity gradually increases. The thickness and the maximum shear wave velocity of the lid structure both vary significantly from cell to cell. We assume that the lid can be interpreted as the seismological lithosphere, and attempt to constrain its thickness by consideration of the velocities in the lower section of the lid. Other studies of the seismological lithosphere have approached this problem in various ways. Gung et al. (2003) describe the lithosphere as a region of seismic velocities 1.5–2 per cent higher than the average of global S-wave tomography models. Frederiksen et al. (2001) use the contour of 2 per cent deviation from the global average velocity to define the base of the lithosphere. In contrast, Debayle & Kennett (2000) and Priestley & Debayle (2003) use the strongest negative gradients in shear wave velocity as an indicator for the base of the lithosphere. In this study, the base of the
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Velocity–depth profiles for model cells shown in Fig. 8, aligned according to latitude. The solid black line shows the best forward model, and the dotted and dashed lines show two standard earth models (iasp91 (Kennett & Engdahl 1991) and PREM (Dziewonski & Anderson 1981)) for comparison. The identification number of the model cell is given in the bottom left-hand corner of each figure.</figcaption></figure>

seismological lithosphere is not well defined from an inspection of velocity gradients. If we use the contour of 2 per cent deviation from iasp91 to define the base of the lithosphere, we obtain values of ~160–180 km along the western profile and ~100–160 km along the eastern profile in southern Greenland from the best-fitting forward models.

The maximum shear wave velocity within the lithospheric lid varies both between the profiles and along the profiles. Comparing the velocities with the values presented in standard earth models we model velocity anomalies of 4–12 per cent above PREM, iasp91 and ak135 (Kennett *et al.* 1995) at depths of 80–140 km. The highest velocity anomalies are modelled at cells 252 and 282 on the western profile, and coincide with the intersection of the western profile with the high phase velocity anomaly observed in the phase velocity maps at periods of 50–80 s.

For the northernmost cells of each profile, the velocity–depth sections change in character compared with the cells further to the south. In particular, the lithospheric lid and low-velocity zone are less pronounced. The velocity structure modelled for these cells is

a reflection of the steepening of the dispersion curves (Fig. 9) compared with the more typically ‘shield-like’ form of the dispersion curves for the cells further south.

The final velocity models for the cells in the two south–north profiles in southern Greenland are shown together in Fig. 12. We have used a linear interpolation routine from the GMT software package (Wessel & Smith 1991) to display the models as contour maps of velocity against depth and along-profile. This is done under the assumption that the structures between adjacent cells change smoothly. The velocity scale is chosen such that the interpreted base of the lithosphere coincides with a clear transition in colour (green to yellow) with depth. The display in Fig. 12 allows for an easier comparison between the various velocity–depth models described in the previous section. We have poor constraint on the model from the surface down to a depth of approximately 60 km, therefore this top region is not included in the figure.

The contour plot allows us to compare more clearly the features of the velocity–depth models shown in Fig. 11. In particular, we can see
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Shear wave velocity contour plot formed by interpolation of the models shown in Fig. 11, using GMT (Wessel & Smith 1991) routines: (a) western profile, (b) eastern profile. The positions of the 1-D models are marked as stars, and labelled with the appropriate cell number (see Fig. 8 for locations). The two contour plots are aligned according to latitude, which increases from left to right.</figcaption></figure>

significant changes in lithospheric thickness along the profiles, as well as lateral variations in absolute velocity within the lithospheric lid and the low-velocity zone beneath it.

The western profile (Fig. 12a) reveals a consistent high-velocity lid with a thickness of ~160–180 km overlying a zone of lower velocities. The lowest sublid velocities are found in the southern part of the profile. The strongest feature along the western profile is a distinct high-velocity core within the lithosphere. This anomaly is located at the depth interval from approximately 90 to 150 km, and is centred around 69° N, 49° W, with a lateral extent of about 400 km.

Tectonically, the anomaly is located underneath the northern Nagssugtoqidian Orogen which is a Palaeoproterozoic deformation

belt resulting from a collision between the Archaean North Atlantic Craton to the south, and the Archaean Illulissat Craton to the north at approximately 1860 to 1840 Ma (van Gool *et al.* 2002). The region of highest lithospheric velocities correlates well with the location of a relatively high Bouguer gravity anomaly in southwestern Greenland (e.g. Dahl-Jensen *et al.* 2003).

The eastern profile (Fig. 12b) is markedly different from the western profile, despite a distance between the profiles of only a few hundred kilometres. A faint signal is left of the high-velocity core of the lithosphere, as expected from the phase velocity maps (Fig. 7), and the lithosphere thickens dramatically away from the near-coastal southern part of the profile, which lies close to the southeast Greenland continent-ocean transition zone. As in the western profile, the
zone of lower velocities underneath the lid is most pronounced to the south.

In both profiles the lithosphere appears to thin in the northernmost part at the same time as the zone of lower velocities beneath the lid begins to fade away. Unfortunately the resolution deteriorates to the north, and the inversion of the individual cells becomes less robust. North of the profiles there are not enough crossing paths to provide the resolution necessary to perform an inversion for individual model cells. The pattern of apparent lithospheric thinning in the northern parts of the profiles is not consistent with the boundary between Proterozoic blocks suggested by Dahl-Jensen et al. (2003); instead it appears to cut the profiles in a roughly west-east orientation. We do not see any clear signal either in the phase velocity maps or in the shear wave velocity models that can be attributed to a plume track in southern Greenland. However, our data coverage is insufficient to rule out any plume influence on the Greenland lithosphere.

In northern Greenland, the sparse phase velocity path coverage made it necessary to model two of the shorter (600–900 km inter-station distance) original two-station dispersion curves available for the region, instead of building up new curves from individual cells of the phase velocity maps. The velocity structures modelled from these two curves (Fig. 13) can represent only an approximation for

northern Greenland and do not reflect the lateral variations that are likely to exist. However, the models give us a useful first indication of the average regional structure. The thickness of the lithospheric lid is approximately 160–180 km, although the base of the lid is difficult to identify due to the low velocity-gradient in the depth range 120–200 km of the DAG–NGR model. In central–north Greenland (ALE–NGR) the dispersion data require a thinner low-velocity zone than that modelled for northeastern Greenland (DAG–NGR). Velocities within the lithospheric lid are comparable to the average values for the southern Greenland models, with an anomaly of ~7–8 per cent above the standard earth models.

## 4.5 Anisotropy

All the models presented in this study assume that the upper mantle of Greenland is isotropic in nature. However, several studies in continental shield areas (e.g. Debayle & Kennett 2000; Freybourger et al. 2001; Saltzer 2002) suggest that this is unlikely to be the case. Instead, we might expect complex patterns of anisotropy in the upper 100–150 km, related to past deformation which has been frozen into the lithosphere, and a simpler pattern of anisotropy below ~150 km, related to present-day deformation as a result of plate motion (Debayle & Kennett 2000; Gung et al. 2003). Since the treatment of anisotropy is beyond the scope of this study, we can only speculate on its nature beneath Greenland, but it is useful to examine how the results presented here might change if anisotropy was included in the study.

In Australia, the relationship between mantle heterogeneity and anisotropy has been studied in detail (Debayle 1999; Debayle & Kennett 2000). The patterns of mantle heterogeneity recovered by multimode surface wave tomography across Australia did not alter significantly depending on whether or not azimuthal anisotropy was included in the calculations. However, the amplitudes of the heterogeneities were significantly altered; in an isotropic inversion, the amplitudes were increased by almost 2 per cent at 100 km depth (Debayle & Kennett 2000). Based on this comparison, the most likely effect of anisotropy on the results presented in this study would be to decrease the apparent shear wave velocities in the lithospheric lid.

# 5 CONCLUSIONS

Heterogeneities in the Greenland lithosphere have been mapped out for the first time, using fundamental-mode Rayleigh waves. The coarse sensitivity of the previous surface wave study (Gregersen 1970) was unable to distinguish these details. Isotropic phase velocity maps show phase velocities varying significantly across Greenland, with a large high-velocity anomaly near the west coast, and lower than average phase velocities near the east coast. Based on linear isotropic inversion in selected model cells, we find that the Greenland Shield consists of a seismological lithosphere ~180 km thick or less, slightly thinner than the average for continental shields, overlying a zone of lower velocities, comparable to the values from reference models such as PREM and iasp91. The zone of lower velocities appears to fade in the northern end of the well-resolved area, north of approximately 70°N. The anomalously high shear velocities found in central-southwestern Greenland are likely to originate from an old, stable core of the lithosphere, whereas the lower velocities and thinner lithosphere in the southeast could be plume related.

The amplitudes of the heterogeneities presented in this study are likely to decrease if anisotropy were to be taken into account; however, the patterns of heterogeneity are likely to remain stable,

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Velocity-depth profiles for two two-station paths in northern Greenland. See Fig. 11 for plotting conventions, and Fig. 1 for the station locations.</figcaption></figure>
based on results from Australia (Debayle & Kennett 2000). Al-
though we have been unable to find a clear signal from the Iceland
plume under Greenland, we have produced a first-order model of
the central Greenland lithosphere which can be used as a basis for
further studies.

ACKNOWLEDGMENTS

The GLATIS project is funded by the following institutions: the
Danish Council for Research Policy (grant 9802989), the Bureau of
Minerals Resources under the Government of Greenland, and the
instrument pool and archive facilities at GFZ-Potsdam, Germany.
The NEAT stations are funded by the Natural Environment Research
Council (UK) and the Danish Lithosphere Centre, and are provided
by SEIS-UK. Stations TULE and UPN use instruments provided by
the Geological Survey of Canada. We thank all those who have been
involved in the deployment and maintenance of the seismograph
stations in Greenland. Frederik Horn and Mathilde Bøttger Sørensen
were responsible for receiving and handling much of the data sent
to KMS from Greenland, and NEAT 2001–2002 data were provided
to us by Sylvana Pilidou at the University of Cambridge, UK. We
thank Peter Voss, Stefan Bernstein, James Chalmers and Jeroen
d van der Gool for valuable discussions; also Nikolai Shapiro and an
anonymous reviewer for their comments.

REFERENCES

Ammon, C.J., Randall, G.E. & Zandt, G., 1990. On the non-uniqueness of receiver function inversions., *J. geophys. Res.*, **95**, 15 303–15 318.

Bock, G. et al., 2001. Seismic probing of Fennoscandian lithosphere, *EOS, Trans. Am. geophys. Un.*, **82**, 621 628–629.

Brune, J. & Dorman, J., 1963. Seismic waves and earth structure in the Canadian Shield, *Bull. seism. Soc. Am.*, **53**, 167–210.

Bruneton, M., Pedersen, H.A., Farra, V. & the SVEKALAPKO Seismic Tomography Working Group, 2003. 3-D surface wave tomography in the central Baltic Shield, *Geophys. Res. Abstracts*, **5**, 10 187.

Chalmers, J.A. & Pulvertaft, T.C.R., 2001. Development of the continental margins of the Labrador Sea: a review, in, *Non-volcanic Rifting of Continental Margins: a Comparison of Evidence from Land and Sea*, Geological Society of London Special Publication 187, pp. 77–105, eds Wilson, R.C.L., Whitmarsh, R.B., Taylor, B. & Froitzheim, N., Geological Society of London, London.

Chian, D. & Louden, K., 1992. The structure of Archean-Ketilidian crust along the continental shelf of southwestern Greenland from a seismic refraction profile, *Can. J. Earth Sci.*, **29**, 301–313.

Cotte, N., Pedersen, H.A., Campillo, M., Farra, V. & Cansi, Y., 2000. Off-great-circle propagation of intermediate-period surface waves observed on a dense array in the French Alps, *Geophys. J. Int.*, **142**, 825–840.

Dahl-Jensen, T., Thybo, H., Hopper, J.R. & Rosing, M.T., 1998. Crustal structure at the SE Greenland Margin from wide-angle and normal incidence seismic data, *Tectonophysics*, **288**, 191–198.

Dahl-Jensen, T. et al., 2003. Depth to Moho in Greenland: receiver-function analysis suggests two Proterozoic blocks in Greenland, *Earth planet. Sci. Lett.*, **205**, 379–393.

Darbyshire, F.A., 2003. Crustal structure across the Canadian High Arctic region from teleseismic receiver function analysis, *Geophys. J. Int.*, **152**, 372–391.

Debayle, E., 1999. *SV*-wave azimuthal anisotropy in the Australian upper mantle: preliminary results from automated Rayleigh waveform inversion, *Geophys. J. Int.*, **137**, 747–754.

Debayle, E. & Kennett, B.L.N., 2000. The Australian continental upper man-
tle: Structure and deformation inferred from surface waves, *J. geophys.
Res.*, **105**, 25 423–25 450.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**, 297–356.

Ekström, G., Tromp, J., & Larson, E.W.F., 1997. Measurement and global
models of surface wave propagation, *J. geophys. Res.*, **102**, 8137–8158.

Franklin, J.N., 1970. Well-posed stochastic extensions of ill-posed linear problems, *J. Math. Anal. Appl.*, **31**, 682–716.

Frederiksen, A.W., Bostock, M.G. & Cassidy, J.F., 2001. S-wave velocity structure of the Canadian upper mantle, *Phys. Earth planet. Inter.*, **124**, 175–191.

Freybourger, M., Gaherty, J.B., Jordan, T.H. & the Kaapvaal Seismic Group, 2001. Structure of the Kaapvaal craton from surface waves, *Geophys. Res. Lett.*, **28**, 2489–2492.

Gomberg, J.S., Priestley, K.F., Masters, T.G. & Brune, J.N., 1988. The structure of the crust and upper mantle of northern Mexico, *Geophys J.*, **94**, 1–20.

Grand, S.P., 1994. Mantle shear structure beneath the Americas and sur-
rounding oceans, *J. geophys. Res.*, **99**, 11 591–11 621.

Gregersen, S., 1970. Surface wave dispersion and crust structure in Green-
land, *Geophys. J. R. astr. Soc.*, **22**, 22–39.

Gregersen, S., 1982. Seismicity and observations of Lg wave attenuation in Greenland, *Tectonophysics*, **89**, 77–93.

Gregersen, S., Clausen, C. & Dahl-Jensen, T., 1988. Crust and upper mantle structure in Greenland, in, *Recent Seismological Investigations in Europe, Proceedings of the 19th General Assembly of the ESC*, pp. 467–469, Nauka, Moscow.

Gung, Y., Panning, M. & Romanowicz, B., 2003. Global anisotropy and the thickness of continents, *Nature*, **422**, 707–711.

Herrmann, R.B. & Ammon, C.J., 2002. *Computer Programs in Seismology: Surface Waves, Receiver Functions and Crustal Structure*, version 3.20. Saint Louis University, St Louis, MO.

James, D.E., Fouch, M.J., VanDecar, J.C., van der Lee, S. & Kaapvaal Seismic Group, 2001. Tectospheric structure beneath southern Africa, *Geophys. Res. Lett.*, **28**, 2485–2488.

Kennett, B.L.N. & Engdahl, E.R., 1991. Traveltimes for global earthquake location and phase identification, *Geophys. J. Int.*, **105**, 429–465.

Kennett, B.L.N., Engdahl, E.R. & Buland, R., 1995. Constraints on seismic velocities in the Earth from travel times, *Geophys. J. Int.*, **122**, 108–124.

Knopoff, L., 1972. Observations and inversion of surface wave dispersion, *Tectonophysics*, **13**, 497–519.

Larsen, H.C. & Saunders, A.D., 1998. Tectonism and volcanism at the SE Greenland rifted margin: a record of plume impact and later continental rupture, *Proc. Ocean Drilling Program Sci. Results*, **152**, 503–533.

Larsen, T.B., Yuen, D.A., & Storey, M., 1999. Ultrafast mantle plumes and implications of flood basalt volcanism in the Northern Atlantic Region, *Tectonophysics*, **311**, 31–43.

Lawver, L.A. & Müller, R.D., 1994. Iceland hotspot track, *Geology*, **22**, 311–314.

Levshin, A., Ratnikova, L. & Berger, J., 1992. Peculiarities of surface-wave propagation across central Eurasia, *Bull. seism. Soc. Am.*, **82**, 2464–2493.

Li, X.D. & Romanowicz, B., 1996. Global mantle shear velocity model developed using nonlinear asymptotic coupling theory, *J. geophys Res.*, **101**, 22 245–22 272.

Masters, G., Johnson, S., Laske, G. & Bolton, H., 1996. A shear-velocity model of the mantle, *Phil. Trans. R. Soc. London*, **354**, 1385–1411.

Megnin, C. & Romanowicz, B., 2000. The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms, *Geophys. J. Int.*, **143**, 709–728.

Pedersen, H.A., Campillo, M. & Balling, N., 1994. Changes in the litho-
spheric structure across the Sorgenfrei-Tornquist Zone inferred from dis-
persion of Rayleigh waves, *Earth planet. Sci. Lett.*, **128**, 37–46.

Priestley, K.F. & Brune, J.N., 1978. Surface waves and the structure of the Great Basin of Nevada and western Utah, *J geophys. Res.*, **83**, 2265–2272.

Priestley, K. & Debayle, E., 2003. Seismic evidence for a moderately thick lithosphere beneath the Siberian Platform, *Geophys. Res. Lett.*, **30**, 1118, doi:10.1029/2002GL015931.

Ritsema, J. & Allen, R.M., 2003. The elusive mantle plume, *Earth planet. Sci. Lett.*, **207**, 1–12.

Ritsema, J., van Heijst, H.J. & Woodhouse, J.H., 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland, *Science*, **286**, 1925–1928.
Ritzwoller, M.H. & Lavely, E.M., 1995. Three-dimensional seismic models of the Earth's mantle, *Rev. Geophys.*, **33**, 1–66.

Saltzer, R.L., 2002. Upper mantle structure of the Kaapvaal craton from surface wave analysis—a second look, *Geophys. Res. Lett.*, **29**, doi:10.1029/2001GL013702.

Saunders, P., Priestley, K. & Taymaz, T., 1998. Variations in the crustal structure beneath western Turkey, *Geophys. J. Int.*, **134**, 373–389.

Schlindwein, V. & Jokat, W., 1999. Structure of the continental crust of northern east Greenland from integrated geophysical studies. *J. geophys. Res.*, **104**, 15 227–15 245.

Schwab, F.A. & Knopoff, L., 1972. Fast surface wave and free mode computations, in, *Methods in Computational Physics 11, Surface Waves and Earth Oscillations*, pp. 87–180, ed. Bolt B.A., Academic Press, New York.

Shapiro, N.M. & Ritzwoller, M.H., 2002. Monte-Carlo inversion for a global shear-velocity model of the crust and upper mantle, *Geophys. J. Int.*, **151**, 88–105.

Simons, F.J., van der Hilst, R.D., Montagner, J.-P. & Zielhuis, A., 2002. Multimode Rayleigh wave inversion for heterogeneity and azimuthal anisotropy of the Australian upper mantle, *Geophys. J. Int.*, **151**, 738–754.

Su, W.J. & Dziewonski, A.M., 1994. Degree 12 model of shear velocity in the mantle, *J. geophys. Res.*, **99**, 6945–6980.

Tegner, C., Duncan, R.A., Bernstein, S., Brooks, C.K., Bird, D.K. & Storey, M., 1998. 40Ar-39Ar geochronology of Tertiary mafic intrusions along the East Greenland rifted margin: relation to flood basalts and the Iceland hotspot track, *Earth planet. Sci. Lett.*, **156**, 75–88.

VanDecar, J.C., James, D.E., & Assumpcao, M., 1995. Seismic evidence for a fossil mantle plume beneath South America and implications for plate driving forces, *Nature*, **378**, 25–31.

van der Hilst, R.D., Kennett, B.L.N., Christie, D. & Grant, J., 1994. SKIPPY: mobile broad-band arrays to study the seismic structure of the lithosphere and mantle beneath Australia, *EOS, Trans. Am. geophys. Un.*, **75**, 177–181.

van der Lee, S. & Nolet, G., 1997. Upper mantle S-velocity structure of North America, *J. geophys. Res.*, **102**, 22 815–22 838.

van Gool, J.A.M., Conelly, J.N., Marker, M. & Mengel, F.C., 2002. The Nagssugtoqidian Orogen of West Greenland: tectonic evolution and regional correlations from a West Greenland perspective, *Can. J. Earth Sci.*, **39**, 665–686.

Wessel, P. & Smith, W.H.F., 1991. Free software helps map and display data, *EOS, Trans. Am. geophys. Un.*, **72**, 441.

White, R.S. & McKenzie, D., 1989. Magmatism at rift zones: the generation of volcanic continental margins and flood basalts, *J. geophys. Res.*, **94**, 7685–7729.

Woodhouse, J.H., 1974. Surface waves in a laterally varying layered structure, *Geophys. J. R. astr. Soc.*, **37**, 461–490.

## APPENDIX: INVERSION AND FORWARD MODELLING—TESTS

When using a linearized inversion scheme like that of the surf96 software (Herrmann & Ammon 2002) it is important to take into account the dependence of the final results on the form of the starting

model, particularly in the regions of the model which are not well constrained by the phase velocity data. We parametrized our model as a stack of layers from 0 to 600 km depth; while direct constraint was only possible in the ~70–300 km depth range, the form of the velocity structure at greater depth had some influence on the dispersion curve at periods of 100 s and more. In addition, the large depth range decreased the likelihood of edge effects in the inversion results.

Of possible concern was the effect of crustal structure. The gradual structure we chose for the starting models reflected the average Moho depth across Greenland, and the fact that large velocity contrasts can give rise to artefacts in the linearized inversion in depth ranges where the data provides little constraint. Although the Moho depth in 20 locations across Greenland has been modelled (Dahl-Jensen et al. 2003), we have very little constraint on velocities within the crust. Some estimates of crustal velocities may be made by comparison with models of the northern Canadian Shield, but it is not clear that this is an accurate reflection of the structure of Greenland.

We tested the inversion using a variety of crustal models, such as the gradational model used in the main study, a model based on the structure of Baffin Island (Canadian Shield; Darbyshire 2003) and the model of Gregersen (1970). In addition, we explored the effect of placing an ice layer 1–3 km thick into the models. The inversion results confirmed that the data set provides no constraint above ~60–80 km depth, and that the choice of crustal model, including the presence of ice, has negligible effect on the mantle models.

We used forward modelling tests to demonstrate that the major features of the velocity models recovered from the inversions were robust, and necessary to match the dispersion curve data. An illustration of the tests is given in Fig. A1, using data from model cell 222 in southwest Greenland (see Fig. 8 for location). The figure illustrates the requirement for three major features in the model:

(1) A high-velocity ‘lid’ structure, controlled by the high phase velocities in the 40–65 s period range.

(2) A low-velocity zone in the 180–270 km depth range, controlled by the relatively flat section of the dispersion curve in the 65–85 s period range.

(3) A positive velocity gradient from 250 km depth to the base of the model, controlled by the gradient of the dispersion curve in the 90–160 s period range.

The removal of any of these major features resulted in a significant loss of fit to the dispersion curve data. In addition, visual inspection of the difference between the data and the synthetic dispersion curves calculated from the iasp91 and PREM models also illustrates the deviation of the Greenland mantle structure from that of standard earth models.

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>

**Figure A1.** Representative forward modelling tests showing the effect on the phase velocity curve of individual features in the velocity-depth model. Two standard earth models, iasp91 (Kennett & Engdahl 1991) and PREM (Dziewonski & Anderson 1981), are shown in each part as a reference: (a) Solid black line, best-fitting forward model; (b) solid black line, high-velocity zone removed; solid grey line, low-velocity zone removed; (c) solid black line, high-velocity zone and low-velocity zone averaged out to a depth of 320 km; solid grey line, structure below 280 km removed.