R A D A R I M A G I N G G L A C I O - V O L C A N I C S T R A T I G R A P H Y : M T . W R A N G E L L , A L A S K A by G U Y M A T T H E W C R O S S B.A.(Geophysics), The State University of New York, 1982 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S Department of Geophysics and Astronomy We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A Apr i l 1987 © G u y Mat thew Cross, 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Geophysics and Astronomy The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date February 24, 1987 DE-6(3/81) ii A B S T R A C T An airborne radar survey was conducted over the ice-filled volcanic caldera at Mt. Wrangell, Alaska. Research reported here involves computer processing and interpre-tation of radio-reflection data acquired over 21 traverses of the summit. In addition to describing useful data enhancement techniques, a dynamic programming approach is introduced for topographically controlled data positioning and spatial correction. Interpretation focusses upon a well defined radio-stratigraphy attributed to high acid-ity horizons deposited at the ice surface during periods of elevated volcanic activity. A comparative analysis of layer character indicates that echoes from the caldera floor are not continuously detected because of anomalously high signal absorption. Con-sequently, results impose a lower limit upon maximum ice thickness. A numerical interpretation scheme, incorporating both glaciological measurements and empirical relations governing the behaviour of firn and ice, is developed to aid interpretation of the glacio-volcanic stratigraphy. Preliminary modelling yields a speculative vol-canic record that roughly matches the known eruption sequence at Mt. Wrangell and suggests a significant extension of the volcanic history. i i i With great admiration for my parents. iv TABLE OF CONTENTS ABSTRACT ii LIST OF TABLES vi LIST OF FIGURES vii ACKNOWLEDGMENTS ix CHAPTER I INTRODUCTION 1 1.1 Mt. Wrangell 1 1.2 Glacier-Volcano Interaction 3 1.3 Radioglaciology 5 CHAPTER II DATA ACQUISITION AND PROCESSING 7 2.1 Data Acquisition 7 2.2 Data Enhancement 8 2.2.1 Background 8 2.2.2 Median Residual Filtering 9 2.2.3 Average Stacking 9 2.2.4 Median Despiking • • ., • 12 2.2.5 Conclusion 12 2.3 Data Positioning and Spatial Corrections 13 2.3.1 Introduction 13 2.3.2 Flight Trajectory and Vertical Distortion 15 2.3.3 Ground Speed Variation and Lateral Distortion 18 2.3.4 Interpretation of the Warping Function 28 2.3.5 Conclusion . 28 2.4 Interpretation Model and Time-Depth Conversion 30 2.4.1 Time-Depth Conversion 30 2.4.2 Interpretation Model 31 2.4.3 Conclusion 33 CHAPTER 111 ICE THICKNESS AND SIGNAL ABSORPTION . . . 35 3.1 Background 35 3.2 Analysis of Deep Reflections 36 3.3 Signal Absorption Rates 44 3.4 Estimation of Volcanic Impurity Concentration 49 3.5 Frequency-Selective Scattering Effects 53 3.6 Conclusion 54 CHAPTER IV RADIO-STRATIGRAPHY AND VOLCANIC HISTORY 55 4.1 Internal Reflection Mechanisms 55 4.2 Practical Estimation of Reflection Coefficients .58 V 4.3 Composite Reflections and the Volcanic Record 59 4.4 Speculation on Volcanic History 67 4.4.1 Determination of Depositional Age 67 4.4.2 Analysis of Surface Flow Divergence 71 4.4.3 Conclusion: A Speculative Volcanic Record 72 REFERENCES 77 APPENDIX GLACIER MOTION MEASUREMENTS . . . . . . . . 81 vi L I S T O F T A B L E S 3.1. Average normalized power variances .45 3.2. Normalized power variances observed in polar studies 45 3.3a. Linear coefficients for HF concentration 7 < 0.01 .50 3.3b. Linear coefficients for HF concentration 0.01 < 7 < 0.035 50 3.4. Acid concentrations in polar ice sheets 51 4.1. Known eruption record for Mt. Wrangell 76 A.l. Horizontal motion components; 1965-66 83 A.2. Horizontal motion components; 1975-76 85 vii LIST OF FIGURES 1.1. Project location map 2 1.2. Evolution of the summit caldera at Mt. Wrangell 3 1.3. 1975 photograph of the summit caldera . 4 1.4. Traverse-line map 6 2.1. Median residual filtering 10 2.2. Average stacking 11 2.3. Median despiking . 13 2.4. Cumulative enhancement effect 14 2.5. Formation of the reference topographic profile 16 2.6. Formation of the data relief profile 17 2.7. Compensation for linear flight trajectory 19 2.8. Vertical rescaling 20 2.9. Generation of the coherency map 22 2.10. Dynamic programming constraints 23 2.11. Spatial warping function 26 2.12. Application of the spatial warping function 27 2.13. Effect of spatial corrections . 29 2.14. Effect of time-depth conversion 34 3.1. Preliminary interpretation of 1976 radio-echoes 37 3.2. Preanalysis of 1978 radio-echoes 38 3.3. Deep reflection events 39 3.4. Sampling strategy for power variance analysis 41 viii 3.5. Power profiles with normalized variance 42 3.6. Normalized power distributions 43 3.7. Absorption rate with frequency-temperature dependence 47 3.8. Absorption rate at 60 and 840 MHz 48 3.9. High frequency conductivity vs. HF concentration 52 4.1. Absorption rate due to thermally activated loss mechanism 60 4.2. Depositional age profile 65 4.3. Spatially corrected depth-section over traverse B l 68 4.4. Evidence for composite reflections 69 4.5. Horizontal surface flow rate and divergence 73 4.6. Speculative volcanic record 75 A.l. 1965-66 motion survey map 82 A.2. 1975-76 motion survey map 84 " ix A C K N O W L E D G E M E N T S Fi rs t and foremost, 1 am grateful to my supervisor, G . K . C . Clarke. His thoughtful guidance and enduring patience are most appreciated. I also thank C S . Benson of the Geophysical Institute, University of A laska, and B . T . Prager, B . B . Narod and W.S .P . Leaney of U B C for fruitful discussions and advice. Final ly, I express my appreciation to the entire faculty, staff and graduate student body of the Department of Geophysics and Astronomy, U B C , for providing an interesting and pleasant environment. Research has been funded by the Nat ional Sciences and Engineering Research Counc i l of Canada. C H A P T E R I 1 I N T R O D U C T I O N 1.1 Mt. Wrangell The Wrangell Mountains are located in the south-eastern interior of Alaska and roughly parallel the coastal Chugach range (Figure l . l ) . Bo th ranges are westward-extending branches of the St. El ias Mounta ins, which are located in the south-western extreme of the Yukon Territory. Owing to their geographical location and high ele-vat ion, the Wrangell Mountains receive abundant snowfall that maintains thorough glaciat ion of the range. Despite an icy outward appearance, the W7rangell moun-tains constitute an active link in the circum-Pacif ic chain of andesitic volcanism. M t . Wrangel l , itself, is a broad shield volcano, the summit of which rises 4317 m (14,163 ft) above sea-level. The summit is characterized by a geologically recent collapse caldera, w i th three post-calderan craters situated along the northern r im (Figures 1.2 and 1.3). Al though recent volcanic activity has mainly occurred within the marginal craters, high geothermal heat flux is also expected at the caldera floor. Nonetheless, glacier ice completely fills the roughly 4 x 6 k m oblong caldera and overflows into Long Glacier at a breach in the southern r im. In addit ion to posing a significant flood hazard, this unique situation provides an ideal opportunity to study glacier-volcano interaction. Figure 1.1. Project locat ion, depicting the Wrangell Mountains in relation to associated ranges. F i g u r e 1.2. Postulated evolution of the summit caldera at M t . Wrangell;. (a and b) Erupt ions from the pre-calderan summit deplete the magma chamber below, (c) Format ion of the caldera by collapse into depleted magma chamber. (d)i Cool ing and further subsidence of the magma chamber give rise to" the present ice mantled condit ion (adapted from Press and Siever, 1974). 1.2 G l a c i e r - V o l c a n o I n t e r a c t i o n Researchers from the Geophysical Institute of the University of A laska measured' surface elevation wi th in the summit caldera in 1965, 1975, 1976 and 1978. W i t h the exception of sites in close proximity to the anomalously active Nor th Crater , the ele-vat ion was found to remain constant, within 1 m, over the 13 year period (Bingham, 1967; Mo tyka , 1983). In addi t ion, a mean surface temperature of -21 °C was re-ported wi th a suggestion that seasonal melting at the surface is negligible. On the basis of these observations, Benson and Motyka (1978) concluded that net annual ice-equivalent accumulat ion at the surface of the caldera is roughly balanced by a F i g u r e 1.3. 1975 photograph of the summit caldera, denoting the location of the three post-calderan craters and Long Glacier . Numbers identify survey control features. Dashed lines indicate survey lines. combinat ion of geothermally fuelled basal melting and outflux of ice v ia Long Glacier . The approximate conservation of ice-equivalent precipi tat ion to and f rom the caldera glacier maintains the surface in a quasi steady-state condit ion and has prompted a "giant calorimeter" analogy. By determining the ice-equivalent surface accumulat ion rate and an approximate ice flux at the head of Long Glacier , a rough estimate of the overall basal melt rate can be obtained. Assuming that the present accumulat ion rate 5 has persisted for the last 2000 years, Benson and Motyka (1978) estimated an upper l imit of 140 H F U (5.88 W / m 2 ) on the equi l ibr ium heat flux at the caldera floor. More precise estimates require detailed knowledge of ice thickness and flow velocity wi th in the caldera, especially in the vic ini ty of Long Glacier. 1.3 Radioglaciology A programme of ground-based radar probing experiments was init iated by the Geophysical Institute in 1976 (Benson, 1982). In addit ion to detecting internal re-flections from wi th in the ice mass, localized echo-soundings suggested a maximum ice thickness approaching 900 m in the central region of the caldera. Detection of internal layer echoes indicated the existence of a well-preserved glacio-volcanic stratigraphy, and supported a proposed ice-coring project that would furnish a combined record of cl imatic var iat ion and volcanic activity. In Ap r i l of 1982, the U B C 840 M H z ice radar system was employed along 21 airborne traverses of the caldera (Figure 1.4). The reconnaissance was intended to guide the selection of core dri l l ing sites and pro-vide continuous caldera geometry profiles for confirmation of previous echo-soundings. A l though reflections from the glacier bed were restricted by anomalously high signal absorpt ion rates, a well defined radio- stratigraphy was recorded, and is attr ibuted to high-acidity horizons deposited at the ice surface during periods of elevated volcanic activity. P r io r to presenting interpretive results, however, the data processing phase of the project is described, including signal enhancement techniques and a computer-based method for posit ioning airborne data. 6 F i g u r e 1.4. Fl ight- l ine map. Fl ight- l ine designation is indicated at the head of each traverse. CHAPTER II 7 DATA ACQUISITION AND PROCESSING 2.1 D a t a A c q u i s i t i o n The U B C radio-echo sounder (Narod and Clarke, 1980; 1983) operates at a centre frequency of 840 M H z . A bandwidth of 40 M H z requires a Nyquist sampling fre-quency of 80 M H z which, in turn, implies a sampling interval of 12.5 ns. To avoid the required sampl ing rate, the U B C system is designed for data acquisit ion based upon a strobing technique. Rather than sample a single returned waveform at the necessary rate, we acquire a 1024-point data trace over 1024 consecutive returns. A sampling t ime base incrementally adjusts the acquisition t ime and generates a trigger pulse for the transmitter. A pulse repeti t ion rate of 10 kHz provides an effective sampling in-terval of 10 ns. Al though small-scale variations are produced by interference effects, the average structure of the waveform is assumed to remain consistent over the 1024 returns required to acquire a sampled version. Pract ical considerations support this assumption. The system antenna is characterized by an E-plane beam width of 18° (3 d B down) parallel to flight heading. A t a ground clearance of 200 m, the radar i l luminates a surface region wi th a dimension of approximately 60 m along the flight-track. In comparison, at a ground speed of 60 m/s (roughly 120 m i / h r ) , the aircraft proceeds approximately 6 m along the traverse during the roughly 100 ms required to fully sample the returned waveform. Clearly, the data redundancy is more than sufficient to justify the strobing assumption. 8 2.2 Da ta Enhancement 2.2.1 Background Histor ical ly, radioglaciological data have been recorded photographically. The re-turned waveform is displayed for recording purposes on an oscilloscope in either A -scope or Z-scope formats. A l though the Z-scope format can be used to generate intensity-modulated t ime sections for interpretation, photographic recording restricts the appl icat ion of post-acquisit ion enhancement and display techniques. The U B C system output , on the other hand, is recorded on analogue magnetic tape and sub-sequently input to an analogue/digital converter. The resulting digital data format allows the use of seismic data enhancement techniques (Sherrif and Geldart , 1983; Prager, 1983). A simple and effective enhancement scheme is employed to improve the signal to noise (S /N) rat io. Elements of this scheme include median-based filtering techniques and lateral trace averaging (stacking). Median filters are part icular ly well suited for removal of spurious noise from signals that are characteristically monotonic over n-length segments or constant over a duration of (n + l ) / 2 samples (Evans, 1982; Nodes and Gal lagher, 1982). A combinat ion of power fading in the returned signal and receiver design causes typical ice-sounding radar data to possess such a character. For this reason, median-based niters are given a prominent role in the applied enhancement scheme. Subsequent sections describe individual procedures in order of application. 9 2.2.2 Median Residual Filtering The power envelope of a radio-echo is typically characterized by a sharp onset to maximum amplitude. Subsequently, the power decays monotonically at a fading rate controlled by surface roughness. Relatively low amplitude reflection events are commonly superimposed upon the fading tail of major echoes, especially upon the surface reflection. Median residual filtering removes gross structure caused by large scale-fading and emphasizes less apparent echoes. The task is accomplished by sub-tracting a median filtered version from the raw data trace. Although for practical purposes the median removes low frequency trends, the residual is not strictly a "high passed" result. Unlike the moving average, which effectively smoothes sharp signal variations, the running median filter (RMF) is non-linear and preserves straight-edged step features while events of relatively short duration are suppressed. As a result, the median residual is characterized by less pronounced side lobes and greater resolution of small scale events (Figure 2.1). Residual filtering is a preparatory procedure for subsequent median despiking (Sec-tion 2.2.4). Claerbout and Muir (1973) suggested that the desired RMF window length may require compensation for effective contraction due to regional signal gradients. Removal of gross structure by residual filtering allows subsequent application of an optimum median despiker. 2.2.3 Average Stacking Average stacking (Robinson, 1970; Prager, 1983) replaces each data trace by a spatially-averaged version estimated over n adjacent traces. The stacked trace is 10 RAW DATA TRACE 0 400 time interval MEDIAN SMOOTHED VERSION time interval RESIDUAL DATA TRACE 0 400 800 time interval Figure 2.1. Median residual filtering, (a) Raw data trace, (b) Smoothed data trace obtained by application of a 51-point running median filter, (c) Residua] data trace formed by subtracting broad scale structure from the raw data trace. 11 RESIDUAL DATA TRACE I 3 5-O 0 400 time interval 800 7-FOLD STACKED TRACE ^0.8 ti 0.4 c3 0.0 — \ (1 1/ 1 \ i i i i • i ( b ) 0 400 time interval 800 F i g u r e 2.2. Average stacking, (a) Median residual data trace (see previous figure), (b) Average stacked trace formed at the midpoint of a 7-trace window. formed at the midpoint of the n-length window. A n underlying assumption for opti-m u m stacking requires that the n adjacent traces be identical except for addit ive noise. For purely random noise, an n-fold average stack improves the S / N ratio by a factor of y/n. In Section 2.1 it was demonstrated that a ground clearance of 200 m and a ground speed of 60 m/s results in roughly 10-fold data redundancy. Nonetheless, the stack fold must be selected wi th regard for a trade-off between S / N ratio and smearing of coherent energy. A conservative 7-fold stack is applied in the present scheme (Figure 2.2). 12 2.2.4 Median Despiking Median residual fi ltering (Section 2.2.2) attempts to retain only gross structure in in i t ia l f i l tering, so that events of shorter duration are emphasized in the residual trace. Wh i le enhancing desirable features, however, the method also preserves spurious noise impulses. Al though stacking reduces random noise appreciably, a running median despiker is applied to further improve the S / N ratio. Appl icat ion of an n-length R M F preserves features wi th durat ion exceeding approximately (n — l ) / 2 samples. Events of shorter duration are suppressed. In comparison to the median residual filter, the median despiker operates over a relatively brief window defined by the duration of the shortest non-noise event. According to Evans (1982), the opt imum R M F window length is given by n = 2{m/s - 3 /2) , (2.1) where m represents the duration of the shortest non-noise event and s is the sampling interval. The U B C system is characterized by a sounding pulse durat ion of 50 ns, including a rise time of 18 ns and a fall time of 28 ns. The effective sampl ing interval is 10 ns. Assuming that the min imum non-noise event durat ion is roughly equivalent to the total pulse durat ion, we obtain an opt imum R M F window length of 7 samples (Figure 2.3). 2.2.5 Conclusion The cumulative effect of data enhancement procedures is i l lustrated by Figure 2.4. The simple methods described above yield impressive enhancement of the raw data. Median-based filtering techniques have proved to be well suited for enhancement of 13 7-FOLD STACKED TRACE 0 400 800 time interval MEDIAN FILTERED TRACE 0 400 800 time interval F i g u r e 2.3. Med ian despiking. (a) 7-fold stacked residual trace (see previous figure), (b) Fi l tered trace obtained by application of a 7-point median despiker. ice-sounding radar data and should find future application in radioglaciology. Possible improvements upon these techniques might include non-uniform weighting (Claerbout and M u i r , 1973) and iterative filter application (Fi tch and others, 1985). 2.3 D a t a P o s i t i o n i n g a n d S p a t i a l C o r r e c t i o n 2.3.1 I n t r o d u c t i o n Sound deductions concerning reflector geometry and ice dynamics require spatially accurate data. A l though airborne traverses are controlled by maintaining a constant flight-heading between obvious topographic features, regional and local variations in (a) (b) Figure 2.4. Cummulat ive enhancement effect, (a) Raw time section, (b) Equivalent t ime section after enhancement processing. Figures 2.1 - 2.3 i l lus-trate the sequence of processing steps leading from the raw time section (a) to the enhanced result (b). 15 al t i tude and ground speed produce spatial distortions in the recorded data. Conse-quently, the surface reflection profile is a spatially distorted record of terrain relief along the traverse. By constraining the distorted profile to match the actual topography, spat ial sampl ing errors are compensated and individual data traces are approximately located in northing, easting and elevation coordinates. A computat ional procedure for profile alignment is described and i l lustrated below. 2.3.2 Flight Trajectory and Vertical Distortion The alignment procedure begins wi th preparation of reference (topographic) and data (surface reflection) profiles. The reference profile is obtained directly from the topographic base map. Assuming accurate orientation of the traverse line, intersec-tions of the flight-track wi th topographic contours are digit ized, providing an elevation profile constrained in northing and easting coordinates (Figure 2.5). A tensioned cubic spline is fit to the resulting profile and resampled at equal intervals. The sampling interval is adjusted so that the reference topographic profile and the surface reflection profile are of equal length. The surface reflection profile is obtained as the inverse profile of first arrival t ime along the traverse (Figure 2.6). A moving average operator is applied to reduce receiver chatter over the resulting profile. If alt i tude and ground speed remain constant over the traverse, reference and data profiles should be identical. Comparison of the resulting profiles establishes the existence of pronounced spatial distortions in the data profile. Most notably, the regional relief gradient is reversed. Th is effect is explained in terms of the gross flight trajectory. Since the surface reflection profile represents a continuous record of sounding range, a simple linear flight trajectory, departing from constant alt itude, 16 (a) ELEVATION PROFILE 1000 2000 3000 easting coordinate SPUN ED TOPOGRAPHIC PROFILE 4 120 ? 4 100 4080 4060 4040 d p 4020 Cu — J CD 4000 3980 (c) 1000 2000 3000 easting coordinate Figure 2.5. Format ion of the reference topographic profile, (a) Intersections between traverse line A l (dotted) and topographic contours are located in horizontal coordinates, (b) Traverse line-contour intersections form the basic elevation profile, (c) Tensioned cubic spline interpolation produces a naturally smoothed profile. RAW SURFACE REFLECTION PROFILE 0 200 400 600 trace SMOOTHED REFLECTION PROFILE 0 200 400 600 trace Figure 2.6. Format ion of the data relief profile, (a) T ime section depict ing raw surface reflection, (b) Surface reflection profile obtained directly as the inverse profile of first arr ival t ime, (c) Appl icat ion of moving average operator produces smoothed relief profile. 18 is capable of producing the observed trend reversal. Al though the actual trajectory is l ikely more complex, the regional gradient is removed from both reference and data profiles as a first approximation (Figure 2.7). The resulting residual data profile retains second-order vertical distortions caused by local alt i tude variations. These local fluctuations are compensated to a large extent by piece-wise linear amplitude rescaling between obvious control features that are mutually apparent in the residual profiles (Figure 2.8). Large-scale topographic extrema serve as control features; in the case of M t . Wrangel l : the caldera r im and central crest. A t this stage, remaining dissimilarities between the residual profiles are assumed to result from regional and localized ground speed variations. 2.3.3 Ground Speed Variation and Lateral Distortion Deviations from constant ground speed result in relative compressions and/or ex-tensions of the data profile. These lateral distortions are corrected by exploit ing a dynamic programming approach to waveform alignment (Myers and Rabiner , 1981; Anderson and Gaby, 1983). The vertical rescaling procedure, described in the previ-ous section, effectively involves the application of a user supplied matching function in vert ical profile coordinates. Dynamic waveform alignment determines a completely analogous warping function in horizontal coordinates by means of a two-stage proce-dure. The init ial stage of the waveform alignment procedure involves a statist ical com-parison of the vertically corrected residual profiles. A coherency map is generated by calculat ing a normalized coherency measure, S(r,d), over a two-dimensional grid. S(r, d) is formed by computing the semblance (Neidell and Taner, 1971) over /-length REFERENCE TOPOGRAPHIC PROFILE DATA RELIEF PROFILE NORMALIZED RESIDUAL TOPOGRAPHY NORMALIZED RESIDUAL SURFACE RELIEF F i g u r e 2.7. Compensation for gross linear flight trajectory. A strong linear flight gradient over the actual topography, (a), produces an apparent reversal of the relief trend recorded, (b). The trend reversal is compensated by removing the regional gradient from reference and data profiles. Corresponding residual profiles are obtained, (c) and (d). 20 REFERENCE TOPOGRAPHY / DATA RELIEF 0 200 400 600 trace REFERENCE / VERTICALLY RE SCALED DATA 0 200 400 600 trace F i g u r e 2.8. Vert ical rescaling. (a) Comparison of normalized residual pro-files, (b) The residual data profile is linearly rescaled relative to the reference profile between control features. Dashed and solid profiles represent data and reference profiles respectively. 21 sliding windows with midpoints r and d representing consecutive reference and data profile samples respectively (Figure 2.9). The coherency map is, thus, defined in terms of local semblance as -m+i) S(r,d) = - , (2.2) 12[R(r-m+i)+D{d-m+i)) for r, d = m + l,m + 2 , N — m, where TV is the profile length, m = ((1/2) - 1/2) is the window width and i ? ( r _ m + t ) and D^-m+t) represent reference and data profile amplitudes at samples (r — m+i) and (d—m+t) respectively. Semblance may be viewed as the ratio of combined energy to the sum of individual energies within the windows of comparison. An acute sensitivity to local amplitude diversity makes semblance a well suited criterion of profile likeness. A second stage in the alignment process involves an augmented path-finding pro-cedure for determination of the lateral warping function [VTr(t),VV^i)], with i = 1,2,n . The procedure optimizes a piece-wise linear path over the coherency map by locally maximizing the semblance level subject to dynamic programming constraints (Figure 2.10). In addition to imposing physically grounded restrictions upon allowable feature association, programming constraints significantly reduce necessary computa-tions. Constraints are applied at both local and global levels and include: 22 GENERATION OF THE COHERENCY MAP oo ^ o co o d d o c j I I 100 200 300 400 500 600 normalized residual topography 700 Figure 2.0. Generat ion of the coherency map. Residual reference and data profiles are compared over sample windows having midpoints r and d respec-tively. The coherency measure is calculated and entered at midpoint coor-dinates. The resulting coherency distr ibut ion is characterized by a well de-fined non-linear ridge that roughly determines the opt imum warping funct ion. G loba l path constraints restrict necessary computations. DYNAMIC PATH FINDING/GLOBAL PATH CONSTRAINTS INSET/LOCAL PATH CONSTRAINTS 0 tOO ZOO 300 400 500 600 700 160 170 180 190 200 210 220 normalized residual topography normalized residual topography Figure 2.10. Dynamic programming constraints. Constraint boundaries are indicated by solid lines. Shaded areas define allowed search regions on both global and local levels. 24 E n d p o i n t C o n s t r a i n t : The data relief profile and reference topographic profile are assumed to span an identical traverse segment. Despite spatial distortion of the data profile wi th in the interval, endpoints are required to coincide. The endpoint constraint init iates and concludes the path-finding procedure by stipulating that \Wr(i),Wi{i)] = [i,i], i= 1,2, . . . , ' m + l , n . (2.4) L o c a l S l o p e C o n s t r a i n t : Fundamentally, the local slope constraint (LSC) re-quires that the path be extended in a monotonically non-decreasing fashion. More stringent restrictions are subjective and related to l imitations upon aircraft accelera-t ion and max imum expected distort ion. The L S C is defined as where L S C > 1 . A n addit ional penalty-based constraint may be imposed in conjunction wi th the L S C . In effect, such a constraint involves the appl icat ion of a local weighting template in favor of linear path extension. The local semblance level is scaled according to a weighting function of the form S ' ( r d ) - l S W S M ' (2 6) b V>d> - \S-WS(r,d), <5>0' where S — [d — W^{i — 1) ) / ( r — Wr{i — 1)) is the local slope, it; is a positive exponent governing penalty severity and S'(r, d) is the scaled semblance. 25 L o c a l I n t e r v a l C o n s t r a i n t : The local interval constraint (LIC) imposes a max-imum l imitat ion upon step-size by requiring that TTC > fWAi)-Wr{i-l), t = m + 2 , m - | - 3 , . . . , n - l , , -\Wd{i)-Wd{i-l), i = m + 2 , m + 3 , . . . , n - 1 ' K ' where L I C > 1 . In other words, the L I C restricts the range over which the procedure is allowed to "look ahead". Selection of the L I C involves a trade-off between accuracy and execution time. G l o b a l S l o p e C o n s t r a i n t : The global slope constraint (GSC) is simply an ex-tension of the L S C applied upon init iat ion of the procedure; 1 WA(I\ < TT-ji < G S C > i = m + 2 ,m + 3 , . . . , n - 1 , (2.8) G S C - Wr(i) where G S C = L S C > 1 . Al though the G S C does not directly influence the path-f inding procedure, the constraint significantly reduces required computations. G l o b a l R a n g e C o n s t r a i n t : The global range constraint ( G R C ) restricts the max imum cumulative interval between data and reference profile samples required to coincide by the resulting warping funct ion. The G R C reflects an estimate of max imum overall distort ion by requiring that Wr(i) - Wd(i) < G R C , i = m + 2 ,m + 3 , . . . , n - l , (2.9) where G R C > 0 . 26 Dynamic path-f inding yields the non-linear warping function depicted in Figure 2.11. Appl icat ion of the function affects relative compressions and extensions required to compensate spatial sampling errors and align the residual profiles (Figure 2.12). A final transformation involves the reversal of linear regression and normalization operations to obtain the corrected surface reflection profile. S P A T I A L W A R P I N G F U N C T I O N 0 200 400 600 reference trace Figure 2.11. Spat ial warping function determined by dynamic path-f inding, subject to global and local constraints. 27 PRE-WARP PROFILE ALIGNMENT CD 3 ft. • t i CO POST-WARP PROFILE ALIGNMENT 0.8 0.4 cu 3 ti ti o. •ti to CU -0.4 ti g "3 4000 3980 1000 2000 easting coordinate TOPOGRAPHIC ALIGNMENT 3000 1000 2000 easting coordinate 3000 Figure 2.12. Appl icat ion of spatial warping function, (a) Profi le alignment prior to warping, (b) Profile alignment resulting upon application of the spatial warping function, (c) Profi le alignment in topographic coordinates. Dashed and solid profiles represent data and reference profiles respectively. 28 2.3.4 Interpretation of the Warping Function Real iz ing that the surface reflection profile and reference topographic profile are temporal ly and spatially sampled sequences respectively, the warping function is cor-rectly interpreted as a time-distance relation. The derivative of the function is, thus, a continuous record of apparent ground speed along the traverse. Quali tat ively, the interpretation roughly reflects the manner in which traverses were pi loted. Quant i -tatively, the apparent ground speed is marginally unreasonable over certain intervals. Ar t i f ic ia l inflations of the actual ground speed are the consequence of compensating residual alt i tude variations that are incorrectly handled as ground speed fluctuations. Despite gradient removal and ampli tude rescaling procedures, small scale vertical dis-tortions persist and pose a minor source of error in the waveform alignment procedure. 2.3.5 Conclusion In addi t ion to rectifying spatial distortions, the alignment process automatically positions indiv idual data traces in northing, easting and elevation coordinates. The end result is a spatial ly accurate t ime section (Figure 2.13). Comparison of original and corrected versions emphasizes the necessity for spatial correction prior to interpre-tat ion. A l though the procedure is presently adequate, a continuous record of aircraft al t i tude should replace linear regression and ampli tude rescaling operations. Precision alt imeters are available for moderate cost. F i g u r e 2 .13 . Effect of spatial corrections, (a) Enhanced t ime section, (b) Equivalent section following appl icat ion of the warping funct ion for compen-sation of spat ial sampl ing errors. 30 2.4 A n Interpretation Model and Time-Depth Conversion 2.4.1 Time-Depth Conversion Al though the surface reflection profile is spatial ly corrected by way of the dy-namic waveform alignment procedure, the trai l ing data traces remain a function of t ime rather than depth. A final processing step involves t ime-depth conversion in the vert ical coordinate to produce true cross-sectional records for interpretation. Conver-sion requires knowledge of the vert ical variation in electromagnetic wave velocity with depth. Electromagnetic wave velocity is given by where v0 represents the propagation velocity in free space and the relative dielectric permit t iv i ty , e, is modelled as a function of density, p, according to an empir ical relation Depth dependent density variation is turn described by a densification model proposed by Benson (1962) for snow/f irn compact ion on the Greenland Ice Sheet. According to Benson's model , (2.10) due to Rob in and others (1969); e = (1 + 8.5 x 10~ 4 p) 2 - (2.11) p p* (2.12) 31 where p is the depth dependent ice pressure, pj is the density of ice and m 1 ?m2 and p*, the cr i t ical pressure, are empir ical constants. The two-way transit t ime, T2, for propagation to a given depth, 2, is determined by integrating - / = — • 2.13 The preceding equations are coupled wi th addit ional differential equations, glacio-logical measurements and empir ical relations in a unified interpretation scheme devel-oped to model depth dependent variat ion of interrelated physical properties wi th in the caldera glacier. 2.4.2 Interpretation Model The numerical interpretation model (Clarke and Benson, 1982; Clarke and others, 1987) yields the surface-normal depth dependent variation in temperature T, heat flux q, ice pressure p, density p, surface-normal component of flow velocity u>, downslope component of flow velocity u, depositional age ta, two-way propagation time T2 and two-way dielectric propagation loss P2. In addit ion to Equations 2.12 and 2.13, the system of coupled linear differential equations includes: dp — — pg cos a (2.16) du dz 32 ^ = - ^ - A (2.17) dz p dz " 2 5 o e X P ( ^ f ) ( i r p t a n a ) n ( 2 " 1 8 ) dta _ _1 dz w dr2 _ 2^/1 (2-19) (2.20) dz v0 = 2D, (2.21) d P 2 dz where K represents the thermal conductivity, c the specific heat capacity, Q the creep act ivat ion energy, R the universal gas constant, F a shape factor, a surface slope, n the flow law exponent, g gravitat ional acceleration, D the propagation loss rate and A = (du/dx + dv/dy) the horizontal flow divergence. The coupled differential equations are integrated, together wi th Equat ion 2.11 and the following addit ional property relations, using the Runge-Ku t ta method: c{T) = 2115.343 + 7.7929(7/ - TQ) (2.22) K(p,T) = 2pKI(T)/(3pI - p) (2.23) K^T) = 2.1725 - 3.403 x 10 _ 3(r - T0) + 9.085 x 10 _ 5 ( T - T0)2 (2.24) D(T) = D0exp{-E/RT), (2.25) 33 Ol 0 < z< z, z > z3 s (2.26) w i th T0 representing the surface temperature, Kj the thermal conductivi ty of ice, D0 a propagation loss constant, E the activation energy for dielectric loss, A 0 the horizontal flow divergence determined at the surface and za a stagnation depth, below which horizontal divergence is assumed negligible. 2.4.3 Conclusion The effect of t ime-depth conversion is i l lustrated in Figure 2.14. Ready transfor-mat ion from t ime to depth section is another distinct advantage of the digital recording format over photographic recording. Various elements of the interpretation model are explained in greater detail where relevant in subsequent chapters. Figure 2.14. Effect of t ime-depth conversion, (a) Spat ia l ly corrected time section, (b) Spat ial ly corrected depth section. C H A P T E R III 35 I C E T H I C K N E S S A N D S I G N A L A B S O R P T I O N 3.1 Background Ground-based radar probing experiments were conducted at M t . Wrangell in 1976, 1977 and 1978 by the Geophysical Institute of the University of Alaska (Mackeith, un-publ ished; Mo tyka and others, 1980; Benson, 1982). Spot measurements of ice thick-ness were made in 1976 using a 5 M H z echo-sounder developed at the Geophysical Institute. Radio-echoes were consistently received from internal layers at approxi-mately 90 and 200 m, and from the glacier bed to a depth approaching 450 m (Figure 3.1). Reflections from the bed were recorded near margins of the caldera only, suggest-ing that insufficient penetration prevented detection of echoes from greater depths. In 1977 the system antennas were modified to enhance sensitivity, and soundings were repeated in the central region of the caldera. A reflection recorded at a two-way t ime of 14.5 ps implies an ice thickness of 1250 m. Owing to poor directionality of the antennas, it was not ascertained whether the return was from the glacier bed, or surrounding terrain. The Scott Polar Research Institute's (SPRI) Mark IV equipment (Smi th , 1971) was employed for similar ice thickness measurements in 1978. It was an-t ic ipated that improved directionality offered by the system would eliminate interfering reflections originating above the ice surface. Bed reflections were detected from depths approaching 900 m and consistent reflections from average depths of 225 and 350 m were at tr ibuted to internal layering. The maximum depth-sounding range achieved 36 wi th the Mark IV system suggested appreciably higher signal absorption rates than expected (Benson, 1982). A n elevated volcanic impuri ty concentration was postulated to explain the anomalous absorption. Re-analysis of reflections recorded in 1978 suggests an absorption loss rate of ap-proximately 2.8 dB/100 m (Figure 3.2). In addit ion, the linear trend observed for reflections attr ibuted to the glacier bed implies a power reflection coefficient ( P R C ) of -34 d B . In comparison, Robin and others (1969) suggested a min imum reflection coefficient of -22 dB for quartz-based materials. Al though Smith (1971) concluded that bed roughness can effectively lower the P R C by as much as 10 d B , the above estimates of basal P R C and bulk loss rate may well be underestimated. 3.2 Analysis of Deep Reflections The present investigation indicates the existence of three prominent reflecting hori-zons that are continuously detected over the caldera. Al though reflections are inter-mittently received from above and below the major horizons, echoes are not recorded from depths exceeding 350 m. Results suggest that high signal absorption rates re-strict detection of deeper horizons and, thus, that the lowermost continuous horizon should be interpreted as an internal layer rather than the glacier bed. The following lines of evidence support such an interpretation. Over much of the recorded section, semi-continuous reflection events are detected below the horizon in question (Figure 3.3). Both isolated and stratif ied reflections are observed and found to extend laterally for distances approaching 750 m. Because geometrical considerations rule out mult iple reflections as the source of deep-lying events, ice thickness must be at least equivalent to the max imum echo range. In 37 A A" »*• oao no K M I H • t*o 9** >»o me taoo irte k t T i « » RADIO ECHO SOUNDING DATA. MOUNT WRANGELL CALDERA Figure 3.1. Preliminary interpretation of radio-echoes recorded in 1976. Vec-tors along survey lines indicate horizontal flow velocities measured during the same field season. 38 RADIO-ECHOES (1978) -80 -40 0 PRC with absorption (dB) Figure 3.2. Radio-echoes recorded in 1978. The P R C , corrected for geo-metrical spreading, is displayed as a function of depth. Isolated clusters are attr ibuted to internal layers. The remaining points exceeding 350 m depth are presumed to represent reflections from the glacier bed. Linear regression of these points suggests an absorption rate of 2.8 d B / l O O m . 39 APPROXIMATE DISTANCE (m) F i g u r e 3.3. Por t ion of depth-section recorded over traverse-line A l , i l lustrat-ing prominent reflection events below major horizons. addi t ion, stat ist ical analysis of the power returned f rom major horizons demonstrates that the physical character of the lowermost horizon closely resembles that of overlying internal layers. Roughness characteristics of reflecting interfaces have been interpreted on the basis of returned pulse shape and lateral variation in peak returned power (Harr ison, 1972; Oswald , 1975; Nea l , 1976; M i l l a r , 1981a). Oswald (1975) examined the distr ibut ion of relative returned power for suites of closely spaced soundings of polar glacier beds. The distr ibut ions were interpreted as functions of sounding range and reflector roughness in terms of three dominant mechanisms for variat ion in returned power: 40 Interference: Where the sounding range is significantly greater than the radius of curvature for irregularities at the interface, numerous reflections are received simul-taneously. Random interference of reflections results in rapid amplitude fluctuation at the receiver. F o c u s s i n g : Where the sounding range is roughly equivalent to the radius of curva-ture for irregularities at the interface, convergence and divergence of reflected radiation result in relatively smooth ampli tude variation at the receiver (Harr ison, 1971). P u l s e s p r e a d i n g : Where the sounding range is wi th in the radius of curvature for deviations from a planar interface, relatively few reflections are received simultane-ously. Ampl i tude variation at the receiver is related to the mean slope of irregularities. To i l lustrate the present analysis, normalized relative power distributions are ob-tained for the three major horizons over a 200 trace range, representing a distance of approximately 600 m along traverse A10 (Figures 3.4, 3.5). The sample range is restricted by requiring that all horizons remain distinct over the interval. In addit ion, an interface between the outer flank of the caldera and the covering ice-pack is sampled for comparison. Examinat ion of the resulting distributions (Figure 3.6) indicates that power returned from the three internal horizons tends to be distr ibuted symmetrical ly about the mean; typically the result where pulse spreading effects dominate. Similar results are obtained for independent samples and suggest that the somewhat bimodal nature exhibited in the case of lowest horizon in Figure 3.6 is a feature of the particular sample only. In contrast, the power returned from the ice/bedrock interface is char-acterized by an asymmetric distr ibut ion, with deviations from the mean exceeding 20 41 APPROXIMATE DISTANCE (m) Figure 3 .4. Portion of time-section recorded over night-line A10. IHl, IH2 and IH3 denote major internal horizons. BH denotes the ice/bedrock inter-face. Arrows along the horizontal axis indicate ranges over which samples of the P R C were obtained for statistical analysis. 42 POWER PROFILES: BED HORIZON INTERNAL HORIZON 1 INTERNAL HORIZON 2 INTERNAL HORIZON 3 Figure 3.5. Power profiles obtained over sample ranges depicted in Figure 3.4. Vp is the normalized power variance. 43 BED HORIZON INTERNAL HORIZON 1 1.0 j? 0.8 g CD Cr ? 0.6 K ti CD 0.4 I 0.2 0.0 II .1 I L__l I I I I L -20 -10 0 10 20 relative power (dB) 1.0 $0.8 CD ? 0.6 < ti CD 5 0.4 "3 I 0-2 0.0 J I I I I I I I L -20 -10 0 10 20 relative power (dB) INTERNAL HORIZON 2 INTERNAL HORIZON 3 -20 -10 0 10 20 relative power (dB) 1.0 j? 0.8 CD o £ 0.6 ti CD 0.4 -3 0.2 0.0 i i i i i i i i i -20 -10 0 10 20 relative power (dB) Figure 3.6. Normalized power distributions obtained over sample ranges de-picted in Figure 3.4. 44 d B . Oswald (1975) demonstrated that such a distr ibution is expected in cases where interference is the dominant mechanism responsible for lateral variat ion in returned power. Further evidence establishing the lowermost horizon as an internal layer is provided by comparison of the normalized power variance, where P represents returned power and angular brackets denote the spatial mean value. In every case studied, the normalized variance calculated for the lowermost horizon is found to agree roughly with the value for overlying internal layers (see Figure 3.5). Average normalized power variances for the three internal horizons are given in Table 3.1. In comparison, the average normalized variance for power returned from the ice/bedrock interface is approximately an order of magnitude larger. Average variance levels established for the three internal horizons are comparable to typical values tabulated by Neal (1977). Neal 's values for an ice/bedrock interface exceed the average found here by another order of magnitude (Table 3.2). 3.3 Signal Absorption Rates Assuming that the maximum ice-sounding range is l imited by signal absorpt ion, a rough estimate of the propagation loss rate is obtained using a modified form of the radar equation (Narod and Clarke, 1983), (3.1) (Pf (P) 2 y->2\2n2 y2 Pr = a, 2A,2 P*exP(~AC)exp(-2ADri), (3.2) 45 I N T E R F A C E Interna] Horizon 1 0.0012 Internal Horizon 2 0.0015 Internal Horizon 3 0.0028 Ice/Bedrock 0.035 T a b l e 3.1. Average normalized power variances, Vp, for prominent internal horizons and ice/bedrock interface at M t . Wrangel l . Values are determined over an average trace interval of 200 samples. I N T E R F A C E S A M P L E S Internal Layer 360 0.007 Internal Layer 360 0.004 Ice/Bedrock 800 0.82 Ice/Bedrock 380 0.62 T a b l e 3 .2. Typ ica l normalized power variances observed in polar studies (from Neal , 1977). where G is the forward antenna gain; A is the wavelength in air; RA is the ampl i -tude reflection coefficient for a plain reflecting target; T 2 is the power transmission coefficient at the air / ice interface (taken to be 1.0); Pt is the transmitted power and r' = ( r f l + - ) (3.3) is the range adjusted for refraction, where n is the refractive index of ice and ra and r t are the range in air and ice respectively. The recorded power, P'r = P r e x p ( - A C ) , (3.4) is the power received at the antenna, Pr, subsequently attenuated by internal losses, C , in decibels. D denotes the one-way propagation rate in decibels, and the conversion 46 factor, A = ( l o g e 1 0 ) / l 0 , is chosen such that 10 l o g 1 0 e x p ( A D ) = D. A t maximum range, the assumption ra << r, is made and we obtain 10 l o g 1 0 r2max + 2 (D) r m a x = 20 l o g 1 0 G + 10 l o g 1 0 ) V r / mai + 2 0 l o g 1 0 XnT7 8?r + 1 0 l o g 1 0 J # + , (3.5) where bracketed variables denote decibel values. A n estimate of the average one-way propagation loss rate in dB /100 m is, thus, given by D = 200r r 2 ( G ) + ( i>) - | -20 log 1 0 XnT2 8TT 1 0 1 o g 1 0 r m a , - ( C ) (3.6) where P = (Pt/P'T)max is the system performance. Supplying the relevant system pa-rameters and assuming an internal layer reflectivity of -45 d B , the effective propagation loss rate for a max imum range of 350 m is approximately 3.9 d B / 1 0 0 m. A similar calculat ion yields a comparable loss rate of roughly 3.4 d B / 1 0 0 m, corresponding to the max imum range of 900 m achieved with the S P R I 60 M H z equipment in 1978. A bed reflection coefficient of -30 d B is assumed for the latter estimate. Al though use of a higher basal P R C improves the agreement between independent estimates, re-analysis of echoes recorded during the 1978 experiment suggests a reflection coefficient less than -30 d B . Recal l that the loss rate estimated directly from the power-depth distr ibut ion is 2.8 d B / 1 0 0 m (Figure 3.2). 100kHz IMHz 10MHz 100MHz IGHz F R E Q U E N C Y F i g u r e 3.7. Absorpt ion per 100 m in ice, as a function of frequency, with temperature as parameter-(from Evans and Smi th , 1969). Measurements at 100 K H z are due to Paren (1970) for samples from Camp Century, Greenland. Higher frequency measurements are due to Westphal (1963) for samples from Tuto Tunnel , Greenland. 48 Figure 3.8. Absorption per 100 m in ice, at 60 and 840 MHz, as a function of temperature. Values were extracted from a digitized version of Figure 3.7. 49 Evans and Smith (1969) published a graphical compilation of laboratory measure-ments of dielectric absorption rates in polar ice (Figure 3.7). Absorpt ion per 100 m path-length is displayed as a function of frequency over a range of ice temperature between -60 and -1 °C. For the present study, relevant absorption/temperature rela-t ionships are extracted at 60 and 840 M H z (Figure 3.8). Assuming Benson (1982) is correct in postulat ing an ice thickness of 900 m, an approximate vertical temperature dist r ibut ion is calculated v ia the interpretation model. Corresponding predictions of the bulk loss rate at 60 and 840 M H z are obtained according to the respective absorption-temperature relationships for polar ice. In both cases, the predicted loss rate, 1.8 d B / 1 0 0 m at 60 M H z and 2.3 dB /100 m at 840 M H z , is more than 1.5 d B below the estimate based upon maximum sounding range. 3.4 E s t i m a t i o n o f V o l c a n i c I m p u r i t y C o n c e n t r a t i o n Accord ing to Walford (1968), dielectric absorption in ice at frequencies between 0.1 and 100 M H z is due to dielectric relaxation. Deviations from the relaxation spec-t rum above a few hundred megahertz are attr ibuted to the low frequency extreme of the infra-red absorption spectrum. Walford concluded that observed variations from the theoretical spectra are probably due to the effects of ionic impurit ies upon the dielectric relaxation mechanism. It is expected that high volcanic impuri ty levels are responsible for the difference between absorption rates at M t . Wrangell and those observed for polar ice. A l though no measurements of dielectric absorption exist for laboratory ice incorporat ing volcanic impurit ies, M i l l a r (1981a; 1982) suggests that ice doped wi th hydroflouric acid (HF) is the nearest approximation. Campl in and others 50 (1978) analyzed the dielectric behavior of HF-doped ice and determined empirical re-lationships between H F impurity concentration and the high frequency conductivity, (Too, over a range of temperature. M i l l a r (1981a) presented linear approximations of the form CTQO = ^1 + H, where 7 represents the H F concentration in moles per cubic metre (Figure 3.9). Appropr iate coefficients for a given temperature and concentration range are supplied in Tables 3.3a,b. T E M P E R A T U R E (°C) A ( x l O - 3 ) B ( x l O - 6 ) -10 1.85 17.0 -20 1.75 6.0 -30 1.55 2.0 -40 1.15 1.0 -50 0.75 1.0 T a b l e 3.3a. Linear coefficients for 7 < 0.01 (see Figure 3.9). T E M P E R A T U R E (°C) A ( x l O " 3 ) B ( x l O " 6 ) -10 2.95 6.0 -20 1.95 4.0 -30 1.09 6.5 -40 0.58 7.0 -50 0.25 6.0 T a b l e 3.3b. L inear coefficients for 0.01 < 7 < 0.035 (see Figure 3.9). Upper l imits on the bulk impurity concentration at M t . Wrangell are obtained from max imum sounding range estimates of the propagation loss rate. The high frequency conduct iv i ty is given by fjoo = 2ire0(.'f tan 6, (3.7) where £ 0 is the permit t iv i ty of free space, e' is the permitt iv i ty of ice and / tancS is the loss parameter. It is assumed that e' = 3.17 ± 0.07 (Evans, 1965) and the loss parameter is / t a n 6 = — r, (3.8) 2TT 4.343(6') = where D is the loss rate in decibels, c is the electromagnetic wave velocity in a vacuo and 4.343 = 1 0 l o g 1 0 e . The interpretation model yields a mean temperature of ap-proximately -10°C for an ice column 900 m thick. Using the appropriate coefficients, estimates of bulk volcanic impuri ty concentration are 0.010 and 0.012 m o l / m 3 for soundings at 60 and 840 M H z respectively. These bulk levels are comparable to layer concentrations measured in polar ice cores, corresponding to major volcanic eruptions (Table 3.4). E R U P T I O N H+] (mo l /m 3 ) Thera (1390 B C ) 0.005 E ld ja (934 A D ) 0.020 Hek la (1104) 0.006 Lak i (1783) 0.021 Tambora (1815) 0.008 Ka tma i (1912) 0.005 T a b l e 3.4. Observed layer concentration of acids in polar ice sheets for known volcanic eruptions (From Mi l l a r , 1981b). 52 HF CONCENTRATION ( m d n f 3 ) Figure 3.9. High frequency conductivity, as a bilinear function of HF concen-tration, with temperature as parameter (from Millar, 1981a). 53 3.5 Frequency-Selective Scattering Effects A minor discrepancy between estimates of volcanic impurity concentration may reflect increased propagation losses due to scattering at the higher frequency. Although the part icular mechanisms responsible are uncertain, a comparison of the relative magnitude of effects at 60 and 840 M H z is furnished by a hypothetical example. Smith and Evans (1972) modelled scattering from spherical ice inclusions within densified firn, according to a Rayleigh dependence upon frequency. Fract ional power loss, per unit path length, due to scattering from m inclusions of radius b is where A c is the wavelength in a vacuo and e.i and c 2 are the relative permitt ivit ies of ice and densified firn respectvely. The relative permit t iv i ty of densified firn is modelled according to an empirical relation given by Robin and others (1969), where p 2 is the firn density. Adapt ing a simple model for spherical ice lenses (p\ = 917.0 k g / m 3 , £ i = 3.17) of radius 50 m m separated by 250 m m wi th in a densified f irn (p2 = 550.0 k g / m 3 , e2 = 2.15), the predicted loss rate is approximately 0.5/A^ d B / 1 0 0 m. Considering radio propagation through numerous th in scattering layers of combined thickness 25 m, scattering losses at 840 M H z amount to 7.5 d B . Taken over a total ice thickness of 900 m, the total loss is represented by a bulk loss rate of 0.008 d B / 1 0 0 m. In comparison, the scattering loss rate at 5 M H z is negligible. (3.9) e 2 = (1 + 8.5 x IO-V2) 2, (3.10) 54 3.6 C o n c l u s i o n High signal absorption rates, due pr imari ly to volcanic impurit ies, restrict the present interpretation of ice thickness over the central region of the caldera. Statistical analysis of the peak power returned from reflectors at max imum range imposes a lower l imi t of 350 m upon maximum ice thickness. Al though depth-soundings to 900 m with the S P R I Mark IV equipment are not confirmed, comparable estimates of the bulk loss rate at 60 and 840 M H z support the earlier findings. Reflections from the glacier bed are also uncertain in marginal regions of the caldera. Prel iminary analysis, based upon statistics of peak returned power, suggests that dipping reflection events are mainly attr ibutable to internal layers. In addit ion, the presence of distinct diffraction hyperbolae may indicate that the glacier bed descends sharply at the margins, causing reflections from the interface to be cri t ical ly refracted. As a result, returns from the bed would not be detected above the ice surface. Estimates of max imum ice thickness reinforce such a conclusion by imply ing that the glacier bed must descend more steeply than observed reflections would indicate. Despite suspected l imitations on the interpretation of ice thickness, continuous depth-sectional coverage provides useful insight into caldera geometry and glacier vol-cano interaction. It is expected that, to a large extent, caldera structure may be inferred from the topography of internal horizons. Where this is not entirely the case, variations between the topography of internal horizons and the present ice surface might be related to heat flow distr ibut ion at the glacier bed. F inal ly , internal horizons are isochronal surfaces and suggest numerous interpretations of ice dynamics. 55 CHAPTER IV RADIO-STRATIGRAPHY AND VOLCANIC HISTORY 4.1 Internal Reflection Mechanisms The electrical character of dielectrics, including glacier ice, is described in terms of two fundamental properties: a complex relative magnetic permeabil i ty and a complex relative permit t iv i ty. Because ice is non-magnetic, the permeabil i ty is simply 1. The complex relative permit t iv i ty is given by where the real part , e', is simply referred to as the relative permitt iv i ty. Electromag-netic waves propagating through a dielectric are subject to part ia l internal reflection by inhomogeneities wi th in the medium. The power reflection coefficient (PRC) asso-ciated wi th a dielectric discontinuity can be expressed in terms of a complex intrinsic admit tance, For contrast ing intrinsic admittances Tj and T2, and assuming normal incidence upon a smooth interface, the P R C is given by * t , • 11 e = t + it , (4.1) (4.2) ( T I - T $ ) 2 ( T ; + T ; ) (4.3) 56 where RA denotes the amplitude reflection coefficient (Paren, 1981). Numerous physical mechanisms have been proposed to explain internal layer echoes. These mechanisms mainly involve relative permit t iv i ty contrasts caused by systematic fluctuations in the physical properties of ice, including density, crystall ine anisotropy and bubble geometry (Robin and others, 1969; Harr ison, 1973; Paren and Rob in , 1975; C lough , 1977; Ackley and Kel iher, 1979). Debris layers, containing dust and ash, have been discounted as a possible source mechanism on account of insufficient permitt ivi ty contrast. Al though correlation between radio-echo soundings and ice core analyses support the previous mechanisms (Gudmandsen, 1975), there is no such evidence to suggest the importance of debris layers. Density fluctuation is generally accepted as the most satisfactory explanation for internal reflections from above 1000 m. The predominant mechanism responsible for echoes originating from greater depths is less certain. Because the crystall ine structure of ice is intrinsically anisotropic, Harrison (1973) suggested that an anisotropic relative permit t iv i ty might play an important role in the origin of deep reflections. Ackley and Kel iher (1979) postulated that layered variat ion in bubble geometry, caused by deep-seated deformational processes, could enhance the P R C for coincident density contrasts. Paren and Robin (1975) investigated the importance of ionic impurity concentration and the distr ibut ion of impuri ty defects wi th in the crystalline lattice. Unl ike previously discussed mechanisms, ionic impurit ies affect the imaginary portion of the complex relative permitt ivi ty, e", which characterizes internal energy absorption due to the diffusion of lattice defects. Dielectric absorption is commonly described in terms of the loss tangent, e" tan (5 = — , (4.4) 57 w i th the complex relative permit t iv i ty defined as e* = t ' ( l - tan 6). (4.5) Because dielectric absorption is directly related to temperature, echoes resulting from variations in the loss tangent are expected to strengthen wi th depth. Various atmospheric impurit ies enter the glacier system by means of condensation and precipitat ion. Stratif ied variations in impuri ty content wi th in the ice mass, thus, reflect temporal fluctuations in the concentration of atmospheric contaminants. A com-binat ion of studies, involving trace element analysis (Boutron and Lor ius, 1979) and conduct iv i ty measurements on melted ice cores from the Greenland Ice Sheet (Ham-mer, 1977; Hammer and others, 1980), have established acids of volcanic origin as the only common impurit ies that display variabi l i ty in concentration on the scale required to produce observed reflections. Violent volcanic eruptions introduce large volumes of H 2 S and S O 2 into the troposphere and stratosphere where subsequent chemical inter-action wi th dominant atmospheric gases produces sulfate aerosols, part icularly sulfuric acid ( H 2 S 0 4 ) . As a result of atmospheric mix ing processes and ensuing "wash-out", global precipitat ion carries an elevated concentration of these volcanic impurit ies. A c -cording to Hammer (1977), precipitat ion deposited on the Greenland Ice Sheet has been found to possess an enhanced electrical conductivity for up to three years after the volcanic event of origin. Comparat ive studies of radio-echo soundings and ice core acidity profiles (Hammer, 1980; M i l l a r , 1981a; 1981b; 1982) have confirmed a suspected cause-effect relationship between horizons of elevated volcanic impuri ty concentration and the radio-stratigraphy. 58 The depositional emplacement of volcanic impurit ies implies an isochronal inter-pretat ion of internal layering. A number of interesting glaciological problems have been approached from this perspective, including ice sheet stabil i ty (Whi l lans, 1976) and absolute dating of polar ice cores (Hammer and others, 1978). Al though for present purposes internal layering at M t . Wrangell is attr ibuted purely to loss tangent var iat ions, prominent horizons might well be due to a combined mechanism involv-ing contrasts in both loss tangent and relative permitt ivi ty. For large-scale eruptive events, it is likely that major surface melting might accompany the deposition of vol-canic acids, causing pronounced density layering in addit ion to increased impurity content. 4.2 P r a c t i c a l E s t i m a t i o n o f R e f l e c t i o n Coe f f i c i en t s Power reflection coefficients are calculated, as a function of range, according to a modif ied form of the "radar equation" presented in Section 3.3 (Equat ion 3.2), P , = GW&AT2 p t e x p ( - ^ C ) e x p ( - 2 ^ D r t ) , (4.6) 647rr where G represents the forward antenna gain, A the wavelength in air, RA the am-pl i tude reflection coefficient for a plain reflecting target, T 2 the power transmission coefficient (taken to be 1.0), Pt the transmitted power and r' = ( r a + r , /n) is the range adjusted for refraction, where n is the refractive index of ice and ra and r t denote the range in air and ice respectively. The recorded power, P'r = P r e x p ( — AC), is the power received at the antenna, Pr, susequently attenuated by internal system losses, C, in decibels. Recal l ing that A = ( l og e 10 ) / l 0 is a convenient conversion factor and D is 59 the propagation loss rate, in d B / m , the power reflection coefficient, (R) = 101og 1 0 R 2A, is given by (R) = (P) - 2 ( G ) - 20 log 10 \XnT2 + (Dn) + (C), (4.7) 87rr' where bracketed variables denote decibel values and P = Pt/P'r is the relative re-ceived power. Relative power is determined by cal ibrat ing the recorded signal over the dynamic range, relative to system performance. Propagat ion losses are modelled, v ia the interpretation model , according to a simple relationship governing a thermally activated loss mechanism, where T represents temperature, R is the universal gas constant and E, the activation energy, and D0 are empirical constants. The constants are determined by fitting the exponential relation to measurements of temperature dependent loss rate at 840 M H z (Westphal , 1963), subject to minimizat ion of squared error (Figure 4.1). Two-way propagation loss, P2, follows by integrating 4.3 C o m p o s i t e R e f l e c t i o n s a n d t h e V o l c a n i c R e c o r d Whi le the radio-stratigraphy recorded at M t . Wrangell is at t r ibuted to past vol-canic act ivi ty, it is emphasized that the reflectivity record is expected to reflect gradual trends, as well as intermittent eruptive events. Due to atmospheric mix ing processes, Antarc t ic and Greenland Ice Sheet stratigraphies record only major volcanic events (4.8) dP2 = 2D. (4.9) dz 60 A B S O R P T I O N L O S S R A T E 5.5 -§ 5.0 § 4.5 4.0 3.5 CD O II CD d CO oo o 3. 2. 2. 30 -20 -10 0 temperature (C) Figure 4.1. Absorption loss per 100 m at 840 MHz as a function of tem-perature. The loss rate is modelled according to a simple relationship for a thermally activated loss mechanism. The model distribution is depicted in relation to actual values extracted from a digitized version of Figure 3.7. 61 of global scale. The M t . Wrangell stratigraphy, on the other hand, could quite con-ceivably reflect minute variations in activity, occurr ing over a matter of days. Whi le strat i f icat ion on such a scale is not resolved by radio-echo sounding, it is important to consider the possible effects of such layering. It has been suggested (Harrison, 1973; Clough,1977; M i l l a r , 1981) that in situations wherein layer thickness and separation are comparable to the spatial pulse length, interference effects give rise to composite reflections wi th power levels exceeding those due to isolated interfaces. The reflec-t ion coefficient for a composite reflection is statistically dependent upon average layer properties over the pulse length. Radio-reflections are caused by contrasts in the complex relative permitt ivi ty, Ae* = A e ' - t c 'A ( tan 6), (4.10) where e' is the relative permit t iv i ty and A e " , A e ' and A(tancS) denote variation in the complex relative permitt iv i ty, relative permit t iv i ty and loss tangent respectively. Harr ison (1973) derived a general expression for the amplitude reflection coefficient due to fluctuation in relative permit t iv i ty, A(tan<5) = 0, and considered three models for depth dependent variat ion, e' = e\ + A e ' ( z ) , where e\ is the mean and A e ' << e^: 1. Isolated thin layers at spacings greater than the pulse length. 2. A random variation with depth, where A e ' is statistically stationary over the depth of the medium (average properties remain constant). 62 3. A random variation with depth, where A e ' is statistically stationary over dis-tances of the same order as the pulse length (average properties smoothly variable on a scale large wi th respect to pulse length). For an isolated thin layer of thickness /, Harr ison obtained an expression previously given by Robin and others (1969); R — \RA\2 — T r 2 / 2 / A e ' \2 (4.11) where R is the power reflection coefficient, RA is the ampli tude reflection coefficient and A m is the wavelength in the medium. The corresponding expression for stationary random variations is R = Ae^ exp (4.12) where pm is the spatial pulse width in the medium, km is the wavenumber in the med ium, angular brackets denote the ensemble average and the characteristic layer spacing, /, is modelled as the correlation length of A e ' , according to a Gaussian auto-correlation function of the form Pe($) = exp I2 + oo / Ae'{z)Ae'(z + $)dz — oo + oo / (Ae>(z))*dz — oo (4.13) 63 Harr ison demonstrated that at any instance, power received at the antenna is the result of a trade-off between isolated and composite effects. According to Equation 4.12, the apparent spacing of reflections caused by composite effects tends to resemble either the pulse w id th , p m , or the actual layer spacing, / , whichever is larger. In addi t ion, it is found that for / << pm the composite reflection coefficient exceeds that for an equivalent isolated contrast by approximately the number of layers within a pulse wid th , pm/l-Of part icular interest here is the final model, which constrains A t ' to be statistically stat ionary on the scale of a pulse length, while allowing for smooth variat ion in average properties on a broader scale. Al though Harrison focussed upon depth dependent th inning of layers caused by densification and deformation, this effect is given only secondary consideration in the present interpretation. Instead, gradual variations in stratif ied layer properties are mainly attr ibuted to periodic trends in volcanic activity at M t . Wrangel l . It is assumed that depth-dependent variabil i ty in complex relative permit t iv i ty is due to gradual fluctuations in the pulse-width-averaged loss tangent. M e a n trends in dielectric absorption are, in turn, attr ibuted to average variations in volcanic impuri ty concentration. As a result, fluctuations in mean complex relative permit t iv i ty are given by A(e* (z ) ) = - » V A ( ( t a n « ) ( * ) ) , (4.14) w i th A ( ( tan«)(z)) 2ireae'f ' (4.15) 64 where eD is the permit t iv i ty of free space, / is the sounding frequency and A (e"(z)), A ((tan S)(z)) and A(oOQ(z)) denote fluctuations in the depth dependent ensemble average of complex relative permit t iv i ty, loss tangent and high frequency conductiv-ity. Relat ive permit t iv i ty, e', is assumed constant and the high frequency conductivity is related to volcanic impuri ty concentration as discussed in Section 3.4. Although variat ions in ionic impuri ty content are assumed to be a direct result of temporal fluc-tuations in volcanic activity, it is necessary to realize that various post-depositional processes influence the glacio-volcanic stratigraphy. Consequently, to investigate the proposed dependence upon volcanic activity, it is useful to examine the variabil ity in reflector strength as a function of depositional age, rather than depth. The relation-ship between depositional age and stratigraphic depth is determined by modell ing the downward transport of surface-deposited precipi tat ion (details are given in a subse-quent section). A l though this vertical coordinate transformation helps to illustrate our interpretat ion, the effects of post-depositional th inning upon the radio-stratigraphy are neither faithful ly nor completely counteracted. A spatially-averaged profile of the power reflection coefficient as a funct ion of depo-sit ional age is depicted in Figure 4.2. Neglecting high frequency fluctuations, a gradual periodic trend in returned power is evident. It is suggested that this observation sup-ports the proposed relationship between depth-dependent variation in complex relative permit t iv i ty and large scale trends in volcanic activity. Intermittent eruptive events, of relatively short durat ion, give rise to isolated layer effects superimposed upon the broad trend in returned power. Because post-deposit ional processes result in a gradual th inning of layers, isolated layer effects are most pronounced wi th in the near surface region. On the other hand, since composite reflection strength exceeds that due to 65 P O W E R - D E P T H D E P O S I T I O N A L A G E P R O F I L E P R O F I L E - 6 0 - 3 0 0 - 6 0 - 3 0 0 P R C ( d B ) P R C ( d B ) Figure 4.2. Spatial ly averaged power-depth profile with corresponding depo-sit ional age profile. Negative values imply depth below surface level and age relative to zero surface age. Posit ive power reflection coefficients are due to approximations inherent in corrections for geometrical spreading and absorp-t ion. 66 an equivalent isolated layer by roughly the number of layers within a pulse width, the average P R C is expected to increase with depth. A simultaneous increase in lateral layer continuity indicates the growing prominence of composite effects. Al though local inhomogeneities can occur in the deposition of volcanic impurit ies, the average con-centration is expected to be laterally homogeneous. Clough (1977) emphasized that, in addit ion to a dependence upon depth-averaged properties, the composite reflection coefficient is a function of pulse width-averaged properties over an individual hori-zon. As a result, in comparison to isolated layer reflections, composite reflections are expected to display relatively greater lateral continuity. M i l l a r (1981a) examined lateral continuity and structure of reflecting horizons recorded over the Greenland Ice Sheet. Apparent radio-echo separation was observed to increase with local thickening of the ice sheet. Where local layer separation exceeded the pulse length, an addit ional echo gradually emerged midway between the original horizons. Internal layering recorded near the southern r im of the caldera exhibits sim-i lar effects where the radio-stratigraphy undergoes pronounced down-warping (Figure 4.3). Two factors could explain the observed warping. F i rs t , by impeding downslope flow below a stagnation depth, roughly coinciding wi th ice thickness at the head of Long Glacier , the southern r im could cause the glacier fabric to buckle. Outflow-ing ice provides the necessary dr iving mechanism by concentrating shear stress below the stagnation depth. Secondly, anomalous volcanic heat flux at the caldera floor could cause a relatively rapid subsidence of glacier strata in the part icular region. It is interesting to note that recent volcanic act ivi ty has concentrated along the op-posite r im , to the north. Consequently, the question of anomalous heat flux below the southern r im hinges upon whether surface activity is related to a relative surplus 67 or depletion of the volcanic heat source. Regardless of the explanation for observed warping, the radio-stratigraphy is strongly influenced by the deformation. Figure 4.4 traces the progression of apparent reflections from a simple f lat- lying configuration, detected wi th in the central region of the caldera, to that observed near the southern r im . Close inspection of layer separation and continuity reveals a complex interplay between neighbouring horizons that appears to confirm the composite nature of deep reflections. Interference effects are analogous to those described by M i l l a r and suggest an average thickening of glacial strata wi th in the synclinal structure. Because apparent layer separation is inconsistent and does not appear to be directly controlled by pulse length, it is possible that the composite reflection mechanism involves the interference of mult iple internal reflections. Al though speculation, along these lines requires further investigation, the foregoing analysis cautions against naive interpretation of deep-lying reflection horizons ih terms of volcanic history. 4.4 Speculation on Volcanic History 4.4.1 Determination of Depositional Age Al though a permanent record of volcanic activity is preserved within the glacio-volcanic stratigraphy at M t . Wrangel l , several post-depositional processes influence the relationship between depositional age and stratigraphic level. Clearly, there ex-ists a fundamental connection between the age of a part icular layer and the rate at which that horizon proceeds vertically downward, relative to the active glacier sur-face. The vertical transport of surface deposited precipitat ion involves densification, NORTH SOUTH 0 1000 2000 3000 APPROXIMATE DISTANCE (m) Figure 4.3. Spatially corrected depth section recorded over traverse Bl. A pronounced down-warping of the radio-stratigraphy is evident near the south-ern end of the traverse. 69 F i g u r e 4.4. Selected interval f rom depth section recorded over traverse B l , depict ing stratigraphic evidence for composite reflections at depth. N l , N2 and N3 denote a sequence of relatively continuous reflecting horizons, as detected w i th in the central region of the caldera. Gradual thickening of deposit ional layers toward the southern r im produces dramatic variat ions in the recorded interference pattern. The same stratigraphic section that gives rise to N l , N2 and N 3 , produces a stratif ied sequence of apparent layers, S i , S2, S3, S4, S5 and S6, at the southern extreme of the interval. The transit ion involves appearance, disappearance, separation and coalescence of apparent layers. 70 lateral spreading and basal melting. From conservation of mass under steady-state conditions, dpjdt = 0, we obtain w^ + p i k k = 0, (4.16) dz where w represents the surface normal component of flow velocity, p denotes depth stratified density and ikk = dVk/dXk — A + dw/dz is the normal strain rate. The horizontal flow divergence, A = du/dx + dv/dy, is determined at the surface, A 0, from ice motion measurements. It is assumed that the two-dimensional divergence remains constant to a stagnation depth, zs, roughly corresponding to ice thickness at the head of Long Glacier. At depths exceeding the stagnation depth, lateral divergence is assumed negligible. The depth rate of change in surface normal transport velocity is, thus, given by ^ = ^ - A , (4.17) dz p dz where = f A 0, 0 \0, z < z < zs The coupled differential equation describing depth dependent densification was pre-sented in section 2.4.1 (Equation 2.12). Depositional age, ta, as a function of surface normal depth, follows by integrating 71 Physical parameters influencing the resulting depositional age distribution include: glacier thickness, horizontal surface flow divergence and rates of surface accumulation and basal melting. Critical sensitivity to the surface flow divergence requires a well controlled estimate obtained by numerical analysis of ice motion measurements. 4.4.2 Analysis of Horizontal Surface Flow Divergence Precise motion surveys were conducted at the surface of the caldera glacier in 1965-66 and 1975-76 (Appendix 1). We describe the distribution of surface flow velocity in terms of third degree polynomial distributions: Vx(x, y) = cxl + cx2x + cx3y + cx4xy + c l 5 x 2 + c l 6 y 2 + cx7x2y + cx8xy2 + cxQx3 + cxXoy 3 (4.19) and Vy(x,y) = c y i + cy2x + Cy3y + Cy4xy + cy5x2 + cyGy2 + cy7x2y + cy8xy2 + cyQx3 + cyl0y3. 4.20) Coefficients are determined by fitting the polynomial surfaces to localized measure-ments of horizontal velocity components, Vx and Vy, subject to minimization of squared error. A second degree polynomial describing horizontal surface flow divergence follows 72 as A 0(x,y) dVx dVy dx dy ( c l 2 + Cy3) + ( 2 c z 5 + Figure 4.5. Two-dimensional distributions for surface flow rate and diver-gence, (a) Summit caldera map, indicat ing region over which the distr ibu-tions are calculated, (b) Surface flow speed (range: 0.104 cm/day to 18.589 cm/day ; contour interval: 0.528 cm/day) , (c) Surface flow divergence (range: - 4 . 9 6 2 x l O - 2 / y r to 1 . 7 3 8 x l 0 " 2 / y r ; contour interval: 1 . 9 1 4 x l O " 3 / y r ) . 74 pulse rise-time. Al though the predicted sequence is found to vary significantly with lo-cality, the record presented here is relatively well constrained. The part icular profile is determined by averaging results obtained over the central region of the caldera, along traverse A l . By design, traverse A l coincides with motion survey transect C — C (Appendix 1) so that flow dynamics are well determined. The prel iminary volcanic record extends to approximately A D 1650. If broad-scale trends in returned power are attr ibutable to gradual periodic fluctuations in volcanic activity, the sequence suggests a periodicity of roughly 100 years. Assuming this conjecture to be plausible, a gradual increase in activity is expected at M t . Wrangel l , w i th a cl imax around A D 2025. Whi le uncertainty in short term events is substantial, it is interesting to compare the predicted near surface record wi th the known eruptive history (Table 4.1). F inal ly , although the glacio-volcanic record preserved wi th in the caldera is pr imari ly attr ibuted to fluctuations in volcanic activity at M t . Wrangel l , we do not rule out the possible influence of addit ional volcanic events in the Pacific northwest. 7 5 S P E C U L A T I V E V O L C A N I C R E C O R D 2 0 0 0 1 9 5 0 1 9 0 0 ^ 1 8 5 0 CO v2 1 8 0 0 § 1 7 5 0 1 7 0 0 1 6 5 0 1982 e !~-^ 1961 ^ 2 1952 7> 1946 1939 Z> 1929 1916 J 1898 -] 1837 — S 1790 1730 J) 1667 I I \ 1644 1 1 1 1 - 6 0 - 3 0 0 P R C ( d B ) Figure 4.6. Speculative volcanic record as a function of reflected power. Posit ive reflection coefficients are due to approximations inherent in corrections for geometrical spreading and absorption. 76 E V E N T R E M A R K S 1930 A n unspectactacular eruption from the West Crater. 1921 Most spectacular recent eruption. Probably from the northeast flank. 1912 Inferred from newspaper reports in Cordova and Valdez. 1899 Triggered by earthquake at Yakutat Bay. T a b l e 4.1. Known eruption record for M t . Wrangel l , A laska. 77 R E F E R E N C E S C I T E D Ackley, S .F . , and T.E.Ke l iher 1979. Ice sheet internal radio- echo reflections and associated physical property changes with depth, Journal of Geophysical Research, 84, 5675-5680. Anderson, K . R . , and J .E .Gaby 1983. Dynamic waveform matching, Information Sci-ences, 31, 221-242. Benson, C S . 1962. Stratigraphic studies in the snow and firn of the Greenland Ice Sheet, SIP RE Research Report 70, 93pp. Benson, C S . 1982. Glaciological-volcanological research on M t . Wrangel l , A laska, Fi-nal technical report on NSF grant EAR 77-15166, Geophysical Institute, University of Alaska. Benson, C . S . , R .J .Motyka 1978. Glacier-volcano interactions on M t . Wrangel l , Alaska, Annual Report, 1977-78, Geophysical Institute, University of Alaska, 1-25. Bingham, D . K . 1967. Ice motion and heat flow studies on M t . Wrangel l , Alaska, M.S . thesis, University of A laska, 117p. Bout ron . C , and C L o r i u s 1979. Trace elements in Antarct ic snows since 1914, Nature, 277, 551-554. Camp l i n , G . C , J . W . G l e n , and J .G.Paren 1978. Theroret ical models for interpreting the dielectric behaviour of HF-doped ice, Journal of Glaciology, 21, 123-142. Claerbout , J . F . , and F .Mu i r 1973. Robust modell ing wi th erratic data, Geophysics, 35, 826-844. C larke, G . K . C , and C.S.Benson 1983. A n interpretation model for the glacier in the M t . Wrangell caldera, A laska, Internal report, Department of Geophysics and Astronomy, University of Br i t i sh Columbia , Vancouver, Canada, 19p. C larke, G . K . C , G .M.Cross , and C.S.Benson 1987. Ai rborne U H F radar measure-ments of caldera geometry and volcanic history, M t . Wrangel l , A laska, Annals of Glaciology, In press. C lough, J . W . 1977. Radio-echo sounding: reflections from internal layers in ice sheets, Journal of Glaciology, 18, 3-14. Evans, J . R . 1982. Running median filters and a general despiker, Bulletin of the Seismological Society of America, 72, 331-338. Evans, S. 1965. Dielectric properties of ice and snow: a review, Journal of Glaciology, 5, 773-792. Evans. S., and B .M.E .Smi th 1969. A radio-echo equipment for depth sounding in polar ice sheets, Journal of Scientific Instruments (Journal of Physics E), Series 2(2), 131-136. 78 F i t c h , E . P . , E .J .Coy le , and N.C.Gal lagher 1984. Root properties and convergence rates of median filters, IEEE Transactions on Acoustics, Speech and Signal Processing, 5, 739-746. Gudmandsen, P. 1975. Layer echoes in polar ice sheets, Journal of Glaciology, 15, 95-101. Hammer , C U . 1977. Past volcanism revealed by Greenland Ice Sheet impuit ies, Na -ture, 270, 482-486. Hammer , C U . 1980. Acid i ty of polar ice cores in relation to absolute dat ing, pat volcanism, and radio-echoes, Journal of Glaciology, 25, 359-372. Hammer , C . U . , H.B.Clausen, and W.Dansgaard 1980. Greenland Ice Sheet evidence of post-glacial volcanism and its cl imatic impact, Nature, 288, 230-235. Hammer , C . U . , H.B.Clausen, W.Dansgaard, N.Gundestrup, S.J.Johnsen, and N.Reeh 1978. Dat ing Greenland ice cores by flow models, isotopes, volcanic debris, and continental dust, Journal of Glaciology, 20, 3-26. Harr ison, C . H . 1971. Radio-echo sounding: focussing effects in wavy strata, Geophys-ical Journal of the Royal Astronomical Society, 24, 383-400. Harr ison, C . H . 1972. Radio propagation effects in glaciers, P h . D . thesis, University of Cambr idge, 193p. Harr ison, C . H . 1973. Radio-echo sounding of horizontal layers in ice, Journal of Glaciology, 12, 383-397. M a c K e i t h P., Unpubl ished field notes for 1976-1978 seasons at M t . Wrangell , 1978. Mo tyka , R . J . 1983. Volcanological research on M t . Wrangell caldera, Alaska ut i l iz ing the natural glacier ice as a calorimeter, P h . D . thesis, University of Alaska. Mo tyka , R . J . , P .MacKe i th , and C S . B e n s o n 1980. M t . Wrangell: ut i l izat ion of glacier ice to measure heat flow and infer thermal regime, EOS, 61, 66. M i l l a r , D . H . M . 1981a. Radio-echo layering in polar ice sheets, P h . D . thesis, University of Cambr idge, 177p. M i l l a r , D . H . M . 1981b. Radio-echo layering in polar ice sheets and past volcanic activ-ity, Nature, 292, 441-443. M i l l a r , D . H . M . 1982. Acid i ty layers in ice sheets from radio-echo sounding, Annals of Glaciology, 3, 199-203. Myers , C .S . , and L.R.Rabiner 1981. A comparative study of several dynamic time-warping algoritms for connected-word recognition, Bell System Technology Jour-nal, 60, 1389-1409. Na rod , B . B . , and G . K . C . C l a r k e 1980. Ai rborne U H F radio-echo sounding of three Yukon glaciers, Journal of Glaciology, 25, 23-31. 79 Narod , B . B . , and G . K . C . C l a r k e 1983. U H F radar system for airborne ice thickness surveys, Canadian Journal of Earth Sciences, 20, 1073-1086. Nea l , C S . 1976. Radio-echo power profi l ing, Journal of Glaciology,' 17, 527-530. Nea l , C S . 1977. Radio-echo studies of the Ross Ice Shelf, P h . D . thesis, University of Cambr idge. Neidel l , N.S. , and M.T.Taner 1971. Semblance and other coherency measures for mult ichannel data, Geophysics, 36, 482-497. Nodes, T . A . , and N.C.Gal lagher 1982. Med ian filters: some modifications and their properties, IEEE Transactions on Acoustics, Speech and Signal Processing, 5, 739-746. Oswald, G . K . A . 1975. Radio-echo sounding of polar glacier beds, P h . D . thesis, Un i -versity of Cambridge, 134p. Paren , J . G . 1970. Dielectric properties of ice, P h . D . thesis, University of Cambridge. Paren , J . G . 1981. Reflection coefficient at a dielectric interface, correspondence, Jour-nal of Glaciology, 27, 203-204. Paren , J . G . , and G.de Q.Rob in 1975. Internal reflections in polar ice sheets, Journal of Glaciology, 14, 251-259. Prager, B .T . 1993. Digi ta l signal procesing of U H F radio-echo sounding data from Ellesmere Island, M.Sc. thesis, University of Br i t ish Columbia, 88p. Press, F. , and R.Siever 1974. Earth, W . H . Freeman and Company, San Francisco, U .S .A . , 945p. R o b i n , G.de Q. , S.Evans, and J .T.Ba i ley 1969. Interpretation of radio-echo sounding in polar ice sheets, Philosophical Transactions of the Royal Society, Series A, 265, 437-505. Rob inson, J . C . 1970. Statist ical ly optimal stacking of seismic data, Geophysics, 35, 436-446. Sheriff, R . E . , and L.P.Geldhar t 1983. Exploration Seismology: Data Processing and Interpretation, Cambridge University Press, Cambridge, U .K . , 221p. Smi th , B . M . E . 1971. Radio-echo studies of glaciers, P h . D . thesis, University of C a m -bridge, 157p. Smi th , B . M . E . , and S.Evans 1972. Radio-echo sounding: absorption and scattering by water inclusions and ice lenses, Journal of Glaciology, 11, 133-146. Wal ford, M . E . R . 1968. F ie ld measurements of dielectric absorption in antarctic ice and snow at very high frequencies, Journal of Glaciology, 7, 89-94. Westphal , J . A . 1963. Personal communication in Evans, 1965. 80 Whi l lans , I .M. 1976. Radio-echo layers and the recent stabi l i ty of the West Antarct ic Ice Sheet, Nature, 264 , 152-155. GENERAL REFERENCES Bogordsky, V . V . , C.R.Bent ley, and P.E.Gudmandsen 1985. Radioglaciology, D. Reidel Publ ish ing Company, Dordrecht, Hol land, 254p. Paterson, W . S . B . 1981. The Physics of Glaciers, Pergamon Press, Oxford , U .K . , 380p. Skolnik, M.I . 1962. Introduction to Radar Systems, M c G r a w - H i l l Book Company, New York , U .S .A . , 648p. APPENDIX 81 GLACIER MOTION MEASUREMENTS Figures A . l and A .2 illustrate stake positions surveyed for ice motion measure-ments in 1965-66 (Bingham, 1967) and 1975-76 (Motyka, 1983) respectively. Corre-sponding data uti l ized for the present analysis of horizontal surface flow divergence are tabulated in Tables A . l and A .2 . 82 Figure A.l. M a p of the summit caldera, depicting original location and horizontal velocity vectors for ice motion stakes surveyed over 1965-66 (from B ingham, 1967). 83 S T A K E E A S T I N G (m) N O R T H I N G (m) I N T E R V A L (days) S3 1322.47 -1207.24 1322.39 -1207.72 56.15 S5 2485.58 -434.72 2485.91 -436.40 56.13 S6 1911.51 -823.54 1913.06 -825.12 56.13 S7 1236.02 -2499.05 1236.38 -2498.26 27.58 S9 1038.03 -2249.52 1039.40 -2248.28 56.15 S16 1037.47 -1618.17 1039.01 -1617.91 56.13 S19 1367.99 -1041.48 1369.70 -1042.26 56.13 S20 1771.40 -692.40 1772.90 -693.81 56.13 S43 1574.26 -1411.84 1575.10 -1411.95 23.58 S44 1979.83 -1166.60 1980.61 -1167.13 23.58 S45 2399.31 -912.32 2399.69 -913.35 23.58 S46 2853.10 -637.20 2853.07 -638.25 23.58 T a b l e A . l . Horizontal motion components; 1965-66 survey (from Bingham, 1967). 84 Figure A.2. Map of the summit caldera, depicting original location and horizontal velocity vectors for ice motion stakes surveyed over 1975-76 (from Motyka, 1983). 85 S T A K E E A S T I N G (m) NORTHING (m) INTERVAL (days) 23 1592.25 266.95 1592.71 264.43 307.0 24 1909.11 -896.43 1911.34 -898.12 71.5 24 aa 1696.84 -1030.71 1699.04 -1032.29 71.5 24bb 2046.99 -1107.67 2049.58 -1109.77 71.5 25 2333.09 -629.14 2334.28 -630.80 71.5 26 2760.73 -359.79 2761.24 -361.32 71.5 27 1483.63 -1163.97 1485.94 -1165.28 71.5 28 1059.63 -1431.80 1060.98 -1431.59 71.5 28a 846.19 -1565.89 847.87 -1565.73 71.5 29 633.21 -1699.71 634.52 -1699.68 71.5 30 1633.96 -475.57 1635.28 -475.67 71.5 31 1358.78 -56.26 1359.99 -57.17 71.5 32 2185.75 -1320.26 2188.93 -1322.96 71.5 33 2463.37 -1745.85 2468.95 -1750.36 71.5 33aa 2252.93 -1885.48 2258.45 -1889.18 73.0 34 2884.00 -1465.74 2887.36 -1470.17 73.0 35 3307.68 -1184.14 3308.72 -1187.47 73.0 * (continued on following page) 86 S T A K E E A S T I N G (m) N O R T H I N G (m) I N T E R V A L (days) 36 2042.99 -2024.97 2047.73 -2027.25 83.0 36a 1831.60 -2165.08 1835.71 -2166.38 83.0 37 1619.22 -2306.82 1622.65 -2306.70 83.0 Table A.2. Horizotal mot ion components; 1975-76 survey (from Motyka , 1983).