Seismic investigations on Gornergletscher

0 downloads 0 Views 54MB Size Report
From the S-wave velocity used for the synthetic seismograms .... 8000 years ago) were due to huge jökulhlaups on the Laurentide ice sheet (Clarke et al. (2004);. 9 ..... The processing and interpretation of the data acquired in these ...... For such investigations on icequakes the knowledge the Q-factor of ice is a key issue.
Seismic investigations on Gornergletscher Diploma thesis March 2007

Author: Valentin Gischig ETH Zurich, Institute of Geophysics

Supervisors: Fabian Walter, VAW Dr. Nicholas Deichmann, Swiss Seismological Service PD Dr. Hansruedi Maurer, Institute for Geophysics

2

Contents Abstract

5

Zusammenfassung

7

1 Introduction 1.1 The scope of this study

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9 10

2 The field campaigns 2004 to 2006 on Gornergletscher 2.1 Site description . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 The Gornergletscher . . . . . . . . . . . . . . . . 2.1.2 The Gornersee . . . . . . . . . . . . . . . . . . . 2.2 Previous field campaigns . . . . . . . . . . . . . . . . . . 2.3 The field campaigns 2004 to 2006 . . . . . . . . . . . . . 2.3.1 Glaciological methods . . . . . . . . . . . . . . . 2.3.2 Geophysical methods . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

13 13 13 13 15 15 16 17

3 The seismic velocity of glacier ice 3.1 Theory . . . . . . . . . . . . . . . . . . 3.1.1 Seismic waves and rays in ice . 3.1.2 Seismic refraction tomography 3.2 Refraction tomography profiles . . . . 3.2.1 Data acquisition . . . . . . . . 3.2.2 Tomography . . . . . . . . . . . 3.2.3 Error estimation . . . . . . . . 3.2.4 Results . . . . . . . . . . . . . 3.2.5 Interpretation . . . . . . . . . . 3.3 Borehole seismics . . . . . . . . . . . . 3.3.1 Data acquisition . . . . . . . . 3.3.2 Coordinates . . . . . . . . . . . 3.3.3 Seismograms . . . . . . . . . . 3.3.4 Seismic velocity estimation . . 3.3.5 Tomography . . . . . . . . . . . 3.3.6 Error estimation . . . . . . . . 3.3.7 Results . . . . . . . . . . . . . 3.3.8 Interpretation . . . . . . . . . . 3.4 Conclusion . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

19 19 19 23 26 26 29 29 30 34 37 37 38 40 40 40 42 43 43 45

3

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

4

CONTENTS

4 Seismic attenuation in glacier ice 4.1 Introduction . . . . . . . . . . . . . . . 4.2 Theory (from Aki & Richards (1980)) 4.2.1 Definitions . . . . . . . . . . . 4.2.2 The problem of causality . . . . 4.2.3 Q-factors from literature . . . . 4.3 The pulse width method . . . . . . . . 4.4 The refraction seismic profiles . . . . . 4.4.1 A two layer model . . . . . . . 4.5 Deep icequakes . . . . . . . . . . . . . 4.6 Conclusions . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

49 49 49 49 50 51 53 54 56 60 61

5 Synthetic seismograms 5.1 The reflectivity method . . . . 5.2 Procedure . . . . . . . . . . . . 5.3 Results . . . . . . . . . . . . . . 5.4 Interpretation . . . . . . . . . . 5.4.1 Crevasses . . . . . . . . 5.4.2 Velocity and attenuation 5.4.3 Anisotropy . . . . . . . 5.5 Conclusion . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . models . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

65 65 67 69 72 72 73 76 77

6 On 6.1 6.2 6.3 6.4

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

79 79 83 84 88

7 Conclusion and Outlook 7.1 Overview of the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91 91 92 93

Acknowledgements Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95 95

the scaling of deep icequakes Source parameters . . . . . . . Data . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . Interpretation . . . . . . . . . .

A Synthetic seismograms

. . . .

. . . .

. . . .

. . . .

101

Abstract Two refraction seismic profiles as well as several borehole surveys were used to determine the seismic velocity and attenuation structure of glacier ice in the ablation area of the Gornergletscher (Valais, Switzerland). The data of the refraction seismic profiles was used to compute P-wave tomography models of the uppermost 20 m of the glacier parallel and perpendicular to a local crevasse pattern. Furthermore, a layer model of the QP -factor was calculated with the pulsewidth method by Gladwin & Stacey (1974), which uses an empirical relationship between the increase of the pulse-width of a signal and the traveltime. Using the borehole survey data a 1Dmodel of the P-wave velocity down to a depth of 120m was calculated and its temporal variation during progression of melt season was investigated. These seismic velocity and attenuation models are required to determine the locations and the seismic moment tensors of deep icequakes in further studies. Such deep icequakes accompany the annual outburst flood (jökulhlaup), which is caused by the sudden drainage of the ice-dammed Gornersee via sub- and intraglacial channels. The investigation of these deep icequakes should give an insight into processes related to these jökulhlaups. The P-wave tomography models showed an increase of the seismic velocity with depth and low velocity zones, that were associated with crevasses. The 1D-models derived from these tomography models gave a P-wave velocity of 3100±100 m/s (perpendicular to the crevasses) and 3420±80 m/s (parallel to the crevasses) at the surface and 3700-3800 m/s at 20 m depth. From the borehole survey it was found, that for the velocity below 20 m a decrease of around 40m/s to a depth of 120 m is likely. Variations of the seismic velocity in the uppermost 20 m near the borehole are likely to be an effect of the rising lake level. For QP a value of around 25 was determined for the uppermost 7m of the ice. In general, the values measured perpendicular to the crevasses were slightly lower than those measured parallel, which can be considered as the attenuating effect of the crevasses. For the ice below a QP -factor ranging from 40 to 60 was estimated. In order to investigate the reliability of these velocity and attenuation models and to test the applicability of 1D-models in further studies, synthetic seismograms were calculated with the reflectivity method proposed by Mueller (1985). The measured seismograms of the profile parallel to the crevasses could be reproduced well with the 1D-model from the tomography. For the profile perpendicular to the crevasses it was found, that crevasses have a strong but rather local effect on the seismic velocity. Between the crevasses the velocity does not vary much from the one parallel to the crevasses. From the S-wave velocity used for the synthetic seismograms a VP /VS -ratio of 1.92±0.04 was calculated. QP was around 20 in the uppermost 10m for both profiles and 40-60 below. The strength of the seismic sources A0 was investigated for a cluster of icequakes located near the glacier bed. These icequakes are likely to originate from tensile cracks. The pulse-width 1/3 τ1/2 apparently does not scale with A0 which was expected for tensile cracks with constant 5

6

ABSTRACT

stress drop and slip velocity (Boatwright (1980)). Instead, stress drop and slip velocity increase with A0 . Furthermore, some evidence was found, that stronger icequakes might be more likely to occur at higher lake levels. This could be explained with the impact of the lake volume on the stress field near the lake.

Zusammenfassung Zwei refraktionsseismische Profile, sowie mehrere bohrlochseismische Untersuchungen wurden benutzt, um die seismische Geschwindigkeits- und Abminderungsstruktur des Gletschereises im Ablationsgebiet des Gornergletschers (Wallis, Schweiz) zu bestimmen. Mit den Daten der Refraktionsseismikprofile wurden P-Wellen-TomographieModelle für die obersten 20 m des Gletschers parallel und senkrecht zu den Spalten berechnet. Mit der Methode von Gladwin & Stacey (1974), welche einen empirischen Zusammenhang zwischen der Pulsbreite eines Signals und dessen Laufzeit benutzt, wurde ein Schichtmodell für den QP -Faktor berechnet. Die Bohrlochmessungen dienten dazu, ein 1D-Modell der P-Wellengeschwindigkeit bis zu einer Tiefe von 120 m zu bestimmen und deren zeitliche Variation wegen der zunehmenden Schmelze im Eis zu untersuchen. Die genaue Kenntnis der seismischen Geschwindigkeit und Abminderung des Gletschereises ist notwendig, um in späteren Studien die genauen Hypozentren und Momenten Tensoren von tiefen Eisbeben berechnen zu können. Solche tiefen Eisbeben treten im Zusammenhang mit den jährlichen Gletscherwasserausbrüchen (Jökulhlaups) auf, welche durch das plötzliche Abfliessen des eisgedämmten Gornersees durch intra- und subglazialen Kanäle verursacht wird. Untersuchungen an solchen Eisbeben können Hinweise auf die Vorgänge, die mit solchen Jökulhlaups einhergehen, liefern. Die P-Wellen-Tomographie Modelle zeigen eine Zunahme der Geschwindigkeit mit der Tiefe, und Zonen mit niedriger Geschwindigkeit, die mit Spalten assoziiert sind. Die 1D-Modelle, die aus den Tomographie-Modellen berechnet wurden, ergaben eine P-Wellen Geschwindigkeit an der Oberfläche von 3100±100 m/s (senkrecht zu den Spalten) und 3420±80 m/s (parallel zu den Spalten), sowie 3700-3800 m/s in 20 m Tiefe. Aus einer Bohrlochseismik-Messung in der Nähe des Sees ging hervor, dass möglicherweise die Geschwindigkeit unterhalb von 20 m um ca. 40 m/s abnimmt bis in eine Tiefe von von 120 m. Geschwindigkeitsvariationen in den obersten 20 m in der Nähe des Bohrlochs sind möglicherweise ein Effekt des steigenden Seespiegels. Für QP in den obersten 7 m konnte ein Wert von ungefähr 25 berechnet werden. Generell waren die Werte senkrecht zu den Spalten tiefer als parallel zu den Spalten, was auf den abmindernden Effekt der Spalten auf seismische Wellen zurückzuführen ist. In der Tiefe ist QP höher und liegt im Bereich zwischen 40 und 60. Um die Verlässlichkeit diese Geschwindigkeits- und Abminderungsmodellen zu untersuchen wurden synthetische Seismogramme berechnet, wobei die Reflektivitätsmethode von Mueller (1985) benutzt wurde. Die parallel zu den Spalten gemessenen Seismogramme konnten mit dem 1D-Modell aus der Tomographie gut reproduziert werden. Bei dem Profil senkrecht zu den Spalten konnte festgestellt werden, dass die Spalten einen starken, aber sehr lokalen Effekt auf die seismische Geschwindigkeit haben, während die Geschwindigkeit zwischen den Spalten kaum von derjenigen parallel zu den Spalten abweicht. Mit den S-Wellen Geschwindigkeiten, welche für die synthetischen Seismogramme verwendet wurden, konnte ein VP /VS -Verhältnis von 1.92±0.04 berechnet werden. Für QP wurden für beide Profile Werte von rund 20 in den obersten 10 m 7

8

ZUSAMMENFASSUNG

und 40-60 in grösseren Tiefen verwendet. Die Stärke A0 von seismische Quellen wurde für ein Eisbeben-Cluster in der Nähe des Gletscherbetts, welches wahrscheinlich von einem Dehnungsbruch stammen, untersucht. Es kon1/3 nte festgestellt werden, dass die Pulsbreite τ1/2 nicht mit A0 skaliert, wie es für einen Dehnungsbruch mit konstantem Spannungsabfall und konstanter Gleitgeschwindigkeit zu erwarten wäre (Boatwright (1980)). Stattdessen, nimmt der Spannungsabfall und die Gleitgeschwindigkeit zu mit A0 . Ausserdem konnte Hinweise gefunden werden, dass stärkere Eisbeben wahrscheinlicher sind bei höherem Seespiegel. Dies könnte auf den Einfluss des Seevolumens auf das Spannungsfeld in der Nähe des Sees hinweisen.

Chapter 1

Introduction Similar to earthquakes in the rigid earth, fracturing of ice can radiate seismic energy. In the same way as seismologists use signals from earthquakes to investigate processes and structures in the earth’s interior, these so-called icequakes can be used to obtain a better understanding of processes near the glacier bed, which are in general not well understood. Glaciers have been subject to seismic investigation since the 1950’s. Initially, mainly active seismic measurements were used to map the depth of the glacier bed or to measure physical properties of ice, such as the seismic velocity (e.g. Roethlisberger (1972); Dewart (1970); Thyssen & Ahmad (1969)). Neave & Savage (1970) were the first to detect icequakes and they associated them with the opening of crevasses. The absence of icequakes deeper than a few tens of meters was consistent with the model of brittle ice existing only at the shallow depth of glaciers. Weaver & Malone (1979) observed low frequency seismic signals caused by processes in the glacier on Cascade volcanoes covered with ice caps. The signal character of these events indicated, that they were associated with stick-slip motions on the glacier bed rather than with crevassing or volcanic activity. Investigations on Unteraargletscher (Deichmann, 1977) or Gornergletscher (Aschwanden, 1992) showed that a few of the recorded icequakes may also origin at depths beneath the crevassed layer. In recent years passive seismic surveys were carried out on many different kinds of glaciers. In general, four processes are proposed to cause seismicity on glaciers: surface crevassing (e.g. Neave & Savage (1970)), calving of icebergs (e.g. Qamar (1988), O’Neel et al. (2006)), stick-slip-motion near the bed (e.g. Wolf & Davies (1986))and hydraulic transients (e.g. St. Lawrence & Qamar (1979)). During three field campaigns in the summers of 2004, 2005 and 2006, various glaciological and geophysical methods were applied on the Gornergletscher, including the operation of seismic networks to detect icequakes. The aim was to investigate the annual outburst flood of the Gornersee near Zermatt, Switzerland. Every spring melt water accumulates at the confluence of the Gornergletscher and the Grenzgletscher thus forming the ice-dammed Gornersee. In the following summer the lake drains completely via subglacial and intraglacial channels, usually within a very short time period of a few days. Such outburst floods (jökulhlaups) are documented for all glacierized areas of the world (e.g. Raymond et al. (2003); Gudmundsson et al. (1997)) and present in general the highest and most far-reaching glacial risk. In the past they caused severe damage of the infrastructure built close to glaciers (e.g. in Zermatt, Raymond et al. (2003)). It could be shown that even abrupt climate changes during the last deglaciation (around 120008000 years ago) were due to huge jökulhlaups on the Laurentide ice sheet (Clarke et al. (2004); 9

10

CHAPTER 1. INTRODUCTION

Broeker & Denton (1990)). The sudden drainage of a large ice-dammed lake, covering an area of several states in the USA and Canada, forced the thermo-haline circulation in the Atlantic ocean to stop and initiated the cold climate of the Younger Dryas. Earlier studies on Gornergletscher showed that the lake drainage can be accompanied by icequakes near the glacier bed. During the recent field campaigns seismic networks were installed with the aim to detect such deep icequakes. Locating the hypocenters and computing moment tensor solutions can give an insight into the tension regime and processes acting during the opening of the draining channels. To increase the precision of the determination of these hypocenter properties an accurate model of the seismic body wave velocity structure and the seismic attenuation (Q-factor) of the glacial ice is required. In previous studies constant values for these parameters were used for the whole glacier ice (homogeneous halfspace). Yet Aschwanden (1992) showed that there is a gradient of the seismic velocity in the ice. It is also possible that there are variations of the seismic velocity due to anisotropy, pore fluid pressure and temperature variations in the ice. In the field campaign of 2005 test explosions in boreholes and at the surface confirmed that the seismic P-wave velocity is lower near the surface than at depth.

1.1

The scope of this study

To improve the seismic velocity model, active seismic surveys were carried out in 2006. Two refraction seismic profiles each with a length of 240 m were shot, both perpendicular and parallel to the strike of the crevasses. Furthermore, vertical profiles in boreholes with a maximum depth of around 120 m were shot 8 times from the top of the boreholes with intervals of around one week in between. In the present work this data is used to compute 2D-tomography models (refraction profiles) as well as 1D-tomography models (boreholes) of the P-wave velocity using the inversion code by Lanz et al. (1998). This provides information about: • the gradient of the P-wave velocity in the uppermost ice layer (refraction profiles), • the gradient of the P-wave velocity down to 120 m depth (borehole profiles) • anisotropy of the seismic velocity • variations due to changes of the pore fluid pressure. To investigate the seismic attenuation of the glacier ice, the "pulse-width method" proposed by Gladwin & Stacey (1974) is applied to the data from the active seismic surveys as well as from deep icequakes located outside the seismic array and recorded in 2006. The aim is to investigate the depth dependence and anisotropy of the Q-factor near the surface as well as the Q-factor at greater depths. Only the attenuation of the P-waves is considered here. The models for the P-wave velocity and the Q-factor are tested by comparing the seismograms of the refraction seismic survey with synthetic seismograms . The synthetic seismograms are calculated with the reflectivity method (Mueller (1985)), which uses a horizontal stack of layers to model the subsurface. This gives also information on the quality of a 1D-model approximation of the glacier ice. Finally, a cluster of deep icequakes located near the glacier bed is examined. The strength, the pulse-width, the amplitude and the released seismic energy of these events are investigated

1.1. THE SCOPE OF THIS STUDY

11

and related to each other. This gives an insight of how different source properties, such as stress drop, rupture velocity, etc., depend on the strength of the seismic source. The text is structured as follows. In chapter 2 an overview of the field area and some general properties of the Gornergletscher are given. Furthermore, the results of previous field campaigns and the results found within the Gorner project are summarized. In chapter 3 the seismic velocity structure of the glacier ice, as calculated from the refraction seismic profiles and the borehole seismic experiments, is presented. In chapter 4 the seismic attenuation structure of the glacier ice is investigated. The results of chapter 3 and 4 are tested using synthetic seismograms, results of which are presented in chapter 5. In chapter 6 a cluster of deep icequakes is investigated with a focus on the scaling of different source parameters. An overall conclusion and outlook is given in the chapter 7.

12

CHAPTER 1. INTRODUCTION

Chapter 2

The field campaigns 2004 to 2006 on Gornergletscher 2.1 2.1.1

Site description The Gornergletscher

The Gornergletscher is located in the southern Valais, Switzerland, near the village of Zermatt. With a length of 14 km and an area of 57.5 km2 (in 2003) it is the second largest ice mass of the Alps after the Aletschgletscher. Its accumulation area is on the Colle Gnifetti in the Monte Rosa massive which reaches an altitude of around 4600 m. With its accumulation area lying above 3000 m it has the highest equilibrium line of all Alpine glaciers. Its maximum thickness is around 400 m. The Gornergletscher is joined by the tributaries Grenzgletscher, Zwillingsgletscher, Schwärzegletscher, Breithorngletscher, Triftjigletscher and Unterer Theodulgletscher. Together they form a large valley glacier with the glacier snout at 2150 m above sea level. Because of the location in an extremely shadowy valley the snout is almost stable at present. During the first half of the 19th century this large glacier system advanced by about 600 m, until it reached its maximal length at the end of the little ice age in 1859. Since then the glacier has retreated about 3 km. Some of its tributaries (Triftjigletscher and Unterer Theodulgletscher) lost their connection to the main stream. The Gornergletscher is one of the few polythermal glaciers of the Alps. Polythermal glaciers consist of ice at its melting point (so-called temperate ice) as well as of ice colder than the melting point (cold ice). The cold ice of the Gornergletscher is contributed mainly by the Grenzgletscher and is accumulated above 4000 m. According to Haeberli (1976) the temperatures of the ice above 4200 m at a depth of 15-20 m is around -10 ◦ C. At 2600 m temperatures as low as -4◦ C were measured in the 1970’s. Since the cold ice formed in the absence of liquid water, it is more porous, less dense, and has a higher albedo than temperate ice (Haeberli (2002)). The polythermal and temperate parts of the glacier may be distinguishable on aerial images.

2.1.2

The Gornersee

The site for our investigations is located on the Grenzgletscher, a few hundred meters west of the Gornersee, and about one hundred meters south of the medial moraine formed by the Gornergletscher and the Grenzgletscher. This area lies at a height of around 2500m and has only minor topography. The Gornersee is formed each spring at the confluence of Gornergletscher and 13

CHAPTER 2. THE FIELD CAMPAIGNS 2004 TO 2006 ON GORNERGLETSCHER 14

Figure 2.1: Map of the Gornergletscher and its tributaries (From Atlas der Schweiz (2004))

2.2. PREVIOUS FIELD CAMPAIGNS

15

Grenzgletscher by melt water from the Monte Rosa area. At the time of maximal water level the lake contains 2 to 4 mio.m3 of water and covers an area of around 0.25 km2 . It is partly covered by icebergs. Around July the lake drains via subglacial channels underneath or in the ice within a short time interval of a few days. In former days these so-called outburst floods (jökulhlaups) caused significant property damage in Zermatt. Nowadays the water is pumped into the Zmutt storage lake by the power company Grande Dixence. During high discharges the water can be stored temporarily in this way.

2.2

Previous field campaigns

Several studies have been conducted on the Gornergletscher since around 1950. They were mainly motivated by the hydropower company Grande Dixence and by the hazard posed by the annual outburst flood for Zermatt. Other fields of interest were the subglacial drainage system, the temperature field in the ice, water pressure and icequakes. The results from studies which are relevant to the present work are summarized here: • Suesstrunk & Knecht (1949) The Grande Dixence power company supported a survey to determine the topography of the glacier bed using active seismic techniques. They determined the maximum ice thickness to be around 400 m at the confluence area. Their deduced bed-topography was used in most of the following studies. They estimated the seismic P-wave velocity to be 3600 m/s to 4000 m/s . • Roethlisberger (1972) developed a theory explaining the evolution of intra- and subglacial channels. This theory was applied to the Gornergletscher for the first time. • Haeberli (1976) In a borehole on the Grenzgletscher at 2600 m the temperatures were around -2 to -3 ◦ C in a depth of 50 m and 180 m below the surface. This proved that the glacier contains areas with ice below the melting point ("cold ice") and thus is polythermal. • Klingéle & Kahle (1977) tried to determine the height of the glacier bed using gravimetry. They found a smaller ice thickness than Suesstrunk & Knecht (1949). • Aschwanden (1992) processed seismic data acquired during a field campaign in 1979. Diurnal fluctuations of the seismicity and a maximum of seismicity after the outburst of Gornersee were found. Almost all events were associated with opening of crevasses. Only two events were found to originate from depths near the glacier bed. Furthermore, the seismic velocity was found to increase with depth near the surface (around 3710 m/s at depth and 3640 m/s at the surface). Aschwanden (1992) assumed that this could influence the precision of the locations of the events. It was suggested, that a more accurate velocity model instead of a homogeneous halfspace would give a better insight into the seismicity inside the glacier.

2.3

The field campaigns 2004 to 2006

In the years 2004 to 2006 a comprehensive set of glaciological and geophysical methods were applied to investigate the surface ice flow field, the water pressure inside the ice, the temperature field, the sub- and englacial water conduits, lake level and seismicity of the glacier. The aim of this project on the Gornergletscher is to better understand the processes triggering these annual

16

CHAPTER 2. THE FIELD CAMPAIGNS 2004 TO 2006 ON GORNERGLETSCHER

jökulhlaups. Therefore, the state of the glacier before, during and after the drainage has to be investigated. A better basic knowledge of glacier dynamics during jökulhlaups is essential for predictions of such events and for further glaciological research in general. Most methods were carried out during all three years to get a sufficient data and to refine results gained in previous years. In the following the applied methods as well as results achieved so far are summarized.

2.3.1

Glaciological methods

• Bed topography: The topography of the glacier bed in the study area was mapped with a dense Ground Penetrating Radar (GPR) survey. • Borehole temperatures: Boreholes were drilled with hot water technique in order to install instruments to measure water pressure, vertical strain and ice temperature. The temperature was found to decrease from 0 ◦ C at the surface to -0.2 ◦ C at 200 m depth, which are higher temperatures than those measured by Haeberli (1976). Thus the temperature effects of the cold ice seem to be negligible (e.g. for numerical models or for the seismic velocity). • Climate: An automatic climate station installed at the margin of the glacier measured air temperature, humidity and precipitation. This data serves as input for climate and mass balance models. • Ice deformation: More than 30 stakes were installed and their positions measured every hour with an automatic theodolite and differential GPS. Due to the presence of water at the glacier bed, the surface of the glacier rises during the outburst of the lake. Both, the uplift of the glacier as a result of the water pressure (buoyancy) and vertical strain of the ice due to horizontal compressing stress applied by the water of the lake causes this vertical motion. The vertical strain was measured via the length change of nearby boreholes. The uplift and flow patterns during the outburst showed that the ice might behave viscoelastically rather than purely viscously. Elastic deformation during the first half of the outburst in 2004 relaxed during the second half of the outburst (elastic rebound). • Monitoring of the lake: The lake level was measured with a water pressure transducer in the lake, which was calibrated with GPS measurements. At the same time the lake was monitored with an automatic digital camera installed above the Gornergletscher. At the snout of the glacier, hydrographs installed by the Grande Dixence company allowed a detailed analysis of the water outflow. The hydraulic conditions inside and beneath the glacier were investigated with dye tracing experiments conducted during the outburst. Interestingly, three different types of lake drainages were observed in all three years, which suggest that three distinct triggering mechanisms might have caused the outburst. In 2004 the lake began to fill on May 15 and reached a maximum volume of 4 · 106 m3 . It started to drain on the surface before it began to drain sub- and englacially on July 2nd. In 2005 the formation of the lake started on May 12 and reached only 1.2 · 106 m3 (water level 18 m lower than 2004). It drained subglacially (without surficial outflow) between June 10 and June 15. The differences in the hydrographs of the lake outflow of these two years suggest that the outburst was triggered differently. In 2004 the lake outflow discharge rose linearly, which can be explained with the initialization of the outburst by flotation of the ice block damming the lake. In contrast, in 2005 the lake outflow discharge rose exponentially,

2.3. THE FIELD CAMPAIGNS 2004 TO 2006

17

which is a so-called classical "slow" jökulhlaup where sub- and englacial channels open via pressure-melting (Nye (1976)). Both kinds of jökulhlaups were observed in the last tens of years (Huss (2005)) as well. However, the hydrographs at the gauging station of both years looked very similar. In 2006 again the lake performed in a different way. It started to fill late due to rather cold spring temperatures. The outflow between July 5 and around July 20 was unusual since the lake did not drain subglacially. Instead it filled until the lake level reached a nearby moulin and formed an almost 50m deep channel in the ice during the drainage. • Theoretical studies: Several theoretical studies are made too understand the physical behavior during the outburst floods. This includes numerical modeling of the subglacial ice conduit hydraulics (Spring & Hutter (1982)), of the subglacial water flow patterns (Flowers & Clarke (1999)) and of the ice flow (Blatter (1995); Gudmundsson (1999)).

2.3.2

Geophysical methods

• Passive seismics: The seismic activity of the glacier was monitored using seismic arrays of three component seismometers. In 2004 a seismic array of 13 surface seismometers (Lennartz, 1 Hz) and one borehole seismometer (28 Hz, at 100 m depth) were installed close to the lake (Figure 2.2). Seven seismometers each were connected to one of two GeodesT M , where the data were digitized with a sampling frequency of 1000 Hz and stored on a notebook. The GeodesT M were also responsible for synchroniz In both 2004 and 2005, between several hundred and a thousand events were recorded each day, most of which were surface events associated with the opening of crevasses. However, a significant fraction (more than one hundred) was located at greater depths. Most deep events belong to clusters of 10-40 icequakes. Several events occurred at intermediate depths. A comparison between the diurnal changes in the seismic activity and borehole water pressure close to the seismic array shows that there might be a strong correlation between water pressure at the glacier bed and seismicity. The locations as well as the source mechanisms of the deep events give insights into the brittle processes during the jökulhlaup, such as extent, orientation, source type (shear cracks or tensile cracks) etc. of the seismic sources. • Active seismics: So far, the locating algorithm models the glacier as a homogeneous full space using a constant seismic velocity. To test whether this assumption is justified several test explosions both near the surface and in boreholes were recorded in 2005. The velocity estimations resulting from that data, led to the conclusion that there exists a low velocity layer of 20 to 40 m thickness near the surface (Gischig (2005)). For this reason in 2006 active seismic experiments (i.e. two refraction seismic profiles and a borehole seismic experiment) were performed in order to determine the seismic velocity and attenuation structure of the ice in more detail. The processing and interpretation of the data acquired in these experiments are the main subject of this study and will be described in the following chapters.

18

CHAPTER 2. THE FIELD CAMPAIGNS 2004 TO 2006 ON GORNERGLETSCHER

Figure 2.2: Sketch of Gornergletscher and the positions of the seismic arrays. In 2004 and 2006 the arrays were located near the lake. In 2005 the array was installed further downstream. At the gauging station the water outflow from glacier is monitored.

Chapter 3

The seismic velocity of glacier ice 3.1 3.1.1

Theory Seismic waves and rays in ice

Physical description of seismic waves and seismic velocity (Aki & Richards (1980)) Seismic waves occur in every medium which responds elastically to applied stress to a certain degree. This phenomenon can be derived using the most general form of Newton’s law that balances linear momentum, surface forces and body forces in a small volume of a medium: ρ

∂ 2 ui ∂ = τji + fi , ∂t2 ∂xj

(3.1)

where ρ is the density of the medium, ui is the ith component of the position vector u, τji represents the stress tensor and fi is the corresponding component of the body force. In a perfectly elastic and isotropic medium Hooke’s law in the form τij = λkk δij + 2µij can applied, where the deformation tensor is ij = 1/2(∂ui /∂xj + ∂uj /∂xi ) and λ and µ are the Lamé constants. Inserting this in the above equation yields : ∂2u = (λ + 2µ)∇(∇ · u) − µ∇ × (∇ × u) + f (3.2) ∂t2 Expressing the displacement u with the potentials Ψ and Φ, u = ∇Φ + ∇ × Ψ, and neglecting all body forces f this equation can be separated after some simplifications into two wave equations. This shows that two kinds of seismic body waves exist in an elastic medium: ρ

∂2 (λ + 2µ) 2 Φ= ∇ Φ 2 ∂t ρ

(3.3)

∂2 µ Ψ = ∇2 Ψ (3.4) 2 ∂t ρ The first of these equation describes the propagation of P-waves, the second S-waves. P-waves are longitudinal waves with displacement in the direction of propagation. In contrast, S-waves are transversal and the displacement vector is perpendicular to the direction of propagation. The waves propagate with a phase velocity s λ + 2µ Vp = (3.5) ρ 19

20

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

in the case of a P-wave and in the case of a S-wave with r µ Vs = ρ

(3.6)

Reflection and refraction A widely used method to describe various seismic phenomena is the ray approximation: The wave is assumed to have an infinite frequency and penetrates the medium as a ray with a velocity Vp resp. Vs . Reflection and refraction of waves at seismic discontinuities can then be described using results known from ray optics. A ray traveling through medium 1 with properties λ1 , µ1 and ρ1 , meeting the boundary to a medium 2 with λ2 , µ2 and ρ2 is then partly reflected, partly refracted according to Snell’s law (Figure 3.1; V stands for both Vp and Vs ) sin α1 sin β1 sin β2 = = , V1 V1 V2

(3.7)

where V1 and V2 are seismic velocities of medium 1 and 2, α is the angle of incidence and β1 and β2 are the angles of the reflected and of the refracted ray, respectively. If V2 > V1 then β2 is larger than α. This means that there exists a critical angle of incidence αc for which β2 becomes 90◦ . The refracted wave then travels along the boundary between the media with the velocity of the "faster" medium 2 and continuously radiates energy back into medium 1. This is a so-called inhomogeneous wave (also head- or mintrop-wave). In a medium with a velocity that increases monotonically with depth, rays coming from a source at the surface are refracted in a manner that they are bent upwards and travel back to the surface. Such rays, as well as critically refracted rays, carry information about the velocity of the medium they have penetrated, i.e. their arrival times are influenced by the seismic velocity structure of the medium. This is used in refraction seismic techniques to image the velocity distribution in the subsurface.

Figure 3.1: Refraction and reflection of a ray at the boundary between medium 1 and medium 2. In red a critically refracted ray.

P- and T- dependence of seismic velocity in ice In glacier ice the elastic parameters λ and µ and the density ρ are functions of temperature, isostatic pressure, porosity as well as water and debris content. The seismic velocities are thus also

3.1. THEORY

21

dependent on these ice properties via equations 3.5 and 3.6. Unfortunately, these relationships between the velocity and those ice properties are difficult to determine and laboratory measurements on ice are often not applicable to glacier ice. Several investigations on different glaciers were made to establish empirical laws describing these effects. Many of them are reviewed by Roethlisberger (1972). Density: The density of glacier ice varies due to the change of porosity and of the material contained within the pores (air in cold ice or a mixture of water and air in temperate ice, Haeberli (2002)). According to Weber & Sandstrom (1960) the influence of density on the seismic velocity dominates all other effects. Of many empirical formulas given in the literature some are listed here: Weber & Sandstrom (1960): dVp m/s = 7.76 dρ 0.001kg/m3

(3.8)

Robin (1958) in Kohnen (1972): ρ = 0.059 + 0.221Vp =⇒ Kohnen (1972):

dVp m/s = 4.47 dρ 0.001kg/m3 ρE

ρ(z) = 1+



VE −V (z) 2.25

1.22

(3.9)

(3.10)

Here ρE is the density of non-porous ice, which is 0.917kg/m3 , VE is the velocity at this density (i.e. for clean ice) and V (z) is the velocity at depth z. This value is not determined uniquely in the literature. Water content: In temperate glaciers the temperature is always at the pressure melting point. Thus there is always some water content which influences the seismic velocity. In Thyssen (1966) the following formula is mentioned: 1 w 1−w = + Vp Vwater VE

(3.11)

which relates the water content w (volume fraction) to the seismic velocity. Accordingly, the seismic velocity decreases by 6m/s per % water content. The problem again is to get an accurate value for the velocity VE in clean ice. Navarro et al. (2005) used the formula proposed by Riznichenko (1949)  Vp = Vwater

100 − w 1+ w

   −1/2 2 (100 − w)ρice (100 − w)ρwater Vwater 1+ 1+ (3.12) 2 wρwater wρice Vice

where w is the water content in percent, to estimate the water content in ice. Though their estimation agreed well with measured values from ice cores, they stated that it might be a little too high. They used Vice = 3800m/s, ρice = 917kg/m3 , Vwater = 1500m/s and ρwater = 1000kg/m3 .

22

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Pressure: There is only very little information about pressure dependency in the literature. However, pressure is expected not to have a large direct effect on seismic velocity, but rather indirectly via the pore volume, which is effected by the isostatic pressure. Temperature: Measurements on single crystals gave gradients of −2.8m/s◦ C for a P-wave parallel to the c-axis and −1.7m/s◦ C normal to the crystallographic c-axis (valid for -20 ◦ C and -5 ◦ C). Similar results were found by Kohnen (1974) in cold ice of Greenland and Antarctica: −2.30 ± 0.17m/s◦ C. A higher gradient was found by Thyssen (1966) using data acquired on an alpine glacier. He gives the formula: Vp = 3750 − 4.5T (Vp in m/s, T in ◦ C). Other values for velocity temperature gradients are listed by Kohnen (1974). The variations of gradient values given in the literature are probably due to the fact that porosity might also vary systematically with the ice temperature. The reason for this is the strong influence the temperature during the formation of the ice has on the porosity. Generally these different gradients show that the effect of the temperature on the seismic velocity will be rather low on Gornergletscher, since minimum temperatures in the investigation area are expected to be only slightly below the pressure melting point. Anisotropy of ice On a scale of single crystals the assumption of isotropy as made in equation 3.2 is not valid anymore. Hooke’s law has to be applied in its most general form: τij = Cijkl kl . Cijkl is a tensor of fourth order and contains 21 independent elements in case of the lowest symmetry possible for a single crystal (i.e. less symmetry elements, such as rotation axes or symmetry planes). Single crystals of ice have a hexagonal symmetry, which means that rotating the crystals by 60◦ around the symmetry axis c leads to a congruent body. Thus all physical properties are different parallel to the c-axis than in the plane perpendicular to the c-axis. The tensor Cijkl has six independent elements in this case. The effects of hexagonal anisotropy in ice on the seismic velocity was studied for instance by Dewart (1970). He showed that the P-wave velocity parallel to the c-axis is 4% higher than normal to the c-axis. Other studies reported by Roethlisberger (1972) give values of around 4030 m/s and 3850 m/s parallel and normal to the c-axis, respectively. However, the assumption of isotropy is still valid in media consisting of an aggregate of single crystals, such as ice, in which the symmetry axes of the crystals are oriented randomly. However, since flow processes influence the structure of glacier ice, one can expect a certain degree of anisotropy that affects the seismic velocity. Four reasons for anisotropic behavior were proposed by several authors (e.g. Roethlisberger (1972) or Dewart (1970)): preferred orientation of the c-axis of the ice crystals, foliation due to layering of bubbly and bubble-free ice, preferred orientation of fractures, such as crevasses and material layering due to snow and firn layers. Since the investigation area of the present work lies in the ablation area, the first three points will be treated in more detail here: Preferred crystal orientation. If the c-axis of crystals in the ice are not oriented randomly, but within a cone with a certain opening angle γ, the ice mass behaves somewhat less anisotropically than in a single crystal. Computed effects for different opening angles amount to a 3.9% increase of the P-wave velocity along the axis of the cone compared to the velocity normal to it for γ = 0◦ , and decreases to 0.4% for γ = 60◦ . A more complete list is given by Roethlisberger (1972). Many different preferred crystal orientations exist in different glaciers. In West antarctic ice

3.1. THEORY

23

sheet vertical c-axes were found in central areas, where strain rates are low and horizontal c-axes near the margin of the ice sheet. In between the c-axes were mainly pointing into flow direction (Bentley (1971)). This shows that preferred c-axis orientations are highly coupled to deformation due to ice flow. Foliation. All kind of metamorphic structures due to deformation in ice are referred to as foliation. The most prominent foliation is the alternating layering of bubble-free and bubbly ice with thicknesses in the range of centi- to decimeters. These layers occur due to annual changes of the snow properties (temperature, water content, etc.) in the accumulation area which are preserved in the ice. Due to ice deformation the layering changes its horizontal orientation and can be oriented vertically or even folded. Trying to measure the c-axis orientation with seismic method, Dewart (1970) found that this kind of foliation can have an even larger effect on anisotropy than the preferred crystal orientation. Unexpectedly, he measured lower velocities in directions parallel to the c-axis in an area where the foliation was nearly perpendicular to the c-axis. Model calculations showed that, at a porosity of 0.7%, 2.8% and 6.2%, the velocity increase parallel to the foliation is 0.2%, 1.9% and 6.2%, respectively. The magnitude of this opposite effect is large enough to mask the crystal orientation effect, which was estimated to be around 2%. The effect is generally larger for small bubble sizes. Fractures. Though fractures can cause anisotropy, it is rather difficult to quantify this effect and to distinguish it from the previous two effects. Since open crevasses with the width in the decimeter range usually contain air, the seismic waves do not travel through them, but are scattered by them. They have a large travel path because they have to travel around the lower ends of the crevasses. Thus in general one can expect lower velocities for rays traveling perpendicular to fractures than for those traveling parallel to them.

3.1.2

Seismic refraction tomography

About tomography As the computing power of modern computers increased, tomography became an efficient and widely used method in many branches of science. It is most popular in medicine, where structures of tissue can be examined by screening it with X-rays. In geophysics tomography is used to investigate the earth’s interior on very different scales, that reach from the shallow subsurface down to the earth’s mantle and core (global tomography). Many different physical properties can be imaged with tomography using e.g. electrical currents, electromagnetic waves, active seismics, teleseismic events, etc. In our case actively generated seismic waves are used to image the seismic velocity distribution in the subsurface of Gornergletscher. Since in many cases the velocity increases with depth, the so-called refraction tomography was developed (e.g. Lanz et al. (1998)), which uses the fact that in such media rays released at the surface are bent upwards back to the surface, where they can be recorded. Such rays are only able to image positive velocity-depth gradients. Seismic tomography as a non-linear inverse problem The traveltime of a seismic signal penetrating a medium bears information about the velocity distribution in this medium. Assuming the raypath is a straight line between source and receiver,

24

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

then the velocity can be easily computed if the medium is sampled sufficiently by many rays. In the case of a simple block model as shown in Figure 3.2 the problem is linear, if we take the reciprocal value of the velocity v, the so-called slowness s.

Figure 3.2: A medium consisting of two blocks with two different seismic velocities v1 and v2 is penetrated by three rays. This medium is sampled sufficiently. Thus the problem is over-determined. The traveltimes can then be expressed as: ti =

2 X

sj lij , i = 1, 2, 3

(3.13)

j=1

whereas the ti is the traveltime of the ith ray, sj the slowness of the j th block and lij the portion of the raypath i within block j. Solving this equations for sj is a linear inverse problem and can be generalized for a medium with m blocks: ti =

m X

lij sj

(3.14)

j=1

or t = Ls

(3.15)

The matrix L again contains the portions of each raypath in each block. If the number of rays traveling through the medium exceeds the number of blocks in the medium and the slowness of each block can be determined uniquely, the problem is called over-determined. The problem can then be solved with the least-squares method, which minimizes the sum of squares of the differences between calculated and measured traveltimes: s = (LT L)−1 LT t

(3.16)

One can also think of cases where a part of the blocks are not sampled by a ray, whereas other blocks are sufficiently sampled. This is the mixed-determined case (mixture of over-determined and under-determined blocks). For those blocks, that are not or insufficiently sampled, so-called a priori information has to be included. There are different ways to do this. One way would be to force the blocks only to deviate from a priori given value if it is absolutely necessary. This is called damping. Another way is to prevent large variations of the values over three neighboring blocks (similar to "second derivative"), unless it is required by the measured traveltimes. This is referred to as smoothing. A priori information is included in the least-squares formula as follows: s = (LT L + pR)−1 LT t.

(3.17)

Here R is the damping or the smoothing matrix and p is the damping factor d or the smoothing factor s, which controls how much the a priori constraints are weighted against the measured data.

3.1. THEORY

25

So far it has been assumed that the raypath is independent of the velocity. As explained in Section 3.1.1 rays are bent due to varying velocities. The raypath and thus the portions of each ray in each block, as expressed in the matrix L, are themselves dependent on the velocity. This makes the problem strongly non-linear. There is no formula that solves this problem directly, and the problem needs to be solved iteratively. In the method used here one starts with a starting model s0 , which is a rough estimate of the slowness distribution in the medium. For this starting model the raypath as well as the traveltimes are computed. Since these traveltimes deviate from the measured data, one has to compute how much the starting model has to be altered in order to get a better fit to the data. Therefore the derivative of the traveltimes with respect to the slowness of each cell is computed: X ∂ti ∆ti = ∆sj (3.18) ∂sj j

Here ∆ti is the difference between the traveltimes calculated for the starting model and the measured ones. The desired deviations from the starting model are ∆sj and can be inverted using the same methods as for the linear problem using the matrix Lij = ∂ti /∂sj . The new model is then s1 = s0 + ∆s, which is used as a starting model for a second iteration. This goes on until the models converge to a model that is stable or until the a measure for deviations of the traveltimes from the data (in our case the rms-value defined as the squareroot of the mean of the squares of all deviations) falls below the error of the measurement. The matrix L has to be computed for each iteration. This requires first the calculation of the raypaths and secondly the computation of the derivatives. This is the so-called forward problem and needs to be solved numerically, which takes most of the computation time for each iteration. The code "inv2d" The C routine "inv2d" was developed by Maurer & Green (1997) to compute 2-dimensional velocity models using refraction seismic data to image buried waste disposals, which in general have lower seismic velocities than the underlying bed rock. It basically solves the above described non-linear inverse problem using both damping and smoothing. These additional constrains are referred to as regularization and are included in the inversion scheme as follows:     L t s (3.19) = D h h and D contain the regularization parameters λ and β; the equation h = Ds is equivalent to     λβs0 λβI = s. (3.20) 0 λ(1 − β)S The matrices I and S are the identity matrix (introduces damping) and the Laplacian smoothing matrix. The parameter λ controls the amount of regularization that is applied (i.e. how strong the constraints are weighted compared to the data), and β controls the relative weights of damping and smoothing. The vector s0 contains the starting model, which corresponds to a 1-D of velocity distribution. This can be obtained from an intuitive estimate of the velocity gradient or by first applying a 1-D-inversion scheme to the data. Having chosen a starting model, the raypath for each source receiver pair has to be computed. The "inv2d" routine contains a modified version of the "brute-force mapping" scheme described by Schneider et al. (1992). This calculates in a very robust way the arrival time of a wave

26

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

penetrating the medium for a dense grid of nodes. In the modified version, raypaths and the derivative of the traveltimes with respect to the slowness of each block can be extracted from these arrival times. This information is then used to build the matrix L. As input for the routine, the following information has to be given: filenames of data file and geometry file (coordinates of sources and receivers), number of iterations, smoothing and damping parameters s and d (used to compute λ and β), number of blocks, maximum model depth, number of nodes within a block, velocity of the surrounding area of the model (e.g. air) and the starting model. The smoothing constraints can be modified in such a way that the horizontal and the vertical smoothing is not equal.

3.2 3.2.1

Refraction tomography profiles Data acquisition

During the field campaign of 2006 on Gornergletscher, a refraction seismic survey was carried out on the 3rd and 4th of June (see Figure 3.3). Two nearly perpendicular profiles were measured, one striking nearly east-west and more or less perpendicular to the crevasses (profile 1), the other striking nearly north-south, parallel to the crevasses (profile 2). The aim was to investigate the seismic velocity distribution in the ice of the shallow subsurface. Further details are shown in table 3.1. A typical shot and the fourier spectrum of a typical trace is shown in Figure 3.4. The dominant signal is the Rayleigh wave. Despite the high quality signal of these signals only the Pwave is focused on in the present work. All traces are shifted with a traveltime corresponding to a seismic velocity of 3800 m/s. It shows that the real velocity is lower than 3800 m/s and is slightly inhomogeneous since the P-wave onsets deviate from the horizontal line. The inhomogeneities are also illustrated in Figure 3.5. All traces with a source-receiver distance of 12.5±1 m and the onsets are shown in Figure 3.5. The traveltimes vary up to 1ms from each other (=8% of the traveltime). This cannot be explained alone with uncertainties in the source-receiver distances or topography effects and have to be due to the inhomogeneities of the velocity near the surface (up to 300 m/s).

Receivers Source Recording length Trigger delay Sampling frequency Profile length Length of geophone spread Number of geophones Geophone spacing Number of shots (2 per shot position) Shot spacing

Profile 1 Profile 2 30 Hz vertical component geophones Pipegun 1s 0.3 s 8000 Hz 240 m 320 m 240 m 2 lines of each 24 geophones 5m 18 24 30 m (once only 18 m)

Table 3.1: Data acquisition parameters.

27 3.2. REFRACTION TOMOGRAPHY PROFILES

Figure 3.3: Aerial picture of the glacier at the confluence of the Gorner- and the Grenzgletscher taken in June 2006. All surveys carried out to measure the seismic velocity are located on the ice of Grenzgletscher which is separated from Gornergletscher by the debris cover (middle moraine) in the upper half of the picture. Profile 2 is nearly parallel to the strike of the crevasses, profile 1 crosses the crevasses. The lake close to the borehole is partly covered with icebergs.

28

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Figure 3.4: A typical common shot gather (shot 1 of profile 1, left) and the fourier spectrum of a typical signal (right). The traces are filtered and shifted in time using a velocity of 3800 m/s. The most prominent wave is a the Rayleigh wave. The P-wave arrival times show small deviations from the arrival time expected from a P-wave traveling with 3800 m/s, which would lie on a horizontal line. These deviations imply that there are inhomogeneities and a gradient in the velocity structure of the shallow ice.

Figure 3.5: All traces with a source-receiver distance within 12.5±1m for both profile 1 and profile 2.

3.2. REFRACTION TOMOGRAPHY PROFILES

3.2.2

29

Tomography

The dominant frequency of the recorded signals was found to be at around 50Hz (Figure 3.4). Thus prior to picking the data was filtered with a causal 4th-order butterworth filter suppressing frequencies below 20Hz and above 400Hz. The filtered data was then used to pick the onsets of the P-waves. This was done with an interactive matlab tool called "picktool.m", developed by Hansruedi Maurer at ETH Zurich. The error of picking was estimated to be at around 0.3ms. With the picked traveltimes the "inv2d" code was run. Since the result of the inversion may be dependent on the starting model and the regularization parameters (non-uniqueness), the following steps were undertaken in order to choose the best end model: • The inversion was run with various different starting models. For all starting models the rms-value converged and remained almost unchanged after the 20th iteration. Similar patterns of velocity variation in the end model appeared independently on the starting model. This shows that the inversion is stable and the end model is independent on the starting model. • To find appropriate regularization parameters the inversion was run for different smoothing and damping parameters. These parameters have to be set as low as possible (not too much weight on a priori information) and as high as necessary (no unrealistic end models). Too low values give rms-values below the picking precision of the traveltimes. In this case the inversion tries to fit noise. This results in patterns in the end models that are too rough or unrealistic (e.g. much higher velocities than 4000 m/s which is impossible for glacier ice). • Since not all areas of the subsurface were sampled equally by rays, some patterns may be a function of the ray coverage. To find out which patterns are caused by ray geometry or by really physical inhomogeneities, the inversion was run using only one half of the original dataset first and then for the other half. Patterns appearing on both models are assumed to be reliable and not an artifact of the ray coverage (see Figure 3.8). • Once the best model was found, outliers were removed from the dataset. These are traveltimes with a difference to the calculated traveltimes of higher than 1 ms or lower then -0.5 ms. Then the tomography was run for this new dataset resulting in the models presented here.

3.2.3

Error estimation

The traveltimes and the coordinates of the shot points and the receiver positions give a certain error, which causes uncertainties in the computed velocity model. Coordinates were determined using differential GPS with a high accuracy of less than 10cm error and thus do not cause significant uncertainties. In contrast the traveltimes were picked with an error of around 0.3 ms. In order to estimate the uncertainty of the velocity model due to this error, Monte Carlo simulations were run: Normally distributed random values with a standard deviation of 0.3 ms were added to each traveltime value and an inversion was run for this altered dataset. Repeating this 50 times gives 50 different velocity models. The standard deviation of these models gives an estimation of the error.

30

3.2.4

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Results

In the following the results of the inversion and the investigations described above are presented. Based on the steps described in Subsection 3.2.2, the velocity models in Figure 3.6A (profile 1) and in Figure 3.7A (profile 2) were decided to be the most reliable ones and were used for further investigations. After removing the outliers, model 1 had an rms-value of 0.34 ms and model 2 one of 0.45 ms. The inversion parameters used for these models are: Parameter number of iterations damping factor d smoothing factor s block size maximal depth velocity gradient in starting model velocity on surface in starting model maximal velocity in starting model

Profile 1 20 0.1 0.5 1.18 m 60 m 25 ms−1 m−1 3000 ms−1 4500 ms−1

Profile 2 20 0.1 0.5 1.60 m 60 m 10 ms−1 m−1 3400 ms−1 m−1 4000 ms−1 m−1

Table 3.2: Inversion parameters Figures 3.6B and 3.7B show the raypaths in the end models. Due to a rather low velocity gradient the rays penetrated only to a depth of around 20 m in both profiles. Only areas that were sampled by rays are interpretable. The starting models of both profiles with the raypaths calculated for these models are shown in Figures 3.6C and 3.7C. They consist of 1D velocity models with a constant gradient and a minimum value at the surface as given in table 3.2. The results of the Monte Carlo simulation are presented in Figures 3.6D and 3.6E, as well as in Figure 3.7d) and 3.7e). As expected, the model calculated by taking the average of the results from 50 inversions (Figures 3.6D and 3.7D) are slightly more smooth than the models from the original datasets. The standard deviation (Figure 3.6E and 3.7E) is greater at the surface, because the same error on traveltimes have a larger effect on short traveltimes than on long traveltimes. The standard deviation lies in a range of 20 m/s and 180 m/s for profile 1 and from 30 m/s to 130 m/s for profile 2. Areas where the standard deviations are large can be considered more sensitive to data errors than other areas. Furthermore, the models calculated from just half of the dataset (i.e. half of the shots) are shown in Figure 3.8. For these models the subsurface consists of only 100 blocks to account for the smaller amount of data used for these models. The inversion code generates the blocks between the outermost shot points, which are not the same for both halves of the dataset. Thus the axes of the related models are not comparable. Nevertheless, these models give insights into the influence of the ray coverage on the velocity patterns: Patterns appearing on both models can be considered as reliable and are not an artifact of the inversion.

3.2. REFRACTION TOMOGRAPHY PROFILES

31

Figure 3.6: Profile 1: A) Result of inversion, B) ray coverage, C) starting model, D) result from Monte Carlo simulation (mean), E) error estimation from Monte Carlo simulation (standart deviation).

32

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Figure 3.7: Profile 2: A) Result of inversion, B) ray coverage, C) starting model, D) result from Monte Carlo simulation (mean), E) error estimation from Monte Carlo simulation (standart deviation).

3.2. REFRACTION TOMOGRAPHY PROFILES

33

Figure 3.8: Profile 1 and 2 calculated with just half of the dataset each. Since the shot points at the beginning and at the end of the profiles are not identical for both halves of the datasets, the profiles do not cover the same area. The pictures are shifted, so that the same area lies just beneath each other.

34

3.2.5

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Interpretation

Crevasses Comparing profile 1 and profile 2 the most remarkable difference are the low velocity patches near the surface in profile 1, which intersects several crevasses. These stains were also visible to a certain degree on the models generated with just half of the dataset and thus cannot be artifacts of the inversion. Though these zones have no narrow steeply dipping shape as one would expect from crevasses, these patterns can be clearly associated with crevasses that are visible on the aerial picture of the glacier (Figure 3.9). Crevasses are characterized by low velocities not only because they are filled with air, but also because rays that cannot penetrate the crevasses have to travel around the crevasses and thus have longer traveltimes. The low velocity zones have velocities around 2800m/s or even lower and have sizes of several meters, though the crevasses are only maximal 50cm wide. This smearing is an artifact of the inversion. In the case of the crevasses in this model the amount of smearing becomes visible and shows the degree of resolution that can be expected form our seismic refraction measurements. Furthermore, it can be stated that crevasses have quite a large impact on the direction of the raypaths and on the shallow velocity structure of the ice. Despite the smearing it can be estimated from the extent of the low velocity zones, that the crevasses reach depths of around 10m. Zones with lower velocities are also present in profile 2 but are less emphasized than in profile 1. These patches could not be associated with crevasses. Since the position of these patches vary in the images calculated with half datasets, it is likely that they are artifacts of the inversion. There is possibly a low velocity layer at the surface, which is smeared along the raypaths due to the smoothing constraints in the inversion.

Figure 3.9: Velocity model of profile 1 and aerial picture showing the crevasses crossing the profile. The crevasses can be clearly associated with low velocity zones in the model. (Red dots: shot positions; black dots: geophones)

3.2. REFRACTION TOMOGRAPHY PROFILES

35

The velocity gradient In many studies in the ablation area of glaciers (e.g. Thyssen & Ahmad (1969)) a velocity gradient in the shallow subsurface down to some tens of meters with velocities as low as 3000m/s at the surface were found. Roethlisberger (1972) proposes the following reasons: Decreasing density due to expansion of the bubbles as the pressure decreases at the surface, internal melting because of solar radiation as well as fractures (crevasses) and small fissures in the ice. Thus this gradient layer can be considered as a complicated zone of varying water and air content and fracture density in the ice. The velocity models above also show that there is a velocity decrease close to the surface, though it is smaller parallel to crevasses than perpendicular. To quantify the velocity increase with depth, a 1D-model was extracted from the tomographic models (Figure 3.10). The aim is to obtain a velocity estimation that can be used for a more accurate location of icequakes and to see how sensitive the velocity is on the direction of crevasses. The 1D-models were extracted by averaging the velocities of all blocks at a certain depth from the surface. The errorbars (thin straight lines) were calculated from the standard deviations (Figures 3.6E and 3.7E)in the same manner as the mean velocities at each depth. The dash-dotted lines are a measure for the variation of the velocities at a certain depth (standard deviation of the velocity values that were averaged at each depth).

Figure 3.10: 1D-models extracted form the 2D-models form the inversion. The dashdotted lines show the variability of the velocities at a certain depth (standart deviation of all velocity at each depth). The thin straight lines are the error estimation from the Monte Carlo simulation (standart deviation from the 50 different models). In profile 1 the velocity increases from 3100 ± 100m/s at the surface to 3690 ± 25m/s at 20m depth. Profile 2 has a velocity of 3420 ± 80m/s at the surface and 3800 ± 25m/s at 22m depth. The rather high values at 20m and 22m depth leads to the assumption that in these depths the maximum velocities of ice are reached (higher values were only measure for single crystals, Section 3.1.1). The effect of the crevasses in profile 1 can be seen in the low velocities and the high variability of velocities close to the surface (dash-dotted line in Figure 3.10). The variability of the velocity in both profiles is also presented in Figure 3.11. All traces were sorted by source-receiver distances of 2.5±1m, 7.5±1m, 12.5±1m, etc. . For each source-receiver class the average and the standart deviation of the traveltimes were calculated. The standart deviation is measure of the variability of the seismic velocity. Figure 3.11 shows, that in profile 2 the variability is lower and nearly constant for all distances, whereas in profile 1 the variability is higher especially at intermediate distances. This can again be interpreted as an effect of the

36

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

crevasses.

Figure 3.11: The standart deviation of the traveltimes at certain source-receiver distances. Blue circles: profile 1, red crosses: profile 2. These values show the variability of the seismic velocity along the profiles.

Anisotropy The velocities of both profiles are compared in Figure 3.12 and show the largest differences at the surface. For clarity only the errorbars, not the variability is included in the plot, which makes no difference considering anisotropy. Profile 1 crossing the crevasses has lower velocities down to 20m, although between 5 and 10m the errorbars mutually overlap the mean velocities and the differences cannot be considered as significant. The relative differences of the velocities (i.e. ∆V /Vprof ile2 ) are 9.4% at the surface and 2.6 % in 20m depth. Since the high differences at the surface cannot be explained by crystal orientation (maximal 3.9 %), nor with layering (maximal 6.2 % in Dewart (1970)), it can be assumed, that the effect of crevassing dominates all possible anisotropy effects. Still these two effects must be present to a certain degree. However, they cannot be quantified since data about the strength and direction of preferred crystal orientation and the strength of layering are not available. Only the orientation of layering is known to dip around 70◦ in the direction 220◦ from north (also visible on the arial picture in Figure 3.3). Thus the strike has an angle of around 35◦ to the strike of profile 1 and 70◦ to profile 2. In layered ice the velocity is expected to be highest parallel to the layering. In our situation the layering should have an increasing effect on the velocities in profile 1 and a lowering effect on profile 2. This is not in agreement to what has been found. The crevassing effect must be superposing this effect and is probably even larger than described above. To clearly seperate the effect of crevasses from the others a similar study has to be done in an area without influence of crevasses.

3.3. BOREHOLE SEISMICS

37

Figure 3.12: Comparison of the 1D-models of profile 1 (blue) and profile 2 (red).

3.3

Borehole seismics

Borehole seismic experiments were conducted in order to obtain a better understanding of the velocity structures gained from the surface seismic survey for the uppermost 20m of the glacier ice. There are two aims of these investigations: First, obtaining information about the velocities beneath the low velocity layer at the surface (i.e. a 1D-model) and second, to investigate the time evolution of the velocity structure during the field campaign caused by the filling of the lake. For the second part the borehole seismic experiment was carried out several times, which makes it possible to compare the seismic velocity on different days.

3.3.1

Data acquisition

The experiment was performed eight times on the dates of 30/05/06, 07/06/06, 10/06/06, 15/06/06, 22/06/06, 28/06/06, 07/07/06 and 25/07/06 (yellow marks in Figure 3.3). A 120m deep borehole located at a few tens of meters from the maximal extension of the Gorner lake was drilled using the hot water drilling technique. A geophone chain, that consisted of four 30Hz three component geophones every 4m, was lowered into the borehole. To get traveltimes along the entire depth of the borehole the geophone chain was first lowered to the bottom of the borehole and then successively lifted 2 m and 16 m until the top of the borehole was reached. In this way traveltime measurements were obtained for every 2 m. (On the first day, 30/05/06, the chain was lifted by 16 m each time, so there are measurements only for every 4 m on this day.) A pipegun was used as source and was triggered in shot holes that were the same for all days. On the dates 07/06/06, 10/06/06 and 15/06/06, two shots were triggered for each depth of the geophone chain in a 30-50 cm deep shothole at 10m distance from the borehole. On the last four days (22/06/06, 28/06/06, 07/07/06 and 25/07/06) two shots each were fired in holes at 10 m and 20 m away from the borehole, which gives four traveltime measurement for each depth and a better ray coverage (crossing rays). The sampling frequency was 8000 Hz (only 4000 Hz on 30/05/06), the recorded time window 2 s and the trigger delay 1s (0.5 s on 25/07/06). In the following calculations the data of the first day is not be considered due to the low quality of the

38

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

dataset.

3.3.2

Coordinates

In contrast to the refraction seismic survey at the surface the geometry was more difficult to determine accurately. The top position of the borehole and the shot positions were measured with GPS only on the first day. The geophone positions within the borehole were determined via the length of the chain cable. Since the geophone chain reached the bottom of the borehole each time (no freezing or other changes within the borehole) and the cable length was nearly constant it can be assumed that the position of the bottom of the borehole was nearly constant (relative to the glacier ice). Markings on the geophone cable every 2 m were used to determine the depth of the geophones. This allows to calculate the coordinates of the geophones at least on the first day precisely. Since on other days again the geophone positions were chosen relative to the bottom of the borehole, the geophone coordinates of the first day can be taken also for the other days. This is sufficiently accurate under the assumption that the ice flow has only a negligible effect on the distance between the shot and geophone positions: Near the surface the horizontal ice flow does not move the shot and geophone position relatively to each other and in depth the relative horizontal motion due to ice deformation has a minor effect on the traveltimes, since the vertical distance is much bigger. In contrast, the change of the vertical coordinates of shot positions has a measurable effect on the raypath (and thus on the velocity) and have to be taken into account. Reasons for a vertical change of the shot positions are the following: • Ablation during the time period of the seismic experiments. • Vertical strain (elongation) of the ice due the increasing horizontal pressure of the lake water caused by the filling of the lake. • Uplift of the glacier ice due to the increasing water pressure (buoyancy forces). It can be assumed that the last reason influences both borehole bottom and shot positions in the same way and does not effect the relative positions. However, the other two reasons can separate shot and receiver positions (vertical strain) or bring them closer to each other (ablation). To estimate these effects the GPS data series recorded at a nearby stake were used (Figure 3.3). This stake was drilled to a certain depth into the glacier ice. A GPS antenna on its top measured the exact position several times per hour and thus also gives the change in vertical position of glacier surface. At the same time the length change of the stake portion above the glacier surface gives the ablation in this area. This was measured once or twice a week. This data is shown in Figure 3.13. The change of the stake foot position was calculated by subtracting the stake length from the coordinates on the top of the stake. The ablation is equivalent to the change of the stake length. As Figure 3.13 shows, at all times the ablation was bigger than the change of the stake foot position by a maximum of 50 cm. This difference includes the vertical strain and the buoyancy and is opposite to the ablation. Unfortunately, the size of the vertical strain and buoyancy effects cannot be distinguished, but they have a similar order of magnitude. The ablation data were used to correct the vertical coordinates of the shot position, while the effects of vertical strain is considered to be negligible.

3.3. BOREHOLE SEISMICS

39

Figure 3.13: Vertical change of the glacier surface relative to the its position on the 25/04/2006.

Figure 3.14: Examples of seismograms recorded at the same depth on different days.

40

3.3.3

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Seismograms

Prior to computing the seismic velocity from the borehole data, the signals recorded at each depth but on different days were compared in order to assess the quality of the data and to observe systematic changes in the arrival times. A typical example is shown in Figure 3.14, where the shot point was 10m away form the borehole. The seismograms have in general not an impulsive shape. The rather ringing character could have been caused by the insufficient coupling of the geophones to the borehole walls. The signals on the same days are usually quite coherent. This is not the case if the different days are compared to each other, since the borehole as well as the shotholes might have altered slightly between the different days. Apparently, a problem occurred on the days 28/06/06 and 07/07/06: At some depths the signal arrived about 1ms earlier for one shot than for the other (which one is earlier is completely random). This time difference cannot be explained with the change of the geophone position or shot position, since this would require a change of shot-receiver-distance of at least 3 m. Since no reasonable explanation was found, it is possible that a defect in the trigger electronics caused the problem. These very uncertain arrival times have not been included in later investigations. Furthermore, the signals systematically arrived up to 1ms earlier in the last days. This can be partly explained by shortening of the shot-receiver-distance due to ablation. Further systematic variations concerning waveforms and arrival times that can be explained by changes of the pore fluid pressure are hardly visible.

3.3.4

Seismic velocity estimation

To find an appropriate starting model for the 1D-tomography a simple 1D-model, which averages the velocity along a straight raypath between source and receiver, was calculated. This means that the velocities were calculated simply by dividing the source-receiver-distance through the traveltime. For this the filtered data (4th order butterworth filter with corner frequencies 200 and 1200Hz) was picked using "picktool.m". The calculated velocities were averaged within layers of 20m size (Figure 3.15). Two systematic errors were introduced in this way: First, the estimated velocities may be slightly too slow, since the real raypaths are curved and thus longer than straight paths. Second, in the lower layers the velocities are underestimated due to the lower velocities of the layers lying above. This effect decreases with depth. Though this calculation does not give the real velocities at depth, for the uppermost layer we obtain a good estimate of the mean velocity in this layer. On all days one similarity in the structures in the depth was visible: almost all vertical profiles showed a tendency to decreasing velocities at lowest depths. This cannot be just due to the effect of the low velocity layer at the surface on the layers below, since this effect is weakest in the lowest layers. This information, as well as the velocities in the surface layer are used for the inversion and its interpretation in the next sections.

3.3.5

Tomography

For the 1D-tomography the code "inv2d" was used, which also includes a version adapted for 1D-models. In this 1D mode of the code a detailed layer model with different velocities in each layer can be used as starting model. Only a damping factor d is used to regularize the inversion. Different starting models and dampings were tested to find the best model. Though for all days the data is overdetermined for a layer thickness of at least 4m, the resulting models were unrealistic with alternating layers of very high and very low velocities. Such models usually have a rms-value below the precision of the traveltimes, which means that the inversion tried to fit picking errors. The reason for this is probably the non-optimal geometry of the experiments

3.3. BOREHOLE SEISMICS

41

Figure 3.15: 1D-velocity estimations assuming straight raypaths and constant velocities between each source-receiver-pair. The red dots are the velocities calculated for each sourcereceiver-pair. These velocities were averaged within layers of 20m thickness (black line). In the lower right corner all velocities are summarized in one plot. On the last day scatter of the dots is quite high due to the bad data quality.

42

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

(i.e. rays penetrating the medium from only one direction is not optimal for tomography). By increasing the damping factor this problem can be solved but then the end models depend strongly on the starting models. Thus as much as possible a priori information has to be included into the starting model. A layering model consisting of six layers of 20m each was chosen in order to make the inverse problem highly overdetermined. For the uppermost layer the velocities estimated in Figure 3.15 were used, for the lower layers a constant velocity of 3700m/s. In order to investigate the reliability of the tomography models, the inversion was run for 40 different starting models. Therefore, a normally distributed random value with a standard deviation of 300m/s was added to the starting model, consisting of a surface layer with a velocity of 3500 m/s and 3700 m/s in the other layers. From the 40 different resulting end models only those with an rms-value below 0.4 ms (picking precision) are presented in Figure 3.16. The black line is the average of all models. The features of these mean models can be considered independent of the starting model (Figure 3.16). Again, the models show a decrease in the velocity with depth. The large variations of the velocities show, that the inverse problem has a broad minimum of the rms-value instead of a well defined single minimum.

Figure 3.16: End models of each 80 different random starting models (red) for 07/06/06 and 28/06/06. The black lines are the means of all end models. All end models presented here had an rms-value below 0.4 ms.

3.3.6

Error estimation

In this inverse problem there are two sources for errors: uncertainties in the picked traveltimes and uncertainties in the source and receiver coordinates. The effect of both kind of errors were examined separately. For the traveltime uncertainties a Monte Carlo simulation was run accordingly to the 2D-tomography models. Normally distributed random values (with standart deviation of 0.3 ms for the first six days, 0.4 ms for 07/07/06 and 0.6 ms for 25/07/06) were added to the traveltimes and for these newly generated datasets the inversion was run. Higher values for the standart deviation were chosen because of the higher rms-values at these days (Table 3.3). Repeating this 50 times gives 50 models of which the mean and the standard deviation was computed. For the geometry errors a more rough error estimation was done. The following assumptions were made: the depth of the shothole can be determined with a precision of ±0.1 m, the depth of the geophones with a precision of ±0.3 m and for the vertical coordinate of the

3.3. BOREHOLE SEISMICS

43

glacier surface an error of ±0.2 m due to uncertainties in the estimation of the ablation and vertical strain of the ice was estimated. These uncertainty estimates were used to calculate the error for the source-receiver-distances. To the vertical coordinates of the sources and the receivers a normally distributed random value according to the estimated errors was added. With these altered coordinates the source-receiver-distance was computed. This is again repeated 50 times and of all 50 results the mean and standard deviation is calculated. Dividing the standard deviation through the mean distance gives the relative error. According to the laws of error propagation (∆v/v = ∆s/s) the relative error is the same for the velocity if we assume that as a first approximation the velocity is linearly dependent on the source-receiver distance s. The errors for each layer given in Figure 3.17 are the errors that can be expected for a mean velocity of 3700 m/s. These errors are added to the errors resulting from the Monte Carlo simulation.

Figure 3.17: Errors resulting form the uncertainties of the source and receiver coordinates for a mean velocity of 3700 m/s.

3.3.7

Results

Figure 3.18 shows the velocity models considered to fit the measured data best. The rms-values are shown in Table 3.3. The high rms-value of the last day is due to the bad data quality, which resulted from insufficient coupling to the borehole walls. This made it difficult to pick the arrival times and results in a large scatter of the data (see also Figure 3.15). All models show similar patterns to some degree: a low velocity layer at the top, the maximal velocity just beneath and below, the velocity again decreases slightly. This decrease is quite small and not clearly present on all models, but the fact that almost all models (except 07/06/06 and 10/06/06) exhibit this pattern to some degree makes it reliable. The surface layer shows high temporal variations. The errors reach from 30m/s at the bottom of the borehole to 100m/s near the surface.

3.3.8

Interpretation

Temporal variations near the surface To examine the variations of the velocities between all days the mean velocities of three layers, 0-20m, 20-60m and 60-120m are presented as time series in Figure 3.19. Although the error bar for the 0-20m layer is the largest, significant temporal variations in the velocities are present. For

44

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Figure 3.18: Results from the 1D-tomography of the borehole data. Day 07/06/06 10/06/06 15/06/06 22/06/06 28/06/06 07/07/06 25/07/06

rms-value [ms] 0.15 0.16 0.18 0.22 0.18 0.42 0.77

Table 3.3: The rms-values of the inversions of all days, representing the discrepancy between the measured data and the fitted data.

3.4. CONCLUSION

45

the other layers the velocities remain nearly constant during the time period. The argument that the variations occur only due to the different starting velocities in the inversion for the different days does not hold, since these starting velocities were taken from the velocity estimations in Figure 3.15, where the same variations also occur. The most obvious change in the vicinity of the borehole that could influence the conditions of the ice, is the rising lake level of the nearby Gornersee which was less than 50m away from the borehole at its maximum level. The lake level measured in units of meter of water column with a pressure sensor inserted into the middle of the lake is shown in Figure 3.20. The velocities and the lake level seem to correlate to some degree. When the lake started to fill the velocity in the uppermost 20m is at around 3350±70 m /s and increases to 3640±55 m/s at the time when the lake started to drain and dropped again to 3450±60 m/s. The only outlier is the drop of the velocity on 28/06/06 (Julian day 179). These changes can be explained with water intruding into previously air-filled spaces, such as thin open fractures and open bubbles, as the lake level rises. The water replacing air in the ice would then have an opposite effect then water replacing ice (due to melting). First, because water has a higher seismic velocity than air and secondly waves that had to "surround" air-filled cracks (scattering) can penetrate these cracks if they are filled with water. The fact that the same variations do not occur at greater depth implies that this effect is limited to the crevassed layer, which also supports the assumption of water-filled cracks. Since there are no empirical laws available describing this effect, it is difficult to quantify this. Seismic velocity at depth Although the resolution of these 1D-models is not very high, it is reasonable to assume, that the velocity gradient at the surface changes in a depth of around 20 m where the velocity reaches its maximal value. This is supported by the 1D-models extracted from 2D-tomography models, where at this depth a velocity is reached that is unlikely to be exceeded at greater depths. Averaging all velocities at depths below 20 m gives a value of 3754±55 m/s, which is slightly higher than proposed by Aschwanden (1992) (3707±30 m/s) but still lies within the boundaries of uncertainty. This is in agreement with the seismic investigations of 2005 in the nearby area, where velocities in the range from 3720 m/s to 3790 m/s, were found for these depths (Gischig (2005)). As already mentioned all models in Figure 3.18 show a slight decrease with depth. Figure 3.19 also shows that the velocities in 60-120 m might be slightly lower than in 20-60 m though the difference is not significant. However, the fact that this can be stated on all days supports this observation. The mean values are 3780±96 m/s in 20-40 m and 3742±30 m/s in 100-120 m depth. An explanation for this is a slightly higher pore water content in greater depth: a decrease of 38 m/s corresponds with 6.3% increase of pore water content (see section 3.1.1). Since the effect of the temperatures on the velocity are only very small in this area (minimal temperatures measured near the survey area: -0.2 ◦ C), a temperature effect cannot stand alone as an explanation. If the ice is assumed to be slightly below the pressure-melting-point, it is also possible that the transition from cold to temperate ice in a certain depth of the borehole can cause this decrease via changes of water content.

3.4

Conclusion

The active seismic investigations used to determine the seismic P-wave velocity of the ice allows us to detect a roughly 20 m thick layer with a rather complicated velocity structure (e.g. lateral variations due to crevasses). In this layer the velocity increases from values between 3100±100

46

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Figure 3.19: Temporal development of the velocities in depth ranges of 0-20m (black), 20-60m (red) and 60-120m (blue)

Figure 3.20: The temporal change of the lake level in meters above sealevel (blue line). The red lines indicate the days when the borehole surveys were carried out. The lake level rises up into the 20m thick layer near the surface.

3.4. CONCLUSION

47

m/s and 3420±80 m/s at the surface to values between 3690±25 m/s and 3800±25 m/s at the depth of around 20 m. Mostly crevasses, but also internal melting due to solar radiation and increasing porosity have a strong influence on this part of the ice. Crevasses have also an anisotropy effect, which is emphasized strongest at the surface where velocities of 3100 m/s perpendicular and 3420 m/s parallel to the crevasses were measured. The effect of the crevasses at the surface is higher than crystal orientation and layering effects. To distinguish between these effects a similar survey in a non-crevassed area has to be carried out. From Figure 3.6A) it can be assumed, that the crevasses reach down to around 20m. Borehole measurements showed that near the Gorner lake the seismic velocity in the surface layer is influenced by the increasing water level. Water penetrating into small fractures and caverns is assumed to have an increasing effect of up to 300 m/s on the velocity. The shapes of the 1D-profiles lead to the assumption that the velocity can decrease again by around 40 m/s from the maximal value determined just beneath the surface layer down to the values found at the bottom of the borehole. This could be due to increasing pore water content or the transition from cold to temperate ice, but the high error bars make certain conclusions impossible. Since refraction tomography cannot image negative velocity gradients (Lanz et al. (1998)), it is impossible to reach depth below 20 m with refraction tomography with longer geophone spreads, if this velocity decrease with depth exists. More detailed information about the seismic velocity structure can be used to determine physical properties in greater depth, such as the pore water content. The geometry in our survey is probably not optimal to provide more accurate velocity estimations. A highly accurate crosshole survey would increase the precision of the velocity estimations and possible temporal variations can be monitored more reliably.

48

CHAPTER 3. THE SEISMIC VELOCITY OF GLACIER ICE

Chapter 4

Seismic attenuation in glacier ice 4.1

Introduction

In the idealized case of a perfectly elastic medium the total energy of the seismic particle motion (the sum of kinetic energy and potential energy) is constant in the total volume of the medium. In a more realistic medium, there is always some energy dissipated by the medium due to various inelastic processes (e.g. changes in crystal defects, motion on grain boundaries) referred to as "internal friction". This results in attenuation of the amplitude of the seismic signal. In earthquake seismology the seismic attenuation is an important property, since the strength of a seismic source can be calculated accurately only if the amount of energy dissipated by the medium along the raypath is known. In order to characterize the deep icequakes recorded on Gornergletscher (Chapter 6) and to calculate synthetic seismograms (Chapter 5) the seismic attenuation has to be calculated first. This is the aim of the following sections, whereas only the Q-factor of the P-waves is focused on.

4.2 4.2.1

Theory (from Aki & Richards (1980)) Definitions

The most general parameter describing the phenomenon of seismic attenuation is the quality factor Q, which is defined as 1 ∆E =− (4.1) Q 2πE where E is the peak strain energy stored in a certain volume and ∆E the energy loss per cycle due to inelastic attenuation. Since the energy of a seismic signal is related to its amplitude like E ∝ A2 and dE ∝ 2AdA, the Q-factor can also be expressed as: ∆A 1 =− Q πA

(4.2)

In the case of propagating waves we are interested in the attenuation of the amplitude per distance dx. The attenuation ∆A of the amplitude per wavelength λ is then: ∆A = −λ 49

dA dx

(4.3)

50

CHAPTER 4. SEISMIC ATTENUATION IN GLACIER ICE

Using the angular frequency ω and the phase velocity c to express λ = 2πc/ω we obtain dA ω =− A dx 2cQ

(4.4)

This differential equation has the solution: ωx − 2cQ

A(x) = A0 e

= A0 e−α(ω)x

(4.5)

with α(ω) = ω/2cQ. The same result can be obtained if we insert a complex wave number k = kr + iα into the solution of the wave equation: A(x, t) = A0 ei(kx−ωt) = A0 e−αx ei(kr x−ωt)

(4.6)

If the attenuation is linear (i.e. an attenuated pulse can be built from its attenuated Fourier components), then the Q-factor can be measured easily by comparing the Fourier components of the original pulse with the components of the pulse having traveled a certain distance.

4.2.2

The problem of causality

Assuming the attenuation is linear and Q and c are both independent of ω, the impulse response of an attenuated pulse wave δ(t − x/c) can be derived analytically. The Fourier transformation of the delta pulse is Z∞ ωx x x F(δ(t − )) = δ(t − )eiωt dt = ei c (4.7) c c −∞

Applying some attenuation to this pulse in the frequency domain gives the frequency response of an attenuated pulse: ωx p(x, ω) = e−αx ei c (4.8) An inverse fourier transformation gives the desired impulse response of the attenuated pulse wave 1 p(x, t) = 2π

Z∞

e−αx ei

ωx c



(4.9)

−∞

The result of this integral is a pulse symmetric about the time t = x/ce , where ce is the phase velocity in a perfectly elastic medium (see Figure 4.1). The pulse is no longer a δ-impulse but a broadened peak. This can be regarded as a general effect of attenuation, which in that way acts similar to a low-pass filter. The major problem of this result is that the onset of this broadened signal arrives earlier at a certain distance than a signal traveling in a perfectly elastic medium. This unphysical behavior violates the elementary assumption of causality. Furthermore, the symmetric shape and the too long rise time (definition below) of the attenuated pulse contradicts the behavior of measured signals. By allowing for frequency dependence of c = c(ω) (dispersion) and Q = Q(ω) more physical results can be obtained. Thus the complex wave number k = kr + iα can be written in the form k=

ω + iα(ω) c(ω)

(4.10)

4.2. THEORY (FROM ?)

51

If c is frequency dependent and the pulse is required to be causal, then it can be shown that ω ω = + H(α(ω)). c(ω) c∞

(4.11)

Here c∞ is the phase velocity c(ω → ∞) and H(α(ω)) is the Hilbert transform of α. Since ω/c(ω) = 2Qα(ω) equation 4.11 becomes 2Qα(ω) =

ω + H(α(ω)) c∞

(4.12)

This equation has no solution if Q is constant, which proves, that also Q must be frequency dependent. Azimi et al. (1968) proposed the following relation: α(ω) =

α0 ω 1 + α1 ω

(4.13)

which gives 1 1 2α0 ln(α1 ω) = − c(ω) c∞ π(1 − α12 ω 2 )

(4.14)

If α1 ω