Application of radar and seismic methods for the investigation of temperate glaciers

Journal of Applied Geophysics 57 (2005) 193 – 211 www.elsevier.com/locate/jappgeo

Application of radar and seismic methods for the investigation of temperate glaciers Francisco J. Navarroa,T, Yuri Ya. Macheretb, Beatriz Benjumeac a

Departamento de Matema´tica Aplicada, ETSI de Telecomunicacio´n, Universidad Polite´cnica de Madrid, Ciudad Universitaria, 28040 Madrid, Spain b Department of Glaciology, Institute of Geography, Russian Academy of Sciences, 29, Staromonetny, Moscow 109017, Russia c Departamento de Geodina`mica i Geofı´sica, Facultat de Geologia, Universitat de Barcelona, Martı´ i Franque`s, s/n, 08028 Barcelona, Spain Received 9 February 2004; accepted 23 November 2004

Abstract The capabilities of seismic and radar methods for the study of ice sheets have been analysed by other authors in the past. The joint use of both techniques has allowed the comparison of information, such as ice thickness, retrieved from both sources. Though these methods, specially the radar sounding, have also been widely used for the study of polythermal and temperate glaciers, the literature lacks joint analysis of their use for the study of temperate glaciers, where physical processes absent in the cold ice masses come into play. We have used seismic and radar methods collected at Johnsons Glacier, a temperate ice mass located in Livingston Island (Antarctica), to show the glaciological information that can be retrieved from such data. The aspects considered include the determination of ice thickness, the retrieval of information concerning the internal structure of the glacier (distinction between accumulation and ablation zones, determination of the depth of the firn–ice transition, detection of buried crevasses), the estimation of physical parameters such as seismic and radio wave velocities, water content in temperate ice and firn density, and the use of radar and seismic data to infer the presence of water channels at the ice–bed interface and to determine the nature of subglacial sediments. D 2005 Elsevier B.V. All rights reserved. Keywords: Seismic reflection; Seismic refraction; Ground penetrating radar; Glacier; Temperate

1. Introduction The propagation of elastic and electromagnetic waves through a given medium depends on very T Corresponding author. Tel.: +34 913367284; fax: +34 913367289. E-mail address: [email protected] (F.J. Navarro). 0926-9851/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2004.11.002

different properties of the material, and thus seismic and radar methods can be jointly exploited to retrieve complementary structural information. The seismic methods have a widespread use for geophysical exploration. The use of radar techniques is also becoming generalised for particular geophysical applications, the main constraint being the limited penetration of the electromagnetic waves through

194

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

most subsurface materials. The shorter wavelengths employed in the radar techniques, however, imply a higher resolution. For this reason, the radar methods are becoming increasingly popular as surface based geophysical methods to study shallow subsurface geology, say down to 40 m below the surface, where the accuracy of the present seismic methods is not always adequate (Daniels, 2004). When dealing with glaciers and ice sheets the penetration is no longer a constraint of the radar methods, as the low conductivity of glacier ice, and hence its low dielectric losses, make it an ideal geological material for the propagation of electromagnetic waves. Penetrations up to several hundred, even a few thousand, metres have been achieved in glaciers and ice sheets (Plewes and Hubbard, 2001). The interest, however, is not restricted to glaciology, as, in many cases, analogues can be found in other branches of geophysics, as will be discussed later. Focusing in glaciological applications, the radar and seismic methods are the most widely used methods for determining the ice thickness of glaciers and ice sheets. Seismic methods have been shown to be an effective tool to determine the basal topography and ice thickness in cold and temperate glaciers (Crary, 1963; Levato et al., 1999; Benjumea and Teixido´, 2001). Seismic surveys allow the investigation of the conditions at the ice–bedrock interface (Nolan and Echelmeyer, 1999), the nature of bedrock material and whether or not bed deformation is occurring (Smith, 1997); they also provide a means for obtaining density–depth profiles in polar ice sheets (Robin, 1958; Bentley, 1975). Radar airborne and ground-based sounding methods have also been used extensively for the retrieval of information concerning ice-thickness, bed conditions, internal layering and physical parameters of the ice in cold, temperate and polythermal glaciers (see, e.g., the reviews by Gogineni et al., 1998; Plewes and Hubbard, 2001). Low-frequency (say, below 15 MHz) monopulse ground-penetrating radars (GPR) have been found to be the most appropriate radar tool for the study of temperate glaciers, since they overcome the strong scattering caused by water inclusions at higher frequencies (Watts and England, 1976). The joint use of both techniques has usually been addressed for comparison of ice thickness in cold glaciers and polar ice sheets (e.g. Weber and

Andrieux, 1970; Drewry, 1975), as well as for calculating the depth-averaged densities of ice shelves (Doake, 1984). Another application has been the use of the density–depth curves in the firn zone, obtained from seismic measurements, to calculate the radio wave velocities needed to convert reflection times to depth and to detect buried crevasses and their direction (Retzlaff and Bentley, 1993; Clarke and Bentley, 1994). In this paper, seismic and radar data collected along coincident profiles are used to analyse the suitability of both methods for the study of ice thickness and physical parameters, bedrock topography, internal structure and subglacial conditions of a temperate glacier in Antarctica (Johnsons glacier). The aim is to show the capabilities of both techniques for the study of temperate glaciers and to illustrate how they complement to each other in accomplishing this task. The study, however, does not intend to be exhaustive, in the sense of discussing the capabilities of both methods in general, but is limited to a case study for a particular glacier and constrained by the instruments used for the field work and the available field data sets. Neither a full glaciological discussion of Johnsons glacier is undertaken in this paper; rather, we have used an approach of discussion of methods and presentation of corresponding examples. The selection of Johnsons glacier was dictated by both logistic reasons (due to its proximity to the Spanish Antarctic station Juan Carlos I) and scientific aspects, as there is a quite large amount of glaciological data available and the glacier dynamics is reasonably well understood. Most of the methodologies discussed in this paper can be extrapolated to other fields of applied geophysics. Just to give a few examples, the technique of using back-scattered seismic surface waves to detect buried crevasses could be used for the detection of fractures in salt mines; a snow-filled crevasse could be considered as an analogue of a dike; the study of water content from radar and seismic velocities is just an application of two-component mixture formulae to the analysis of physical properties of porous media; the use of GPR for water content estimate can be applied to groundwater studies of sandy sedimentary environments of high electrical resistivity; and so on.

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

2. Geographical setting Johnsons Glacier is a tidewater glacier located in Livingston Island, South Shetland Islands, Antarctica (Fig. 1), at 62840VS, 60830VW. It is delimited by a local ice divide (with altitudes between 200 and 330 m a.s.l.) that defines it as a separate glacier basin within the Hurd Peninsula ice cap. The glacier covers a total area of about 5 km2 and terminates in a 50 m high ice cliff extending 500 m along the coast. The basement consists of sandstones (Johnsons Bay side) and contact metamorphic rocks (Smellie et al., 1995). Thermally, it is a temperate glacier, as revealed by temperature and density profiles measured at boreholes (Furdada et al., 1999; data from M. Pourchet and J.M. Casas reported in Ximenis, 2001). This means that the ice mass is at pressure melting point throughout, except at a thin upper layer (10–30 m) in the accumulation area, in contrast to cold glaciers, which are at temperatures below melting point throughout, and polythermal glaciers, made of an upper cold layer overlaying a temperate layer. The flow lines in the northern part are shorter and have larger slopes (e.g., ca. 108 for L4 in Fig. 1) than those in the southern part (e.g., ca. 68 for L3 in Fig. 1). The confluence of the northern and southern flows results in a folded and highly fractured terminal zone (Ximenis, 2001). Ice surface velocities increase downstream from the ice divide, reaching values about 40 m year
1 near the glacier terminus (Ximenis, 2001). Accumulation and ablation rates show a large spatial and temporal (yearly) variability, with maximum accumulation rates about 1 m we year
1 and maximum ablation rates up
4 m we year
1 measured during the last eight years (Ximenis, 2001; unpublished data from F. Navarro). The equilibrium line (which separates the zones of net accumulation and net ablation along the year) altitude is about 150 m in the northern area and 180–260 m in the southern area (Ximenis, 2001). The shallowest part of the accumulation zone is characterized by a complex pattern of alternate firn and regelation ice layers, associated to episodes of intense summer surface melting and subsequent percolation and refreezing. Thus, the firn–ice transition is not clearly determined, though it can be considered to occur, near the divides, at depths about 22 m, though sometimes as low as 10 m, according

195

to density data from ice cores (Furdada et al., 1999; data from M. Pourchet and J.M. Casas reported in Ximenis, 2001).

3. Data acquisition and processing Our data set consists of seismic and radar data collected during different austral summer campaigns, primarily intended for ice thickness measurements. The seismic data correspond to 1996–1997 and 1997– 1998 campaigns. They include both refraction and reflection profiles, with a total length of 4.5 km (Benjumea and Teixido´, 2001). Radar profiles with a total length of 22 km (4.5 of them coincident with those of the previous seismic lines) were performed during the 1999–2000 campaign and complemented by some 5 km collected during the 2000–2001 campaign. Both seismic and radar profiles along and transverse to the glacier flow lines are available. 3.1. Seismic data The seismic data set consists of both refraction (L1, L7, L8 and L9 in Fig. 1) and reflection (L3 and L4 in Fig. 1). L3 and L4 approximately follow glacier flow lines; most of L1 is transverse to flow directions; L7, L8 and L9 are smaller profiles near divide areas. Shot gun SISSY and low-energy explosives (noise makers) were used for seismic data acquisition. Both sources were shot at 1 m depth. Different geometries were used for the data acquisition of refraction and reflection seismic profiles. For the refraction profiles, 48 geophones were placed with a spacing of 5 m. Five shots were fired along each spread, located with a nearest offset of 50 m at both sides of the spread, two at each end and one placed in the midpoint. A single spread was used for L7, L8 and L9, while L1 profile is the result of 8 consecutive spreads, shown in Fig. 1. Data acquisition for the reflection profiles was performed using 24 channels. The selected shot interval was of 10 m and the nearest shot-geophone distance (offset) was 30 m. For both refraction and reflection profiles, single vertical 40 Hz geophone stations were deployed and a 48-channel digital Seismograph (BISON 4098) was used for acquiring the seismic data with a sample rate of 0.1 ms and a record length of 500 ms.

196

Johnsons

a

B.A.E.

Reina Sofía

L7 L3

A

L1 L1

ANTARCTIC PENINSULA

L9 3048000

c b

Refugio

Y UTM (m)

South Shetland Islands

Hurd Peninsula

L8

Johnsons Glacier

BH1

c

L3 637000

634000

X UTM (m) Fig. 1. (a) Location map of Livingston Island (South Shetland, Antarctica). (b) Location of Hurd Peninsula and Johnsons Glacier within Livingston Island. (c) Layout of the coincident seismic and radar profiles (labelled as Ln, n being a number; L1 and L3 are labelled at the beginning and end of the line) and radar-only profiles (unlabelled straight lines) on Johnsons Glacier. The solid dots in L1 mark the position of the boundaries between the different sections (refraction profiles) that make up this profile. The solid square dots labelled BH1 and BH2 correspond to shallow borehole locations. The light grey areas are ice-free zones and the dark grey represents the sea. The dotted line shows the approximate location of the ice divide. The ice surface topography (contour interval 10 m) is also shown in this map.

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

BH2

CHILE

Livingston Island

L4

3049500

ARGENTINA

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

3.2. Radar data

Fig. 2 shows raw and filtered data acquired for the profile L8 as an example of the quality of seismic data. Several reflections can be observed in both raw and filtered data: R1, interpreted as the reflection from the bottom of the glacier ice; R2, probably resulting from the presence of a thin layer of sediments or till underneath the ice (see more detailed information in the subglacial conditions section); and R1*, which corresponds to the P-SV wave from the ice–bed interface. R1* is present at large offsets. The processing of reflection profiles consisted of a conventional common mid point (CMP) method, where the main difficulty was the elimination of high amplitude surface waves. Refraction data were also used to obtain stacked sections. These are characterized by a low fold due to the refraction geometry used. Nevertheless, good signal-to-noise ratio allowed the retrieval of structural information. For all data sets, band-pass filtering, automatic gain control, spectral balancing, velocity analysis, stacking, and post-stack filtering where required, constituted the main stages of the processing scheme. As final step, approximate depths were calculated using the stacking velocities. 235 m 0

197

The radar data were acquired using a lowfrequency monopulse ice-penetrating radar VIRL-2, specially designed not only for temperate glaciers but also useful for cold/polythermal glaciers (Vasilenko et al., 2002). The VIRL-2 equipment consists of transmitter, receiver and digital recording system (DRS). The antennae are resistively loaded half-wave dipoles of 5.8 m length. The transmitter generates pulses of 25–30 ns, with a peak power of 1.5 kW, at a pulse repetition frequency of 20 kHz. The transmitting system has a centre frequency of 15 MHz; by centre frequency we mean the frequency of the peak of the power spectrum of the radiated pulse. The receiver has a logarithmic amplifier with 100 MHz bandwidth and 80 dB input dynamic range. Synchronisation between transmitter and receiver is accomplished by a dedicated radio channel. The DRS sampling interval is 5 ns and 4082 samples are recorded for each waveform. For profiling, transmitter and receiver were placed on separate sledges towed by a snowmobile. Recording was made, using a common-offset geometry, at 1 s time interval between records, keeping a constant

50 m 0

1 2

50

50

2

4

100

4 5

150

3

4

Time (ms)

Time (ms)

R1

3

R2

100

R1* 150 M 200

200

a

b

Fig. 2. Examples of seismic data acquired at L8 profile, with shots located at the centre of the spread (left side of figures) and at 50 m offset from the last geophone (right side of figures), as indicated by the asterisks in panel a. (a) Raw data showing the main arrivals: (1) direct wave; (2) refracted waves; (3) reflection from the ice–bed interface; (4) shear and surface waves; (5) back-scattered energy of shear and surface waves from ice crevasses. (b) Filtered data (Butterworth filter 90–125–300–400 Hz) showing: (R1) ice-bottom reflection; (R2) reflection from the subice material; (R1*) P-SV converted wave; (M) first multiple of R1.

198

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

speed of about 3 m s
1 or, alternatively, a record for each odometer wheel rotation, the latter implying a distance between adjacent records of 1.55 m. The DRS controlled the radar transmitter and logged navigation data from GPS and the odometer wheel in real time. In the 1999–2000 campaign, the transmitting– receiving antennae were arranged parallel to each other at a distance of 4.6 m and transverse to the profile direction, while in the 2000–2001 campaign they were arranged collinearly, following the profile direction, and at a distance of about 7 m between antenna centres. In addition to the different antennae layout, some changes were introduced in the radar equipment used during the second campaign (VIRL2A). These included changes in the transmitter electronics, resulting in a better input–output impedance coupling between transmitter and antenna, and the use of dipole antennae not fully loaded (loading South 0

500

coefficient 0.25) in order to obtain more energy efficiency. These changes are manifested in the radar sections shown in Fig. 3. The uppermost zone obscured in the radar records (except for very strong reflections) is about 40 m in panel (a), while only about 15 m in panel (b). This is due to both the larger distance between antennae and the better coupling transmitter–antenna. A more clear bottom reflection is observed in (b), particularly in its central part, because the diffraction hyperbolae are clearly enhanced in (a) and, consequently, the smaller amount of energy transmitted to the lower levels results in a reduced contrast of the ice–bed interface. The enhancement of the energy from the diffraction hyperbolae can be attributed to: (1) the directionality of the radar antennae with respect to the flow direction and the crevasse orientation (Nobes, 1999), as the antennae are collinear with the flow direction and cross the crevasses perpendicularly in (b), and (2) an increased

Distance (m) 1000

North 1500

0.5

40

1

80

1.5

120

2

160

a

2.5 0

Time (µs)

200 500

0

1000

1500 0

0.5

40

1

80

1.5

120

2

160

2.5

Approx. depth (m)

0

b

Approx. depth (m)

Time (µs)

0

200

Fig. 3. Sample of radar data collected, at the same location (profile L3), in the 1999–2000 (a) and the 2000–2001 (b) campaigns. L3 approximately follows a glacier flowline; in (a), the antennae were arranged parallel to each other and transverse to the profile direction; in (b), they were arranged collinearly and following the profile direction.

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

water content in the glacier, consistent with meteorological data available. In Nobes experiments, the internal reflectivity is significantly less and the basal ice-rock reflection is clearer for the axially orientated radar antennae, in agreement with our results. According to Nobes, the cause of the directional dependence of the radar response is not the preferential orientation of ice crystals but the preferential orientation of structural features within the ice relative to the polarization of the radar signal. The processing of radar data included DC correction, amplitude scaling, band-pass filtering, and migration where required. The conversion to depth was done using our best estimate for the radio wave velocity, as will be discussed later.

4. Physical parameters of ice The seismic and radio wave velocities in ice are fundamental magnitudes, in the sense that other important parameters, such as ice thickness, water content or density are derived from them. Thereby, the highest attention should be paid to their accurate measurement. 4.1. Radio wave velocity The radio wave velocities were calculated based on the travel time to diffractors observed in the radar records. Only symmetric diffraction hyperbolae were used, which correspond either to point or to horizontal linear diffractors crossed by the radar profile (Clarke and Bentley, 1994). To measure the shape of asymmetric hyperbolae, which typically correspond to inclined linear diffractors such as water-filled channels, we must detect and measure them on multiple profiles to correctly estimate the dip and strike of the linear body and the radio wave velocity (Macheret, 2000). The accuracy of the radio wave velocity estimate could be improved by using field methods such as the common mid point or the wide angle reflection (Davis and Annan, 1989). However, the application of these methods in our case was not possible because the radio wave synchronisation between transmitter and receiver limited the separation between antennae to a distance of about 10–15 m.

199

At Johnsons glacier, we counted on a data set of 42 clearly identified symmetric diffraction hyperbolae, corresponding to diffractor bodies between 28 and 162 m depth. The associated velocities ranged from 157 to 174 m As
1, with an average of 164.9F4.2 m As
1. These are usual values for temperate glaciers, except those above 168 m As
1, which can be attributed to the effect of the snow and firn layers. This is confirmed by the fact that, with just a couple of exceptions, such high velocities corresponded to relatively shallow diffractors in the accumulation area. 4.2. Seismic wave velocity Different approaches were employed to obtain seismic velocity information, mainly depending on the configuration used for data acquisition. For the eastern and western sectors of refraction profile L1, where first arrivals from the ice–bed contact are observed, we used the routine refraction delay-time method (Pakiser and Black, 1957). The curve fitting method (Hunter, 1971) was applied to obtain velocity–depth profiles for the uppermost 50 m of the refraction profiles (L1, L7, L8 and L9). Reflection travel times were used for obtaining root-mean-square velocity from surface to bed in refraction profiles lacking information of first arrival from the bed. Three different techniques were employed depending on the data quality and bed topography: the T 2
X 2 technique, the T
Dt method and the calculation of the semblance coefficient. These methods were used for L7 and L8 profiles. Source-generated noise (surface waves, refracted waves) precluded the detection of the reflection in the L9 profile, where the glacier bed is shallow. For estimating the seismic velocities between the glacier surface and the bed for the reflection profiles (L3 and L4), stacking velocities were used. At Johnsons glacier, the seismic velocity shows a spatial variability larger than that of radio wave velocity. Profiles L1, L3, L7 and L8 show velocities ranging from 3146 to 3718 m s
1, with estimated errors within 20–192 m s
1 (i.e., relative errors in velocity between 1% and 5%). L4 reflection profile is characterized by extremely low seismic velocities (2800–3000 m s
1), which could be attributed to the

200

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

combination of a thick snow–firn layer and presence of crevasses. 4.3. Water content in temperate ice The estimate of water content W in temperate ice from two independent data sources
seismic and radar wave velocities
was the subject of a previous paper by the authors (Benjumea et al., 2003). Therefore, we will limit here to present the main ideas and results. The water content in temperate ice can be estimated from the two-component dielectric mixture formulae by Looyenga (1965) for ice with spherical water inclusions  W  1=3 1=3 1=3 e1=3 ¼ ei þ ew
ei ð1Þ 100 and the relation Vr=c/e 1/2, where e, e i and e w are the relative dielectric permittivities of the mixture, solid ice and water, and Vr and c=300 m As
1 are the radio wave velocities in glacier ice and air, respectively. This relationship takes implicitly into account the density of the mixture, as W for temperate ice is taken as (1
h) 100=W, h=q/q i, q and q i being the densities of the mixture and of solid ice, respectively (Macheret et al., 1993), where the assumption is made that the ice is water saturated (keep in mind that Eq. (1) is a two-component mixture formula). When using Eq. (1) to estimate the water content of temperate ice, values VrN168 m As
1 should be discarded, as they are associated either to dry glacier ice or to dry or wet firn (Macheret et al., 1993). Discarding such velocities from our original Johnsons glacier velocity data set, 29 diffraction hyperbolae remained, giving velocities within 157–167 m As
1 and an average of 162.8F2.9 m As
1. The associated water contents, obtained from Eq. (1) using e i=3.19 and e w=86 at melting point (Macheret and Glazovsky, 2000), range from 0.2% to 2.3%, with an average of 1.1F0.6%. If we eliminate diffraction hyperbolae far from the seismic profiles, a data set of 12 hyperbolae remains, giving velocities within 157–166 m As
1 and an average of 162.3F3.0 m As
1, with associated water contents ranging from 0.4% to 2.3% and averaging 1.2F0.6%. In a similar way to that used for radar data, Benjumea et al. (2003) proposed to estimate the

water content in temperate ice from the seismic velocity data using the equation, due to Riznichenko (1949), describing the dependence between both variables for a two-component heterogeneous elastic medium:    100
W ð100
W Þqi 1þ Vs ¼ Vw 1 þ W W qw  
1=2 2 ð100
W Þqw Vw  1þ ; ð2Þ W qi Vi2 where q w and q i are the densities of water and solid ice at 0 8C, and V s and Vi are the seismic wave velocities in these media, respectively. Only velocities higher than 3500 m s
1 were considered for the water content estimation at Johnsons glacier. Lower velocities could be attributed to the presence of dry ice, a thick snow–firn layer or crevasses. With such restriction, the seismic velocities at Johnsons glacier profiles range from 3530 to 3720 m s
1, with estimated errors within 20–192 m s
1. The corresponding water content estimates obtained from Eq. (2), using Vi=3800 m s
1 (Kohnen, 1974), V w=1500 m s
1, q w=1000 kg m
3 and q i=917 kg m
3, range from 0.9% to 3.2%, with an average of 2.2F0.9%. The error estimates for W from radar and seismic velocity data were obtained applying error propagation to the two-degree polynomial curves approximating the set of points obtained from Eqs. (1) and (2), respectively, for W values from 1% to 10%, in steps of 1.0. Our estimates of water content are close to the values of 1–1.7% measured directly from ice cores in temperate glaciers (Raymond and Harrison, 1975) and also to those estimated from common mid point and borehole radar measurements in temperate glaciers (Macheret et al., 1993; Murray et al., 2000). Though our water content estimates obtained from radar and seismic data show a rather good agreement, the estimates derived from seismic data are generally higher than those computed from radar data. The reason is that seismic velocity increases with density, while radio wave velocity decreases with permittivity, which is, in turn, an increasing function of density. Thereby, the water content estimates obtained from seismic and radio wave velocities should be understood as upper and lower limits, respectively, to the actual values.

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

4.4. Firn and ice density Density versus depth in the accumulation area can be calculated from the seismic velocity–depth curves using the equation obtained by Robin (1958). Fig. 4

provides an example for the seismic refraction profiles L1, L7, L8 and L9. The seismic velocities were obtained using the curve fitting method (Hunter, 1971). Both seismic velocity and density as functions of depth are shown. L7 and L8 are located in the

Seismic velocity (m s-1)

Depth (m)

1500 0

2000

2500

3000

3500

Seismic velocity (m s-1) 4000

1500 0

10

10

20

20

30

30

L8

2000

2500

3000

3500

4000

L9 BH1 40

40 400

500

600

700

800

400

900

500

600

Density (kg m-3) 1500 0

2000

2500

3000

3500

700

800

900

Density (kg m-3)

Seismic velocity (m s-1)

Depth (m)

201

Seismic velocity (m s-1) 4000

1500

2000

2500

500

600

3000

3500

4000

0

10

10

20

20

30

30

L7 BH2

L1

40

40 400

500

600

700

800

Density (kg m-3)

900

400

700

800

900

Density (kg m-3)

Fig. 4. Seismic velocity–depth and density–depth curves obtained from curve fitting of the first arrivals of the refraction profiles L1, L7, L8 and L9. For L1 profile, the curve corresponds to its third sector starting from the west (see Fig. 1). The error bars are for computed density values. For L7 and L9, the densities measured at neighbouring boreholes are also shown (square dots).

202

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

accumulation zone, the sector of L1 corresponds to the ablation zone and L9 is in the transition between the accumulation and ablation areas. For L7 and L9, density values measured at neighbouring shallow boreholes (see their location in Fig. 1c) are also shown. Densities computed from seismic velocities show a rather good agreement with those measured from ice cores, though, for L7, the observed ones are systematically above the computed ones; the most significant discrepancies, however, correspond to regelation layers. Also notice in Fig. 1c that the borehole locations are not so close to the seismic lines (this is especially true for BH1-L9). Density versus depth curves allow the identification of the firn–ice transition, as will be discussed in a later section. Moreover, firn layers behave as low seismic velocity layers and taking them into account is fundamental for an accurate determination of ice thickness and water content from seismic data. As a sample, for profile L8 the correction for such low velocity layer with respect to the full column averaged seismic velocity implies a change in the water content estimate, using seismic velocity data, from 2.4% to 1.0%. Density–depth curves obtained from both seismic velocity data and ice cores have also been used by other authors as an indirect way to determine the radio wave velocity profile, as a function of density, in the firn layer of ice sheets (Retzlaff and Bentley, 1993). They use another Looyenga equation, in terms of density instead of water content and valid only for dry firn. Seismic and radar data have also been combined to estimate the depth–averaged density of ice shelves (Doake, 1984). Again, however, we should stress that these applications are restricted to cold ice masses.

5. Ice thickness and bedrock topography For the equipment used, the vertical resolution provided by radar data was about twice that of seismic data. Considering a seismic wave velocity in ice of 3550–3630 m s
1 and a 150 Hz reflection seismic signal, we get a wavelength of approximately 24 m. On the other hand, considering a radar wave velocity in temperate ice of 165 m As
1 and a dominant frequency for the reflection radar signal of 15 MHz, we get a wavelength of 11 m. Taking as vertical resolution a quarter of the wavelength, seismic and

radar data resolutions in ice are about 6 and 3 m, respectively. In spite of such relatively small difference between radar and seismic wavelengths, diffractions are far more common in radar sections than in seismic sections. This is due to the presence of water inclusions in temperate glaciers (interstitial, intrusions, water lenses), which have a strong effect upon the radio wave propagation in temperate ice and not so strong for seismic waves. Because of this, the seismic method is a useful complement of the radar method for studying the ice thickness and bedrock topography of temperate glaciers. As an example, Fig. 5 shows, for profile L3, the stacked seismic section (panel c) and the radar section (panel d). The seismic section depicts a clear bed reflection; diffractions from this interface can also be observed. A cluster of internal diffractions, as well as bed diffractions, can be observed in the central part of the radar profile. The internal diffractions are attributed to the presence of water inclusions, such as water lenses, or water channels crossed by the radar profiles. These are rather usual features of temperate glaciers such as Johnsons glacier, where intense summer melting occurs at the surface. Diffractions at the ice–bed interface are more abundant in radar profiles than in seismic sections. This could be due to the different wavelengths of seismic and radar waves as compared to the size of the diffractors. However, even if diffractors of intermediate size between both wavelengths were involved, seismic and radar records could show a different bed diffraction pattern, as seismic and radio waves show a different response (stronger contrast for the latter) to phenomena such as the presence of water channels at the ice–bed interface. The comparison of seismic and radar data shows that both methods provide valuable information concerning ice thickness and bedrock topography. Ice thickness values for Johnsons glacier obtained with both methods do agree in general (differences b10 m). In some areas, radar data suggest a bedrock deeper (differences N10 m) than that determined from seismic methods. This is most likely due to the effect of the firn layer, which behaves as a low velocity layer for the seismic waves while as a high velocity layer for the radio waves. Its effect on the travel times could not be thoroughly taken into account because of the

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

South

Distance (m)

0

500

203

North 1000

1500

2000

1500

2000

Time (ms)

0

20

a 0 50

b Distance (m) 0

500

1000

Time (ms)

40 80

50

120 160

100

Approx. depth (m)

0

0

200

c

150

Distance (m) 0

500

1000

1500

0.5

40

1

80

1.5

120

2

160

2.5

d A

Approx. depth (m)

0

0

Time (µs)

2000

200

B

Fig. 5. Results for survey profile L3: (a) 30 m common-offset profile and (b) common-receiver stacked section. The arrows mark the position of the crevasses along this profile, shown by first arrival delays in (a) and back scattered energy in (b). (c) Seismic stacked section and (d) radar section. Diffraction hyperbolae denoted A and B in zoomed areas of panel d are discussed in the text.

204

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

unavailability of information on thickness of the firn layer for some of the profiles in the accumulation zone. Because of both instrumental constraints and the geometry used for radar and seismic data acquisition, basement shallower than 30–35 m for seismic data could not be resolved in the present study; in the case of radar data, this zone extended to 40 m depth for the 1999–2000 campaign data, and about 15 m for the 2000–2001 campaign data, because of both improvements in the equipment electronics and in the antennae layout, as previously discussed.

6. Glacier structure The wavelengths of our seismic and radar signals do not allow to identify any layering within the glacier firn or ice. Nevertheless, useful information concerning glacier structure can be extracted from seismic and radar data, as discussed below. 6.1. Crevasse detection Both seismic and radar techniques are appropriate for recognizing buried crevasses. Seismic data is restricted to the identification of buried near-surface crevasses, and the feasibility of their recognition depends on their size and orientation. Radar data allows the estimation of the geometry and position of buried crevasses, as well as other englacial bodies, such as water inclusions (Macheret, 2000). Concerning seismic data, they can be used for the detection of near-vertical discontinuities within the ice. Back scattered shear and surface waves energy, manifested as events whose arrival times decrease with the source–receiver distance (e.g., 5 in Fig. 2a), indicates the presence of these discontinuities, interpreted as crevasses or closed-up crevasses. In order to graphically illustrate this effect and to locate these discontinuities along the continuous reflection profile L3, an F-K filter has been applied to remove the forward energy. Linear-Moveout (LMO) allows the alignment of the back-scattered surface waves to be done at the receiver domain. Stacking of common receiver data has been done to enhance this energy. The comparison of this common-receiver stacked section and the 30 m common-offset display (Fig. 5,

panels a and b) shows the good correlation between the position of these discontinuities (front of the backscattered energy) and the areas where time delay in the first arrival can be observed (pointed by arrows on panel a). In the radar records, the crevasses are usually manifested in the form of stacks of diffraction hyperbolae extending from near the surface to depths down to 40 m or even deeper. Very clear examples of such stacks can be found, e.g., in Fig. 3 of Moore et al. (1999). Other possible patterns are discussed in Clarke and Bentley (1994), who analyse a full set of hyperbola patterns and associated interpretations as crevasses, faults, sharp foldings and sagging snow bridges. Notice, however, that these studies are based on radar records from cold and polythermal glaciers, respectively. Interpretation is rather more difficult in the case of temperate glaciers, where hyperbolae arising from such kind of features mix with those arising from water-related features (water channels or water pockets), either full of water or void, and those associated to refreezing of percolating meltwater, such as ice lenses (though these water related features are also present in polythermal glaciers, their uppermost layer is made of cold ice, so allowing an easier interpretation of crevasses in radar records). Our radar section for profile L3 (Fig. 5d) provides an example of such difficulty, which is further complicated because the uppermost 15 m are not resolved in our radar records. In this radar section, large clusters of diffractions are present at many locations. Though most of the features interpreted as crevasses from panels a and b correspond to diffractions in panel d, some others do not have a clear counterpart in panel d. Fortunately, the information provided by the seismic records can help in the interpretation of some of the hyperbolae. For instance, Clarke and Bentley (1994) point out that freshly opened, straight-edged, smoothwalled crevasses with burial depths of less than 1 m show a single strong hyperbola at depth (bottom of crevasse) but almost no surface or wall signature. This could be the case for the faint hyperbolae labelled A and B on the zoomed areas of panel d, taking into account the depth of their apex at about 30–40 m and that they correspond to clear crevasse locations identified from panels a and b (the arrows showing their horizontal positions have been replicated on panel d for easy identification).

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

A large amount of hyperbolae from panel d do not correspond to crevasses identified from panels a and b. Almost all of them are likely associated to water inclusions (in particular, all the deep ones, but also many others), though we should remark that some crevasses that write a signature on the radar records could have no signature on the seismic records because of their size or orientation. The point to stress here is that the joint use of seismic and radar methods can help in the identification of buried crevasses, which is undoubtedly a highly interesting subject, as it has remarkable applications, such as the estimation of the stagnation time of ice streams from the depth of buried crevasses (Retzlaff and Bentley, 1993). Further comments on crevasse detection from diffraction hyperbolae, including the determination of the crossing angle from the shape of the hyperbola, can be found in Clarke and Bentley (1994) and Macheret (2000). Another way of determining the crossing angle is possible if parallel profiles crossing the crevasse are available. The procedure is illustrated in Fig. 6. The crossing angle a can be determined from a=tan
1 (d/D), where d is the distance between the parallel profiles and D is the horizontal distance between the apexes of the diffraction hyperbola in the radar sections (D is shown in the figure for hyperbola H1). According to this formula, the crossing angles for the crevasses corresponding to H1, H2 and H3 are 338, 398 and 448, respectively. 6.2. Distinction between accumulation and ablation zones Seismic data offer two ways to distinguish between accumulation and ablation zones. One is based on the analysis of the horizontal variation of the seismic refraction velocities. This can be measured using the differences in first-arrival travel time between adjacent records (Lawton, 1989). This procedure is adequate for the acquisition geometry of the reflection profiles (end-on shot configuration), since it does not require reciprocal shooting. An example is shown in Fig. 7c, which displays the variation of seismic velocity along profile L3. The high-velocity anomaly observed at the beginning of the profile corresponds to a zone where bedrock is close to the surface,

205

resulting in refraction arrivals from this unit. This anomaly is followed by a sector characterized by low seismic velocity, which is interpreted as the accumulation zone. In the last part of the profile, the velocity in ice varies between 3300 and 3800 m s
1. Values as low as 3100 m s
1 in this part have been interpreted as an indicator of a crevassed area. The average value for the near-surface ice velocity in the ablation area is 3525F198 m s
1. The second method is based on the analysis of the changes in the characteristics of the seismic surface waves, considering that their seismic signatures depend on the properties of the area of the glacier where they were acquired. In the accumulation area, surface waves are dispersive (i.e., their velocities depend on frequency) due to the presence of a velocity gradient near the surface (snow, firn), while data acquired on the ablation area is characterized by non-dispersive waves with velocities higher than those for the accumulation area. This can be observed in Fig. 7, where the shot gather corresponding to the accumulation area (panel a; from profile L7) shows a surface wave train that spreads out with increasing offset, very different of what is shown in the shot gather for the ablation zone (panel b; from profile L1). This fact is further illustrated in Fig. 6 from Benjumea and Teixido´ (2001), where plots of velocity versus frequency are presented. The possibility of distinguishing between accumulation and ablation zones from radar data is strongly dependent on the frequency of the radar used. High frequency radars (say, centre frequency higher than 100 MHz) offer a resolution that allows a detailed analysis of the structure of the shallowest zone of the glacier, allowing to identify the snow and firn layers, the transition between firn and ice, and so on. This, however, depends on the thermal character of the glacier. In cold and polythermal glaciers, the transition between firn and cold ice is easily detected (e.g., Figs. 3 and 5 of Moore et al., 1999; Kohler, 2002) as the firn shows a layered structure that sharply contrasts with the underlying cold ice, which shows almost no diffractions. In temperate glaciers, or the areas of polythermal glaciers where firn overlies warm ice, the identification is not so easy, because warm ice usually shows much scattering due to englacial water bodies, that can be also present

206

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

N

100 m

S

0

H1

H2

H3

Time (µs)

0.5

1

1.5

a 2

D 0

H1

H2

H3

Time (µs)

0.5

1

1.5

b 2 Fig. 6. Radar sections for two parallel profiles crossing several crevasses. Their corresponding diffraction hyperbolae are labelled as H1, H2, H3. (a) is a portion of L3 profile, while (b) is parallel to it at a distance of 10 m.

within the firn layers. In this case, the identification is made based on the layered structure of firn (e.g., Figs. 2 and 5 of Moore et al., 1999; Figs. 5, 6 and 10 of Arcone, 2002). For our Johnsons glacier radar data, however, two reasons prevent the distinction between the accumulation and ablation zones: first, the resolution is low because of the low frequency of the radar; second, the radar signal, as mentioned earlier, is obscured for the uppermost 40 m (1999–2000 data) or the uppermost 15 m (2000–2001 data), precluding the identification of the snow–firn layers.

6.3. Depth of firn–ice transition An evident application of the density–depth curves obtained from seismic velocity–depth curves in the accumulation zone, discussed earlier, is the identification of the depth of the firn–ice transition. A look at the density–depth curves in Fig. 4 illustrates it. The shape of the curve for L7 is consistent with a curve for the accumulation zone: density smoothly increases from values about 460 kg m
3 near the surface to densities typical of ice (830 kg m
3) at 20 m depth. L8 curve is similar, though shows a sharp increase in

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

0

60

60

Time (ms)

Time (ms)

235 m 0

120

180 4900

Accumulation zone Bedrock

4700

235 m

120

180

a

207

b Ablation zone

4500

Velocity (m s-1)

4300 4100 3900 3700 3500 3300 3100 2900 2700

c

2500 4685

4765

4854

4952

5031

5089

5152

X UTM (m) Fig. 7. Sample shot gathers corresponding to the accumulation (a) and ablation (b) areas, collected at profiles L7 and L1, respectively; (c) nearsurface seismic velocity variation along profile L3 obtained from differences in first-arrival travel times between adjacent records. Accumulation and ablation zones can be distinguished based on both the dispersion characteristics of the shot gathers and the velocity values.

density in the uppermost 2 m and a firn layer thinner (about 15 m thick) than that of L7. These depths for the firn–ice transition are consistent with those measured from ice cores at divide locations and shown in Fig. 4 for boreholes BH1 and BH2 overlapped with the density–depth curves for their closest seismic lines (L9 and L7, respectively). Another borehole at a neighbouring location outside Johnsons glacier also shows a firn–ice transition in the depth range 15–20 m. Density–depth profile for L1, on the contrary, is characteristic of ablation zone, showing just a thin (b1 m) snow cover immediately followed by ice. L9 profile is transitional between the

other two types, reaching the density typical of ice between 5 and 10 m depth after a very sharp increase in density from surface to 5 m depth. As mentioned earlier, our radar data do not allow to resolve the snow–firn structure, so preventing the detection of the depth of the firn–ice transition.

7. Subglacial conditions Seismic data allow the determination of the seismic velocity in the underlying bedrock. The routine refraction method (Pakiser and Black, 1957), based

208

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

on a layered model, allowed us to compute seismic velocities for bedrock in areas where the latter is close to the glacier surface. These values range from 4200 m s
1 to 4300 m s
1. In addition, the seismic method can provide information about the acoustic bed properties that are critical in the investigation of glacier dynamics. The use of wavelet energies of both ice–bed interface reflection and its first multiple allows the calculation of the seismic reflection coefficient at bed (Smith,

1997). Both primary and multiple reflections can only be observed in seismic profile L8 (Fig. 2a). Following Smith (1997), we determined the seismic reflection coefficient at the bed using an average trace resulting from four traces with an incidence angle less than 88. The absorption coefficient was calculated from the internal friction Q
1 obtained by Clee et al. (1969) for a temperate glacier. Assuming a frequency of 150 Hz and a velocity of 3630 m s
1, the absorption coefficient is 210
3. The resulting seismic reflection

Northwest

Southeast

Time (ms)

0

50

BR 100

CW X UTM=5890 Y UTM=7863

a X UTM=6190 Y UTM=7652

0

Time (µs)

0.5 1

1.5

ID ID BR

2

BD

Power/Sample(mV2/sample)

2.5

b IRP

0.1

0.001

BRP

1E-005

1E-007 280 m

c

Fig. 8. (a) Seismic stacked section for L8 profile. C.W.=seismic stacked P-SV converted waves, B.R.=bed reflections. (b) Radar section for L8 profile. I.D.=diffraction within the ice, B.R.=bed reflections, B.D.=bed diffractions. (c) Relative power level of Internal Reflection Power (IRP) and Bed Reflection Power (BRP) for L8 profile, calculated from radar data.

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

coefficient is F0.17F0.006. The sign of this coefficient is given by the polarity observed in the ice–bed reflection. Since the direct and reflected waves show the same polarity, we assume a positive polarity for the ice–bed reflection. Considering a seismic velocity in ice of 3630 m s
1 near the glacier bed and a density of 917 kg m
3, the acoustic impedance for ice Z i would be 3.33106 kg m
2 s
1, which would lead to an acoustic impedance of the bed material Z b of 4.5106 kg m
2 s
1. According to data from Smith (1997), Atre and Bentley (1993) and references therein, this acoustic impedance would correspond to a porosity of the subglacial material b0.30, that could be interpreted as lodged till or poorly lithified sedimentary rock (non-deforming sediments) underneath this section of the glacier. The presence of these subglacial sediments is confirmed by the stacked section for L8 (Fig. 8a), where several reflectors below the glacier bottom clearly appear. The earliest and strongest reflection (labelled BR) is interpreted as the reflection from the ice–subglacial material interface. This interface is also clearly defined in the radar section for L8 (Fig. 8b). The southern part of the profile shows a zone of near-surface multiple diffractions whereas the northern part is more transparent. Only two clear diffractions are observed in the latter part. Radar data can be used to study the space distribution of the bed reflection power (BRP), which can be related to changes in basal properties (Gades et al., 2000). From the radar data collected at profile L8, we have calculated the Internal Reflection Power (IRP—Fig. 8c), defined as one-half the sum-ofsquared amplitudes (after correcting for the effects of the logarithmic amplifier of the receiver) divided by the number of samples within a time window starting after the direct wave and ending before the bed reflection arrival times (Gades et al., 2000). This value is assumed to be an index of the power attenuated/ scattered in the ice. For the Bed Reflection Power (BRP) calculation, the time window is centred on the reflected pulse and has a length of 0.07 As, wide enough to include a full wavelength in ice. Fig. 8c shows the variation of BRP and IRP along the profile. Since the ice thickness variation is not high, the variations of BRP along the profile cannot be due to changes in geometrical spreading from trace to trace. The increment of IRP towards SE is related to the

209

strong scattering observed from surface down to approximately 1 As (Fig. 8b). Correspondingly, there is a decreasing trend in the BRP variation towards SE. However, the smaller scale variations in BRP cannot be explained by the IRP variations. The NW end is characterized by a high value of BRP followed by a significant drop (by 25 dB) in this value. The high BRP in the NW end is likely due to the presence of subglacial water in this part of the profile, which corresponds to a zone where the glacier bed suddenly flattens after a very steep slope.

8. Conclusions Radar and seismic methods have been extensively used for the study of ice thickness, internal structure, physical parameters of ice and bed conditions of glaciers and ice sheets. Radar profiling provides an easier and faster rate of data acquisition as compared to seismic profiling (recording from airplanes or snowmobiles versus manual deployment of shot points and seismic receivers), thus allowing a more detailed study of the ice masses. The radar method should thus be preferred. In temperate glaciers, however, the water content is high (especially during the melt season) and the radar sounding can be effectively complemented by seismic profiles in order to get a better understanding of the water content distribution within the ice and, thereby, an improved control on the estimates of ice thickness, especially in the accumulation zone. A precise knowledge of water content distribution is also important for studying the dynamics of temperate glaciers, because of the dependence of viscosity upon the water content (Duval, 1977). Additionally, we have shown that seismic methods complement to radar methods in the interpretation of field data. In particular: 1) Seismic data have been used for the detection of near-vertical discontinuities within the ice: back scattered shear and surface waves energy indicates the presence of these discontinuities, interpreted as crevasses or closed-up crevasses. These are manifested in the radar records as near-surface diffractions and, correspondingly, a lower level of bedrock reflections.

210

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211

2) The bed reflection power calculated from radar data can be used for inferring the presence of water at the glacier bed. Some of the radio wave diffractions on the ice–bed interface could also be attributed to the presence of water channels at this interface. 3) Density–depth curves retrieved from seismic velocity–depth profiles have been employed to determine the depth of the firn–ice transition. In our case, this interface could not be detected by radar methods because of the low frequency of the radar employed. 4) The near-surface seismic velocity of refracted waves and the seismic signature of surface waves have been used as a tool for the distinction between accumulation and ablation zones. Again, in our case this distinction was not possible using radar methods due to the low frequency of the radar. 5) The acoustic impedance calculated from seismic data has allowed to estimate the nature of subglacial sediments and, in particular, whether or not they are deformable. As our study has made evident, the fact that seismic and radar waves sample different properties of the medium make the estimates from both methods independent, so that in general they do not reduce the error of each other. This could be interpreted as not getting an added benefit from their joint use. However, the main benefit arises when the medium properties or conditions make one of the methods efficient while the other is not (for instance, measuring ice thickness or detecting the firn–ice transition), or when the use of both methods provides bounds for the physical property being measured, as is the case for water content estimates. In addition to the latter, in our case study the joint use of both techniques showed its greatest benefit for crevasse detection.

Acknowledgments This work was funded by project ANT99-0963 from the Spanish Antarctic Research Program, project REN2002-03199/ANT from the Spanish Natural Resources Research Program, grant 99-05-65551 from the Russian Fund of Basic Research and funding from the RED Research Program.

References Arcone, S.A., 2002. Airborne-radar stratigraphy and electrical structure of temperate firn: Bagley Ice Field, Alaska, USA. J. Glaciol. 48 (161), 317 – 334. Atre, S.R., Bentley, C.R., 1993. Laterally varying basal conditions beneath Ice Streams B and C, West Antarctica. J. Glaciol. 39 (133), 507 – 514. Benjumea, B., Teixido´, T., 2001. Seismic reflection constrains on the glacial dynamics of Johnsons Glacier, Antarctica. J. Appl. Geophys. 46 (1), 31 – 44. Benjumea, B., Macheret, Yu.Ya., Navarro, F., Teixido´, T., 2003. Estimation of water content in a temperate glacier from radar and seismic sounding data. Ann. Glaciol. 37, 317 – 324. Bentley, C.R., 1975. Advances in geophysical exploration of ice sheets and glaciers. J. Glaciol. 15 (73), 113 – 135. Clarke, T.S., Bentley, C.R., 1994. High-resolution radar on Ice Stream B2, Antarctica: measurements of electromagnetic wave speed in firn and strain history from buried crevasses. Ann. Glaciol. 20, 153 – 159. Clee, T.E., Savage, J.C., Neave, K.G., 1969. Internal friction in ice near its melting point. J. Geophys. Res. 74 (4), 973 – 980. Crary, A.P., 1963. Results of United States traverses in East Antarctica, 1958–1961. IGY Glaciological Report 7. Daniels, D.J., 2004. Ground Penetrating Radar, 2nd edition. IEE, London, UK. 726 pp. Davis, J.L., Annan, A.P., 1989. Ground penetrating radar for high resolution mapping of soil and rock stratigraphy. Geophys. Prospect. 37, 531 – 551. Doake, C.S.M., 1984. Ice-shelf densities from a comparison of radio echo and seismic soundings. Ann. Glaciol. 5, 47 – 50. Drewry, D.J., 1975. Comparison of electromagnetic and seismicgravity ice thickness measurements in East Antarctica. J. Glaciol. 15 (73), 137 – 150. Duval, P., 1977. The role of the water content on the creep rate of polycrystalline ice. Int. Assoc. Hydrol. Sci. 118, 29 – 33. Furdada, G., Pourchet, M., Vilaplana, J.M., 1999. Characterization of Johnsons Glacier (Livingston island, Antarctica) by means of shallow ice cores and their tephra and 137Cs contents. Acta Geol. Hisp. 34 (4), 391 – 401. Gades, T.M., Raymond, C.F., Conway, H., Jacobel, R.W., 2000. Bed properties of Siple Dome and adjacent ice streams, West Antarctica, inferred from radio-echo sounding measurements. J. Glaciol. 46 (152), 88 – 94. Gogineni, S., Chuah, T., Allen, C., Jezek, K., Moore, R.K., 1998. An improved coherent radar depth sounder. J. Glaciol. 44 (148), 659 – 669. Hunter, J.A., 1971. A computer method to obtain the velocity–depth function from seismic refraction data. Report of Activities Part B. Pap.-Geol. Surv. Can. 71-1B, 40 – 48. Kohler, J., 2002. Polythermal glacier firn and ice stratigraphy imaged with ground-penetrating radar. Eos Trans.Am. Geophys. Union 83 (47) (Fall Meet. Suppl., Abstract C62A-0906). Kohnen, H., 1974. The temperature dependence of seismic waves in ice. J. Glaciol. 13 (67), 144 – 147.

F.J. Navarro et al. / Journal of Applied Geophysics 57 (2005) 193–211 Lawton, D.C., 1989. Computation of refraction static corrections using first-break traveltime differences. Geophysics 54 (10), 1289 – 1296. Levato, L., Veronese, L., Lozej, A., Tabacco, E., 1999. Seismic image of the ice–bedrock contact at the Lobbia glacier, Adamello Massif, Italy. J. Appl. Geophys. 42 (1), 55 – 63. Looyenga, H., 1965. Dielectric constants of heterogeneous mixture. Physica 31 (3), 401 – 406. Macheret, Yu.Ya., 2000. Otsenka soderzhaniya vody v lednikakh po giperbolicheskim otrazheniyam (Estimation of water content in glaciers from hyperbolic reflections). Mater. Glyatsiologicheskikh Issledovaniy [Data of Glaciological Studies] 89, 3 – 10. Macheret, Yu.Ya., Glazovsky, A.F., 2000. Estimation of absolute water content in Spitsbergen glaciers from radar sounding data. Polar Res. 19 (2), 205 – 216. Macheret, Yu.Ya., Moskalevsky, M.Yu., Vasilenko, E.V., 1993. Velocity of radio waves in glaciers as an indicator of their hydrothermal state, structure and regime. J. Glaciol. 39 (132), 373 – 384. Moore, J.C., P7lli, A., Ludwig, F., Blatter, H., Jania, J., Gadek, B., Glowacki, P., Mochnacki, D., Isaksson, E., 1999. Highresolution hydrothermal structure of Hansbreen, Spitsbergen, mapped by ground-penetrating radar. J. Glaciol. 45 (151), 524 – 532. Murray, T., Stuart, G.W., Fry, M., Gamble, N.H., Crabtree, M.D., 2000. Englacial water distribution in a temperate glacier from surface and borehole radar velocity analysis. J. Glaciol. 46 (154), 389 – 398. Nobes, D.C., 1999. The directional dependence of the ground penetrating radar response on the accumulation zones of temperate Alpine glaciers. First Break 17 (7), 249 – 259. Nolan, M., Echelmeyer, K., 1999. Seismic detection of transient changes beneath Black Rapids Glacier, Alaska, U.S.A.: I. Techniques and observations. J. Glaciol. 45 (149), 119 – 131.

211

Pakiser, L.C., Black, R.A., 1957. Exploring for ancient channels with the refraction seismograph. Geophysics 22 (1), 32 – 47. Plewes, L.A., Hubbard, B., 2001. A review of the use of radio-echo sounding in glaciology. Prog. Phys. Geogr. 25 (2), 203 – 236. Raymond, C.F., Harrison, W.D., 1975. Some observations on the behavior of liquid and gas phases in temperate glacier ice. J. Glaciol. 14 (71), 213 – 234. Retzlaff, R., Bentley, C.R., 1993. Timing of stagnation of Ice Stream C, West Antarctica, from short-pulse radar studies of buried surface crevasses. J. Glaciol. 39 (133), 553 – 561. Riznichenko, Yu.V., 1949. O rasprostranenii seysmicheskikh voln v diskretnykh y geterogennykh sredakh ([On propagation of seismic waves in discrete and heterogeneous media]). Izv. Acad. Sci. USSR Ser. Geogr. Geophys. 13 (2), 115 – 128. Robin, G. de Q., 1958. Seismic shooting and related investigations. Glaciology III, Norwegian–British–Swedish Antarctic Expedition 1949–1952. Sci. Results 5, 1 – 134. Smellie, J.L., Liesa, M., Mun˜oz, J.A., Sa`bat, F., Palla`s, R., Willan, R.C.R., 1995. Lithostratigraphy of volcanic and sedimentary sequences in central Livingston island, South Shetland Islands. Antarct. Sci. 7, 99 – 113. Smith, A.M., 1997. Basal conditions on Rutford Ice Stream, West Antarctica, from seismic observations. J. Geophys. Res. 102 (B1), 543 – 552. Vasilenko, E.V., Sokolov, V.A., Macheret, Yu.Ya., Glazovsky, A.F., Cuadrado, M.L., Navarro, F.J., 2002. A digital recording system for radioglaciological studies. Bull.-R. Soc. N. Z. 35, 611 – 618. Watts, R.D., England, A.W., 1976. Radio-echo sounding of temperate glaciers: ice properties and sounder design criteria. J. Glaciol. 17 (75), 39 – 48. Weber, J.R., Andrieux, P., 1970. Radar sounding on the Penney ice cap, Baffin Island. J. Glaciol. 9 (55), 49 – 54. Ximenis, L., 2001. Dina`mica de la Glacera Johnsons (Livingston, Shetland del Sud, Anta`rtida). (PhD thesis, University of Barcelona. 194 pp.)