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 # +