REMOTE SENSING, MORPHOLOGIC ANALYSIS, AND ANALOGUE MODELING OF LAVA CHANNEL NETWORKS IN HAWAI‘I by HANNAH ROSE DIETTERICH A DISSERTATION Presented to the Department of Geological Sciences and the Graduate School of the University of Oregon in partial fulfillment of the requirements for the degree of Doctor of Philosophy June 2014 ii DISSERTATION APPROVAL PAGE Student: Hannah Rose Dietterich Title: Remote Sensing, Morphologic Analysis, and Analogue Modeling of Lava Channel Networks in Hawai‘i This dissertation has been accepted and approved in partial fulfillment of the requirements for the Doctor of Philosophy degree in the Department of Geological Sciences by: Katharine V. Cashman Chairperson Alan W. Rempel Core Member David A. Schmidt Core Member Mark A. Fonstad Institutional Representative and Kimberly Andrews Espy Vice President for Research and Innovation; Dean of the Graduate School Original approval signatures are on file with the University of Oregon Graduate School. Degree awarded June 2014 iii © 2014 Hannah Rose Dietterich iv DISSERTATION ABSTRACT Hannah Rose Dietterich Doctor of Philosophy Department of Geological Sciences June 2014 Title: Remote Sensing, Morphologic Analysis, and Analogue Modeling of Lava Channel Networks in Hawai‘i Lava flows are common at volcanoes around the world and on other terrestrial planets, but their behavior is not fully understood. In Hawai‘i, advances in remote sensing are offering new insights into lava flow emplacement. In this dissertation, I develop new techniques using satellite-based synthetic aperature radar, aerial photographs, and airborne lidar to produce three-dimensional high-resolution maps of lava flows from data collected before, during, and after emplacement. These new datasets highlight complex lava channel networks within these flows, which are not incorporated into current predictive or probabilistic lava flow models yet may affect flow behavior. I investigate the origin and influence of these channel networks through morphologic analysis of underlying topography, network topology, and flow morphology and volume. Channel network geometries range from distributary systems dominated by flow branching around local obstacles to tributary systems constricted by topography. I find that flow branching occurs where the flow thins over steeper slopes and that the degree of flow branching, network connectivity, and longevity of flow segments all influence the final flow morphology. Furthermore, because channel networks govern the distribution of lava supply within a flow, changes in the channel topology can dramatically alter the effective v volumetric flux in any one branch, which affects both flow length and advance rate. Specifically, branching will slow and shorten flows, while merging can accelerate and lengthen them. To test these observations from historic eruptions and morphologic analysis, I use analogue experiments to simulate the interaction of a lava flow with a topographic obstacle and determine the conditions under which the flow branches and the effects of the bifurcation on flow advance rate. These experiments support the earlier results but also demonstrate the importance of flow dynamics and obstacle morphology on governing when flows may overtop obstacles. Consideration of channel networks is thus important for predicting lava flow behavior and mitigating flow hazards with diversion barriers. One video of Kīlauea lava flow activity from 2003–2010 accompanies this dissertation as a supplemental file. This dissertation includes both previously published and unpublished co-authored material. vi CURRICULUM VITAE NAME OF AUTHOR: Hannah Rose Dietterich GRADUATE AND UNDERGRADUATE SCHOOLS ATTENDED: University of Oregon, Eugene Pomona College, Claremont, CA DEGREES AWARDED: Doctor of Philosophy, Geological Sciences, 2014, University of Oregon Bachelor of Arts, Geology, 2009, Pomona College AREAS OF SPECIAL INTEREST: Physical Volcanology PROFESSIONAL EXPERIENCE: Graduate Teaching and Research Fellow, Department of Geological Sciences, University of Oregon, 2009 - 2010, 2013 - 2014 Student Researcher, Department of Geosciences, Oregon State University, 2007 – 2009 Teaching Assistant, Department of Geology, Pomona College, 2008 GRANTS, AWARDS, AND HONORS: Conference Support, IAVCEI General Assembly, 2013 Graduate Research Fellowship, National Science Foundation, 2010-2013 Staples Fellowship, Department of Geological Sciences, University of Oregon, 2009, 2011, and 2013 Conference Support, American Geophysical Union, 2012 Johnston Scholarship, Department of Geological Sciences, University of Oregon, 2011 and 2012 vii Donald B. McIntyre - H. Stanton Hill Geology Award, Department of Geology, Pomona College, 2009 Elected to Phi Beta Kappa, 2009 Isabel F. Smith and Donald H. Zenger Award, Department of Geology, Pomona College, 2008 Elected to Sigma Xi, 2008 Outstanding Student Paper Award, American Geophysical Union, VGP Section, 2008 PUBLICATIONS: Dietterich, H. R., and K. V. Cashman (in press), Channel networks within lava flows: Formation, evolution, and implications for flow behavior, J. Geophys. Res., doi:10.1002/2014JF003103. Dietterich, H. R., S. A. Soule, K. V. Cashman, and B. H. Mackey (in press), Lava Flows in 3D: Using airborne lidar and pre-eruptive topography to evaluate lava flow surface morphology and thickness in Hawai‘i, in Hawaiian Volcanoes, From Source to Surface, Geophys. Monogr. Ser., AGU, Washington, D. C. Cashman, K. V., S. A. Soule, B. H. Mackey, N. I. Deligne, N. Deardorff, and H. R. Dietterich (2013), How lava flows: New insights from applications of lidar technologies to lava flow studies, Geosphere, 9, 1664-1680, doi: 10.1130/GES00706.1. Dietterich, H. R., M. P. Poland, D. A. Schmidt, K. V. Cashman, D. R. Sherrod, and A. T. Espinosa (2012), Tracking lava flow emplacement on the east rift zone of Kīlauea, Hawai‘i, with synthetic aperture radar coherence, Geochem. Geophy. Geosy., 13, Q05001, doi:10.1029/2011GC004016. Dietterich, H.R., S. L. de Silva, and G. Salas (2010), Sulfur yield of the 1600 eruption of Huaynaputina, Peru: Contributions from magmatic, fluid-phase, and hydrothermal sulfur, J. Volcanol. Geoth. Res., 197, 303-312, doi:10.1016/j.jvolgeores.2010.01.003. Peterson, M. G., H. R. Dietterich, B. Lachenbruch (2007), Do Douglas-fir branches and roots have juvenile wood? Wood Fiber Sci., 39, 651-660. viii ACKNOWLEDGMENTS This dissertation has benefited from the contributions of many people. I would first like to thank my advisor, Kathy Cashman, for offering guidance and support throughout my PhD. She has been an amazing mentor, and her scientific curiosity, creativity, enthusiasm, and generosity have inspired me as I begin my own scientific career. I would also like to thank my committee, including David Schmidt, Alan Rempel, and Mark Fonstad, for their feedback, advice, and contributions to this dissertation. My dissertation is a product of much collaboration, and includes the contributions of many co-authors apart from Kathy and myself. Chapter II is co-authored by Mike Poland, David Schmidt, Dave Sherrod, and Arkin Tapia Espinosa, Chapter III is co- authored by Adam Soule and Ben Mackey, and Chapter V is co-authored with Alison Rust. However, co-authors do not represent all of the people that have contributed to this dissertation. The USGS Hawaiian Volcano Observatory staff has shared a tremendous amount of data and knowledge, and offered feedback on all of this work. University of Bristol students, faculty, and staff have assisted with my experiments, rheology measurements, and equipment fabrication. Although molten basalt experiments largely did not make it into this dissertation, I have to thank Jeff Karson and Bob Wysocki at the Syracuse University Lava Project, as well as Einat Lev. Mark Fonstad and James Dietrich have contributed equipment, expertise, and Structure from Motion data processing to both the molten basalt experiments and field 3D mapping in Hawai‘i. Three of these chapters have undergone peer-review and been improved by many thoughtful and helpful reviews from journal editors and reviewers, including Scott Rowland, Matt Patrick, Jim Kauahikaua, and an anonymous reviewer for Chapter II, ix Rebecca Carey, Francesco Mazzarini, and an anonymous reviewer for Chapter III, and Alexander Densmore, Jon Pelletier, Christopher Hamilton, Matthew Patrick, and an anonymous reviewer for Chapter IV. Members of the Cashman lab, both at University of Oregon and University of Bristol, have been a core group of friends and colleagues for which I am incredibly grateful. At University of Oregon we have a fantastic community of graduate students who will surely be lifelong friends and colleagues. The faculty and office staff have also been incredible and always available for help and conversation. At University of Bristol I’m especially grateful to the postgraduate students, faculty, and staff for welcoming me and welcoming me back each time I’ve visited. Finally, I thank my family and friends for their support over the years. My parents introduced me to the wonders of the natural world and to science as a pursuit and career. Thanks also to my brother Noah for entertaining me, challenging me, and always being there. Special thanks, as well, to my extended family and friends for putting up with all of my geology chatter while being so encouraging. This work was largely funded by the National Science Foundation through a Graduate Research Fellowship under grant DGE-0829517, along with research grants awarded to Kathy Cashman, NSF EAR 0738894 and 1250554. The AXA Research Fund and a Royal Society University Research Fellowship for Alison Rust funded work performed while at the University of Bristol. The University of Oregon Department of Geological Sciences and the Graduate School supplied travel grants and equipment awards. I have also received travel grants from the American Geophysical Union and the International Association of Volcanology and Chemistry of the Earth’s Interior. x For my family xi TABLE OF CONTENTS Chapter Page I. INTRODUCTION ................................................................................................. 1 II. TRACKING LAVA FLOW EMPLACEMENT ON THE EAST RIFT ZONE OF KĪLAUEA, HAWAI‘I, WITH SYNTHETIC APERTURE RADAR COHERENCE........................................................................................................ 7 1. Introduction........................................................................................................ 7 2. Methodology ...................................................................................................... 12 2.1. SAR Data Sources and Processing ........................................................... 12 2.2. Image Analysis.......................................................................................... 14 2.2.1. Vegetation Removal......................................................................... 15 2.2.2. Image Thresholding ......................................................................... 18 2.2.3. Time Series Assembly ..................................................................... 19 2.2.4. Flow Identification........................................................................... 20 3. Results................................................................................................................ 21 4. Discussion .......................................................................................................... 25 4.1. Evaluation of Results ................................................................................ 25 4.1.1. Errors in SAR-Based Mapping ........................................................ 26 4.1.2. New Views From SAR-Based Mapping.......................................... 28 4.2. Analysis and Applications ........................................................................ 32 5. Conclusions........................................................................................................ 38 6. Bridge................................................................................................................. 40 III. LAVA FLOWS IN 3D: USING AIRBORNE LIDAR AND PRE-ERUPTIVE TOPOGRAPHY TO EVALUATE LAVA FLOW SURFACE MORPHOLOGY AND THICKNESS IN HAWAI‘I ............................................ 41 xii Chapter Page 1. Introduction........................................................................................................ 41 1.1. Characterizing Lava Flow Surfaces .......................................................... 42 1.2. Characterizing Lava Flow Morphology.................................................... 44 1.3. A Brief Review of the Target Lava Flows................................................ 45 2. Methods.............................................................................................................. 47 2.1. Data Sources ............................................................................................. 47 2.1.1. Lidar................................................................................................. 48 2.1.2. Photogrammetry............................................................................... 50 2.2. Derived Datasets ....................................................................................... 52 2.2.1. Lidar Surface Properties .................................................................. 52 2.2.2. Surface Properties Classification ..................................................... 56 2.2.3. Flow Thickness ................................................................................ 57 3. Results and Discussion ...................................................................................... 57 3.1. Intensity Analysis...................................................................................... 58 3.2. Surface Roughness Analysis..................................................................... 61 3.3. Classification of Lidar Surfaces................................................................ 64 3.4. Flow Thickness Analysis .......................................................................... 67 3.4.1. Flow Volume ................................................................................... 68 3.4.2. Thickness of Lava Flow Features .................................................... 73 3.4.3. Influence of Pre-Eruptive Topography on Flow Thickness............. 77 4. Conclusions........................................................................................................ 80 5. Bridge................................................................................................................. 82 xiii Chapter Page IV. CHANNEL NETWORKS WITHIN LAVA FLOWS: FORMATION, EVOLUTION, AND IMPLICATIONS FOR FLOW BEHAVIOR...................... 84 1. Introduction........................................................................................................ 84 2. Background ........................................................................................................ 87 2.1. Controls on Lava Flow Emplacement....................................................... 89 2.2. Morphologic Observations of Lava Flows ............................................... 90 3. Methods.............................................................................................................. 92 3.1. Quantifying Channel Networks ................................................................ 94 3.1.1. Braiding Index ................................................................................. 94 3.1.2. Network Measures ........................................................................... 95 3.1.2.1. Ordering .................................................................................. 96 3.1.2.2. Brain Connectivity Toolbox ................................................... 99 3.2. Quantifying Flow and Channel Morphology............................................ 100 3.2.1. Morphology at Flow Cross-Sections ............................................... 101 3.2.2. Morphology of Channel Segments .................................................. 101 4. Results................................................................................................................ 103 5. Discussion .......................................................................................................... 106 5.1. The Formation of Channel Networks ....................................................... 107 5.2. The Influence of the Channel Network on Flow Morphology and Evolution .................................................................................................. 110 5.3. The Influence of the Channel Network on Flow Emplacement Behavior ................................................................................................... 115 5.3.1. Influence on Lava Flow Length....................................................... 115 xiv Chapter Page 5.3.2. Influence on Advance Rate.............................................................. 117 6. Implications........................................................................................................ 118 6.1. Flow Prediction......................................................................................... 119 6.2. Flow Mitigation ........................................................................................ 121 6.3. Wider Implications of Channel Networks for Lava Flow Emplacement ............................................................................................ 124 7. Bridge................................................................................................................. 126 V. EXPERIMENTAL MODELING OF LAVA FLOW INTERACTIONS WITH TOPOGRAPHIC OBSTACLES............................................................................ 128 1. Introduction........................................................................................................ 128 2. Lava Flow Rheology and Theory ...................................................................... 132 3. Methods.............................................................................................................. 135 3.1. Experimental Apparatus and Procedure.................................................... 137 3.2. Experiment Series ..................................................................................... 138 4. Results................................................................................................................ 141 5. Discussion .......................................................................................................... 148 5.1. Viscous Flow without an Obstacle ........................................................... 149 5.2. Analysis of Overtopping Experiments...................................................... 151 5.3. The Influence of Flow Velocity................................................................ 153 5.4. Influence of Obstacle Shape and Size....................................................... 154 5.5. Interpretations from the Stationary Wave Analysis.................................. 157 5.6. Qualitative Summary of Observations and Rationalizations .................... 166 6. Conclusions........................................................................................................ 168 xv Chapter Page VI. CONCLUSIONS ................................................................................................... 171 1. Dissertation Summary........................................................................................ 171 2. Future Research Directions................................................................................ 174 REFERENCES CITED................................................................................................ 180 SUPPLEMENTAL FILES VIDEO: KĪLAUEA FLOW MAPS 2003–2010 xvi LIST OF FIGURES Figure Page CHAPTER II 1. Location map of the study area on the east rift zone of Kīlauea, Hawai‘i............. 13 2. Flow chart of the SCM method.............................................................................. 15 3. Schematic demonstration of how the time series inversion captures flow activity over a single pixel ..................................................................................... 17 4. All 210 flow maps in the time series displayed in four panels over the SRTM DEM....................................................................................................................... 22 5. SCM flow map from January 2003 to January 2010, SCM flow map from July 2003 to October 2010, USGS flow map from 1983 to January 2010, USGS flow map highlighting the 2007 to 2010 flow field and the March to October 2010 flow ........................................................................................................................ 23 6. SCM map of the July 2007 flow extent as of October 17, 2007 compared to the ALOS coherence image of the July 2007 flow field from July 16 to October 16, 2007.................................................................................................... 27 7. SCM flow map within the boundary of the USGS map ........................................ 29 8. View of the pali and coastal plain in a coherence image shows breakouts over the pali, the path of a lava tube .............................................................................. 30 9. Oblique view of the SCM flow map and a FLIR image ........................................ 31 10. Extent of the Thanksgiving Eve Breakout flow displayed with time .................... 33 11. The length of time each pixel is decorrelated after initial emplacement matches well with surveyed flow thickness contours .......................................................... 34 12. Comparison of flow thicknesses and the average time pixels are decorrelated within each thickness contour................................................................................ 36 13. Envisat interferogram spanning October 21 to November 25, 2009 ..................... 39 xvii Figure Page CHAPTER III 1. Location maps of the study area ............................................................................ 46 2. Lidar intensity map with lidar hillshade overlay ................................................... 49 3. Assembly of the pre-eruptive DEM and thickness map ........................................ 51 4. Demonstration of errors and corrections of the pre-eruptive DEM....................... 53 5. Example lidar surface properties maps .................................................................. 54 6. Map of normalized lidar intensities on Mauna Loa NERZ with mapped geologic units outlined.......................................................................................................... 59 7. Intensity increases with decreasing elevation within individual flows.................. 60 8. Age-intensity relationships for lava flows on the NERZ of Mauna Loa ............... 62 9. Derived roughness map and down-flow roughness for the December 1974 lava flow on the SW Rift Zone of Kīlauea .................................................................... 64 10. All-layer classification results for eight classes, four classes, and two classes ..... 66 11. Surface roughness ratio classification results for eight classes ............................. 68 12. Flow thickness map with lidar hillshade of the Mauna Loa 1984 flow................. 69 13. Volume distribution among flow lobes.................................................................. 71 14. Distribution of volume with distance along the flow............................................. 74 15. Flow thickness maps of the stable channel of Flow 1 and the distal portion of Flow 2 .................................................................................................................... 75 16. Decrease in flux down the main channel of Flow 1/1A with distance from the vent as recorded on April 3–4, 1984...................................................................... 77 17. Comparison of lava flow thickness to underlying slope at a variety of length scales ...................................................................................................................... 79 xviii Figure Page CHAPTER IV 1. Locations, lidar hillshade maps, and traced channel networks for the Mauna Loa 1984 flow and the Kīlauea December 1974 flow ........................................... 86 2. Flow thickness map of the ML84 flow.................................................................. 92 3. Map of the channel network of the ML84 flow with cross-section lines color-coded by braiding index ............................................................................... 95 4. Schematic of the different channel ordering schemes used ................................... 98 5. Map of estimated flux order and Shreve order for the KDec74 flow .................... 99 6. Map of the normalized betweenness centrality measure with classified nodes..... 100 7. Example of the channel segment morphology measurements............................... 102 8. Comparison of how the different network measures change with distance along flow for distributary and tributary flows................................................................ 105 9. Distance along flow shows a modest correlation with channel width and little relationship with average levee height................................................................... 106 10. Braiding index of Flow 4 of the ML84 flow overlain on underlying slope........... 108 11. Estimated flux order and Shreve order show relationships with channel width.... 111 12. Average levee height and sinuosity compared to Shreve order............................. 112 13. Box and whisker plot of Shreve order correlating with channel longevity ........... 113 14. Longevity also shows a relationship with flow morphology ................................. 114 15. Plot of effusion rate versus flow length for channelized lava flows erupted at Mauna Loa and Kīlauea ......................................................................................... 116 16. Map and advance data from the Pu‘u ‘Ō‘ō-Kupaianaha eruption of Kīlauea ....... 119 xix Figure Page CHAPTER V 1. Experimental viscosities ........................................................................................ 136 2. Schematic drawing of the experimental setup ....................................................... 138 3. Obstacle shapes and sizes used in experiments ..................................................... 140 4. Photographs of stationary waves on the upslope side of the obstacle ................... 146 5. Flow height upslope of two 4 cm triangular obstacles plotted with time .............. 146 6. Flow height measurements for PTFE large cylinders plotted against flow velocity and plastic triangular obstacles plotted against obstacle angle.............................. 147 7. Comparison of the flow heights of stationary waves behind small and large obstacles................................................................................................................. 147 8. Comparison of the flow height upslope of the obstacle as measured behind obstacles with different surface properties ............................................................ 149 9. Evolution of the down-slope and cross-slope dimensions of an example experimental gravity current .................................................................................. 151 10. Comparison of the flow height calculated and measured ...................................... 152 11. Overtopping experiment results compared to flow heights and flow velocities.... 153 12. Plot of flow velocity versus the wave height ......................................................... 154 13. Plot of obstacle angle compared to wave height divided by average velocity0.69 .............................................................................................................. 156 14. Obstacle angle compared to the difference between the flow heights upslope of a large, 30 cm triangular obstacle and a small, 4 cm one of the same angle at the same flux and slope.......................................................................................... 157 15. Annotated picture of flow into and along a 150° large obstacle............................ 160 16. Along obstacle advance in the down-slope and cross-slope directions compared with obstacle angles ............................................................................................... 161 xx Figure Page 17. Experimental data for flow advance along the obstacle compared to modeled advance .................................................................................................................. 162 18. Plot of the flow height predicted by the Jefferys equation compared to the flow height within the stationary wave .......................................................................... 164 19. Flow advance before, along, and after the obstacle compared flow to modeling ................................................................................................................ 165 CHAPTER VI 1. Oblique view of a structure from motion DEM where the channel divides around an obstacle.................................................................................................. 176 2. Example data from a molten basalt experiment..................................................... 178 xxi LIST OF TABLES Table Page CHAPTER III 1. All layer unsupervised classification interpretation............................................... 67 CHAPTER V 1. Physical properties of golden syrup at room temperature (20°C) ......................... 136 2. Summary of experimental parameters and results ................................................. 142 1 CHAPTER I INTRODUCTION Lava flows are the repaving mechanism of the Earth’s surface. They represent the largest eruptions by volume on Earth in the form of flood basalts and make up Earth’s longest mountain range through eruptions along the mid-ocean ridges that circle the globe. However, they are also a frequent hazard on a smaller scale at volcanoes around the world. Although rarely deadly, lava flows can travel quite quickly and destroy everything in their paths. Scientists and emergency managers have therefore sought to develop the ability to monitor, predict, and divert lava flows to protect populations and property. To do this, much effort has been put in to developing techniques for mapping lava flows as they erupt and spread over the landscape, investigating the important parameters and physics of lava flow behavior, and combining both of these to successfully observe, predict, and mitigate the hazards posed by effusive eruptions. In this dissertation, I seek to contribute to these goals, advancing new remote sensing techniques to collect data from active and recent lava flows in Hawai‘i that I then analyze with the aim of improving our understanding of how underlying topography and channel network geometry influence flow emplacement behavior. Lava flow monitoring in Hawai‘i currently consists of weekly (or more frequently during emergencies) field mapping campaigns combined with satellite imagery of the active Kīlauea Volcano, on the island of Hawai‘i. Through these efforts, the Hawaiian 2 Volcano Observatory (HVO) is able to track the advance of flows toward property and infrastructure while watching the overall rate of lava effusion. However, field mapping, especially via helicopter, is time-consuming and costly, while much of the satellite data used to map active flows is often obscured by clouds and has low temporal and spatial resolution. In Chapter II, I outline a new technique to map both the advance and internal activity of Hawaiian lava flows using satellite Synthetic Aperture Radar (SAR) data. This SAR imagery is already collected regularly to measure ground deformation at volcanoes with interferometric SAR (InSAR), a method in which before and after images are compared to measure centimeter-scale changes in the surface elevation. However, these before and after images can also serve to detect lava flow emplacement, based on where the surface has been radically altered (i.e., repaved) between the two passes of the satellite, which I measure with SAR coherence. I use a time series of these data from 2003–2010 to develop an algorithm to process SAR coherence images into flow maps. These maps can then be used to make measurements of flow advance rates, and even estimates of flow volumes. The technique offers new insights from new views of the active flow field at Kīlauea and has now been put into service to aid in flow mapping at HVO. Remote sensing can also supply high-resolution topographic data that, even without the repeated deployment required for flow monitoring, provides unprecedented views of lava flow morphology. In Chapter III, I use aerial photographs and airborne lidar (light detection and ranging, a form of laser altimetry) data to build digital elevation models (DEMs) of the pre- and post-eruptive ground surface of the large lava flow emplaced during the 1984 eruption of Mauna Loa Volcano. From the lidar data it is 3 possible to asses the surface morphology of flow, identifying channels, levees, and the distribution of ‘a‘ā and pāhoehoe within a flow. The lidar intensity (amplitude, similar to an infrared image) also offers information as to the age of the flow surface and the amount of revegetation after emplacement. From the difference between the pre- and post-eruptive surfaces, I can also calculate flow thickness, and create the most detailed flow thickness map yet produced for a Hawaiian lava flow. Analysis of these thickness data yields valuable information about the volume distribution within the flow, with implications for the influence of underlying topography and flow branching. Equipped with these datasets, I undertake an investigation in Chapter IV into how underlying topography may form channel networks, and how those networks in turn affect flow morphology and behavior. For this purpose, I build maps of underlying slope to compare with channel network paths. I find that there are greater numbers of parallel channels on steeper slopes, suggesting that flows must split around obstacles as they thin over increased slopes. These channel networks can also be analyzed using metrics from rivers and neural networks and compared to the channel morphology captured in the lidar. This examination reveals that the most connected, important channels are the longest-lived and narrowest, with overall higher levees. These observations support existing models for the evolution of channels through time. Analysis of a set of historical eruptions suggests further implications of branching and merging channels for flow behavior. My examination shows that branched flows advance slower and have shorter final lengths than flows that are confined by topography. This implies that branching and confinement, and the topography that controls them, need to be consider in models of flow prediction. 4 I can test this field evidence for the influence of underlying topography in Chapter V with analogue experiments that send flows into obstacles. Observations from the previous chapters offer hypotheses as to when a flow must split around an obstacle and thus form a flow branch, or simply overtop it. In Chapter IV, the thickness of the flow appears to be related to whether it can overtop obstacles, while in Chapters III and IV, velocity controlled lava run-up onto obstacles is apparent. I perform a series of experiments with golden syrup, a common analogue for lava in the literature, flowing down an inclined plane where it intersects an obstacle. By varying the flux of syrup and the underlying slope I can control the flow thickness and velocity, and by changing the obstacle I can modify its size and shape. Surface tension mostly prevents thicker flows from overtopping obstacles, but flows are able to easily overtop at high velocity. This is because the syrup hitting the obstacle gets stuck upslope of it and forms a thick stationary wave that sends it over the obstacle top. Obstacles that are larger or are V-shaped with a greater internal angle also increase the height of this wave. The interaction between the flow and the obstacle is thus very complex, but can be described qualitatively as being controlled by the amount of syrup impacting the obstacle, and how easy it is for that syrup to get past the obstacle. The results of these experiments have implications for the scale of topography that can cause flow branching or topographic confinement. This also applies to interactions with the artificial topography used to divert lava flows, so this work further informs the size and shape of diversion barriers necessary to successfully mitigate a lava flow hazard. Collectively, these chapters describe new remote sensing and morphologic analysis techniques that yield new observations and insights into lava flow behavior, 5 which I then test with analogue modeling. This dissertation therefore contributes to the goals of advancing our ability to monitor active lava flows and improving our understanding of how pre-eruptive topography and channel network geometry influence lava flow emplacement, with direct implications for flow prediction and mitigation by diversion barriers. All of these chapters are co-authored by my advisor, Kathy Cashman, who offered guidance and editorial assistance throughout all of this work. Chapter II has four other co-authors, Michael Poland (U.S. Geological Survey, Hawaiian Volcano Observatory) processed all of the SAR imagery, David Schmidt (University of Oregon, now at University of Washington) co-advised the project, and David Sherrod (U.S. Geological Survey, Cascades Volcano Observatory) and Arkin Tapia Espinosa (University of Panama) shared field mapping data. Chapter III also represents a collaboration, and is co- authored by Adam Soule (Woods Hole Oceanographic Institute), who contributed his own work on lidar surface analysis, and Benjamin Mackey (University of Canterbury), who processed aerial photographs to make the pre-eruptive DEM. Chapter V greatly benefitted from the involvement of co-author Alison Rust (University of Bristol), who co-advised the project and helped set up the analogue experiments. Each chapter has reached a different stage of the publication process. Chapter II was published in Geochemistry, Geophysics, Geosystems in 2012. Chapter III is currently in press as a chapter in an American Geophysical Union Geophysical Monograph entitled, Hawaiian Volcanoes, From Source to Surface, with an expected publication date in 2014. Chapter IV is currently in press in the Journal of Geophysical Research: Earth 6 Surface, and Chapter V is in preparation for the Journal of Geophysical Research: Solid Earth. 7 CHAPTER II TRACKING LAVA FLOW EMPLACEMENT ON THE EAST RIFT ZONE OF KĪLAUEA, HAWAI‘I, WITH SYNTHETIC APERTURE RADAR COHERENCE This chapter was published in Geochemistry, Geophysics, Geosystems in May, 2012. I was the lead author on the paper, doing the method development, analysis, and writing. The second author, Michael Poland (U.S. Geological Survey, Hawaiian Volcano Observatory), processed all of the satellite imagery for the work. My co-authors Kathy Cashman (University of Oregon) and David Schmidt (now at University of Washington) served as advisors for this project, helping with method development, as well as manuscript editing. The final two authors, David Sherrod (U.S. Geological Survey, Cascades Volcano Observatory) and Arkin Tapia Espinosa (University of Panama) provided mapping data that served to ground truth our results. 1. Introduction Mapping of lava flows from the Pu‘u ‘Ō‘ō-Kupaianaha eruption on the east rift zone of Kīlauea serves to document the ongoing eruption, while yielding insights into how lava flow fields develop [e.g., Mattox et al., 1993; Kauahikaua et al., 2003]. However, the locations, extents, and timing of flows are challenging to record with traditional field methods because of the large size and inaccessibility of the affected area. Remote sensing offers a means of gathering this information that circumvents these problems. Interferometric synthetic aperture radar (InSAR) is widely used to measure 8 deformation by detecting minute changes in ground surface elevations at points that stay correlated during repeat observations [Massonnet and Feigl, 1998 and references therein]. Correlation of the radar echoes between two SAR acquisitions requires that the radar-reflecting surfaces within the radar footprint remain unchanged, permitting measurements of the line-of-sight change in range with InSAR. The eruption and emplacement of lava, however, changes the scattering properties of the surface, thereby disrupting the coherence of the radar echoes and precluding the application of interferometric techniques. Nevertheless, the abrupt change in surface properties where the ground is repaved by lava allows these new flows to be mapped with SAR coherence images [Zebker et al., 1996]. We extend the application of coherence mapping to extract time-dependent flow maps from SAR coherence images. Our method, which we term SAR Coherence Mapping (SCM), is independent of look angle or satellite path, and it therefore allows synthesis of all available SAR data in a region to maximize the temporal sampling. We demonstrate the technique by examining the advance of lava flows in Hawai‘i from 2003 to 2010, and explore ways that these data can be used to investigate lava flow behavior during and after emplacement. At many volcanoes around the world, lava flows are recurring natural hazards that impact populations and property. Mitigation strategies include both probabilistic approaches to hazard zone mapping [e.g., Wadge et al., 1994; Kauahikaua et al., 1995] and the development of models to address both conditions of lava flow advance [Crisci et al., 1986; Harris and Rowland, 2001] and terrain analysis [Kauahikaua et al., 2003]; the latter allow real-time decision-making in response to evolving eruptions. Observations of how flow fields develop over time are therefore crucial for tracking the path and velocity 9 of flows, and for building an understanding of lava flow emplacement that is necessary for predicting future flow behavior. Flow mapping is the principal means of gathering information on the advance and morphology of active lava flows. Traditionally, active flows are mapped by ground-based field campaigns or airborne visual surveys [e.g., Kauahikaua et al., 2003]. During the current east rift zone eruption of Kīlauea, the boundaries of new flows are recorded approximately weekly with handheld Global Positioning System (GPS) receivers for the whole flow, or more frequently as lava approaches populated areas (T. Orr, personal communication, 2012). Thermal imagery, in the form of forward-looking infrared (FLIR) collected from a helicopter or on the ground, or by a variety of satellites with infrared detectors, is also used to identify new flows. Repeat surveys constrain the aerial coverage and advance rates of active flows. These data are then used with flow thickness data to estimate volumes and average effusion rates. SCM can also be used to map lava flow emplacement and augments traditional mapping techniques in a few important ways. For example, direct access to the flows is not required and the state of the entire flow field is recorded in a single image. SAR also has the ability to see through cloud cover and is not limited to making observations during daylight hours. Additionally, the high spatial resolution and continuity of the satellite images yield more internal detail than GPS points collected along flow margins. Finally, a constellation of SAR satellites can collect images of the same region with greater frequency than the on-average weekly repeat time for field mapping. In real-time applications, ground-based methods are still critical given the lack of operational control 10 in acquiring satellite images. Additionally, ground-based methods can better resolve flow fronts that enter vegetated regions and provide ground truthing of the coherence maps. The ability of SAR coherence to detect changes in the ground surface over time provides the basis for the lava flow mapping technique. SAR satellites travel in repeating orbits and record radar echoes backscattered to the satellite from the surface. Surface change can be measured by comparing these radar echoes from two SAR scenes collected at different times from approximately the same viewpoint. Where the ground surface is identical between the two images, the radar echoes will have high coherence. However, where the surface has been altered between the acquisition dates of the images, it will have low coherence. This comparative property between SAR scenes can be defined mathematically as the correlation coefficient, a spatially averaged measurement of the similarity of the radar signals in phase and amplitude, with values between 0 and 1 [Zebker and Villasenor, 1992]. The correlation coefficients for each pixel combine to create a coherence image. Decorrelation is often associated with vegetation and water, whose surfaces change significantly between SAR acquisitions [Zebker and Villasenor, 1992]. Any variation in the viewpoint of the satellite between the SAR scene pairs also impacts the overall coherence of the images by changing the angle and arrangement of scatterers that the satellite sees. The perpendicular baseline, which is the distance between the satellite locations during the two SAR acquisitions, must therefore be limited for images to be sufficiently correlated [Gatelli et al., 1994]. The length of time between SAR acquisitions further influences the overall coherence, because non-eruptive changes in the ground surface accumulate with time [Lu and Freymueller, 1998]. 11 New lava flows create decorrelation by repaving the ground surface over the extent of a flow, completely changing the ground’s scattering properties between images collected before and after flow emplacement. However, once a flow becomes inactive, its surface remains stable between image acquisitions, yielding high radar coherence. Using these observations, Zebker et al. [1996] showed that coherence images of the active flow field on Kīlauea can be used to create flow maps, because decorrelated new flows are surrounded by highly correlated older flows. Therefore, new flow areas can be mapped by thresholding the image with a correlation value that represents the margin of the new flow [Zebker et al., 1996]. This technique was also used to define flow boundaries from the 1995 Fernandina and 1998 Cerro Azul eruptions in the Galápagos [Rowland et al., 2003], and the areal extent of the 1997 lava flow in the caldera of Okmok volcano [Lu et al., 2000]. Similarly, non-lava sources of surface change, including pyroclastic flows at Soufrière Hills Volcano, Montserrat [Wadge et al., 2002], and the Portuguese Bend landslide in California [Calabro et al., 2010], have been mapped with SAR coherence. In combination with InSAR, SAR coherence has also been used to measure post- emplacement behavior at lava flows on Okmok and Etna volcanoes. The 1997 Okmok lava flow and the 1986–1987, 1989, and 1991–1993 lava flows on Etna remained decorrelated for several years due to continued surface change from cooling, settling, and contraction [Lu et al., 2000; Stevens et al., 2001; Lu et al., 2005]. Once a flow is correlated it is possible to extract a deformation signal with InSAR, which reveals post- emplacement subsidence [Briole et al., 1997; Stevens et al., 2001; Lu et al., 2005]. Our study builds on these previous mapping applications of SAR coherence by creating an algorithm to map successive lava flows from the ongoing Pu‘u ‘Ō‘ō- 12 Kupaianaha eruption of Kīlauea, Hawai‘i. This eruption on the middle east rift zone of the volcano began in 1983 with episodic lava fountaining from the Pu‘u ‘Ō‘ō vent that lasted until 1986. The eruption then continued through a constant effusion phase from the Kupaianaha vent from 1986 to 1992 and a nearly continuous effusion phase from Pu‘u ‘Ō‘ō from 1992 to 2007, before shifting downrift to a vent between Pu‘u ‘Ō‘ō and Kupaianaha called the “D” vent, because it was along Fissure D of the July 21, 2007 eruption, from 2007 to 2011 [Heliker and Mattox, 2003; Poland et al., 2008; Patrick et al., 2011; Patrick and Orr, 2012]. This eruption provides an ideal locality and event sequence to develop a SAR coherence flow mapping technique because of the multitude of flow field SAR scenes collected over time, the lack of snow or other environmental factors that could interfere with the radar detection of flows, and the detailed documentation of the eruption by the U.S. Geological Survey Hawaiian Volcano Observatory (HVO). By integrating a range of archived SAR data we assemble a detailed history of lava flow behavior at Kīlauea over an eight-year period, while developing an image analysis methodology to automate the detection of new flows. Our results provide flow emplacement data that suggest numerous potential applications of this technique to the study of lava flows around the world. 2. Methodology 2.1. SAR Data Sources and Processing We use 211 scenes from six tracks of Envisat SAR data (two ascending and four descending) to capture lava flow activity from January 2003 to October 2010 (Figure 1). Using the GAMMA processing package, we produce hundreds of coherence images 13 between SAR scenes on respective orbital tracks for 35- to 105-day intervals (based on a 35-day repeat time), although the small number of SAR acquisitions early in the time sequence results in some coherence images with longer durations. To maintain overall coherence, we assemble coherence images only from SAR scene pairs with perpendicular baselines of less than 600 m. We maximize spatial resolution by processing the images at high resolution (one look in the range direction) and georectifying them with an oversampled digital elevation model (DEM) produced from the Shuttle Radar Topography Mission (SRTM) [Farr and Kobrick, 2000] to yield a final pixel size of 20 m. Finally, we crop the coherence images from all tracks to the same geographic extent, with each pixel assumed to represent an identical location on the volcano’s surface. Figure 1. Location map of the study area on the east rift zone of Kīlauea, Hawai‘i Island (Big Island), Hawai‘i. Red boxes show the extents of Envisat SAR tracks and frames, while line-of-sight (LOS) vectors indicate the satellite look directions for ascending and descending tracks. 14 The coherence images document the degree of surface change (i.e., growth of vegetation, lava flows, etc.) over the duration of the image. Unlike deformation, which is directional and can only be measured along the satellite line of sight with interferometric SAR, the change in coherence from resurfacing is independent of the satellite viewpoint. This allows us to combine data from multiple tracks, and potentially even multiple SAR satellites, and assemble a detailed time series with a much higher temporal resolution than a single track of any one satellite. 2.2. Image Analysis Our SCM technique combines and analyzes SAR coherence images to extract maps of lava flow activity with time, utilizing a series of image analysis algorithms in MATLAB®. The method maps new lava flows by thresholding the coherence images to identify decorrelated areas. Flow identification is complicated by nonvolcanic causes of decorrelation within the images, which must be masked out, and fluctuations in the overall coherence of the images, which require flexibility in the threshold value. Because coherence images from different tracks overlap in time, we can combine all of the images to create a time series, in which each time series frame spans one period of overlap, called a time series epoch. We then filter these time series frames to remove spurious results and thereby produce a final series of flow maps. The methodology is laid out in Figures 2 and 3, with the flow chart in Figure 2 displaying how an example initial coherence image from 2007 is processed and combined with overlapping images to create a time-dependent flow map. 15 2.2.1. Vegetation Removal The cropped original coherence images display the degree to which the surface has changed between two SAR scene acquisitions. Existing, stable flows have high radar coherence, or correlation values near 1, whereas active lava flows result in regions of very low coherence, with the mean correlation value of mapped flows equal to 0.13 ± 0.07 (1σ) (Figure 2a, Step 1). However, coherence also tends to be low in vegetated areas and over water, requiring these features to be masked in the images. The coastline of the ––––––––––––––––––––––––––––––––– Figure 2 (next page). Flow chart of the SCM method. Image processing proceeds in the following order: (Step 1) A coherence image (a) is produced from a pair of SAR scenes and displays the degree to which the scattering properties of the ground surface in each pixel have changed between May 14 and June 18, 2007. Areas that are unchanged have high coherence, with values near one, while regions with altered surfaces have low coherence, with values near zero. Correlated features including the 1983–2007 flow field, labeled “recent flows,” and the Mauna Ulu flow field, can be identified in the image. Visibly decorrelated areas include vegetation beyond the edge of the flow fields, along with regions within the flow field that represent new activity. The ocean appears completely black in (a) because it is masked by the extent of the DEM used to create the coherence image and represents a region of no data. (Step 2) Vegetated areas within (a) are removed with a vegetation mask to create an image (b) where the remaining decorrelated pixels are potential lava flows. All pixels in the vegetation mask, as well as the ocean, have been made white. The vegetation mask for (a) is built by creating an initial vegetation map, and then removing any pixels that are vegetated initially, but become correlated by the end date of the image (06/18/2007). (Step 3) Multiple coherence images, each with its vegetation masked, are thresholded to select the most decorrelated pixels. The threshold applied to each image is a function of the perpendicular baseline, duration, and time of year of the images. For example, applying a threshold of 0.376 to (b) produces (d), which is a binary image where pixels with values below the threshold are colored black. The offset timelines in the center show the temporal locations of images (c) and (d). (Step 4) These two thresholded images (c and d) combine with sixteen others that overlap the 05/14/07–05/30/07 time series epoch to produce the time series frame (e; see Figure 3). Only pixels that are decorrelated in all eighteen images appear as gray or red in the time series frame. (Step 5) By selecting only decorrelated areas of the time series frame that fall within the flow field and filtering out any regions that are within the scale of the noise (less than 0.008 km2 in area) or represent sieve-like partially vegetated regions (area to perimeter ratios of less than 1.4, determined by trial and error), we isolate the final flow map for the 5/14/07–5/30/07 epoch, shown in red in (e). 16 DEM used to process the coherence images masks the ocean. Vegetation could be masked with an independent map of vegetation at the start of the time series. Alternatively, we can create a vegetation map from the SAR data by assuming that these features stay decorrelated throughout the whole time series unless covered by lava. We determine the coherence levels of vegetated areas by calculating the average coherence over a 0.1 km2 vegetated area in each image, and then use this to define a threshold to identify vegetated regions. We assemble an initial vegetation mask from pixels with values below this threshold in the first 26 coherence images (through November 2004). This captures persistently decorrelated areas that represent the extent of the vegetation at the beginning of the time sequence. 04/25/2007−05/30/2007 155.1°W 155°W 19 .3 °N 19 .4 °N 155.1°W 155°W 19 .3 °N 19 .4 °N 05/14/2007−06/18/2007 05/14/2007−06/18/2007 155.1°W 155°W 19 .3 °N 19 .4 °N 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 05/14/2007−06/18/2007 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Ocean Vegetation Recent flows New flow 155.1°W 155°W 19 .3 °N 19 .4 °N Mauna Ulu Flows C A B E D 04/25/07 05/30/07 05/14/07 06/18/07 Image C Image D Time series epochs Image E{+ 16 other overlapping images 05/14/2007−05/30/2007 155.1°W 155°W 19 .3 °N 19 .4 °N Time series frameThreshold = 0.376 Threshold = 0.376 2. Mask vegetation 3. Apply thresholds 4. Combine into time series 5. Filter and identify flows 1. Process images 17 Figure 3. Schematic demonstration of how the time series inversion captures flow activity over a single pixel. Two hypothetical lava flows over the pixel are shown in the first timeline. These cause decorrelation in the coherence images (shown as timeline segments) that span their emplacement in three simulated satellite tracks. The tracks each record the same activity, but are offset in time, allowing us to combine them to produce a time series with greater temporal resolution, divided into epochs. Each epoch spans the time between each SAR scene, regardless of its track. The pixel is defined as decorrelated during the epoch only if it is decorrelated in all images that overlap the epoch in time. In this way we are able to document the flow activity with great accuracy and better resolution than a single satellite track. As flows penetrate the vegetation at the margins of the flow field, the vegetation mask must shrink accordingly. Once a flow enters vegetation and becomes inactive, the flow becomes correlated. Therefore, if pixels that are part of the vegetation mask become correlated, a lava flow must have been emplaced in the previous time step. The algorithm accounts for this by assembling a dynamic vegetation mask that is modified through time, removing pixels from the mask for a given image that are correlated in any of the coherence images that overlap or immediately follow it in time (Figure 2b, Step 2). This Lava flow activity: Coherence images: time Correlated Decorrelated Track 1 Track 2 Track 3 Image 1 Image 2 Image 3 Image 4 Recorded flow activity: Time series: 1 2 3 4 5 6 7 8 9 10 11 12 13Epoch 18 allows pixels that become correlated in the next time step to be identified as lava flows in the current time step. 2.2.2. Image Thresholding With the vegetation and ocean masked, the decorrelated areas within the images should represent lava flows. These flows are identified by applying a threshold to the newly masked images to create a binary image, with the pixels of interest below a threshold correlation assigned values of 1, and all others assigned a value of 0 (Figure 2, Step 3). However, because each of the coherence images was created from a slightly different orbital geometry and surface condition, the threshold correlation value that marks the margin of a flow varies from image to image. We therefore calculate a unique threshold for each masked correlation image based on the perpendicular baseline, duration, and time of year of the image. We compute the threshold by measuring the average value of four 0.1 km2 areas of existing flows at widespread locations within the image and fitting curves to these average correlation values plotted against baseline, duration, and the start date of the image. We find that none of these variables exerts a strong influence on the coherence of the images, in contrast to locations such as Alaska that have a strong seasonal signal [Lu and Freymueller, 1998]. The effects of elevation, rainfall, and slope were also tested but revealed little influence. Only the pali, a steeply sloping terrain where lava flows have covered a large fault scarp, has poor coherence, and only at long baselines. Through this process, each image is assigned a unique threshold value based on the set of calibration variables, with an average threshold value of 0.36 for the entire suite 19 of images. For each masked coherence image, the corresponding binary images highlight the decorrelated features that do not fall within the vegetation mask (Figures 2c and 2d). 2.2.3. Time Series Assembly The thresholded images identify resurfacing of the ground surface within 35 to 350-day windows between SAR acquisitions on each track. Although the thresholded images can only be produced from SAR scene pairs within the same track, the tracks are offset in time, making it possible to create a new time series of images that represents the areas that are decorrelated where images from any track overlap in time (Figure 2e, Step 4). Figure 3 demonstrates schematically how this time series inversion produces an image between each consecutive pair of SAR acquisition dates, which we call time series epochs. The combination of the tracks into a time series is based on the principle that every thresholded image that spans the period of the emplacement of a lava flow must be decorrelated over the extent of the flow. Therefore, only pixels that are decorrelated in all of the thresholded images that overlap with a given time series epoch are mapped as decorrelated within the epoch. Because this removes pixels that are not consistently decorrelated over a time series epoch, it acts to filter out noise from areas with correlation values barely below the threshold. Importantly, it allows us to record lava flow activity at higher temporal resolution than a single track would permit. However, lava flows that have entered (decorrelated) vegetation are identified only when they become correlated. For this reason, the extent to which flows have advanced into vegetation can be mapped only retroactively, and may therefore overlap in time with images that are fully correlated. To account for this difference, pixels that become correlated between the start 20 and end of each time epoch (i.e. those pixels that represent flows that entered vegetation) are added to the time series image for that epoch to create the final flow map. 2.2.4. Flow Identification Despite these image processing techniques, each time series frame shows some level of scatter in the signal that does not correspond to the active flows during the epoch (shown as gray pixels in Figure 2e). One significant source of error is decorrelation in some time series frames due to the steep slopes of the pali. By using the DEM to mask out parts of the image where the slope is greater than 18°, the worst of this noise is eliminated. We then assume that during an epoch, new lava will occupy a contiguous region within the flow field, as defined by a simplified polygon. A filter based on the areas and perimeters of these regions removes local noise and the sieve-like areas of partial vegetation that are occasionally decorrelated. These filters allow us to map regions that fit the previously defined criteria, have an area greater than 0.008 km2 (20 pixels), and an area to perimeter ratio of greater than 1.4, as lava flows. Perimeters calculated in this analysis represent the distances between the centers of adjoining edge pixels, while areas include the full pixel extent. Threshold values for area and the ratio of area to perimeter were determined from trial and error to remove known noise without eliminating known lava flows. This image analysis algorithm produces a flow map such as that shown in red in Figure 2e. 21 3. Results Each of the 210 flow maps displays the activity associated with Kīlauea’s east rift zone eruption over an epoch; epochs range in duration from 1 to 137 days, with an average of 13.4 days. When viewed chronologically, the maps show the progression of the eruption from January 2003 to October 2010 (video available in the supplemental files included in this dissertation). The maps can be grouped by the major episodes as documented by HVO (Figure 4). Here flows are colored according to time sequence, with blue being the oldest, and red the youngest, of each time period. We then compare our composite flow maps with composite flow maps from HVO (Figure 5). The flows active over the first year and a half of the time sequence, from January 20, 2003 to July 21, 2004 (Figure 4a), belong to episode 55 of the Pu‘u ‘Ō‘ō-Kupaianaha eruption, which began in 1997 and lasted until 2007. During this part of episode 55, flows traveled down the western margin of the flow field (blue flows) and then found a different path to the ocean in 2004 (red flow). Rootless shields [Kauahikaua et al., 2003] over the lava tube were active during the entire period and are delineated by the zone of decorrelation to the south and east of Pu‘u ‘Ō‘ō. Older rootless shields formed in 2002 appear as an hourglass-shaped region of decorrelated pixels southeast of Pu‘u ‘Ō‘ō. The widespread blue regions of decorrelation (Figure 4a) are due to noise and not included in the following analysis. Immediately to the southeast of Pu‘u ‘Ō‘ō, the red flows represent the earliest stage of the next period of activity. 22 Figure 4. All 210 flow maps in the time series are displayed in four panels over the SRTM DEM [Farr and Kobrick, 2000]. Within each, the oldest maps are shown in dark blue grading up to the newest maps in maroon, with the number of days since the start of the panel shown on the adjoining color scale. (a) Middle part of episode 55. The large decorrelated regions in light blue are due to noise. The hourglass shape to the east of the main flows reveals persistent activity at several rootless shields that became quiet in August 2002. (b) The latter part of episode 55. Continued activity and persistent decorrelation of the flows emplaced during the interval shown in part A can be seen in the dark blue areas to the west of the main flows. (c) The July 21, 2007 flow and the Thanksgiving Eve Breakout (TEB) of episode 58. The July 21, 2007 D-vent in the southwest part of the July 2007 flow field and breakouts along the tube are apparent. (d) Continued breakouts from the TEB tube produce the March to October 2010 flows along the eastern margin of the TEB flow field. The extensive decorrelation in the final maroon map is noise. (e) A timeline of the study period displays the temporal distribution of SAR scenes. Note the lower density of scenes toward the beginning of the time sequence. 23 Figure 5. (a) SCM flow map summarizing flow field activity from January 2003 to January 2010. The colors show the different eruption episodes as defined by the USGS. The extents of Figures 6, 7, and 8 are shown as black outlines. (b) SCM flow map of flow field activity from July 2003 to October 2010, displaying combined pre-2007 flows, the July 2007 to January 2010 flow field, and the March to October 2010 flow. (c) USGS flow map of activity from 1983 to January 2010, corresponding with (a). (d) USGS flow map highlighting the 2007 to 2010 flow field and the March to October 2010 flow, for comparison with (b). The period from July 21, 2004 to June 17, 2007, corresponding to the last three years of episode 55, was marked by activity within a flow field and set of tube systems distinct from those of previous flows (Figure 4b). Lava flows generally migrated from west to east during this period and created many different ocean entries. West of these 2004–2007 flows, the narrow 2004 flow is still evident as a chain of blue pixels in the flow maps. The rootless shields above the tube from the previous period are also 24 decorrelated and remain so until late 2006. During July 17–19, 2007, an intrusion and small eruption occurred uprift of Pu‘u ‘Ō‘ō—Episode 56 of the long-term east rift zone activity. After an 11-day pause, lava reappeared at Pu‘u ‘Ō‘ō on July 1 and began refilling the cone (Episode 57) [Poland et al., 2008]. Neither of these episodes are specifically identifiable in the flow maps. On July 21, 2007, activity shifted from Pu‘u ‘Ō‘ō to a series of fissures in the northeast part of the flow field to mark the start of episode 58 of the Pu‘u ‘Ō‘ō- Kupaianaha eruption [Kauahikaua, 2007]. Within a few days the eruption localized at a single vent—the July 21, 2007 D vent—and began feeding lava to the northeast, eventually forming a 2-km-long perched lava channel [Patrick et al., 2011]. Northeast- directed flows lasted until November 21, 2007, as shown by the dark blue flows in Figure 4c. At this time, the flow direction changed when a breakout at the D-vent initiated the Thanksgiving Eve Breakout (TEB) flow, which formed a new lava tube to the ocean [Patrick and Orr, 2012]. Breakouts from this tube occurred throughout 2008–2010. The source of the TEB flow, and the rootless shield complex that formed along the lava tube to the southeast, remained decorrelated through 2009. Breakouts along the tube yield a chain of fan-shaped regions of decorrelation, but the tube itself is correlated, and therefore not mapped with SCM. The advance of some of the breakouts to the ocean can be seen in the progressive color change along the individual flows. The last nine months of the time sequence are dominated by breakouts from the TEB lava tube between January 22 and October 13, 2010 (Figure 4d). Early in this period, a breakout that initiated in January traveled south from the vent but failed to make it to the coast (shown in blue). On March 12, a larger breakout started feeding a flow to 25 the southeast. Once established, the tube system of the March 2010 breakout transported lava to the ocean on several separate occasions between April 2010 and until the end of the time series. This flow is highlighted in Figure 5b for comparison with the flow map produced by HVO. 4. Discussion 4.1. Evaluation of Results To assess the accuracy of the maps produced by SCM we compare them to the maps created by HVO. The coherence maps (Figures 5a and 5b) correspond well to the USGS maps produced using both aerial and field surveys (Figures 5c and 5d). In detail, however, there are some differences between the SAR-generated and field-generated flow maps. First, in Figures 5a and 5c, the northwestern-most lobe of the 2004–2007 flows is not fully captured by the SAR because these flows were emplaced on top of the rootless shields from the 2002–2004 flows, which had not yet become correlated. Second, a comparison of Figures 5b and 5d shows differences in the detail and extent of the March–October 2010 flow over the pali. This is caused by the difference in spatial resolution between the SAR and the USGS data. Satellite imagery is spatially continuous, whereas field mapping is often limited to the flow margins. The internal structure of the active flows, including the kīpuka, is therefore imaged more completely in the SCM maps. We examine this point in more detail below. It is also important to note that Figures 5a and 5c are not directly comparable because of the absence in the SCM map of flows emplaced before 2003, when the coherence time series begins. 26 4.1.1. Errors in SAR-Based Mapping Early in the time series, there are many fewer SAR acquisitions, reducing the number of overlapping images and the temporal resolution (Figure 4e). This led to problems with some of the early flow maps, which are quite noisy and therefore were not included in the summary maps (blue regions in Figure 4a). Limited image acquisition meant that each map was produced using only one coherence image. The SCM technique works best when many images overlap a time epoch and can therefore be combined to filter out noise from images with poor coherence or decorrelated regions due to steep topography. Another source of error relates to difficulties in generating SAR-based flow maps in vegetated regions, as illustrated by the mismatch of maps showing the northeastern July 2007 flow (Figures 5a and 5c). This highlights a major weakness of SCM. Although vegetated pixels can be mapped as lava once they become correlated, many of these pixels do not appear in the final maps, because the flow does not become correlated all at once; instead, parts of the flow remain decorrelated for longer periods than others. When this happens, our flow identification algorithm filters out parts of these flows that are small in size and spatially isolated. Our ability to map flows into vegetation therefore depends on how long the flow stays decorrelated. Flows that travel into vegetation can only be mapped after they become correlated, which may be months after emplacement. This prevents mapping flows into vegetation in near-real time. We illustrate this problem by comparing an SCM flow map from July to October 17, 2007 with a flow map from a field survey that we conducted on October 16, 2007 (Figure 6a). This comparison shows that the coherence- 27 derived flow map is confined primarily to the previously non-vegetated areas, the exception being flow lobes within the vegetation that had become correlated by October. Flow lobes that remained decorrelated for long time periods (up to 2.5 years) are not fully mapped because of the piecewise manner in which the flow becomes inactive and correlated. Figure 6. (a) The SCM map of the July 2007 flow extent as of October 17, 2007 (pink) is overlain on a pre-eruptive satellite image, and compared to the flow extent mapped on October 16 (yellow outline). The SCM technique cannot detect the flow where it enters vegetation, except where it has become inactive and correlated (older kīpuka in the center and white arrows). Even once the entire July 2007 flow becomes correlated, SCM fails to capture the full extent (compare to Figure 11 SCM map extent). (b) ALOS coherence image of the July 2007 flow field from July 16 to October 16, 2007. The longer wavelength (L-band; 23.6 cm) of the ALOS satellite yields higher coherence of the vegetation, allowing the discernment of the full extent of the July 2007 into the vegetated regions. Incorporation of ALOS data into this technique has the potential to solve the problem of mapping flows into vegetation. Possible solutions to this vegetation problem include independent measurements of flow extent from thermal imagery or other mapping to document the true flow margins. Alternatively, the addition of SAR imagery from another satellite with a longer wavelength has enormous potential to both solve the vegetation problem and enhance the temporal detail of the time series. We demonstrate this potential with an image from the 28 ALOS satellite, which uses a wavelength of 23.6 cm (L-band), compared to the 5.6 cm wavelength of the C-band Envisat satellite. The longer wavelength makes it less sensitive to minor changes in the ground surface and increases the coherence of vegetated areas [Zebker et al., 1997]. For this reason, the decorrelation signal of the active lava flows is more easily recognized (Figure 6b). The longer wavelength, and therefore lower sensitivity, also yield a quicker return of the flow to correlation. 4.1.2. New Views From SAR-Based Mapping The SCM technique contributes new insight into flow emplacement processes in that it goes beyond mapping the surface flow margins that have been meticulously documented by the USGS. This is because decorrelated pixels not only indicate new surface flows, but also lava tubes, ocean entries, and eruptive vents, which are areas of strong deformation and numerous small breakouts. An example of internal detail within the margin of a flow being captured by SCM is demonstrated by a comparison of a SAR- based flow map to USGS tube and flow margin maps from May to June 2007 (Figure 7). The SAR-based flow map shows extensive activity near the ocean entry over its 19-day duration. The origin of this decorrelation can be determined by examining a short- wavelength infrared (SWIR) image from the ALI instrument on the EO-1 satellite that captures thermal activity on the flow field at an instant in time. Here hot spots in the SWIR image (yellow dots in Figure 7) lie directly on top of the edges of the SAR- mapped flows near the coast, in addition to areas on the tube. This correspondence confirms that much of the extensive area of decorrelation near the ocean entry is caused by numerous breakouts. 29 Figure 7. The SCM flow map from May 30 to June 18, 2007 (red) lies within the boundary of the USGS map from June 13, 2007 (pink) (T. Orr, personal communication, 2012). The SAR displays activity over the lava tube (black) and the ocean entry, where severe deformation and breakouts can cause decorrelation. Yellow regions display the hottest pixels from a short-wavelength infrared satellite image collected by EO-1 ALI on June 3, 2007. These hot areas represent active breakouts and correspond well to locations of SAR decorrelation. Lava tubes, which can be difficult to map in the field without the aid of thermal imaging, can also be observed in some of the coherence images. Once a tube is established below a roof of inactive flows, it becomes sufficiently correlated to be omitted from the final flow maps. However, active lava tubes can often be identified as a moderately decorrelated path through the flow field in the original coherence images (Figure 8). The reduced coherence of the tube is likely a product of deformation and disturbances of the surface over the tube, which may occur during inflation of filled tubes or by changes in flux through open tubes [e.g., Kauahikaua et al., 1998; 2003; Orr, 2010]. 30 Figure 8. View of the pali and coastal plain in a coherence image from 07/28/2009 to 10/06/2009. The fully decorrelated area in the center (1) shows breakouts over the pali near the Royal Gardens subdivision. The thin line of moderately decorrelated pixels leading from these breakouts to the coast (2) reveals the path of a lava tube. The route matches well with USGS mapping of the tube from the summer of 2009 (T. Orr, personal communication, 2012). Finally, coherence maps highlight breakouts and active regions contained within an active flow field, which are details that are not generally included in USGS maps. For example, the SAR- and field-based maps of the March to October 2010 flow (Figures 5b and 5d) are similar in outline but quite different in internal flow structure. To evaluate the source of this difference, we compare a single (3-day) coherence flow map with an (instantaneous) oblique FLIR image from early May 2010 (Figure 9) (M. Patrick, personal communication, 2012). This comparison shows that each individual decorrelation region of the March to October 2010 flow corresponds to an active or recently active breakout, as shown by hotter surface temperatures (warmer colors) of the FLIR image. Other regions of the flow are both cooler and correlated, despite an active 31 lava tube beneath the surface. The flow map does not extend all the way to the 2010 coastline, which has been built by new lava deltas since 2000, when the DEM used to create the coherence images was created. Figure 9. (a) Oblique view of the SCM flow map (looking toward the northwest) from May 4 to May 7, 2010 draped over topography and a satellite image. (b) FLIR image collected on May 7, 2010 from a helicopter looking upslope to the northwest (M. Patrick, personal communication, 2012). The colors show approximate surface temperature, with the hottest colors representing new breakouts of lava. Locations 1, 2, and 3 mark common features in the images. Decorrelated pixels (red) correspond to new and recent breakouts, while cooler, existing parts of the flow appear correlated. The different coastlines of the SCM map and 2010 imagery show the growth of the island into the ocean, because the SCM map is limited to the extent of the 2000 coastline. In summary, we find that SAR-based flow maps correspond well with field-based maps except where flows enter vegetated regions (or form new lava deltas, e.g., Figure 9). Problems in mapping lava entry into vegetated regions can be alleviated by using longer wavelength (L-band) SAR images. We also show that SAR-based flow maps provide detailed views of flow emplacement, because they highlight active regions inboard of flow margins that are often inaccessible from the ground. These include active lava tubes and surface breakouts. SCM thus complements both traditional flow mapping techniques and more recent thermal imaging strategies, which are limited by cloud cover 32 (satellite-based thermal imaging) and spatial extent (handheld or helicopter-mounted FLIR cameras). 4.2. Analysis and Applications The ability to map lava flows at high spatial and temporal resolution has numerous potential applications to the study and monitoring of volcanoes. We have demonstrated that flow maps generated by SCM offer insights into planform flow morphology, which is useful for characterizing flow emplacement behavior. Flow lengths, widths, and areas can be measured as on any map and time sequences provide information on flow advance rates. For example, the TEB event produced a flow that reached the ocean on March 5, 2008. Its progress can be documented with the SCM maps from the time of its breakout (constrained by SCM between November 21 and November 25, 2007) until it reached the ocean (constrained by SCM between March 5 and March 13, 2008). Initial average rates of flow advance were 2.4 m/hr; flow advance then accelerated to 13.6 m/hr (Figure 10). These data agree with estimates from field mapping (Figure 10 inset) (M. Patrick, personal communication, 2012), except in the length of the western branch of the flow. Near its final extent, this branch traveled through vegetation and narrowed, preventing SCM from identifying its full length. Flow advance rates can be calculated with SCM only for flows that are emplaced over more than one time epoch; thus, they are time-averaged rates applicable to flows with durations of days or months, rather than hours. The final length of the TEB flow as of July 7, 2008 is 12 km, with an average width of 550 m, and a total area, including the July 2007 flow extent, of 16 km2. This total area is less than the 21 km2 mapped in the field (T. Orr, personal 33 communication, 2012) because the SCM technique missed regions where the flows went into vegetation. Figure 10. The extent of the Thanksgiving Eve Breakout flow, which began November 21, 2007, is displayed with time. From this map we can calculate the advance rate of the flow, shown in the inset plot of the length of the flow with time. Our data mostly compare well with field-based advance data plotted as a black line (M. Patrick, personal communication, 2012), with the exception of the distal extent of the western branch, plotted in green. The end of this flow is lost with SCM due to its passage through vegetation and narrow morphology. The flow accelerates as it advances over the pali. SCM lava flow maps also provide important information on other aspects of lava flow emplacement. Particularly intriguing is the wide variation in the time over which different parts of a flow remain decorrelated after initial emplacement. This variation is particularly apparent in the northeastern flows that were emplaced in 2007 [Patrick et al., 34 2011], where decorrelation durations range from 13 to 903 days. When examined in map view, it is clear that the center of the flow remained decorrelated much longer than the edges (Figure 11). Figure 11. The length of time each pixel in the July 2007 flow is decorrelated after initial emplacement matches well with the surveyed flow thickness contours from October 16, 2007. The thickest part of the flow, 35 m, is decorrelated an average of 375 days and corresponds to the perched lava channel [Patrick et al., 2011], whereas the thinnest, 2-m- thick flow margins are decorrelated for an average of 94 days. The area of the flow just south of the thickest section is decorrelated the longest, up to 903 days, because it is the source of the TEB flow in November 2007 (July 21, 2007 D-vent), where activity persisted long after the rest of the July–November 2007 flow field became correlated. Where the flow entered vegetation and was recovered by the flow-mapping algorithm, the decorrelation duration is a minimum, because it only represents the duration of the epoch prior to when the flow became correlated. This yields the yellow holes within the orange region in the center of the flow. To test the extent to which this decorrelation pattern reflects variations in flow thickness, we overlay our flow thickness data onto the SAR decorrelation map (Figure 35 11). The thickness data were acquired in mid-October 2007 by kinematic GPS using Leica SR520 receivers carried on numerous traverses while sampling at 1 Hz frequency. Base station data came from nearby continuous GPS stations. The kinematic data (x, y, z) are accurate to within a few centimeters. Each traverse started and ended a few meters off the new lava flows so that altitudes could be correlated with SAR digital elevation data acquired in 2005. Flow thickness correlates reasonably well with decorrelation time (Figures 11 and 12), with the thickest area of the flow showing an average decorrelation duration of 375 days, compared to 94 days within the thinnest regions. The thickest part of the flow corresponds to the site of the perched channel, which localized lava accumulation in the area and remained consistently active until the TEB breakout [Patrick et al., 2011]. The location of maximum decorrelation duration lies just south of the thickest portion, because the D-vent and rootless shields led to persistent activity in the area long after the July 2007 flow became inactive [Patrick and Orr, 2012]. Our data underestimate decorrelation durations for the parts of the flow that entered the vegetation because of the problems with flow identification in this area that are outlined above. To test the relationship between decorrelation duration and flow thickness, we calculate the average number of days decorrelated within each contour; when this average is compared to the measured thickness of the flow, we find a linear relationship (r2=0.97; Figure 12). Here the minimum duration of decorrelation (94 days) reflects the average time over which the flow was emplaced. Continued decorrelation after flow emplacement is caused by changes in the flow surface after emplacement. Thus the observed relationship between decorrelation duration and flow thickness reflects both post- 36 Figure 12. Comparison of the surveyed flow thicknesses of the July 2007 flow and the average time pixels are decorrelated within each thickness contour from Figure 11. Each point represents the mean length of time a pixel is decorrelated in the areas between each of the thickness contours, while the errors bars show the standard deviation. Error also arises from the uncertain timing of activity due to the duration of the epochs, but it is minor compared to the range in time decorrelated within each contour. The horizontal black line marks the end of emplacement in November 2007. Decorrelation prior to this represents flow emplacement, while continued decorrelation is a product of further change within the flow, including cooling, settling, and subsidence of the flow [Stevens et al., 2001]. Some parts of the flow were emplaced over a period up to 112 days, but the flow edges were correlated after an average of 94 days. emplacement processes (such as deflation) and the response of the SAR coherence to these processes. The duration of decorrelation is in part a function of the wavelength used by the satellite, owing to the magnitude of change recognized at increasingly short wavelength. This effect becomes evident when comparing the decorrelation duration of the July 2007 lava flow as viewed using two different SAR satellites. In the ALOS imagery, which uses a wavelength of 23.6 cm, most of the July 2007 flow is correlated after only ~90 days. In contrast, the same flow is decorrelated for ~600 days in the 5.6 cm wavelength Envisat imagery. 0 5 10 15 20 25 30 35 400 50 100 150 200 250 300 350 400 450 500 Emplacement time Post- emplacement Flow thickness (m) D ay s de co rre la te d r2 = 0.97 37 The linear relation between flow decorrelation and flow thickness defined by the Envisat images of the July 2007 lava flow can potentially be used to estimate flow thicknesses for any flow in the time series and, when combined with flow areas and flow advance rates, can provide estimates of both flow volumes and average effusion rates. For example, using data for decorrelation durations ≤ 400 days, we calculate a volume of 62×106 m3 ± 9×106 m3 for the July 2007 flow. We can compare this to the volume determined from the field survey, calculated by standard GIS methods. First, the GPS data were contoured at 2-m intervals and gridded to derive a digital elevation model (DEM) depicting the upper surface of the new lava field. Then, the DEM of the preexisting topographic surface was subtracted to calculate the residual volume, which corresponds to the bulk-rock volume of the lava flows mapped, about 58×106 m3 emplaced between July 21 and October 16, 2007 (87 days). We estimate the error of this volume to be about 10 percent. The SAR- and field-based values correspond well, although the comparison is approximate because of underestimates of flow extent in the SAR coherence map (because of vegetation) and continued activity in the flow field until the end of November, which increased the flow thickness beyond the October measurements. However, it points to the possibility of using remote sensing data to obtain accurate measurements of flow thickness (volume) data that are otherwise time-consuming to obtain from ground-based surveys. A similar relationship between thickness and duration of decorrelation has been observed in lava flows at Etna, Italy [Stevens et al., 2001]. The origin of the decorrelation signal may lie either in cooling-related contraction [Briole et al., 1997; Stevens et al., 38 2001; Lu et al., 2005] or changes in the roughness of the lava surface (and scattering properties) because of weathering [e.g., Kahle et al., 1995; Mazzarrini et al., 2007]. If a lava flow remains decorrelated after emplacement because it is subsiding faster than the fringe rate for the satellite phase [Stevens et al., 2001], it will become correlated when the subsidence rate slows. Once the flow is correlated, a differential interferogram can be used to measure this subsidence rate. We have generated a differential interferogram that shows deformation of the 2007 lava flow over a 35-day period in 2009 (Figure 13). In this image, the thickest part of the flow is subsiding fastest (shown as darker blue). Moreover, nearly two years after the TEB event terminated the July–November 2007 flow, the subsidence was occurring at a rate of 1.4 mm/day. This is calculated from the range change measured at a near-vertical incidence angle of 18.8° and assuming only vertical motion. A correlation between thickness and subsidence rate has also been observed at Okmok and Etna volcanoes [Stevens et al., 2001; Lu et al., 2005]. Together these data show that the combination of coherence maps and interferograms has the potential to document post-emplacement processes in lava flows. 5. Conclusions We find that SAR Coherence Mapping (SCM) provides a superb means of documenting lava flow activity in Hawai‘i. By combining eight years of data from many satellite tracks, it is possible to map volcanic activity on the east rift zone of Kīlauea at high spatial and temporal resolution. This remote sensing technique has some advantages over traditional field mapping, including a large spatial extent, high spatial resolution, automation, and the ability to detect activity apart from surface flows. Disadvantages 39 Figure 13. Unwrapped Envisat interferogram from track 322 spanning October 21 to November 25, 2009 shows negative range change from the satellite in centimeters, with blue regions moving away from the satellite. Gray regions contain no data, corresponding to decorrelated areas, while the outline shows the mapped extent on the flow as of November 25, 2009 (T. Orr, personal communication, 2012). The contour interval for the topographic basemap is 50 meters. The interferogram reveals subsidence of up to 5 cm over the duration of the image assuming all motion is vertical, with 1.4 mm/day at the center of the July 2007 flow two years after emplacement. Some of the range change reflects a topographic component because the DEM was not updated with the morphology of the recent flow. However, this component is likely small (< 0.27 cm over the duration of the interferogram) given the small perpendicular baseline of 11.6 m and a maximum increase in elevation of 40 m with the July 2007 flow emplacement [Massonnet and Feigl, 1998]. include occasional low temporal resolution, the ambiguity of what decorrelated regions represent, persistent decorrelation masking new activity, and problems with mapping flows that enter vegetation. Overall, this method has the potential to be a valuable tool for monitoring remote volcanoes and to provide a wealth of flow emplacement data for exploring the behavior of lava flows. Future applications and directions of the SCM method are numerous. This technique can be used to map any form of surface change, such as the extent of ash falls, 40 floods, landslides, and pyroclastic flows. Furthermore, SAR data from any track of any satellite can be combined to produce the largest possible dataset with the highest temporal resolution. We hope that this tool can continue to be developed and facilitate the investigation of many geologic phenomena, as well as lava flow emplacement and volcanic processes. 6. Bridge In this chapter, I developed techniques to assemble a time series of planform maps of active lava flow fields from satellite-based SAR data. These maps allowed me to calculate flow advance rates through time, detect intra-flow activity, and map lava tube pathways. I also discovered the potential to collect three-dimensional data from the duration of post-emplacement decorrelation in the SAR imagery. This work developed an automated active flow mapping system that is now in service at the U.S. Geological Survey Hawaiian Volcano Observatory to aid in flow mapping. In the next chapter, I will build on this study to create methodologies to use airborne topographic data to map the extents, ages, and morphologic features of lava flows. These data are of higher resolution than the SAR data in Chapter II, but do not cover the development of active flows. Instead there are only two time steps, before and after the emplacement of a large ‘a‘ā flow. However, these are snapshots of topography, and offer detailed information about the surface morphology and volume distribution of the flow. 41 CHAPTER III LAVA FLOWS IN 3D: USING AIRBORNE LIDAR AND PRE-ERUPTIVE TOPOGRAPHY TO EVALUATE LAVA FLOW SURFACE MORPHOLOGY AND THICKNESS IN HAWAI‘I This chapter is in press as a book chapter in the American Geophysical Union Monograph, Hawaiian Volcanoes, From Source to Surface that is set to be published in 2014. I was the primary author, doing the bulk of synthesizing and interpreting data and writing the manuscript. Co-author Adam Soule (Woods Hole Oceanographic Institute) contributed related work on lidar analysis for this chapter, including the study of revegetation of the Mauna Loa flows and the analysis of pāhoehoe and ‘a‘ā distribution within the Kīlauea December 1974 flow. Kathy Cashman (University of Oregon) advised the project and helped with the manuscript and Benjamin Mackey (University of Canterbury) processed the pre-eruptive aerial photographs to create a pre-eruptive digital elevation model. 1. Introduction Mafic lava flows represent a persistent and widespread form of volcanic activity that affects communities around the world. Central to lava flow hazard assessment and mitigation are the construction of probabilistic flow hazard maps [e.g., Favalli et al., 2005; 2009a,b; 2012; Tarquini and Favalli, 2010] and development of tools for real-time prediction of flow paths, flow advance rates, and final flow lengths. Both hazard 42 assessment and predictive modeling require high resolution digital topography, as well as a fundamental understanding of the ways in which advancing lava flows interact with the Earth’s surface. For this reason, volcanologists have embraced new methods of imaging the landscape, focusing particularly on airborne lidar because of its combination of high spatial resolution and extensive spatial coverage [e.g. Mazzarini et al., 2005; Pyle and Elliot, 2006; Coltelli et al., 2007; Ventura and Vilardo, 2008; Marsella et al., 2009; Favalli et al., 2010; Deardorff and Cashman, 2012; Tarquini et al., 2012]. New multi- temporal imaging of active lava flows is now providing remarkable spatial and temporal coverage of the flow emplacement process [e.g., Marsella et al., 2009; Favalli et al., 2010]. Lidar data also provide important information on lava flow surfaces in the form of both intensity data and surface roughness. The combination of high-resolution topography and quantifiable surface characteristics has the potential to transform our understanding of lava flow emplacement. That this potential has not yet been reached attests to the need for both new analytical tools and improved models of lava flow behavior. We begin with an overview of previous work on the description of lava flow surfaces and the use of lidar to investigate flow morphology before introducing the work presented here. 1.1. Characterizing Lava Flow Surfaces Lava flow surface characteristics have been examined primarily in the context of discriminating flows of either different ages or flow type [e.g., Mazzarini et al., 2007; Morris et al., 2008; Bretar et al., 2013]. One measure of flow surface properties is derived from the intensity of the lidar return signal. A study of lava flows at Mt. Etna 43 [Mazzarini et al., 2007] shows that the lidar intensity is controlled primarily by surface reflectance and by the distance between the receiver and the imaged surface. When normalized to a standard elevation to correct for receiver distance (plane height in the case of airborne lidar), lidar data show intensity variations that reflect both flow surface roughness and flow surface characteristics such as surface weathering and vegetation coverage. For this reason, lidar intensity maps are difficult to interpret without additional information (such as flow age, vegetation cover, or an additional measure of surface roughness). Studies of lava flow surface roughness have been developed primarily by the planetary community for mapping the surfaces of the terrestrial planets. The most straightforward measure of surface roughness is the root mean square (RMS) height within a cell [usually taken as a square matrix of pixels; e.g., Shepard et al., 1995; Campbell and Shepard, 1996; Morris et al., 2008; Bretar et al., 2013]. The RMS height is strongly dependent on cell size and the scale of surface features and measures only variations in vertical elevation. It is also affected by slope, and thus requires detrending of elevations prior to calculation [eg., Campbell and Shepard, 1996]. Horizontal roughness can be assessed by determining the relationship between RMS height deviation and a correlation length scale [Morris et al., 2008; Bretar et al., 2013]. Periodic flow surface features, such as folds, can be measured by fast Fourier transform (FFT) of linear profiles perpendicular to the feature of interest [Pyle and Elliott, 2006]. Surface analysis has also been used to highlight abrupt topographic offsets [Morris et al., 2008] and to discriminate among different flow surface morphologies (e.g., tephra, smooth pāhoehoe, slabby pāhoehoe and ‘a‘ā) [Bretar et al., 2013]. 44 1.2. Characterizing Lava Flow Morphology Volcanological applications of lidar data have focused primarily on developing high resolution DEMs for probabilistic lava inundation models [e.g., Favalli et al., 2005; 2009a,b; 2012] and for morphometric analysis of (mostly recent) lava flows [e.g., Mazzarini et al., 2005; Pyle and Elliot, 2006; Coltelli et al., 2007; Neri et al., 2008; Ventura and Vilardo, 2008; Marsella et al., 2009; Favalli et al., 2010; Deardorff and Cashman, 2012; Tarquini et al., 2012]. Morphometric analysis includes assessment of along-flow changes in flow characteristics, planform analysis of flow margins, and detailed studies of specific flow features. Many of these studies have had, as a goal, increasing the accuracy of flow volume estimates and effusion rates, and, where detailed flow chronologies exist, changes in effusion rate with time. Another goal has been to identify topographic controls on both flow morphology and flow emplacement conditions. These studies have been supplemented by photogrammetric studies [Applegarth et al., 2010; Marsella et al., 2012] and field measurements of flow cross- sections [Branca et al., 2013]. Key to the analysis of topographic controls on lava flow emplacement is an accurate pre-eruption DEM or other information about the landscape traversed by the flow. Pre-eruption lidar data are rare, although current lidar mapping campaigns, particularly at volcanoes such as Stromboli and Etna in Italy, will provide reference topography for future lava flows. For this reason, the resolution of thickness maps constructed using post-eruption lidar is constrained by the resolution of the pre-eruptive DEM. Resulting maps have been used to measure total flow volumes and average effusion rates [e.g., Coltelli et al., 2007; Marsella et al., 2009] and to evaluate the effect 45 of topography on the development of individual lava channels [e.g., Favalli et al., 2010; Tarquini et al., 2012]. 1.3. A Brief Review of the Target Lava Flows The studies reviewed above have developed methodologies for lidar processing and demonstrated the potential of lidar for investigations of lava flow emplacement. None of these studies, however, has used flow thickness data to examine mass distribution and dynamics of multi-lobed flows, nor have morphometric and surface analysis characteristics been combined to address along-flow changes in both flux and flow characteristics. Here we do both, with our focus on two recent (December 1974 and March–April 1984) and well-characterized lava flows from Kīlauea and Mauna Loa volcanoes, Hawai‘i (Figure 1). The December 1974 lava flow erupted along a series of right-stepping NE-SW fissures southwest of Kīlauea caldera (Figure 1b). Fire fountaining migrated from the NE to the SW end of the fissure system over the course of the ~6 h eruption. The average volumetric effusion rate was ~ 270 m3/s and the total erupted volume was ~5.9×106 m3 [Lockwood et al., 1999]. Abutment against the Koa‘e fault zone forced the lava into a single channel along the last 9 km of the total 12 km flow length. At the same time, the flow surface changes from pāhoehoe at the vent to ‘a‘ā at the distal end, a consequence of lava cooling and crystallization within the active channel [Soule et al., 2004]. The March–April 1984 eruption of Mauna Loa volcano occurred on the northeast rift zone (Figure 1c). The eruption lasted three weeks and produced a large (~1.98×108 m3) lava flow. The main fissures that opened on March 26 initially fed four main flow 46 Figure 1. Location maps of the study area. (a) Hillshade map of the Island of Hawai‘i with the extents of the lava flows included in this study shown. (b) Map of the Kīlauea December 1974 flow in solid black, with the extent of the lidar in the black outline and the Koa‘e fault zone marked by white dashed lines. (c) Map of the Mauna Loa 1984 flow in solid black overlain on the Wolfe and Morris [1996] geologic map with the ages of surrounding flows. The lidar extent is shown in the black outline and the major branches of flow labeled after Lipman and Banks [1987] (Flows 1–4, along with Flow 1A and Flow 1B). 155°0'0"W156°0'0"W 20 °0 '0 "N 19 °0 '0 "N 155°15'0"W155°20'0"W155°25'0"W 19 °3 5' 0" N 19 °3 5' 0" N 155°20'0"W 19 °2 0' 0" N 0 5025 km a b c b c 3,000-10,000 1,500-3,000 750-1,500 200-750 0-200 1984 flow Flow ages (ybp) 0 42 km Faults 1974 flow 0 21 km 1A Flow 11B 2 3 4 Kīlauea Ma un a L oa Mauna Kea Kohala Hualālai 47 lobes (Flows 1–4) with an effusion rate of ~275 m3/s. Flows 2–4 stagnated within two days. Flow 1 advanced 27 km from the vents within the first few days of the eruption. A blockage of the Flow 1 channel and subsequent overflow on March 29 created Flow 1A, which had traveled nearly as far by April 5 (Q ~ 100 m3/s). Reduced flux and further blockages produced successive shorter overflows throughout the remainder of the eruption, including Flow 1B on April 5. The velocity, temperature, and channel dimension measurements along the length of this flow were closely monitored because of concern that it could reach the city of Hilo. The intense monitoring effort produced an unprecedented set of observations along a large active ‘a‘ā channel that greatly improved understanding of channel construction and flow emplacement [Lipman and Banks, 1987]. 2. Methods 2.1. Data Sources We use two types of airborne remote sensing to quantify lava flow morphology for both lava flows. Airborne laser swath mapping (i.e., ALSM or airborne lidar) provides topographic data of the present-day flow surfaces along with intensities of backscattered laser pulses. Photogrammetry allows us to view grayscale imagery from orthorectified aerial photographs and to construct topography of the pre-eruptive surfaces. Data products derived from these two data sources provide detailed views of the extents, properties, and internal structures of both lava flows. 48 2.1.1. Lidar The ALSM datasets were collected during a survey in 2009 by the National Center for Airborne Laser Mapping (www.ncalm.org) (Figure 1). This paper describes surveys over a portion of the Northeast Rift Zone (NERZ) of Mauna Loa and the December 1974 lava flow from the Southwest Rift Zone (SWRZ) of Kīlauea Volcano conducted with an Optech GEMINI Airborne Laser Terrain Mapper (specifications for this instrument at www.optech.ca). These and two additional surveys of the summit regions of Mauna Loa and Kīlauea Volcanoes, as well as operational details, are available at www.opentopography.org. ALSM creates a GPS-rectified x-y-z point cloud of topography by sending laser pulses from an aircraft to the ground surface and recording the time-of-flight of the reflections. The accuracy of individual lidar returns is determined by the size of the laser spot when it reaches the ground. At our survey altitude of ~1800 m, vertical accuracy is 5–30 cm and horizontal accuracy is ~30 cm. The survey collected ~5 returns per square meter, and the point clouds were gridded to produce a digital elevation model (DEM) of the ground surface topography. We did not filter laser returns from vegetation before gridding, because of the minimal vegetation cover. For this reason, the DEM does not fully represent the ground surface in the few areas where pioneering plants are present. We removed spurious points manually before using a linear weighted mean within a 3 m x 3 m window to generate DEMs with 1-m horizontal resolution, a significant improvement over existing DEMs for these regions. In addition to topographic data, ALSM provides surface imagery in the form of the amplitude of the laser returns. The near-infrared laser is more effectively backscattered by some surfaces than others, so that a map of lidar “intensity” reveals 49 differences in the near-infrared color of different features on the surface (Figure 2). The strength of the laser reflection is strongly dependent on the aircraft altitude, and thus requires altitude corrections to produce a normalized map of lidar intensity [Mazzarini et al., 2007]. In addition, laser intensity may be affected by differences in atmospheric transmission, e.g. due to humidity or presence of volcanic gases. This was the case for the Mauna Loa 1984 flow where there were persistent clouds due to the orographic effect at lower elevations. We corrected for atmospheric transmission effects by shifting normalized intensities on parallel flight lines so that features covered by both flight lines had consistent intensity values [Favalli et al., 2009a]. Figure 2. Annotated zoom of the lidar intensity map with lidar hillshade overlay to demonstrate the intensity differences between surfaces of different ages, cones, flows, and channels. Ages from Wolfe and Morris [1996]. The morphometric properties of this region are shown in Figure 5 and classified in Figures 9 and 10. 50 2.1.2. Photogrammetry To visualize topography prior to flow emplacement, we used aerial photographs taken in 1977 by the U.S. Geological Survey (USGS) to construct a DEM of the pre- eruptive surface of the Mauna Loa 1984 flow (Figure 3). The aerial photographs were high-resolution (14 µm scan, ~1 m ground resolution) scanned photographs acquired from the USGS; the accompanying camera calibration report allowed us to constrain the internal orientation of the camera. A set of seven images provided stereo-pair coverage across the whole flow path (Figure 3a). Both the DEM and accompanying orthophotos were generated using Leica Photogrammetry Suite (LPS) software. We identified 50 suitable ground control points (GCP’s) on unfiltered shaded-relief lidar maps, and used the 1984 flow outline mapped from newly collected lidar elevation and intensity data to ensure that the GCP’s were outside the area covered by the 1984 flow. GCP’s were typically distinctive small-scale patterns of vegetation (shrubs), or structural features in the older (pre-1984) lava flows. The lidar DEM was used for GCP elevation control. We then identified the corresponding points in the aerial photographs using the LPS graphical interface. We identified points common to multiple photographs, with at least 3 control points between overlapping photos, and augmented this set with the LPS automated tie point generation feature. The resulting triangulation produced a 10-meter resolution DEM with horizontal and vertical root mean square error (RMSE) of <2 m in all directions (Figure 3b). After the photogrammetric DEM construction, comparison to the lidar DEM revealed that the entire dataset was tilted. This is likely a product of minor misregistration across the whole DEM, which we corrected by fitting a planar ramp to the offset and 51 Figure 3. Local view of the assembly of the pre-eruptive DEM and thickness map for the Mauna Loa 1984 flow. (a) Orthorectified aerial photograph from a USGS survey in 1977 with an outline of the 1984 lava flow. (b) Hillshade of the 10-meter resolution photogrammetric DEM with the 1984 lava flow outline. (c) Hillshade of the 1-meter resolution lidar DEM with the 1984 lava flow extent shown in blue. (d) Thickness map produced by DEM-differencing (c minus b) within the extent of the lava flow. subtracting it from the DEM [Hollingsworth et al., 2013]. The photogrammetric DEM captures the upper surface (canopy) of vegetation, so where forest has been cleared for development or buried by new lava, the difference between the lidar and the photogrammetric DEM is strongly negative (Figure 4a and b). We corrected for this by identifying a band of well-fit pixels around the negative excursions and interpolating a planar surface between them, replacing the artificially high pixels in the photogrammetric DEM (Figure 4c). This does not fully remove the problem but greatly reduces the offsets 52 (Figure 4d). Other positive excursions from pixels resulted from clouds in the lidar and other errors. We corrected these and produced a final DEM by removing any pixels with values outside ± 3σ. We then used the difference between the lidar DEM and the photogrammetric DEM outside of the Mauna Loa 1984 lava flow to estimate the photogrammetric DEM error. The RMSE (vertical) of these areas is 2.36 m. In an unvegetated region near the Mauna Loa 1984 vent, the RMSE is 1.03 m, while a region within a cleared part of forest has an RMSE of 4.46 m. 2.2. Derived Datasets 2.2.1. Lidar Surface Properties Lidar elevation data form the basis of the morphologic analysis presented here. Additionally, we convert DEMs into maps of individual surface properties to highlight small-scale flow features and characteristics. For this purpose we use derived products of the lidar DEM that utilize neighborhoods of pixels (3x3 pixel region with the target pixel in the center) to calculate properties of the surface, such as slope, curvature, or roughness (Figure 5) to assist in feature identification. We discuss a variety of methods to process lidar data that allow extraction of quantitative morphometric information from within, between, and along lava flows. Slope and curvature (the first and second spatial derivatives) are frequent metrics used to investigate surface morphology. Slope (the maximum slope at a pixel in units of degrees) reveals steeper areas within a lava flow, such as channel or flow margins, in comparison to the shallow slopes along a channel or flow lobe [Ventura and Vilardo, 2008]. A slope map of the pre-eruptive surface beneath the lava flow reveals the path of 53 Figure 4. Demonstration of errors and corrections of the pre-eruptive photogrammetric DEM and resulting thickness map due to dense pockets of pre-eruptive vegetation. (a) Uncorrected photogrammetric DEM overlain on a 1977 orthophoto. Artificial “hills” are visible within the heavily vegetated areas. (b) Thickness (DEM-difference) map produced from subtracting the DEM in (a) from the lidar DEM. Where the flow buried vegetation, the “hill” artifacts produce strongly negative thicknesses outlined by a 0 meter thickness contour line. (c) Photogrammetric DEM that has been corrected by outlining the problem regions and fitting a plane to the pixels surrounding them in order to interpolate a new surface within them. (d) Thickness map resulting from the corrected photogrammetric DEM removes much of the error, while some negative regions remain, outlined by a 0 meter thickness contour line. Reoccupation of a pre-eruptive channel by the 1984 flow is also apparent in the comparison between the pre-eruptive imagery and the post-eruptive lidar. steepest descent, which is commonly used to predict flow paths [e.g., Favalli et al., 2005]. We calculate a slope map by measuring the local gradient of a surface fit to a moving 3x3 neighborhood of pixels (Figure 5a) [Horn, 1981]. Curvature provides 54 Figure 5. Example lidar surface properties maps within the extent of the classification analysis: (a) slope map in degrees, (b) total (surface) curvature, (c) surface roughness R parameter from McKean and Roering [2004] where high values represent smooth terrain, (d) natural log of the ratio of the first and second eigenvalues (ln(S1/S2)) showing clustering of vectors normal to the ground surface (larger values are more clustered and therefore smoother), (e) natural log of the ratio of the second and third eigenvalues (ln(S2/S3)) showing linearity of vectors normal to the ground surface (larger values are more linear or girdled), (f) natural log of the ratio of the first and third eigenvalues (ln(S1/S3)) showing combined clustering and lineation. In the roughness maps (c–f), blue is smooth and red is rough. All measures are produced from a resampled 2-meter lidar DEM to match the 2-meter intensity dataset. information on inflection points within a surface and whether a region is concave (positive curvature) or convex (negative curvature). In geomorphology, curvature is used to distinguish river channels (concave) from hillslopes (convex) and to calculate amounts of erosion and accumulation [e.g., Passalacqua et al., 2010]. Curvature of a lava flow 55 surface can highlight features such as levees [Tarquini et al., 2012]. We create a second derivative map of curvature by fitting a fourth-order polynomial surface to a pixel and its immediate neighbors (3x3) with Lagrange polynomials and then calculating the second derivative (Figure 5b) [Zevenbergen and Thorne, 1987]. Surface roughness is also valuable for differentiating features within and between lava flows; when quantified using orientation surfaces, surface roughness can be used to identify linear features, such as such as levees and surface ridges in landslides [McKean and Roering, 2004]. We test the effectiveness of four of these orientation-based measurements of surface roughness for analysis of lava flow surfaces features. All four parameters are defined by calculating unit vectors perpendicular to each pixel within a moving window (a 3x3 neighborhood). The parameter R is the length of a vector sum of all of the unit vectors in the window (Figure 5c). For example, the maximum value of R is 9 for a 3x3 window with nine pixels; this represents a perfectly flat, smooth surface with nine parallel unit vectors. Three other parameters can be calculated from the normalized eigenvalues of an orientation matrix of the nine pixels in the window (S1, S2, S3). The natural logs of ratios of these eigenvalues reflect clustering and alignment of the vectors [Figure 5 in McKean and Roering, 2004]. High values of ln(S1/S2) designate highly clustered vectors and a smooth surface (Figure 5d), while high ln(S2/S3) represents highly lineated or girdled vectors (Figure 5e). Ln(S1/S3) is a combination of both clustering and lineation (Figure 5f). Together these parameters can be used to identify channel levees, flow margins, and surface roughness. 56 2.2.2. Surface Properties Classification Surface property measurements can be used individually to highlight flow surface features, or in concert to map elements of an entire flow surface. We use an unsupervised classification of all surface layers (intensity, slope, curvature, and roughness) to extract flow margins and internal features with the ERDAS Imagine software. Unsupervised classification differentiates classes statistically by finding pixels with similar values across different layers. We specify the number of classes and initialize the class clusters along the principle axis at the start of each classification. The ERDAS software assigns pixels to the nearest cluster, and then calculates a new cluster mean and reassigns pixels iteratively until the pixel class assignments stabilize and each pixel is classified. We apply this technique to a small section of the northeast rift zone of Mauna Loa that contains older flows, an old cone, and a channeled segment of the Mauna Loa 1984 flow (Figures 2 and 5) and use two different sets of these input layers (all of the layers, and only the surface roughness eigenvalue ratio layers) and a range of specified numbers of classes. This small target area allows us to test that capabilities and shortfalls of automated classification of lava flows in a region with a range of features and flow ages. The limited spatial extent reduces the number of individual units, and therefore the number of classes required, and minimizes the effects of varying elevation (vegetation) and along-flow rheological changes that can influence surface characteristics within individual units across larger spatial extents. 57 2.2.3. Flow Thickness We have generated a thickness map for the Mauna Loa 1984 lava flow by calculating the vertical difference between the post-lidar DEM and the pre-eruptive photogrammetric DEM (Figure 3d). The different resolutions of the two DEMs (Figure 3b versus c) required us to resample the 1-meter (horizontal) lidar DEM to a DEM with a 10-meter resolution aligned with the photogrammetric DEM grid. The error in the thickness values can be estimated by propagating the error of both DEMs and taking the square root of the sum of the squares, which yields an approximate flow thickness error of 2.38 m. This assumes a uniform whole-flow error, which is not an ideal assumption, given the higher error of the photogrammetric DEM in vegetated areas. We therefore estimate a maximum error of 4.47 m using the photogrammetric DEM RMSE calculated among partially cleared vegetation, which encapsulates the most negative thickness values within the forest. A minimum error, from the RMSE of unvegetated regions, is 1.07 m. The thickness shows large variations across the flow that correspond to vents, levees, overflows, and the channel network. 3. Results and Discussion Our surface intensity, roughness, classification, and flow thickness results illustrate the types of data that can be derived from lidar and photogrammetric DEMs. The remarkable complexity and detail visible in the flows illustrates the value of using high resolution topographic data to relate flow morphology to the dynamics of flow emplacement. These data, individually and in concert, also provide insight into post- 58 emplacement biological colonization of flows [e.g. Deligne et al., 2012], and assist in mapping flows and flow features. 3.1. Intensity Analysis Intensity of lidar reflections has been used to map the relative (and perhaps absolute) ages of volcanic surfaces [e.g., Mazzarini et al., 2007]. On the Northeast Rift Zone (NERZ) of Mauna Loa, lidar intensities (Figure 2) are most strongly controlled by the abundance of low-lying vegetation (e.g., grasses, lichens, ferns). As vegetation on young lava flows is strongly controlled by elevation [Raich et al., 2000], lidar intensity can discriminate flows by age only when corrected for elevation. The NERZ of Mauna Loa is an ideal location to evaluate lidar intensity-age relationships because the lava flows boundaries and ages are well defined [Wolfe and Morris, 1996], as are the local climatic conditions that control revegetation rates on the flow surfaces [Aplet and Vitousek, 1994; Aplet et al., 1998; Raich et al., 1997]. The intensity map we have generated includes lava flows with ages between 30 and >5000 years [Wolfe and Morris, 1996] and elevations between 3100 m and 2000 m (Figure 6). Variations in lidar intensities are clear and closely match flow boundaries mapped by aerial photographs and field studies. Intensities also vary within flows as a function of elevation (Figure 7); all flows have higher intensities at lower elevations. These data match measurements of changing biological productivity during primary succession on flow surfaces [e.g., Aplet and Vitousek, 1994], which mirrors patterns of decreasing temperature and precipitation with increasing elevation along the NERZ (Figure 7 inset). Primary succession on new lava surfaces initiates with the lichen Stereocaulan vulcani, 59 which can become established on newly erupted lava flows within a few years [e.g., Jackson, 1971]. The productivity of this, and other, early colonizers is amplified by increasing temperature and to a lesser extent rainfall [Anderson-Teixeira et al., 2008; Raich et al., 2000]. Figure 6. Map of normalized lidar intensities on Mauna Loa NERZ with mapped geologic units outlined [Wolfe and Morris, 1996]. This can be compared to Figure 1c for flow ages. We calculate an intensity value for individual flows and flow units (including older flows preserved in kīpukas within younger flows) by determining the mode of the distribution of intensity values from within the mapped flow boundaries. When viewed as 60 Figure 7. Intensity increases with decreasing elevation within individual flows (i.e. of constant age). This increase is consistent with increasing biological productivity due to climatic gradients along the NERZ of Mauna Loa including air temperature and precipitation that are shown in the inset (after Aplet and Vitousek [1994]). The rate of intensity change is greater for older flows, perhaps due to changes in plant types that are present at different stages of succession. Linear fits to the data suggest that the intensity difference between flows of different ages will be small or negligible at high elevations. a whole, the age-intensity relationships of lava flows within the study area show increasing intensity with age [e.g., Mazzarini et al., 2007], albeit with significant scatter (Figure 8). We can fit the mean intensities of each age group using a (natural) logarithmic function, which implies a rapid increase in intensity within the first few hundred years, and slower increase over the next few thousand (Figure 8). We attribute the scatter in the data to the fact that flows within any age group traverse an elevation range along the rift 0 200 400 600 800 1000 0 10 20 30 40 50 60 Elevation below 3100 m N or m al iz ed in te ns ity 1984 200−750 ybp 1500−3000 ybp 6 4 2 0 Pr ec ip ita tio n (cm )30 20 10 0 0100020003000 Te m pe ra tu re (˚C ) Altitude (m.a.s.l.) 61 zone. Using the linear fits to the elevation-intensity data determined for individual flows (Figure 7), we construct additional curves that would reflect age-intensity relationships at constant elevations of 2000 and 3000 m (Figure 8). We note that the fit to the mean intensities of each age group is nearly identical to the age-intensity relationship at 2500 m, the approximate mean elevation within the data. The logarithmic form of age-intensity relationship, and flattened form of that relationship at high elevations, mimic the simulated form of plant productivity during primary succession on windward Mauna Loa slopes [Raich et al., 2000], where initial primary productivity declines with increasing elevation. Our results support the hypothesis that biological colonization is the primary driver of age-intensity variations [Mazzarini et al., 2007] and suggest that absolute ages could be determined from lidar intensity of volcanic surfaces where local climatic and biological colonization patterns are known. Alternatively, lidar intensity could be used to map vegetation variations with elevation in locations where lava flow ages are well constrained. The fidelity of these models will be greatest in the steeper regions of the fit curves (i.e., at young ages), the age range of which will depend on rates of biological productivity (elevation). 3.2. Surface Roughness Analysis We use the December 1974 lava flow of Kīlauea [Lockwood et al., 1999] to test methods of quantifying the fine-scale roughness resolved by the ALSM data. Previous satellite-based measured of roughness [Gaddis et al., 1990] and mapped variations in the proportions of pāhoehoe and ‘a‘ā surface morphologies [Soule et al., 2004] of this flow 62 Figure 8. Age-intensity relationships for lava flows on the NERZ of Mauna Loa show an increase in intensity with age. Lava flow ages are from Wolfe and Morris [1996]. The data can be fit by a logarithmic function defining an initial rapid increase in intensity over a few hundred years that slows over the next several thousand years. The solid line is a fit to the mean of each age group (flows <200 years are grouped with a mean age of 78 years). Significant scatter in the data is due to the fact that intensities are calculated for flows that span a large range in elevation, over which climatic and thus biological productivity during primary succession will be highly variable. Additional fits to the data are shown for altitudes of 2000 and 3000 m by dashed lines, constrained by fits to elevation-intensity relationships within single flows that span large elevation ranges (Figure 7), explain much of the variability within each age group. allow us to compare, directly, the results of the lidar analysis with those both satellite- and ground-based measurements. We measure roughness using the method of vector dispersion from McKean and Roering [2004] as described above. The calculations are made over a 3m x 3m window, thus are tuned to identify roughness on 1–2 meter length scales. We determined 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 10 20 30 40 50 60 70 80 Age (years before present) N or m al iz ed in te ns ity ln fit to group means & 2500 m ln fit to 200 0 m ln fit to 3000 m individ. flow intensity mean of age group int. @ 2000 m.a.s.l. int. @ 2500 m.a.s.l. int. @ 3000 m.a.s.l. 63 roughness along the ~11.5 km flow using ~1 km (along flow) bins in order to resolve changes in morphology that occur over length scales of hundreds of meters. This bin spacing is the same as used in field-based interpretations of flow morphology [e.g., Soule et al., 2004], which facilitates comparison between the two methods. Bin boundaries are dictated by spacing along the presumed flow pathway and orientations that are normal to the dominant flow direction, although flow directions are more variable in the near-vent region (Figure 9a). Given the non-Gaussian distribution, the roughness within individual bins is determined by the mode of the ln(S2/S3) measurements within each bin. Of the four roughness parameters, we find that ln(S2/S3), which identifies lineations, and ln(S1/S3), which combines linear and clustered features, provide the most useful results. Overall, ln(S2/S3) highlights a down-flow sigmoidal shift from low roughness in the area mapped as entirely pāhoehoe to high roughness in the area mapped as entirely ‘a‘ā (Figure 9b), consistent in form with that obtained from field and aerial photographs [Soule et al., 2004]. Relative to that study, however, the inflection points marking the onset and conclusion of the pāhoehoe-to-‘a‘ā transition occur ~ 1 km closer to the fissure vents. This offset might reflect variations in the vertical scale of roughness along the flow from changes in the ‘a‘ā clast size [Gaddis et al., 1990; Soule et al., 2004], or subtle changes in pāhoehoe surface textures that were not captured in the binary pāhoehoe-‘a‘ā analysis of Soule et al. [2004]. Despite these discrepancies, the common sigmoidal form of the mapped pāhoehoe-to-‘a‘ā transition indicates that the ALSM data and roughness metrics can be used to evaluate flow-scale surface morphology as well as the characteristic roughness of end member morphologies. 64 Figure 9. Derived roughness map and down-flow roughness for the December 1974 lava flow on the SW Rift Zone of Kīlauea. (a) Map of roughness and bins. (b) Plot of vector dispersion roughness metric, ln(S2/S3), for bins along the December 1974 flow of Kīlauea. A fit to the relative proportions of pāhoehoe and ‘a‘ā to observations from airphoto and field-based mapping of the flow is shown by the solid line. A modification to the fit line to better fit the roughness measurements is shown by the dashed line. 3.3. Classification of Lidar Surfaces As demonstrated above, lidar elevation and intensity data can be combined to discriminate flows of different ages and surface textures. Here we use all of these input layers to test their collective ability to map a section of the Mauna Loa 1984 flow and surrounding area using simple automated classification (Figure 2). We use two different sets of input layers - all layers, and only the surface roughness eigenvalue ratio layers (Figure 5). Of the specified layers, we find that intensity and slope have the greatest variability between different features in area, and so are most useful for classification. Both the number of layers and the number of assigned classes determine the degree of detail represented by the classification. When we use all of the layers, an eight- class result captures most of the different units and features in the image (Figure 10a; Table 1). This particular classification can be used to identify differences in flow age as 0 2 4 6 8 10 12 1.7 1.8 1.9 2.0 2.1 2.2 Distance (km) ln (S 2/S 3) 2.3 1.0 0.8 0.6 0.4 0.2 0 R el at iv e fra ct io n of p āh oe ho e/ ‘a ‘ā 0.5 1.0 1.5 2.0 2.5 3.0 252000 21 44 00 21 40 00 21 36 00 256000 260000 baSurface roughness (ln(S2/S3)) 65 well as the channels within individual flows. Importantly, half of the classes (1, 3, 6, and 8) are distinct (unmixed) in the features that they include, while the other half range from somewhat mixed to random. Reclassifying the dataset into only four classes yields simpler results and still captures differences in surface age, but retains some mixed classes (Figure 10b; Table 1). As a result, some of the noise is reduced but the channels are not mapped as continuous features because the levees are partially grouped with the rest of the flow. In the four-class analysis, classes three and four each make up half the flow surface and could easily be combined, and yet they are not merged when a three- class analysis is run. Finally, requiring only two classes creates a map of the prehistoric surface (class 1) relative to the historic lava flows (class 2; Figure 10c). Running the unsupervised classification with only the eigenvalue ratios as input layers produces a classification that discriminates among only the flow or cone features and not the flow ages. Allowing multiple classes increases the scatter within flow segments outside of the main channels (Figure 11a). Removal of all classes except the channels (class 1) and the cone (class 8), however, can generate a map that discriminates features from different textures (Figure 11b). The results presented above demonstrate that unsupervised (automated) classification of lidar-derived surface metrics can be used to map distinctive features such as lava flows of different ages, cones, and lava channels. Here we have focused on a small part of the lidar image. Although the same statistical descriptions of each class (the class signature) produced could be used in a supervised classification of a larger area, several of the parameters are sensitive to elevation or distance from the eruptive vent. This sensitivity complicates development of automated classification protocols for the 66 Figure 10. All-layer classification results for (a) eight classes, (b) four classes, and (c) two classes. 67 Table 1. All layer unsupervised classification interpretation Class Description Eight classes 1 Very bright areas with the 1500–3000 year old units 2 The rest of the 1500–3000 year old units 3 Inner levee slopes 4 Scattered pixels within 1942 and 1984 lava flows 5 200–750 year old flow, 1942 channel, and overflows within the 1984 flow 6 Majority of the non-channel 1984 and 1942 flow surface 7 Scattered pixels within 1942 and 1984 lava flows, some concentration in the channels 8 Outer side of levees Four classes 1 1500–3000 year old units 2 200–750 year old flow, the 1942 channel, and the rest of the 1500–3000 year old units 3 Half the surface of the 1942 and 1984 lava flows 4 Half the surface of the 1942 and 1984 lava flows entire flow. Addition of input layers, such as high-resolution color, radar or infrared imagery, could improve classification results. Radar amplitude, in particular, measures centimeter scale surface roughness that could facilitate differentiation of flow ages and fine-scale changes in texture [e.g., Gaddis et al., 1990]. Remote sensing of lava flows at higher spatial, spectral, and temporal resolutions will continue to improve our quantification of flow morphology and understanding of flow emplacement. 3.4. Flow Thickness Analysis Our thickness map of the Mauna Loa 1984 flow provides important new insight into both the overall development of the flow field and the distribution of lava within the flow (Figure 12). Analysis of these data is enhanced by detailed observations, made both during and after the 1984 eruption, related to the evolution of lava channels and spatial and temporal changes in lava rheology [e.g., Lipman and Banks, 1987; Crisp et al., 1994]. 68 Figure 11. Surface roughness ratio classification results for eight classes showing (a) all of the classes, or (b) just the two classes that show channels and the cone. 3.4.1. Flow Volume Although our lidar dataset is limited to the upper portion of the Mauna Loa 1984 flow field, it captures much of the main flow lobe (Flow 1 in Lipman and Banks, [1987]), and all of the southern lobes that were active for the first 48 hours of the eruption (Flows 2–4 in Lipman and Banks [1987]; Figure 1). We calculate a total volume of the lidar- 69 Figure 12. Flow thickness map with lidar hillshade of the Mauna Loa 1984 flow within the extent of the lidar data. Boxes indicate the extents of other figures. imaged 1984 lava flow by summing the flow thickness values of each pixel within the flow and multiplying by pixel area (100 m2). A maximum error can be estimated by multiplying the thickness error by the flow area. A minimum estimate assumes uncorrelated errors and is calculated by dividing the maximum error by the square root of the number of pixels to calculate an error range [Favalli et al., 2010]. The resulting volume of 0.894 (± 0.001–0.489) ×108 m3 accounts for approximately 45% of the ~1.98×108 m3 total volume of the flow estimated by Lockwood et al. [1985]. The area of the flow within the lidar extent is ~67% of the total flow area, which requires the more 70 distal section of the flow to be significantly thicker, on average, than the upper two- thirds. In fact, the volume estimate of Lockwood et al. [1985] would require the distal section to have an average flow thickness of almost eleven meters. If, instead, we apply our average thickness of 4.35±0.01(std error) m to the final third of the flow, we obtain a new total volume of 1.34 (± 0.001–0.733) ×108 m3. This is value is 33% smaller than the estimate made at the time of eruption and represents a minimum volume, as flow thickness typically increases with distance from the vent because of increasing viscosity [Moore, 1987; Lister, 1992]. Although the uncertainty in these estimates could be large because of the error in the pre-eruptive DEM, the high spatial resolution of the DEM differencing improves volume estimates by accounting for local thickness variations. We note that part of the discrepancy in volume estimates may be accounted for by volume reduction due to post-emplacement contraction and settling [e.g., Stevens et al., 2001; Lu et al., 2005; Dietterich et al., 2012]. The thickness map illustrates the inhomogeneous distribution of lava volume within the flow; it also demonstrates that the flow is thickest in regions of major flow branching. To a first order, flow lobes of longer duration, which are subjected to repeated channel overflows and levee development, are also thicker [Lipman and Banks, 1987; Coltelli et al., 2007; Branca et al., 2013]. From this perspective, the Mauna Loa 1984 lava flow can be divided into two parts. The first is the portion of the channel that transported lava for the entire three weeks of the eruption, including Flow 1 and 1A and their subsidiary branches, Flow 1B and what we term Flow 1 north lobe and Flow 1 east lobe (Flow 1 + branches; Figure 13). This part of the flow contains 0.519 (± 0.001– 0.281)×108 m3, or ~58% of the total volume within the lidar image. The other three main 71 flow branches, Flows 2–4, plus the shield developed over the lower vents, account for the remaining 42% of the flow volume within the lidar extent (individual lobe volumes are shown in Figure 13). These southern branches started at the same time as Flow 1, but were active for only the first 48 hours of the eruption. Their combined volume of 0.375 (± 0.001–0.208)×108 m3 highlights the large volume that went into these branches to the south, instead of along the dominant channels of Flow 1, during the early (high effusion rate) phase of the eruption. Figure 13. Volume distribution among flow lobes. Colors reference the inset map of the different major flows. Flows 1 and 1A Flow 1B East lobe North lobe Lower vents Flow 3 Flow 4Flow 2 0 0.5 1 1.5 2 2.5 3 3.5 Volumes by flow lobe Flow 2 Flow 3 Flo w 1 + b ran che sTotal volume: 8.94 x 107 m3 Volume in Flow 1+ 5.19 x 107 m3 Volume in Flows 2-4: 3.75 x 107 m3 Flow lobe Vo lu m e (×1 07 m 3 ) 35.5% 7.74% 9.51% 5.39% 5.56% 14.1% 18.2% 4.09% Flow 4 72 The four primary flow lobes are separated by topography, with each occupying a different lavashed (equivalent to a watershed) delineated by drainage divides [Kauahikaua et al., 2003]. Confinement of individual flow lobes within single lavasheds confirms the use of lavasheds as a simple hazard mapping tool. We can also use the known flow lobe durations (48 hours) to calculate time-averaged discharge rates for the each flow component: the respective average fluxes for Flows 2, 3, and 4, are 73 m3/s, 94 m3/s, and 21 m3/s. These make up 68% of the ~275 m3/s instantaneous effusion rate at the vents measured during the emplacement of these lobes, which requires that Flow 1 account for the remaining 32%. No detailed observations of these southern flow lobes were made at the time of the eruption, but our estimated time-averaged discharge rates are similar to flow rates measured at distal sites of Flow 1 and 1A during the eruption [Lipman and Banks, 1987]. The three-fold variation in flux most likely reflects the initial flux from the fissure segment feeding each lobe, the partitioning of lava between different flow branches, and the actual duration of each flow branch (we have assumed a constant value of 48 hours for all three). Understanding the differences in flux between branches is crucial for determining both rates of flow advance and final flow length; for this reason, collecting syn-eruptive multitemporal topographic data could provide important input data for future eruptions. Branching lobes with different fluxes have also been observed in other locations such as the 1669 flow field at Mt. Etna [Branca et al., 2013]. Another way to view volume distribution is to measure the cumulative volume in the downflow direction, which shows where lava is being stored within the flow field (Figure 14). We have done this by measuring the cumulative volume upslope of generalized flow cross-sections spaced every kilometer along the flow path. We plot 73 separately curves for the entire flow (black), Flow 1 and its associated branches (blue), and Flows 2–4 (shown in red). The constant slope of the whole-flow line shows that the volume within each one-kilometer bin is approximately constant, and the flow (all branches) therefore has a consistent total cross-sectional area. In contrast, the curve showing Flow 1 (blue) documents lava accumulation both near the vent (within the first few kilometers) and where the flow branches extensively and widens at a distance of 12 km from the vent. The red curve emphasizes the importance of Flows 2–4 as major (early) components of the flow volume. Dashed lines show extrapolation of the whole- flow volume distribution curve to the total volume of the flow estimated by Lockwood et al. [1985] and our total calculated from the average flow thickness and total inundated area. A final estimate of the total volume of the flow can be made from our observation that the cumulative volume curve for the whole flow has a constant slope. If this slope is extrapolated with a linear fit, we find a total flow volume of ~1.53×108 m3, which falls between the previous estimates. This value represents an average flow thickness of ~6.3 m over the final third of the flow, which is within the range of distal flow thicknesses measured by Lipman and Banks [1987]. 3.4.2. Thickness of Lava Flow Features Within the Mauna Loa 1984 flow, differences in flow thickness record the details of flow emplacement. The thickest portions of the flow correspond to the eruptive vents (≤ 32 m), flow levees (≤ 17 m), and major channel overflows and junctions where the flow splits (≤ 24 m; Figure 12). Flow thickness also varies with the degree of channel development, described by Lipman and Banks [1987] as the stabilized channel, 74 Figure 14. Distribution of volume with distance along the flow. Different lines show the volume distribution for the whole flow (black), Flow 1 and all of its branches (blue), and Flows 2–4 (red). Dashed lines extend to the total volume reported in Lockwood et al. [1985] (short dashes), our linear fit extrapolation of total volume (dash-dot), and our estimate of total volume calculated by extending the average flow thickness in the first 15 km of the flow to the total flow area (long dashes). The blue line could also be extrapolated to the maximum flow length with a set of dashed lines parallel the whole- flow ones. transitional and dispersed flow zones, and the flow toe. In total, these patterns of flow thickness reveal much about the evolution of the flow field. Locally thick regions of the flow show both the development and failure of lava channels. The main (continuously active) channel of Flow 1, which makes up much of the narrow, central portion of the lidar-imaged flow, is characterized by distinct levees that form the margins of a confined channel that is only 10–20 m wide (Figure 15a). This section of the channel represents the stabilized channel zone, where flow was restricted to the central channel conduit [Lipman and Banks, 1987]. Here the channel is so narrow that Lockwood et al. (1985) Total volume estimate Whole flow Flow 1 Main channels 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Distance along flow (km) Cu m ul at ive v ol um e (×1 08 m 3 ) Only Flow 1+ branches Only Flows 2 - 4 Minimum total volume estimate from average thickness and total area Linear extrapolation 75 the 10-meter resolution of the thickness map is barely enough to capture the thin channel between the thick levees. A profile line across this section of the channel (A-A’) shows pre- and post-eruptive elevations, as well as the flow thickness from the thickness map. At 15 meters height, the levees in this segment are some of the tallest in the whole flow, probably because the longevity of this part of the channel allowed prolonged levee construction [Lipman and Banks, 1987]. Figure 15. Flow thickness maps with lidar hillshade overlay of (a) the stable channel of Flow 1 with a large overflow and (b) the distal portion of Flow 2, showing the broad channel of the transitional channel zone and the flow toe. Profile lines of the main channel, overflow, transitional channel zone, and flow toe are also shown, with the pre- and post-eruptive topography in black, and the flow thickness in green. Channel overflows that caused local flow thickening were common throughout the Mauna Loa 1984 eruption and formed both at locations where the channel was blocked and during times of temporary increases in lava flux (surges). For example, the overflow shown in Figure 15a increased the flow thickness by 1–3 m (see profile B-B’ in 76 comparison to A-A’). Here the thickest part of the flow lies at the downslope edge of the junction where the overflow splits from the main channel and probably reflects the combined effect of the overflow volume, ponding of lava in the channel prior to overflow, and material added by run-up of lava on the central barrier as the flow divided. Locally increased thicknesses from overflows have also been documented in the 2001 and 2006 lava flows at Mt. Etna [Coltelli et al., 2007; Favalli et al., 2010]. Overflows (channel breakouts) redirect much of the lava flowing down the channel. For example, field observations made along the Flow 1/1A channel during a 24- hour period show a 4-fold decrease in the volumetric flux in the channel with distance from the vent (Figure 16) [Lipman and Banks, 1987]. The thickness data show that much of this missing flux was diverted to levees, overflows, and branches. From this perspective, thickness maps provide important information about changes in lava flux along individual flow lobes, information that is key to accurate models of lava flow emplacement (and is not fully accounted for in existing flow models). Flow thickness data can also inform our understanding of the more distal portions of a lava flow. Down flow of the stabilized channel zone is a broad transitional channel zone that leads to a zone of dispersed flow and flow toe with no central channel. These features, originally defined by observations of advancing flows during this eruption, are nicely preserved in Flows 2–4, whose full lengths are captured by the lidar. A map and set of profiles from Flow 2 shows the wide channel and marginal lenses and ridges of the transitional channel zone and the convex shape of the flow toe (Figure 15b). Also evident are multiple small lobes at the flow toe, a feature that is common in lava flows at Mt. 77 Figure 16. Decrease in flux down the main channel of Flow 1/1A with distance from the vent as recorded on April 3–4, 1984. Data from Lipman and Banks [1987]. Etna [Applegarth et al., 2010; Favalli et al., 2010]. The lobe directly downslope from the central channel traveled the furthest, as also seen at Mt. Etna. 3.4.3. Influence of Pre-Eruptive Topography on Flow Thickness Reconstruction of the pre-eruptive surface also allows us to investigate the influence of pre-eruptive topography on flow thickness. Cross-slope profiles (Figure 15) show that local topography affects both the path of the flow and asymmetry in flow thickness. All flow profiles show a dip in the cross-slope pre-eruptive topography, which shows that individual flow lobes followed pre-existing depressions. Interestingly, in many places the flow appears to have traveled at an angle to the slope. This appears to have affected the flow thickness, in that the downslope side of the flow is thicker than the upslope side (e.g., A-A’). The overflow in Figure 15a also occurs on the downslope side. 0 5 10 15 200 200 400 600 800 1000 Distance from vent (km) Vo lu m et ric fl ux (m 3 /s ) Data from April 3-4, 1984 Apparent DRE 78 This suggests that channels traveling across the slope may thicken preferentially, and overtop, in the downslope direction. On a more fundamental level, viscous flow theory predicts that the slope of the pre-existing surface is a major control on flow thickness, such that flows thicken on low slopes and thin on high slopes [e.g., Lister, 1992; Griffiths, 2000]. As described above, when flows traverse irregular topography, the effect of slope may be overprinted by levee formation and overflows. To investigate the influence of slope on flow thickness we track the thickness of the flow in a channel within a single, short-lived lobe (Flow 2). By restricting our analysis to a short-lived flow lobe we remove complications related to levee development and variations in lava flux. We calculate the local slope by measuring the maximum slope of a local surface approximated by a neighborhood of pixels around the pixel of interest [Horn, 1981] (Figure 17). Although this measure of pre-eruptive slope is noisy on the pixel-scale (Figure 17a), the length scale of the slope that controls flow thickness is likely to be larger than the pixel size. To evaluate this hypothesis, we use the Soil Land Inference Model (SoLIM) software (University of Wisconsin-Madison) to change the neighborhood size used for the slope calculation [Zhu et al., 2008]. Correlations between slope and flow thickness are weak (Figure 17b), but peak at a neighborhood size of 70–90 meters (Figure 17c). This corresponds to approximately twice the flow width and suggests an upper limit to the DEM resolution required for accurate flow modeling. One portion of the channel is omitted from the analysis due to a pile-up of lava in the channel that did not drain. This is probably the origin of the overall general scatter in the plot, as different channel segments clearly drained different amounts. 79 Figure 17. Comparison of lava flow thickness, represented along a channel in Flow 2, to underlying slope at a variety of length scales. (a) Plot of slope along the channel showing the effects of using a different neighborhood size. (b) Example linear regression of slope and flow thickness for a 90 m neighborhood slope calculation. (c) Correlation, measured by the r2 of linear fits, as in (b), for a range of neighborhood sizes. Two other parameters important for determining the flow thickness are flux and viscosity. In most eruptions, volumetric flux from the vent decreases with time. Additionally, the complexity of the channel network within the flow causes down flow variations in flux within and between individual flow lobes. Viscosity also increases with distance because of increasing crystallinity and decreasing vesicularity [Moore, 1987; Riker et al., 2009], although in high effusion rate flows from Mauna Loa, substantial viscosity increases occur only after flow transport distances of several kilometers. 80 Our analysis of the flow thickness map highlights variations in flux and volume distribution within flow lobes and between flow branches, variations that will control flow advance rates and final flow lengths [Kauahikaua et al. 2003]. Our data also reveal information about flow routing and conditions under which flows traverse topographic lavasheds. We show that flows travel down depressions, even when they cut across slopes, which can lead to overflows and thickening along the downslope side of the channel. The length scale over which slope influences flow thickness also reveals the necessary length scale, and therefore resolution, for modeling flow thickness changes and resulting flow behavior. We anticipate that increased availability of high resolution DEMs will make routine the creation of high quality flow thickness data that will permit more detailed comparisons of theoretical and actual conditions of lava flow emplacement. 4. Conclusions Our use of lidar and photogrammetric topographic data to investigate the emplacement history of Hawaiian lava flows suggests new methods to quantify flow morphology and surface properties, and to interpret the evolution and style of flow emplacement from 3D flow data. Importantly, high-resolution lidar topographic and intensity data allow unbiased analysis of surface characteristics that can be used to track the biological colonization of lava flows and to date flow surfaces. Surface roughness can also be used to detect morphological changes within lava flows, particularly the transition from ‘a‘ā to pāhoehoe that is a characteristic indicator of down flow rheological changes in Hawaiian lava flows. Together, these surface measures, combined with those of slope and curvature, offer an unbiased and unsupervised means of identifying lava flow 81 features. The results of our classification provide statistics on the surface properties of flows and cones of separate ages and channel features that suggest potential for automated flow mapping. When integrated with pre-eruptive imagery to build a photogrammetric DEM, the lidar data constrain measurements of flow thickness, adding a third dimension to our view of lava flow morphology. Analyses of volume distribution and local flow thickening reveal where lava was stored within the flow, emphasize the importance of side branches, and explain measured flux decreases along the Mauna Loa 1984 channel. Examination of the pre-eruptive surface with respect to the flow route and flow thickness completes our analysis, exposing some of the influence of pre-eruptive topography on flow emplacement and morphology. Taken together, these data illustrate the power of integrating pre- and post-emplacement DEMs with field observations of emplacement chronology, effusion rates, and viscosities; our analysis demonstrates that such integration can provide significant insight into the dynamics of lava flow emplacement. This work expands on new techniques and analyses that have been developed for Italian volcanoes by combining both surface and 3D analysis to address questions related to the emplacement of two Hawaiian lava flows. Detailed observations of Hawaiian lava flows, including work on flow surface morphologies and revegetation studies, provide necessary baseline data for our analysis. Our work also extends previous analyses by pioneering the 3D study of lava flow morphology in Hawai‘i, and at the same time provides a parallel and well-constrained dataset to compare flow emplacement styles at Mt. Etna and Hawai‘i. An important result of our work is demonstrating the effect of pre- existing topography on Hawaiian lava, with associated implications for predictive models 82 of future flow paths. Taken together, our combined analysis of lava flow surface characteristics and volume distribution shows the potential of such data for improving our understanding of lava flow emplacement and informing efforts to mitigate the hazards posed by effusive eruptions. 5. Bridge In Chapter III, I developed tools to extract 3D flow morphology data from lidar and photogrammetric DEMs. Lava flow intensity data provided information about flow age and revegetation. Surface roughness revealed intra-flow morphology, including vents, channels, and the distribution of ‘a‘ā and pāhoehoe. Importantly, flow thickness and volume analysis exposed the highly uneven distribution of lava within the flow field, highlighting thickening where the flow splits and how variations in the underlying slope influenced the local flow thickness. These data demonstrate the large portion of lava that flowed into side branches, and thus how the flow branching geometry influences the lava flux distribution. The next chapter draws directly from this work to further explore this observation of flow branching, as well as merging. In Chapter IV, I seek to investigate what creates the branches and confluences that make up the flow and channel network geometry, and how do they affect the flow. I employ the datasets from Chapter III to examine the relationships between channel network geometry, the pre-eruptive eruptive surface, and the flow and channel segment morphology. Channel maps, pre-eruptive slope maps, and flow thickness data derived from the analysis in Chapter III are all used in Chapter IV to 83 test what controls the formation of channel networks, and how the network geometry in turn affects the flow morphology and behavior. 84 CHAPTER IV CHANNEL NETWORKS WITHIN LAVA FLOWS: FORMATION, EVOLUTION, AND IMPLICATIONS FOR FLOW BEHAVIOR This chapter is in press for the Journal of Geophysical Research: Earth Surface. I was the first author, and the paper presents my data, methods, analysis, interpretations, and writing, with supervision and editorial assistance from my co-author and advisor, Kathy Cashman (University of Oregon). 1. Introduction Lava flows surface much of the Earth, including numerous volcanic regions and the seafloor, but their dynamics remain poorly understood. Hawaiian lava flows erupt as crystal-poor melts and construct channels as they flow, cool, and solidify. In Hawai‘i, lava flows can be classified into two main flow types: pāhoehoe, which advances as a succession of small lobes, and ‘a‘ā, which has a distributed rubbly flow front behind which levees develop and the flow channelizes. ‘A‘ā flows commonly form from multiple eruptive vents, bifurcate at the flow front, generate overflows and breakouts from existing channels (equivalent to river avulsions), and produce flow confluences where channels merge. The processes that produce these networks, and their effect on how fast or far flows travel, are not well known. Prior studies of lava flows have identified eruption rate, underlying slope, and viscosity as key parameters that influence both the rate of flow advance and surface 85 morphology [e.g., Walker, 1973; Cashman et al., 1999; Kauahikaua et al., 2003]. Here we assess the additional importance of an often-overlooked feature of Hawaiian ‘a‘ā flows, which is the complexity of lava channel networks. Specifically we examine controls on network formation and show ways in which the topology of the network can dramatically influence the flow behavior. This work has important implications for interpreting emplacement conditions of past lava flows, including those on other planets, and for constructing effective barriers to protect property from future lava flows. High resolution lidar topography for two historical lava flows – the Mauna Loa 1984 flow (ML84) and the Kīlauea December 1974 flow (KDec74) – reveals the complexity of Hawaiian lava channel networks (Figure 1). These lava flows represent endmembers of a spectrum of possible channel network geometries that ranges from distributary (ML84, Figure 1a) to tributary (KDec74, Figure 1b). The distributary ML84 flow was emplaced over 3 weeks in March–April 1984. The eruption initiated with four main flow branches that had a combined effusion rate of ~275 m3/s (Flows 1 to 4 in Figure 1a). Flows 2 to 4 ceased activity within 48 hours, thus preserving their initial lava channel networks without the overprinting of later overflows and breakouts. Flow 1 continued to advance, reaching a final length of 27 km within four days from the start of the eruption. The channel network within Flow 1 then continued to evolve as blockages caused overflows that evolved into new flow lobes (e.g., Flow 1A, 1B). Activity finally ceased after twenty days [Lockwood et al., 1985; Lipman and Banks, 1987]. The tributary KDec74 flow was emplaced over only ~6 hours. Lava erupted from a series of en echelon fissures that opened from NE to SW (Figure 1b). With an average effusion rate of ~270 m3/s [Lockwood et al., 1999], the lava from the fissure system flowed downhill to the 86 south until it was confined against the Koa‘e fault scarps, which forced the discrete flow branches to merge and redirected the flow to the southwest. This extreme example of topographic confinement created a flow with more confluences than bifurcations, yielding a tributary network similar to that of a river. Figure 1. Locations, lidar hillshade maps (illumination from the northwest), and traced channel networks for two lava flows on the Island of Hawai‘i, (a) the Mauna Loa 1984 (ML84) flow and (b) the Kīlauea December 1974 (KDec74) flow. Other recent lava flows at Mauna Loa and Kīlauea mentioned in this paper are shown in the inset location map, with all flow outlines from Wolfe and Morris [1996]. In (a) the channel network is dominantly distributary, branching to either side of the Mauna Loa 1942 flow shown for reference in green, while in (b) the en echelon fissure vents feed a broad flow field that is confined into a tributary system by the north-facing scarps of the Koa‘e Fault Zone. The presence of channel networks in lava flows has implications for predictive models that are essential for estimating where, how far, and how fast lava flows travel. Predictive models used in Hawai‘i have included maps of topographically confined lava sheds derived from hydrologic analysis [Kauahikaua et al., 2003], flow paths calculated from steepest descent routes over a stochastically perturbed surface (DOWNFLOW) [Favalli et al., 2005], and a thermo-rheological model for channelized lava flow 87 (FLOWGO) [Harris and Rowland, 2001]. However, the steepest descent paths used in the lava shed maps and DOWNFLOW code represent single flow pathways, and the FLOWGO model predicts the behavior of flow through only a single channel. The models therefore do not account for the complex channel networks that are ubiquitous in Hawaiian lava flows. 2. Background Lava flows are similar to flows of water, sediments, and ice in that they are gravity-driven flows that are strongly influenced by topography. They can travel at velocities comparable to debris flows and create their own levees by lateral movement of material at the flow front, forming lateral ridges akin to those that develop in debris flows, earthflows, and pyroclastic flows. However, lava flow rheology is regulated by the bulk viscosity of a crystallizing fluid, while these other geologic flows have rheologies dominated by particle interaction within fluidized granular mixtures or the shearing and basal sliding of colluvium, weathered bedrock, or ice [e.g., Crisp et al., 1994; Iverson, 1997]. Like glaciers, lava flows undergo phase transitions during flow; lava, however, has multiple components (ice has one) and solidification occurs over a wide temperature range (~200°C instead of a single pressure-dependent temperature). The crystal content of a lava flow affects the rheology, somewhat analogous to the particulate content of a debris flow. For this reason, down flow variations in cooling and crystallization provide the dominant control on lava flow emplacement. Because the formation and evolution of lava flow channel networks have not previously been studied, we use observations of networks in other geologic flows to 88 guide our analysis. The tributary and distributary networks that we have identified in the lava flows are common features of river networks. In river systems, the topology of water drainages creates tributary networks made up largely of confluences [e.g., Shreve, 1966]. Flow branching in rivers occurs in depositional environments (i.e., low slopes and deltas), where sedimentation widens and elevates existing channels above their surroundings, leading to avulsions and the formation of new flow branches at lower elevations [e.g., Jerolmack and Mohrig, 2007]. Branching is associated with low flow velocities because of reductions in hydraulic radius and slope [e.g., Makaske et al., 2009], which creates a feedback loop involving sedimentation, velocity reduction, and branching [Pelletier and DeLong, 2004]. These processes are driven by sedimentation, which is not an important process in lava flows. However, the methods used to describe the topology of river networks, and the fundamental controls on branching, including slope and avulsions, can inform our study. In lava flows, branching can be initiated at the flow front or caused by a breakout or overflow from an existing channel, comparable to a river avulsion. Flow front bifurcations that form the primary channel network often record interaction with topographic obstacles. Instabilities such as viscous fingering [e.g., Huppert, 1982] are probably responsible for the scalloped margins of lava flows, but their origins or contribution to larger-scale branching are unknown [Bruno et al., 1992]. Breakouts (levee failures) and overflows produce secondary bifurcations and can be caused by blockages within an existing channel [Lipman and Banks, 1987] or pulses in lava supply to the flow front [Tarquini and Vitturi, 2014]. During initial flow emplacement (e.g., ML84 Flows 2–4 or all of KDec74), we assume that the channel network is produced by flow front 89 processes. Channel networks evolve through time by breakouts and overflows as illustrated by ML84 Flows 1A and 1B [Lipman and Banks, 1987]. Therefore, we examine early flows to learn about the initial formation of channel networks, and channels with prolonged activity to explore network evolution. These processes and their effects on flow length and advance rate are fundamentally regulated by lava flow physics. 2.1. Controls on Lava Flow Emplacement Investigations into the major controls on lava flow emplacement offer a theoretical framework for understanding channel network formation. Lava flow emplacement is controlled primarily by effusion rate (volumetric flux), underlying slope, and rheology (a function of melt composition, temperature, and crystallinity). Effusion rate governs how fast and far flows may advance. Average effusion rates correlate positively with final flow lengths [Walker, 1973]. In Hawai‘i, effusion rate is also a good predictor of flow advance rates, with higher effusion rates feeding faster flows [Kauahikaua et al., 2003]. However, the effusion rate at the vent does not necessarily record the lava supply to the flow front. When a flow branches or merges, the lava flux within those branches must be divided or summed. As both flow length and advance rate are proportional to lava flux (volumetric or mass flow rate), the upslope history of a flow branch influences both its length and advance rate [Macdonald, 1943; Guest et al., 1987; Pinkerton and Wilson, 1994]. From this perspective, understanding the development and evolution of a channel network is critical for anticipating flow behavior. Additional external and internal factors that control flow behavior and morphology include both slope and lava rheology. Flows are faster, narrower, and thinner 90 on steeper slopes [e.g., Lister, 1992; Kerr et al., 2006], and are slower, thicker, and wider for higher viscosity lavas [e.g., Lister, 1992; Lyman et al., 2005; Kerr et al., 2006]. This applies to comparisons between low-viscosity basaltic flows and high-viscosity rhyolitic flows [Walker, 1973], but also for changing rheology along a single flow [Lipman and Banks, 1987; Riker et al., 2009]. During emplacement, heat loss increases the lava viscosity because of both cooling and crystallization [e.g., Crisp et al., 1994]. Flows therefore are expected to slow, thicken, and widen as they travel away from the vent [e.g., Kerr et al., 2006; Kerr and Lyman, 2007]. These basic controls delineate expectations for channel network studies. The initial channel network is controlled by interactions between the flow front and the local topography. If a flow front is able to overtop an obstacle up to the height of the flow thickness [e.g., Favalli et al., 2005], then flow branching should increase as flow height decreases. Flows should thin on high slopes, where the flow accelerates [e.g., Lister, 1992]. When a flow splits, the lava flux supplying the flow front also splits, resulting in new flow branches that should be both slower and narrower [Kerr et al., 2006]. Confluences would have the opposite effects. These theoretical relationships between underlying topography, network geometry, and flow morphology and behavior can be tested with morphologic data and flow observations. 2.2. Morphologic Observations of Lava Flows The effects of flux, slope, and rheology on flow and channel morphology have been documented by observations of recent eruptions in both Hawai‘i and Italy. Measurements made along individual flows show that channel morphology, in particular, 91 responds to changes in flux, slope, and rheology [Lipman and Banks, 1987]. High- resolution lidar datasets collected both during and after flow emplacement provide new information on channel and flow dimensions [e.g., Mazzarini et al., 2005; Pyle and Elliot, 2006; Coltelli et al., 2007; Ventura and Vilardo, 2008; Favalli et al., 2010b; Deardorff and Cashman, 2012; Tarquini et al., 2012]. These data have been used to examine correlations between flow emplacement parameters and morphology, to relate measurements of flow volume and morphology to eruption observations [e.g., Pyle and Elliot, 2006; Coltelli et al., 2007] and to infer the emplacement conditions of unobserved flows [e.g., Deardorff and Cashman, 2012]. Prior work demonstrates that lidar data also offer new opportunities to investigate channel networks, as they permit quantification of both planform topology and 3D flow and channel dimensions [Dietterich et al., in press]. With a pre-eruptive digital elevation model (DEM), whole-flow thickness data and characterization of the pre-eruptive surface is possible. These data have been used to test simple controls on flow emplacement, for instance verifying that pre-eruptive slope is negatively correlated with flow thickness along channels of the ML84 flow [Dietterich et al., in press]. However, these datasets also record potential influences of the channel network on flow morphology that motivate further investigation. Thickness analysis of the ML84 flow has highlighted areas of local thickening at sites of flow bifurcations (Figure 2), showing a relationship between the network structure and flow morphology. Integrating these thickness data into lava volume, illustrates the influence of branching on lava volume distribution between each flow branch. Flow 1 was the longest-lived branch and it has the largest total volume. Interestingly, the data in Figure 2 show a long section in the middle of the flow (outlined 92 in Figure 2) that is quite thin, despite funneling the extensive lava volume down slope. These observations hint at the small and large-scale impacts of channel networks, and prompt us to address the role of channel networks in Hawaiian lava flows. Figure 2. Flow thickness map of the ML84 flow. Inset shows local flow thickening at bifurcations. Highlighted narrow section of Flow 1 demonstrates that regions with low thickness may still carry high lava volumes. Volumes for the four main branches within the lidar extent (Figure 1a) are also shown [Dietterich et al., in press]. 3. Methods We use several different datasets and techniques to quantify Hawaiian lava channel networks and their associated flow and channel morphologies. As described in 93 Dietterich et al. [in press], DEMs of the upper two-thirds of the ML84 flow and the full extent of the KDec74 flow were generated by airborne lidar at a resolution of 1 meter; these data were used to map flow features and to quantify flow surface characteristics. Associated maps of lidar intensity record the infrared “color” of the lidar signal, which is controlled by surface age and roughness [e.g., Mazzarini et al., 2007]. Intensity variations can be used to identify boundaries between flows of different ages, as well as channels within the flow. Mosaicked digital orthophoto quadrangles produced by the National Digital Orthoimagery Program complement the lidar intensity data with color imagery of the flows. To analyze planform and 3D channel morphologies, we first produce maps of flow outlines, levees, and channel centerlines. Flow outlines are assembled by tracing flow margins in the lidar intensity data; we use color orthophotos to help distinguish the flows. Slope and curvature maps of the lidar DEM are used to identify and map the steep margins of levees [e.g., Tarquini et al., 2012]. We then digitize channel centerlines between identified levees, assisted by the contrast in the lidar intensity between the channel and the bounding levee. Levees and associated channel centerlines are not continuous where channels are buried by later activity, but extend along most flow pathways from the vent toward the flow toe. Visible levees are absent, however, in the pāhoehoe-dominated proximal flows, and in distal regions where ‘a‘ā channels transition to dispersed flow morphologies [e.g., Lipman and Banks, 1987]. Lidar DEMs provide high-resolution elevation data of the present day flow surfaces. To determine the role of pre-eruptive topography on channel development, we use stereophotogrammetry to make a 10-meter resolution pre-eruptive DEM using 94 imagery from a USGS aerial photography survey from 1977 [root mean squared error (RMSE) of 2.4 meters from comparison to the lidar DEM off the flow]. Flow thicknesses and volume can then be calculated by subtracting the present-day elevations from the pre- eruptive elevations. We characterize the pre-eruptive slope by averaging over a length scale of 90 meters, the slope length scale that best correlates with flow thickness [Dietterich et al., in press]. 3.1. Quantifying Channel Networks To quantify channel network geometry, we draw on a range of methods from different sources. From geomorphology, we employ methods developed for quantifying the channels in braided rivers [Egozi and Ashmore, 2008], as well as the distributary patterns in deltas [Edmonds and Slingerland, 2007]. From neuroscience, we repurpose network topology tools for measuring brain connectivity [Rubinov and Sporns, 2010; Marra et al., 2013]. All metrics are based solely on geometry with no inherent assumptions about physics, allowing their direct application to lava flows. Each metric acts at a different scale, requires different inputs and reveals different features of the lava channel networks. 3.1.1. Braiding Index The first and simplest metric for quantifying lava channel networks is the braiding index, defined as the number of channels that intersect a given cross-section line [Egozi and Ashmore, 2008]. In distributary flows, like the ML84 flow, we are interested in local changes in channel network geometry. For this reason, we apply the braiding index 95 metric not on the entire flow, but instead within the main flow branches (Flows 1, 1A, 1B, 1 North, 1 East, 2, 3, and 4; Figure 1a). To do this, we select a path along each branch (e.g., a main channel centerline) and use the HEC-GeoRAS software (US Army Corps of Engineers, http://www.hec.usace.army.mil/software/hec-georas/) to create flow cross-sections spaced 100 meters apart along the flow length. We then link the map of channel centerlines to the flow cross-sections in ArcGIS and extract the number of channels intersected by each cross-section intersects (the braiding index; Figure 3). Figure 3. Map of the channel network of the ML84 flow within the lidar extent with flow section cross-section lines color-coded by braiding index. 3.1.2. Network Measures We can also treat the channel system as a true inter-connected, directional 96 network. This allows us to quantify the connectivity between channels, paths through the network, positions of individual channels and junctions within the network. As these measurements require a completely connected network of channels, missing channel segments must be reconnected by hand. The resulting channel map can be converted into a network dataset composed of channel segments, or edges, that split or merge at nodes. Each channel segment has nodes of origin and destination, as well as other properties such as length and slope. We then translate the network into a connectivity matrix, where each row and column is a node, and a channel segment that runs from a row node to a column node is represented by a value of one in the matrix at (row node, column node). Zeros indicate that the nodes do not connect in the direction of row node to column node. Values in the matrix can also represent weighted metrics, such as the length of the channel segment. These lengths can be summed down flow to calculate the distance along flow from the vent to a given node or channel segment. The lengths of each channel segment can also be used to calculate channel sinuosity, a common measure in fluvial geomorphology, by dividing the segment lengths by the straight-line distances between the start and end nodes for each segment. 3.1.2.1. Ordering The order of segments of the lava flow channel system can be used to quantify the overall distributary or tributary form. Tools for ordering have been developed for studying river systems (stream ordering). Used to classify parts of a river network and define the size of a river system, stream ordering assigns values to each river segment in the network, building along the flow as segments merge (Figure 4). Strahler ordering 97 assigns to all initial streams a value of one; this value increments where two streams of the same order merge (Figure 4a) [Strahler, 1957]. Alternatively, Shreve ordering gives initial streams values of one, but the values of merging streams sum along the flow (Figure 4b) [Shreve, 1966]. Lava channel networks can also be distributary (ML84); here the ordering works backward, such that each flow toe has a value of one, and the values increase toward the vent (Figures 4c and d). Lava flow channels also split and reconnect, which requires modification of the ordering algorithms to allow the order to both decrease and increase [Gleyzer et al., 2004]. After modification, Strahler order is interval in nature, while Shreve order utilizes ratios. Ordering of lava channels provides a way to classify the network complexity and to identify the major (backbone) channels within a flow. Distributary systems, such as river deltas, have also been described using the bifurcation order, which is an integer that represents the number of bifurcations upstream of a channel (Figures 4c and g) [Edmonds and Slingerland, 2007]. The order therefore increases along the distributary system. This methodology cannot account for merging channels, but documents flow splitting within a given channel segment. None of these metrics provide an ordering scheme that moves with the flow direction (like the bifurcation order), while also allowing confluences (like the stream ordering schemes). We therefore introduce a final ordering metric, which we call the estimated flux order (EFO; Figures 4 d and h). We calculate the EFO by first distributing a hypothetical unit flux of 1 from the vents to the flow toes. As channels branch along the flow, this number divides evenly among the channel segments. If channels merge, their estimated fluxes are summed. We convert the estimated flux to an ordering scheme as EFO = log2(estimated flux-1). For tributary flows, the estimated flux is the same as the 98 Figure 4. Schematic of the different channel ordering schemes used. Strahler and Shreve order are directional, requiring calculation from the flow toe upstream for distributary flows and from vents downstream for tributary channel networks. Bifurcation order and estimated flux order (EFO) are always tabulated downstream from the vent(s). Estimated flux values used to calculate EFO are shown in parentheses in gray in (d) and (h). Shreve order, producing a direct correlation with EFO (Figure 5). In distributary flows it mimics the bifurcation order, but includes confluences. This ordering metric is designed specifically to describe the distributary nature of most lava channel networks and to allow ratio, rather than solely integer, values. 99 Figure 5. Map of estimated flux order (EFO) and Shreve order for the KDec74 flow. For tributary systems like this, EFO = log2(Shreve order-1) and both show channel dominance based on how many tributaries feed each channel segment. The values correspond well to the channel widths visible in the inset Google Earth satellite image. 3.1.2.2. Brain Connectivity Toolbox Additional analysis of the channel network uses the Brain Connectivity Toolbox [Rubinov and Sporns, 2010], which was designed for neuroscience applications but has also been employed to investigate braided river networks [Marra et al., 2013]. This MATLAB toolbox analyzes the connectivity matrix of the network to calculate the betweenness centrality and node degree, from which we compute the normalized betweenness centrality (NBC) and classify bifurcations and confluences (Figure 6). The betweenness centrality is the proportion of shortest paths through the network that pass through a given channel segment or node [Rubinov and Sporns, 2010]. It is therefore a measure of channel, or node, importance within the network. The betweenness centrality 100 values in our channel network are strongly influenced by the boundaries of the network, because (unlike the brain) the network has start and end points. We therefore normalize betweenness centrality values by the product of the number of nodes along the flow (the potential number of shortest paths that could pass through the node in a directed network) to remove the boundary effects and generate the NBC metric [Marra et al., 2013]. We also measure the node degree quantity, which represents the number of channel segments that enter (in-degree) or leave (out-degree) a given node, and can be used to identify and classify all of the bifurcations (out-degree > 1) and confluences (in-degree > 1) within a network. Figure 6. Map of the normalized betweenness centrality (NBC) measure with classified nodes produced with the Brain Connectivity Toolbox [Rubinov and Sporns, 2010]. Values of NBC are very small because they are from a peripheral flow branch that makes up a small portion of the great number of channel segments within the entire ML84 flow. 3.2. Quantifying Flow and Channel Morphology We quantify our channel network at both the flow segment scale (braiding index) and the channel segment scale (network measures) and therefore measure flow and channel dimensions at the equivalent scales for comparison. We compute flow segment morphology along the cross-sections used for calculating the braiding index and channel segment morphology by averaging measurements made on channel cross-sections within each segment. 101 3.2.1. Morphology at Flow Cross-Sections The braiding index measures the number of parallel channels along a flow section. Possible morphology metrics that could relate to the braiding index include flow width, flow thickness, and underlying slope. To calculate each of these, we use a GIS- based method from Deardorff and Cashman [2012] that buffers each braiding index cross-section line by adding ten meters on either side to create 20-meter-wide measurement boxes oriented perpendicular to the flow. To measure flow width, the total flow area within each box is calculated as the intersection of the flow area polygon (outline) with each box. We divide this area by the width of the box (20 m) to compute an average flow width for each cross-section. The same cross-section boxes are used to calculate the average thickness of the flow and the underlying slope for each braiding index line. We use zonal statistics in ArcGIS to average all thickness or slope values within each box. This requires oversampling the 10-m resolution datasets to 1 m to accommodate the differing resolutions of raster and vector data; it also necessitates iterating the zonal statistics tool for overlapping boxes. This method allows us to calculate all metrics for each braiding index measurement without sensitivity to local excursions at cross-section lines. 3.2.2. Morphology of Channel Segments To characterize the morphology of individual channel segments, we average measurements from cross-section lines along a given channel segment to determine the width, depth, flow thickness at the channel center (channel height), and flow thickness at the channel levees (levee height) (Figure 7). We use the HEC-GeoRAS software to 102 construct cross-sections spaced ten meters apart along each channel centerline segment. We trim or extend lines to intersect corresponding levee traces, and select only lines that intersect levees on both sides of the channel; this excludes channel segments that do not have two levees (including near-vent overprinted regions and more distal segments that lack well-defined channel morphologies). The channel center points are the intersections of each cross-section line with the channel centerline; intersections with the levees yield left and right levee points that provide a measure of the channel width. Channel and levee heights are calculated by subtracting the pre-eruptive elevation from the lidar elevation. Channel depth is the difference between the average levee and channel heights. We average these values across all of the cross-section lines in a given channel segment to compare to the network data. Figure 7. Example of the channel segment morphology measurements from Flow 4 of the ML84 flow. Traced levees (red) and channel centerlines (black) are used to produce numerous flow cross-sections (yellow) and extract location and elevation data from the channel center (green) and levee points (blue). Morphology measures for each cross- section line are shown in the inset cartoon. 103 4. Results All measures of the channel network show that some parts of the flows are more complex than others (Figure 2), and show different patterns in their network geometry and morphology (Figures 4–7). These data allow us to compare network measurements for the contrasting distributary (ML84) and tributary (KDec74) channel networks. By examining each of the network metrics, it is possible to identify those that are most distinctive and therefore valuable for examining the relationships between the network and morphology metrics. We focus on the ML84 flow for morphometric analysis because of its well-developed channels, more complete observations, and complex history that produced both short- and long-lived branches. At a fundamental level, lava flow channel networks are either tributary or distributary based on whether they merge or branch (measured by comparing the total number of confluences versus bifurcations). The contrast between tributary and distributary networks is highlighted by trends in the network metrics with distance along both flows (Figure 8). In the distributary ML84 flow, the Strahler and Shreve orders increase toward the vent where the flow branches originate from a common point (Figures 8a and b, black). Bifurcation order and EFO show the opposite trend, emphasizing the tendency for the flow to split away from the vent (Figure 8c and d, black). In contrast, the KDec74 flow has a tributary form, because of topographic confinement by the Koa‘e fault scarps (Figure 5). This network geometry is characterized by increasing Strahler and Shreve order with distance as tributaries merge, much like a river network (Figures 8a and b, red). Because bifurcation order counts only points where the flow splits and ignores confluences, the bifurcation order also increases, albeit slowly, 104 down the KDec74 flow (Figure 8c, red). EFO, however, which mimics the bifurcation order but also includes confluences, decreases away from the vent (Figure 8d, red). Like the bifurcation order and EFO, the normalized betweenness centrality (NBC) is measured identically whether a network is tributary or distributary in form. It highlights the most important channel segments, but still shows a trend with distance (Figure 8e). In the ML84 flow, the highest NBC values occur near the vent where the dominant channel segments connect to the largest portion of the flow. The values then decrease toward the flow toe (black points). In contrast, the KDec74 flow has the highest NBC values far away from the vent, along the part of the channel network that converges along the fault scarps (red points). Comparing Figures 8b and e, it is apparent that both Shreve order and NBC provide a measure of channel dominance, so they are correlated strongly in the ML84 flow and somewhat less so in the KDec74 flow (Figures 8f and g). From the trends in the network metrics, it is possible to identify those that best differentiate channel network geometry. There are three metrics that show different trends for the different network structures: Shreve order, EFO, and NBC. Because Shreve order and NBC both measure channel importance, we use Shreve order and EFO as our primary network metrics. The ML84 flow, with its well defined ‘a‘ā channels and levees, allows us to relate flow and channel morphology. Average channel width gradually increases with distance from ~15 m near the vents, to up to 160 m distally (Figure 9a). The trend resets after the broad, distal channels of Flows 2–4 end, giving way to the long, narrow section of Flow 1 (outlined in Figure 2; 15–20 m channel width) that feeds wider channels with distance. This observation matches measurements made along long-lived (Flow 1) ML84 channels 105 Figure 8. Comparison of how the different network measures change with distance along flow for distributary (ML84, black, left y-axis) and tributary (KDec74, red, right y-axis) flows. (a) Strahler order, (b) Shreve order, (c) bifurcation order, (d) estimated flux order, (e) normalized betweenness centrality, and (f) and (g) which are comparisons between Shreve order and normalized betweenness centrality, which both measure channel importance, in ML84 and KDec74, respectively. immediately following the eruption (gray squares) [Lipman and Banks, 1987]. Channel widening is an expected consequence of increasing viscosity [e.g., Kerr et al., 2006], 106 which changes from ~102 Pa s at the vents to 107 Pa s at the flow toe [Moore, 1987]. Increasing viscosity should also increase the average levee height, although we see no obvious trend in our data (Figure 9b), where the scatter probably reflects the effect of secondary overflows on individual levee segments. Figure 9. Distance along flow shows a modest correlation with (a) channel width and little relationship with (b) average levee height. Post-emplacement field measurements by Lipman and Banks [1987] are included as the gray squares for validation. The gap in the data show the narrow section of Flow 1 that lies in between the ends of Flows 2–4 and before the breakouts of Flow 1A and 1B (Figure 1a). 5. Discussion The simplest hypothesis about the origin of channel networks is that they result from interactions between advancing lava flows and topographic obstacles. To test this hypothesis, we examine relationships between the network geometry, the local slope and flow morphology. We then view channel networks from the opposite perspective, and evaluate the effects of the channel network on flow behavior. To do this, we compare network metrics to the suite of morphometric measurements, as well as observations of overall flow lengths and advance rates. 107 5.1. The Formation of Channel Networks We can compare the primary channel network (as recorded in the short-lived ML84 Flows 2–4) with underlying topography to determine whether our results support the hypothesis that obstacle interaction is the process driving network formation. Examining individual obstacles is difficult with the coarse (10-meter) DEM resolution available for the pre-eruptive surface, but given an arbitrary distribution of obstacle sizes, we can use underlying slope, flow thickness, and branching data to test our hypothesis. At the flow segment scale, braiding index (channel branching) is positively correlated with slope (Figures 10a and b). Flows are also thinner on higher slopes (Figure 10c), thus our results are consistent with a model of flow branching caused by flows thinning over steep slopes and being forced to split around more obstacles that it cannot overtop. Conversely, thicker flows on shallower slopes may cover the topography more completely. This simple model would suggest that flows should be distributary, as is apparent even in the KDec74 flow, which splits into a number of branches when it escapes from the confining effects of an individual fault scarp segment (Figure 1b). While slope may control interactions with obstacles on the scale of the flow front, large topographic barriers also shape lava flow channel networks. A clear example is the ten-meter high fault scarps that funnel disparate channels together to produce the tributary form of the KDec74 flow (Figures 1b and 5). Here the Koa‘e fault scarps cut across the direct downhill path of the flows emanating from a chain of fissure vents, forcing flows to merge and be redirected to the southwest. Exceptions include flows from the westernmost vents, which do not merge with the main channel, in part because the first major fault scarp does not extend sufficiently far to the west. The relative timing of 108 Figure 10. (a) Braiding index and channel network of Flow 4 of the ML84 flow overlain on underlying slope calculated from the photogrammetric DEM of pre-eruptive topography at a length scale of 90 meters. The map reveals that sections with high braiding index correspond to sections with high slope, suggesting a correlation. (b) Plot showing the positive correlation between slope and braiding index in (a). (c) Plot showing the negative correlation between slope and flow thickness. (d) Flow width positively correlates with braiding index. Data from the whole flow (empty circles) is largely chaotic due to the effects of breakouts and spillovers from the long-lasting Flow 1. Focusing on only Flows 2–4 (gray) and particularly without the most proximal or distal sections where channels are less distinct (using only areas with identifiable levees, black), the correlation between braiding index and flow width becomes quite strong (r2 = 0.77). fissure vent activity was also important: the westernmost vents were the last to be activated and did not erupt sufficient lava volume to send flows to the next fault scarp, where they would have been forced to merge with the main channel. A similar situation exists along much of the 50 km length of the ML1859 flow. This large flow began as a distributary flow that traveled down the northwest slope of Mauna Loa (Figure 1 inset). The flow was confined, however, when it reached the saddle between Hualālai and Mauna Loa volcanoes; the result was a flow that focused into a single channel for the rest of its length [Riker et al., 2009]. 109 Evidence of local topographic confinement can also been seen in flows where networks are primarily distributary. For example, the central, narrow section of ML84 Flow 1 comprises a single channel for about six kilometers (outlined in Figure 2). This segment follows, on the down slope side, a ~3 m high levee that marks the margin of a lava flow emplaced in 1942 (Figure 1a). A 1942 channel levee cuts at a slight angle across the slope, and thus shifts the course of the ML84 flow. It seems likely that topographic confinement also prevented the ML84 flow from splitting extensively in this section. As a consequence, this was one of the longest-lived and highest-flux channel segments [Lipman and Banks, 1987]. Confining features may also act as dividing features, depending on their location and orientation; the vent of the same ML1942 flow functions in part to separate Flow 1 from Flows 2–4 in the ML1984 flow (Figure 1a). Morphologic features suggest other factors that may contribute to primary channel network development. Individual branches of the ML84 flow, including Flows 1 North, 1 East, 1B, and 2, have braiding indices that decrease with distance along the flow (Figure 3). This suggests that the increase in viscosity with flow distance may reduce flow branching. One explanation for this behavior could be that higher viscosities cause the flow to thicken, which might reduce topographic interactions. Another common feature is local flow thickening at bifurcations (Figure 2) [Dietterich et al., in press]. Some thickening may be due to lava ponding within channels prior to breakouts [Lipman and Banks, 1987]. In many cases, however, topographic obstacles are preserved at the point of channel bifurcation and are often plastered on the upslope side with lava that has the morphology of a stationary wave, like lava run-up [e.g., Soule et al., 2004]. The physics of stationary waves in viscous flows is not fully understood [e.g., Baxter et al., 2009], but 110 their presence suggests that flow dynamics, especially velocity, may also influence topographic interactions. The formation of a stationary wave would make it possible for a flow to overtop an obstacle that is higher than the thickness of the flow without the wave. 5.2. The Influence of the Channel Network on Flow Morphology and Evolution Channel networks govern the distribution of lava supply and provide lava pathways that can evolve to create a range of flow and channel segment morphologies. On the scale of flow segments, higher braiding indices correspond to thinner, wider flows (Figure 10). Both thinning and branching are responses of lava to an increase in slope; dividing lava flux among multiple flow branches should further reduce the flow thickness [Lister, 1992]. The overall flow width increases by branching, as shown by the correlation between flow width and braiding index (Figure 10d). This supports the observation that complex flows have wider flow fields [e.g., Guest et al., 1987] and demonstrates that the physics of simple flows is not fully valid for describing the behavior of branched flows. Figure 10d also illustrates the scatter in the flow segment data. There is no clear trend in the braiding index and width data from the entire flow (open circles). However, limiting data to channelized regions of Flows 2–4 (black circles) highlights the positive correlation, because excluding the longer-lived Flow 1, with its many secondary overflows and breakout events, and non-channelized regions (near vent and distal) largely removes the effects of viscosity and channel evolution from our analysis. Within channel segments throughout the ML84 flow field there is also evidence for a decrease in channel width with increased branching. Theoretical analysis shows that 111 channel width should narrow with decreased flux [Kerr et al., 2006], and therefore increased EFO, but capturing the branching effect alone is challenging. Although channel width generally increases with EFO (Figure 11a), both are strongly dependent on the distance from the vent because of increasing flow viscosity with distance (Figures 8d and 9a). We can remove the effect of viscosity by plotting channel width versus EFO while controlling for distance. Within a narrow range of distances (e.g., 11–14 km) we see the expected decrease in channel width with increasing EFO (Figure 11b). Figure 11. Estimated flux order (a) and Shreve order (c) show relationships with channel width in the ML84 flow that mostly reflect the effect of distance (Figures 8 and 9). Using only data from an approximately constant distance along flow, increasing estimated flux order (b) and Shreve order (d) correspond to narrower channels, reflecting the narrowing of channels with increasing flow branching and channel importance. 112 Additionally, there is a connection between measures of channel importance and channel morphology. Channels with high values of both Shreve order and NBC are connected to large portions of the flow. These metrics provide a measure of channel importance, in that a blockage in a high value channel would affect a relatively high percentage of the flow. Channel width decreases with increasing Shreve order (Figure 11c), although Shreve order also varies somewhat with distance along flow (Figure 8b). Subsampling these data to eliminate the effect of distance shows the same trend (Figure 11d). Channel dominance also influences other aspects of channel morphology, with higher values of Shreve order weakly corresponding to increased levee heights and decreasing channel segment sinuosity (Figure 12). Both show a wide range of values at low Shreve order but narrow to a more limited range at higher order. These results send mixed messages because the expectation is that major channels would transport higher fluxes, and would therefore have higher levees, but also be wider [Kerr et al., 2006]. This leads to another important factor controlling channel width, which is the duration of activity with the channel. Figure 12. (a) Average levee height and (b) sinuosity compared to Shreve order. 113 We can expand our analysis to include temporal changes in channel morphology by incorporating the observed longevity of channel segments in the ML84 flow. We apply the length of time major flow branches were active to the channels within them to estimate the duration of activity of each channel segment [as reported in Lipman and Banks, 1987]. Shreve order correlates with channel longevity (Figure 13), which shows that dominant channels were active for longer. Overall, it also appears that longer-lived channels are narrower, and have higher levees, than short-lived channels (Figure 14). Figure 13. Box and whisker plot of Shreve order correlating with channel longevity. The boxes and whiskers denote quartiles, with the median marked in the center and plus symbols indicating outliers beyond ±2.7σ (assuming a normal distribution). The dashed arrow highlights the overall increase in duration of channel activity with Shreve order. The processes behind these relationships can be assessed from observations made during the eruption. Lipman and Banks [1987] described channel evolution as a transition from wide dispersed zones behind the flow toe to narrow levee-bounded stable channel zones driven by progressive solidification inward from the channel margins. Our 114 Figure 14. Longevity also shows a relationship with flow morphology, with (a) channel narrowing and (b) average levee height increasing with channel longevity. The dashed arrows show the interpreted trends. The channel segments active the longest, corresponding to the central, narrow section of Flow 1, have especially narrow channel widths and also smaller levees. observations that longer-lived channels are narrower with higher levees support the interpretation of increased channel stability with time. Similarly channels with higher Shreve order are generally straighter (Figure 12b). This relationship can be explained by observations that compression and drag by lava flowing in the channel triggered levee collapse on the inside banks of meanders, causing individual channel segments to straighten over time [Lipman and Banks, 1987]. If higher Shreve order implies greater channel longevity, this process can explain the weak negative correlation between Shreve order and sinuosity. A variant on these trends is provided by the longest-lived section of the ML84 flow, represented by the central, confined section of Flow 1 (outlined in Figure 2). This channel segment is very narrow (Figure 14a) but also has low average levee heights (Figure 14b). The latter observation suggests that confinement of this flow section may have limited channel overflows, and therefore levee heights. In fact, observations at the time showed extremely efficient transport of lava through this section [Lipman and 115 Banks, 1987]. The narrowness and low channel profile of this channel segment also emphasizes the lack of correlation between channel morphology and total volume of lava transported (Figure 2). In this example, the entire lava volume of Flows 1 and 1A were fed through this narrow channel segment. 5.3. The Influence of the Channel Network on Flow Emplacement Behavior We can use records from these and other Hawaiian channel-fed flows to investigate the impacts of channel networks on flow length and advance rate. Given that bifurcations and confluences serve to divide or merge volumetric flux, we would expect that flows should advance more slowly, and be shorter, if they branch, and the opposite if they merge. We can test this from records of flow emplacement, where we compare (1) the lengths of flows that differ only in their degree of branching or confinement and (2) the advance rates of flows before and after splitting or merging. 5.3.1. Influence on Lava Flow Length Lava flow length is a function of many parameters; to isolate the effect of network topology, we compare flows with similar effusion rates and durations. Past studies have shown that flow length and effusion rate tend to be positively correlated [Walker, 1973; Malin, 1980; Pinkerton and Wilson, 1994], although this relationship is not very robust in Hawai‘i (Figure 15). In fact, over a wide range of effusion rates, the flows appear limited to about 25 km in length [Riker et al., 2009]. An exception is the Mauna Loa 1859 lava flow, which traveled 50 km to the coast. Riker et al. [2009] showed that this flow was erupted at an unusually high temperature (~1200˚C) and they therefore inferred that the 116 great distance reflected more protracted cooling and crystallization. Here we investigate the extent to which branching and confinement may help explain these data. Figure 15. Plot of effusion rate versus flow length for channelized lava flows erupted at Mauna Loa (squares, dark gray region) and Kīlauea (diamonds, light gray region). Regions summarize points from each volcano, while lava flows discussed in the text are shown as labeled points, including two points for the ML84 flow representing Flows 1 and 1A and Flows 2–4 separately (ML84 F1/1A and F2–4). Flows are marked (C) if confined and all flows have durations in hours (h) or days (d) included. The lines of constant duration are from an equivalent plot by Pinkerton and Wilson [1994]. Regions and points produced from data in Wolfe [1988], Heliker et al. [2001], Heliker and Mattox [2003], and Riker et al. [2009]. We start by comparing the ML84 and ML1859 flows, which have similar effusion rates and durations, but different channel topologies (Figure 1). Both flows were erupted at rates ~102 m3/s (~275 m3/s for initial advance of ML84 and ≤391 m3/s for ML1859, which had an average effusion rate of 116–235 m3/s during the ‘a‘ā phase) [Lipman and Banks, 1987; Riker et al., 2009]. Despite similar durations, the ML 1859 flow traveled 117 almost twice as far as ML84 before reaching the ocean (Figure 15). Importantly, the ML1859 flow was also confined, by neighboring Hualālai volcano, to a single channel for much of its length. More generally, we find that confined flows in Hawai‘i are often unexpectedly long. The length of the KDec74 flow, for instance, is surprising, given that its 6-hour duration would be expected to produce a much shorter flow (Figure 15 duration lines) [Pinkerton and Wilson, 1994]. However, as illustrated above, the confinement of the flow by fault scarps focused the KDec74 flow into a single channel, allowed the flow to travel farther down slope [Soule et al., 2004]. Another pair of flows with similar effusion rates offers a final comparison between confined and unconfined flows. Episodes 40 and 43 of the 1983–1986 fountaining episodes at Pu‘u ‘Ō‘ō cone (Kīlauea East Rift Zone) were short-duration eruptions with very different channel network geometries. Both were distributary, but Episode 40 was confined within a topographic low, with a maximum Strahler order of 2, maximum Shreve order of 4, and average braiding index of 1.2. Episode 43 has a maximum Strahler order of 3, maximum Shreve order of 7, and average braiding index of 2.3. With nearly identical average effusion rates of ~230 m3/s, Episode 40 traveled 8.4 km in 14 hours, while the wider, more heavily branched Episode 43 traveled only 5.3 km in 12 hours (Figure 15) [Heliker et al., 2001]. 5.3.2. Influence on Advance Rate Historical accounts of flow advance rates provide information on the effect of branching on flow advance rates. The most complete records can be found for the early fountaining episodes at Pu‘u ‘Ō‘ō cone, when flow front locations were meticulously 118 mapped with time [Wolfe, 1988; Heliker et al., 2001], allowing the velocities before or after flow bifurcations to be measured. Most flow bifurcations are accompanied by a decrease in flow advance rate. For instance, in Episode 11, the NE flow was originally advancing at a rate of 330 m/hr, but then split to create two branches that advanced at rates of 142 m/hr and 123 m/hr (Figure 16a). A compilation of all available pre- and post-bifurcation advance rate pairs for Episodes 1–48 shows that, in most instances, the advance rate drops approximately in half after a flow split (plotting along the 1:2 line in Figure 16b). All except one of the points that fall above the 1:1 line have significant slope increases after the bifurcation, which would have caused individual flow branches to accelerate. In the case of Episode 21, another branch in the flow ceased activity, which would have redirected the lava supply to the new branches. There are few confluences recorded in this dataset that can be used to test whether the flow accelerates when branches merge. In the only well constrained example (Episode 29), the flow advance rate doubles when the two branches merge. This analysis supports the notion that with flow branching comes decreased effective flux and therefore slower flow advance, and visa versa for flow merging. 6. Implications Our analysis of lava flow channel networks and their influence has numerous implications for lava flow hazard assessment and mitigation, and for modeling lava flow emplacement. Understanding the conditions that control channel network formation, and the acute effects of channel networks on flow behavior, offers new insight to improve predictions of lava flow behavior. Analysis of the impacts of splitting flows and their 119 Figure 16. (a) Map and advance data from Episode 11 of the Pu‘u ‘Ō‘ō-Kupaianaha eruption of Kīlauea (November 5–7, 1983) show the slowing of the flow after the bifurcation (redrafted from Wolfe [1988]). (b) Plot of advance rates before and after bifurcations from the maps of Episodes 1–48 (1983–1986) of the Pu‘u ‘Ō‘ō eruption. Symbols represent other factors that may influence the flow advance effect, influences slope increase or simultaneous branch shutoff. One confluence (square) is plotted in reverse (the pre-confluence advance rate as the y-value and the post-confluence advance rate as the x-value). Data come from Episodes 3, 5, 6, 7, 11, 16, 21, 29, 30, and 31 as documented in Wolfe [1988] and Heliker et al. [2001]. interactions with topographic barriers also helps inform lava flow diversion efforts. Finally, networks are present in many different styles of lava flow emplacement, ranging from braided lava tube networks to the complex geometries of flows on other terrestrial planets [e.g., Kauahikaua et al., 2003; Byrne et al., 2013]. Our work therefore can be applied to, and further draw from, this vast array of lava transport networks to expand our understanding of lava flow behavior. 6.1. Flow Prediction We have shown that lava channel networks exert significant control on where, how far, and how fast flows will travel. Distributary networks (ML84 flow) produce flows that are limited in length (typically to ≤25 km on Mauna Loa). However, the 120 opposite is true for tributary flows, which may accelerate at confluences and travel farther than expected (Figure 15). These observations suggest caution when using existing predictive models of channelized flow [e.g., Harris and Rowland, 2001], which do not consider the formation or evolution of complex channel networks. Models of flow through a single channel will yield maximum (fully confined) estimates of the advance rate and flow length. Incorporating lava channel networks into models of channelized flows requires that both the geometry of the channel network (distributary or tributary) and the number and location of branches be used to modify the lava flux along the flow. These modified flux values must be used to predict variations in both flow advance rate and flow length. To capture primary network formation, it is crucial to consider the underlying topography, including the pre-eruptive slope and topographic obstacles that may cause the flow to split or focus. Watershed and steepest descent path analysis can highlight major features that could confine flows (e.g., fault scarps). For instance, a lava shed analysis of the northeast rift zone of Mauna Loa using a 30-meter resolution DEM successfully divides the four main branches of the ML84 flow [Kauahikaua et al., 2003]. However, a higher resolution DEM would be required to capture features such as the confining ML1942 flow margin, which is only about 3 m high. Additional work is required to assess the locations and origins of the secondary overflows and breakouts that also influence the overall network geometry. During an eruption, real-time observations of channel network evolution, such as major breakouts, the shut-off of the vents or branches, and oscillations in flux [e.g., Favalli et al., 2010a], are needed to update any models of advance. 121 Predictions of flow length and advance rate are currently based on effusion rates, slope, and viscosity (including cooling), but do not include the effects of channel networks. As flows break up or become confined, their thermal efficiency will change with the flow dimensions (i.e., channel width) [Harris and Rowland, 2001], so that flow cooling will also be a function of the channel network. An initial step would be to include flux distribution and thermal efficiency from the overall tributary or distributary geometry into an estimate of the prediction uncertainty. Advances in cellular automata or other 3D lava flow models, such as those used at Mt. Etna, Italy [Hidaka et al., 2005; Vicari et al., 2007; Crisci et al., 2010], may make it possible to incorporate multi-channel flow behavior into predictive flow models. 6.2. Flow Mitigation When infrastructure and property have been threatened by lava flows, mitigation has been attempted using forced bifurcations, damming or diversion barriers [e.g., Barberi and Carapezza, 2004]. Most of these intervention strategies are intended to control lava flow branching and interaction with topographic obstacles. Only rarely, however, has the likely response of flows to these endeavors been assessed [e.g., Moore, 1982]. Flow diversion by splitting has been employed at Mauna Loa, Hawai‘i, and Mt. Etna, Italy, with the goal of redirecting lava flux from the flow lobe of concern into a new branch. In Hawai‘i, this has been done via aerial bombing campaigns during the 1935 and 1942 eruptions of Mauna Loa (Figure 1), although neither was particularly successful. In 1935, the eruption ceased (coincidently) shortly after bombing, while during the 1942 122 eruption, a new branch formed but quickly rejoined the original channel down slope [Lockwood and Torgerson, 1980]. At Mt. Etna, ground-based explosives have been used to breach levees and divert lava into artificial channels. This was first done in 1983 with small and temporary effects, but was more successful in 1992, when about two-thirds of the lava in the original channel was diverted [Barberi and Carapezza, 2004]. Diversion barriers have had an imperfect track record in Hawai‘i and Iceland but significant success in Italy. During the 1955 and 1960 eruptions of Kīlauea (Figure 1), barriers built to protect homes were almost all overtopped [Richter et al., 1970]. In 1973, barriers were combined with water-cooling to mitigate the effects of the Heimaey eruption, Iceland [Williams and Moore, 1983]. These mitigation strategies had the intended effect of saving Heimaey’s harbor, but caused diversion of one flow branch into the town. At Mt. Etna, mitigation strategies have included attempts to slow the advance of flows by using obstacles placed orthogonal to the flow front (e.g., 1983), as well as interventions to nudge flows away from infrastructure with oblique obstacles [Barberi and Carapezza, 2004]. These strategies have been generally successful, although earthen barriers have eventually been either overtopped or undermined. Our analysis of the impact of flow bifurcations validates flow diversion by splitting as a form of flow mitigation and offers data to describe the potential effects. First, if the flux is halved by a bifurcation, there can be a similarly halving of the rate of flow advance (Figure 16). Flow branching also limits the overall flow length (Figure 15). While bombing flow levees is challenging, our results suggest that artificial bifurcations could be used to slow and shorten flows, if there is sufficient up-flow area to permit resulting flow widening. From this perspective, upslope-pointing V-shaped diversion 123 barriers designed to both split and divert flows, like the existing barrier at the Mauna Loa Observatory [Moore, 1982], should prove effective. Our observations of flow confinement can also provide insight into the scale and effect of lava diversion barriers. Importantly, even small-scale features can impact flows (~3 m in the case of the ML1942 flow margin). This scale is determined, at least in part, by the slope of the underlying surface through its control on flow velocity and thickness. Diversion by barriers, however, confines flows, which inhibits flow branching and therefore maximizes the efficiency of the flow. For this reason, barriers that capture multiple flow pathways could produce major confluences, with consequent flow acceleration and lengthening (e.g., KDec74). The combination of flow slowing and shortening by bifurcation and flow rerouting by barrier diversion presents a significant toolset for flow mitigation. To take advantage of both effects, it may possible to merge the two into “leaky” barriers that split the flow while diverting it away from populated areas. Critical for such a scheme would be the analysis of limiting conditions to prevent the flow branches from re-joining below the diversion structure. Beyond offering information on flow mitigation strategies, channel network analysis can also help choose the best locations for flow intervention structures. The most important channel segments last longest, and transport the greatest lava volume. Early identification of these segments using network measures such as the Shreve order and betweenness centrality could be used to choose targets for intervention. Disruption of backbone channel segments would divide or divert large portions of the total flux and affect a large percentage of the flow field. For instance, a forced bifurcation of the main channel of the KDec74 down slope of the initial confluences from the eastern fissure 124 vents would have a much more significant impact than a barrier built below the westernmost fissure vent. Improving the success of flow mitigation strategies requires consideration of channel segment importance, as well as other location-dependent parameters, such as the underlying slope or the presence of topographic obstacles that could immediately reverse an intervention (e.g., a fault scarp below a forced bifurcation). 6.3. Wider Implications of Channel Networks for Lava Flow Emplacement Network analysis can inform comparisons of whole-flow behavior beyond the differing evolution of individual flow branches or channels. For example, the braiding index can be measured along an entire flow and then averaged to describe the mean number of parallel channels in the flow field. The maximum Strahler or Shreve values in a network can be used to describe the overall scale and complexity of a flow. The whole- flow network properties can be quantified by examining trends in the slope or directionality of the estimated flux order (EFO), Shreve order, or normalized betweenness centrality (NBC) along the flow. In this way, whole-flow metrics offer another tool to quantify a channel network and thus be able to compare networks to generalized flow parameters or behavior, such as effusion rate or final flow length. Lava transport networks are not unique to high effusion rate ‘a‘ā flows in Hawai‘i. Lava tubes, the main lava transport system within pāhoehoe flows, also have complex networks. Branching and merging within lava tubes has been documented with thermal imagery within flow fields in Hawai‘i [Kauahikaua et al., 2003], and evidenced by multiple simultaneous ocean entries or shatter rings and breakouts over multiple paths through the flow field during the ongoing Pu‘u ‘Ō‘ō-Kupaianaha eruption [Orr, 2010]. 125 The tubes form as pathways through the molten interiors of pāhoehoe lobes that have solidified on the outside but now carry lava flux from the vent to flow toe. The advance of pāhoehoe lobes is erratic and crust-dominated [e.g., Hon et al., 1994], but their routes are also governed by underlying topography [Crown and Baloga, 1999; Hamilton et al., 2013]. Patterns in lava tube networks could therefore record conditions of flow field evolution in the same way as ‘a‘ā channel networks. Finally, channel network analysis can be used to study emplacement of unobserved terrestrial and extraterrestrial lava flows. Morphology and network analysis of prehistoric or planetary flows could be used to infer parameters such as underlying slope, flow rheology, effusion rates, and flow field evolution. Of particular interest are both the form of the network and the location of dominant channels. Tributary and distributary channel networks have also been recognized on many other terrestrial planets and moons, but are often interpreted as fluvial because of their network structure [e.g., Byrne et al., 2013]. This raises the important question of whether there are network structures that are unique to lava flows, or at least distinguishable from river flows. For example, our work in Hawai‘i suggests that distributary networks are common features of lava flows, and thus contrast with the predominant tributary structure of river channel networks. As a result, lava flows decrease in flux along the network, in contrast to river networks, which gain in volume downstream. In summary, lava flow channel network analysis contributes to our ability to predict, influence, and understand lava flow emplacement. We have documented topographic controls on networks and their influence, in turn, on flow behavior. We have also discussed the implications of network analysis for both predictive models and 126 effective lava diversions. By combining analysis of 2D network geometry and 3D flow morphology, we can identify both primary controls on network development, and the influence of the network on flow evolution. Our results demonstrate that the flow morphology preserves signals of both intrinsic (viscosity change) and extrinsic (eruption duration) parameters. From this perspective, similar techniques could be applied to forensic study of older flows. Finally, our ideas about network formation could be tested and improved using 4D data sets that are now being obtained by repeat lidar surveys during effusive eruptions [e.g. Favalli et al., 2010a]. Further work to investigate the physics of interactions with topographic obstacles, the effects of flow composition and rheology, and the network geometry and morphology of lava tube networks, will better document the origins of channel networks and their effects for the broadest spectrum of lava flows. 7. Bridge In Chapter IV, I used pre-eruptive and post-eruptive morphology data, along with channel network maps, to investigate the creation, development, and impacts of lava flow channel networks. This analysis suggested that flow branching is the product of interaction with obstacles. I found that the flow was unable to overtop obstacles at high pre-eruptive slopes when the flow is thinnest. I also observed that flow merging and confinement is possible along extended obstacles (i.e., fault scarps) that cut across multiple flow lobes. From historical records of flow emplacement, the degree of confinement relates to how far or fast a flow may go, with branched flows advancing slower than confined ones and confined flows achieving long lengths. The network 127 geometry itself also controls the flow morphology. The parts of the flow that are more highly connected or spawn more branches than others have channel segments that are the longest-lived and narrowest, with high levees. This chapter therefore describes a few hypotheses related to lava flow interaction with topography that I test experimentally in Chapter V. My analysis in Chapter IV determined that flow thickness controlled whether a flow could overtop an obstacle and therefore not be forced to branch around it. I also observed thickening at flow branches in Chapters III and IV that suggest some local influence of flow run-up onto an obstacle, and therefore velocity influence in the obstacle interaction. I use analogue experiments with flows intersecting obstacles to test how these two parameters, as well as obstacle characteristics such as shape and size, dictate how a flow reacts to an obstacle collision. Watching the behavior of these experimental flows after splitting also offers an experimental investigation to the flow advance rate effects of splitting that I observed in historical records of flow emplacement in Chapter IV. 128 CHAPTER V EXPERIMENTAL MODELING OF LAVA FLOW INTERACTIONS WITH TOPOGRAPHIC OBSTACLES This chapter is in preparation for the Journal of Geophysical Research: Solid Earth. As lead author, I wrote the manuscript and performed all of the experiments, analysis, and interpretation for this chapter. My co-authors, Kathy Cashman (University of Oregon) and Alison Rust (University of Bristol), have helped me with the interpretation of my results, as well as editorial assistance. Alison Rust also aided in setting up the experiments at the University of Bristol Geological Fluid Dynamics Laboratory. 1. Introduction The parameters that control lava flow interactions with underlying topography are not well understood but play a significant role in governing lava flow emplacement, lava flow diversion, and efforts to calculate emplacement conditions. When lava flows encounter topographic obstacles, they may be diverted around them, run up onto them, or overtop them, with different implications for flow behavior. Flow diversion or overtopping controls the flow path, which can produce features like flow branching [e.g., Dietterich and Cashman, in press] or be used to artificially redirect lava flows with manmade diversion barriers [e.g., Richter et al., 1970; Williams and Moore, 1983; Barberi and Carapezza, 2004]. When flows run up onto, but do not overtop obstacles, the 129 run-up height has been used to estimate flow velocities during emplacement based on kinetic to potential energy conversion [Guest et al., 1995; Kauahikaua et al., 2002; Soule et al., 2004]. However, these diversion attempts and velocity calculations are largely theoretical because no study has investigated how lava flows interact with obstacles, or what factors might control the scale or form of topography that a flow may surmount during an interaction. We present hypotheses derived from field evidence for factors that influence obstacle interaction, and then test these with a suite of analogue experiments. Topographic obstacles, including pre-existing cones, flows, and fault scarps, control the topology of lava channels in Hawai‘i by creating interactions that split or confine the flow [Dietterich and Cashman, in press]. These have considerable consequences for flow lengths and advance rates, with confined flows traveling farther and faster than branched ones. Attempts to relate pre-eruptive topography to lava flow channel networks reveal that greater flow branching occurs over higher pre-eruptive slopes, which are associated with flow thinning, suggesting that flow thickness determines whether a flow must split around or will overtop an obstacle [Dietterich et al., in press; Dietterich and Cashman, in press]. The flow thickness parameter has been previously used to evaluate the influence of topography on lava flow paths, as well as the scale of lava flow diversion barriers. Simple steepest descent paths have been employed to calculate lava sheds and potential flow routes in Hawai‘i, but these ignore the potential of flows to surmount topography [e.g., Harris and Rowland, 2001; Kauahikaua et al., 2003]. To rectify this, Favalli et al. [2005] have used a stochastic approach to perturb the underlying surface by a characteristic flow thickness over repeated simulations, therefore allowing a simple 130 steepest descent line to pass through a barrier at the scale of the height of the flow. These simulations do not have a strong physical basis, but make the assumption that flow is able to overtop an obstacle of a size up to the flow thickness. A similar approach has been used to design the scale of diversion barriers required to redirect a flow. When installing the large V-shaped barrier that currently protects the NOAA Mauna Loa Observatory on the north slope of Mauna Loa, the thicknesses of typical Mauna Loa flows were tabulated from historical data and modeled, to support the choice of barrier dimensions [Moore, 1982]. This flow thickness hypothesis is therefore widely used, but has never been fully assessed. There is extensive field evidence that flow velocity may also influence the ability of flows to surmount obstacles. It is a well-observed phenomenon that when a flow encounters an obstacle, such as a rock in a river, the flow can run up the upslope side of the obstacle. This feature can be preserved in coatings of mud in debris flows or as lava plastered onto the front of an obstacle in lava flows [e.g., Pierson, 1985; Guest et al., 1995]. Run-up height, the increase in flow height upslope of the obstacle recorded by the flow veneer, is typically described as the height acquired when the flow converts all of its kinetic energy to potential energy against the obstacle [e.g., Guest et al., 1995]. The run- up height (Δh) is therefore a function of flow velocity (u), as well as gravity (g), This equation is frequently used to estimate flow velocities of unobserved flows [Guest et al., 1995; Kauahikaua et al., 2002; Soule et al., 2004], but has been found to be € Δh = u 2 2g , (1) 131 somewhat unreliable, yielding slower velocities than observed in debris flows [e.g., Pierson et al., 1985]. However, the run-up wave observations suggest that flow velocity could affect the thickness of the flow behind an obstacle, and therefore the likelihood of obstacle overtopping. Theoretical work for industrial purposes on the disturbance of thin film flows by obstacle encounters supports these field-derived parameters. In models of free surface thin film flow, the flow thickens upstream of an obstacle followed by a local dip in the flow surface immediately below the obstacle [Hayes et al., 2000; Baxter et al., 2009]. The shape as it splits around an obstacle can be related to the slope, the obstacle geometry and the Bond number, which is a function of density (ρ), gravity (g), flow thickness (H), and surface tension (σ) [Baxter et al., 2009]. At higher slopes, which create faster flows, and higher Bo, which increases with flow thickness, the wave upslope of the obstacle grows in height. When the flow intersects a cylinder of larger radius, the height of the wave also increases. This modeling is of the disturbance of a uniform, unbounded flow surface without a flow front or width, and includes significant surface tension. It cannot therefore be applied directly to the behavior of lava flows, but it does offer support to our flow height and velocity hypotheses. Studies of glacier response to underlying topography and obstacles in Antarctica reveal similar flow thickening behind the obstacle [Corti et al., 2003; 2008]. They also € Bo = ρgH 2 σ , (2) 132 show slowing of the flow as it approaches the obstacle, and uplift within the flow. This behavior may be analogous to that of lava flows, because both are viscous and have significant thickness, although glaciers flow much more slowly. We can therefore identify factors that may control the behavior of a lava flow intersecting an obstacle. Flow thickness and velocity may play a significant role in obstacle interactions. Obstacle height is an implied parameter, as well, to determine whether a flow of a given thickness will overtop the obstacle or not. It is also likely that the size and shape of the obstacle could have an effect, as observed in thin film flow modeling. There is field evidence for this, as well, in lava flow diversions at Mt. Etna, Italy. During eruptions at Etna in 1992 and 2001, barriers built orthogonally to the flow direction slowed the flow but were eventually overtopped [Barberi and Carapezza, 2004], while barriers built at an angle were able to nudge the flow laterally without easily being overrun [Barberi et al., 2003]. We investigate the roles of each of these parameters with analogue experiments. 2. Lava Flow Rheology and Theory Lava flow behavior and rheology have been studied extensively and allow us to identify the necessary physics to include in analogue modeling. Basaltic lava flows have initial viscosities of ~100 Pa s [e.g., Shaw, 1969], but values may reach as high as 107 Pa s with cooling away from the vent [Moore, 1987]. With flow thicknesses (H) on the order of meters [e.g., Wolfe et al., 1988], velocities (u) up to ~10 m s-1 [e.g., Kauahikaua et al., 2003], and densities (ρ) around 1500 kg m-3 [e.g., Lipman and Banks, 1987], basalt flows therefore have a Reynolds number, defined as 133 of order 10-3–102, indicating laminar flow. Lava flow rheology changes with cooling and crust formation, producing flows with an internal yield strength or crustal yield strength [e.g., Lyman et al., 2005], but the initial spreading behavior of Hawaiian lava flows compares well to that of a Newtonian fluid [e.g., Takagi and Huppert, 2010; Castruccio et al., 2013]. In lava flows, Bo ~105>>1 and surface tension can be ignored, although the crust of a flow may behave similarly [e.g., Hon et al., 1994; Kerr and Lyman, 2007]. Unconfined viscous flow of Newtonian fluids with negligible surface tension is well studied analytically and experimentally [e.g., Huppert, 1982; Lister, 1992], so it is this system we will use to investigate flow interactions with a topographic obstacle. This is a crucial starting point that must be understood before considering further complications such as the development of a crust. Outlining the basic behavior of unconfined viscous flows is fundamental to this goal. Viscous flow emanating from a point source (i.e., a vent) will initially spread axisymmetrically both in the down-slope (X) and cross-slope directions (Y), as defined by where ρ, g, and µ, are still density, gravity, and viscosity, Q is the volumetric flux, α is the slope, and t is time [Lister, 1992]. This flow also has a thickness (H) that scales as € Re = ρuH µ € X =Y ~ ρgQ 3 cosα µ $ % & ' ( ) 1 8 t1 2 , (3) , (4) 134 After a length of time the flow transitions to a long-time slope-dominated regime where the dimensions of the flow have a down-slope extent of a cross-slope extent equivalent to one half of the total flow width and a thickness € H ~ Qµ ρgcosα $ % & ' ( ) 1 4 € T ~ µ 3 cos5α Q ρg( )3 $ % & & ' ( ) ) 1 4 1 sin2α € X ~ ρg( ) 3Q4 sin5α 3µ( )3 cos2α $ % & & ' ( ) ) 1 9 t 7 9 € Y ~ Qcosαsinα # $ % & ' ( 1 3 t1 3 € H ~ Q 2µ3 ρg( )3 cosα sin2α $ % & & ' ( ) ) 1 9 t −1 9 . (5) , (6) , (7) , (8) . (9) 135 These scaling relationships define how eruptive conditions (e.g., Q) and pre-existing topography (e.g., α) affect the flow parameters of thickness and velocity that we want to investigate. 3. Methods Analogue experiments allow the exploration of natural phenomena at the lab scale. We performed a series of unconfined viscous gravity current experiments with an obstacle to investigate the parameters that control when flows will overtop obstacles, the local thickness of a flow upslope of an obstacle, and how obstacle interaction may affect flow behavior. For these experiments, we use golden syrup, a readily available viscous, Newtonian fluid with physical properties described in Table 1. With a room temperature viscosity of ~80 Pa s and a surface tension comparable to that of water and basaltic magma, golden syrup has frequently been used as an analogue for magma and lava in the earth sciences [e.g., Llewellin et al., 2002; Bagdassarov and Pinkerton, 2004; Castruccio et al., 2010; Beckett et al., 2011]. The viscosity of the syrup is an influential parameter, but it is highly temperature dependent. Experimental temperatures are therefore recorded and then reproduced by water bath for viscosity measurements of the syrup after the experiments by rotational rheometer (Haake RV20 with a concentric cylinder geometry in 2011, and Haake MARS with parallel plate and cone and plate geometries in 2012 and 2013, respectively). Golden syrup viscosities, shown in Figure 1, match well with the observed temperature dependence in Beckett et al. [2011]. Lower values are a product of slight dilution of the syrup with water from the pump used to extrude syrup during the 136 experiments (0.7% H2O by volume fits the most affected samples). The piston-style pump contains a leather washer at the bottom that must be soaked in water before use. This introduces a small amount of water to the base of the pump. Syrup samples for viscosity measurements were collected from the diluted base of the pump in 2011, and from the top of the pump in 2012 and 2013. Table 1. Physical properties of golden syrup at room temperature (20°C) Density (ρ) 1443 kg/m3 Viscosity 80 Pa s Surface tension (σ) 0.08 N/ma aLlewellin et al. (2002) Figure 1. Experimental viscosities measured from samples collected from the pump from which syrup is extruded in the experiments. Data compare well with measurements by Beckett et al. [2011], but show some dilution with water. Samples from 2011 were collected at the base of the pump after the end of each experiment, where a water-soaked leather washer sits during the experiment. In 2012 and 2013, syrup samples were taken from the top of the pump before the experiments and are thus less affected. 137 Since these syrup flows are gravity currents, their density is also important to quantify. The density of pure syrup ranges from 1444 kg m-3 at 18°C to 1441 kg m-3 at 22°C (F. Beckett, personal communication, 2011). Even with dilution of the syrup by water, the density of a 0.7% H2O by volume mixture at 22°C is expected to be 1438 kg m-3. We therefore do not consider density as a variable. In our experiments, the syrup flows have Reynolds numbers from 10-5 to 10-4 and are therefore in the laminar regime with viscous forces dominanting behavior. To eliminate surface tension as a significant resisting force (i.e., Bo >> 1), experiments were run at high enough fluxes and low enough slopes to produce flows with H > 6 mm which corresponds to Bo > 6. 3.1. Experimental Apparatus and Procedure Our experimental setup consists of a piston-style pump connected by tubing to a point source in an inclined plane on which an obstacle is mounted (Figure 2). The piston has a leather cap that is soaked in water overnight to ensure a tight seal. In preparation for an experiment, the piston within the pump is lowered and the pump is filled. During an experiment, the piston rises at a constant velocity, forcing the syrup through into the tubing to extrude a constant volumetric flux of syrup through the point source in the inclined plane. The plane is an 80 cm x 100 cm reinforced plastic sheet set at a slope, so that the extruding syrup produces a gravity current that intersects an obstacle attached to the surface 25 cm down slope from the point source. A caliper mounted above the obstacle allows flow thickness to be measured orthogonally to the sloped surface ~2 mm upslope of the obstacle throughout the experiment. An overhead video camera records the 138 2D flow behavior while hand-held photographs capture side-views of the obstacle interaction and otherwise document the experiment. These videos and images are analyzed with Tracker, an open source program that allows point tracking and measurements. Figure 2. Schematic drawing of the experimental setup. Variables include slope (α), flux (Q), and the obstacle. 3.2. Experiment Series Our experiments were divided into three series, each testing a different set of parameters. Series 1 investigated the effects of flux and slope on the ability of a flow to overtop an obstacle of constant shape and height; series 2 investigated the effects of flux and slope on the flow height upslope of an obstacle too tall to overtop, and series 3 investigated the effects of changing obstacle shape and size. We also tested the influence of contact angle on the flow-obstacle interaction for a range of flow and obstacle parameters. These experiments allowed us to determine the influence of both the flow and the obstacle on their interactions. 139 In the first series, we varied the experimental flux and slope and documented the advance of the resulting flow and whether it overtopped or split around an obstacle. The obstacle is a 6.25 mm tall, 5.7 cm diameter hollow polytetrafluoroethylene (PTFE, or Teflon) cylinder with a height that is close to the flow thickness. We performed these experiments at fluxes ranging from 0.5 to 2 mL/s and at slopes of 3 to 15°. An experiment was considered split fully if no syrup covered the top of the obstacle, and to overtop if syrup began flowing on top of the obstacle. Transitional experiments had a small amount of syrup that flowed onto the obstacle but then ceased to advance. No flow thicknesses were measured in these experiments. The overtopping experiments showed that a stationary wave formed upslope of the obstacle and appeared to determine the experiment outcome. To investigate this wave, we repeated the experiments at a range of fluxes and slopes, but with the same PTFE cylindrical obstacle made twice as tall to prevent overtopping. We also set up a caliper to measure the height of the stationary wave. Corresponding control experiments with no obstacle capture the form of flows without obstacle interference. Flux and slope are the parameters that govern flow behavior, but obstacle interaction may also be controlled by the nature of the obstacle itself. To investigate this, we ran a third set of experiments with different obstacle shapes and sizes. For obstacles, we used plastic cylinders and triangular (V-shaped) obstacles, pointing upslope. The original 5.7 cm PTFE cylinder was replaced by 4 cm (small) and 5.7 cm (large) diameter plastic ones (Figure 3). For the triangular obstacles, the lengths of the upslope sides of the isosceles triangle are 4 cm (small) or 30 cm (large), while the angle between them varies from 30° to 150°. We made a 180° obstacle by rotating the 90° obstacle so that the long 140 side faced upslope and created a 5.7 cm (small) or 42.4 cm (large) long wall. We ran these experiments at constant flux and slope, and measured both the height of the stationary wave and the advance rate of the flow. These experiments were repeated at five different flux and slope conditions, which also matched the previous control (i.e., obstacle-free) and cylinder experiments (0.75 mL/s and 10°, 0.5 mL/s and 15°, 0.75 mL/s and 15°, 1 mL/s and 15°, and 1.5 mL/s and 15°). Selected experiments were performed as pairs with the same obstacle shape and flow parameters, but different cylinder diameters or triangle side lengths. Figure 3. Obstacle shapes and sizes used in experiments. Finally, although the syrup flows are thick enough that the effect of surface tension is negligible for their pre-obstacle behavior, the contact angle between the syrup and the obstacle could play a significant role in shaping the stationary wave [e.g., Baxter et al., 2009]. To determine how this might effect our experimental measurements, we ran duplicate experiments with obstacles made of different materials (5.7 cm PTFE cylinder versus 5.7 cm plastic cylinder). With PTFE as the least-wetting material, we captured the fully wetting endmember by smearing a thin layer of syrup on the faces of the plastic obstacles. This produces a pre-wetted surface that reduces the contact angle to zero. Even 141 though the syrup flow fully wetted the obstacle, the flow height was still measured underneath the caliper, ~2 mm upslope. 4. Results The results of all of the experiment series capture a range of parameters that control the ability of a flow to overtop an obstacle and the flow advance behavior (Table 2). In experiments with the short, cylindrical obstacle, overtopping occurs with flows at high slope and flux. We observe significant thickening of syrup upslope of the obstacle, forming a stationary wave that can drive the flow over the top of it (Figure 4). Measuring the height of this wave through time, the flow thickness upslope of the obstacle approaches a steady state value (Figure 5), which we take as the flow height (h) in Table 2. Flow advance and width through time are captured by the overhead video camera. Video analysis provides velocity measurements, and we use the velocity of the flow when it reaches the obstacle as the representative flow front velocity (v) in Table 2. The final columns show the wave height, which is the difference between the flow height of the experiment with an obstacle (hexp) and the flow height of the control experiment (hcontrol), and the average velocity (vavg), representing the average of the velocities of the experiment and its corresponding control experiment. Differences in viscosity between the obstacle and control experiments may produce differences in the flow heights that are not due to the presence of the obstacle, but we do not see this over our viscosity range. 142 Table 2. Summary of experimental parameters and results Flux Q (mL/s) Slope α (°) Obstacle (angle θ) (°) Temperature T (°C) Viscosity µ (Pa s) Overtop?a Flow width w (± 0.20, cm) Flow velocity u (mm/s) Flow height h (mm) Wave heightb Δh (mm) Average velocityc uavg (mm/s) Overtopping experiments 2 11 PTFE O (short) 21.0 34.9 Y 20.08 1.05 ± 0.14d - - - 2 5 PTFE O (short) 21.2 42.5 Y 27.22 0.63 ± 0.10 - - - 2 3 PTFE O (short) 21.5 42.5 Y 32.02 0.47 ± 0.04 - - - 1.5 11 PTFE O (short) 21.0 34.9 Y 17.07 0.94 ± 0.07 - - - 1.5 7 PTFE O (short) 21.0 40.0 Y 23.06 0.68 ± 0.12 - - - 1.5 5 PTFE O (short) 21.2 40.0 Y 26.68 0.53 ± 0.07 - - - 1.5 3 PTFE O (short) ~21e 40.0 N 33.25 0.37 ± 0.06 - - - 1 15 PTFE O (short) 21.4 37.2 Y 13.86 0.97 ± 0.18 - - - 1 11 PTFE O (short) 20.5 45.2 Y 16.83 0.70 ± 0.04 - - - 1 7 PTFE O (short) 21.2 37.2 T 21.85 0.58 ± 0.09 - - - 1 3 PTFE O (short) 21.3 37.2 N 32.34 0.28 ± 0.02 - - - 0.75 15 PTFE O (short) 20.8 43.0 Y 13.45 0.77 ± 0.06 - - - 0.75 11 PTFE O (short) 20.7 45.2 T 15.64 0.65 ± 0.10 - - - 0.75 11 PTFE O (short) 21.5 42.5 N 15.69 0.68 ± 0.02 - - - 0.75 7 PTFE O (short) 20.8 45.2 N 20.93 0.50 ± 0.10 - - - 0.75 3 PTFE O (short) 21.1 37.2 N 30.70 0.27 ± 0.02 - - - 0.5 15 PTFE O (short) 21.0 40.0 N 11.31 0.68 ± 0.06 - - - 0.5 11 PTFE O (short) 20.7 45.2 N 13.15 0.54 ± 0.01 - - - 0.5 7 PTFE O (short) 21.3 40.0 N 18.97 0.37 ± 0.04 - - - Changing flow parameters 2 25 PTFE O 20.5 60.2 12.88 1.80 ± 0.00 16.29 ± 0.09 7.60 ± 0.19 1.85 ± 0.10 2 18 PTFE O 19.0 65.8 15.22 1.64 ± 0.08 15.29 ± 0.35 6.47 ± 0.41 1.57 ± 0.14 1.5 15 PTFE O 21.0 58.3 15.87 1.26 ± 0.08 13.33 ± 0.22 5.43 ± 0.26 1.27 ± 0.10 1 15 PTFE O 19.6 78.1 15.06 0.79 ± 0.30 11.49 ± 0.05 4.36 ± 0.17 0.82 ± 0.30 1 15 PTFE O 19.3 80.0 14.87 0.83 ± 0.08 11.20 ± 0.14 4.07 ± 0.22 0.84 ± 0.10 143 Table 2. continued Flux Q (mL/s) Slope α (°) Obstacle (angle θ) (°) Temperature T (°C) Viscosity µ (Pa s) Overtop?a Flow width w (± 0.20, cm) Flow velocity u (mm/s) Flow height h (mm) Wave heightb Δh (mm) Average velocityc uavg (mm/s) 1 10 PTFE O 19.5 55.7 19.67 0.60 ± 0.10 10.76 ± 0.32 3.37 ± 0.37 0.62 ± 0.12 0.75 15 PTFE O 18.0 101.0 15.19 0.64 ± 0.09 11.43 ± 0.09 4.51 ± 0.12 0.67 ± 0.10 0.75 10 PTFE O 17.5 85.3 18.45 0.52 ± 0.10 9.94 ± 0.42 2.56 ± 0.43 0.51 ± 0.10 0.75 5 PTFE O 19.3 82.5 25.92 0.30 ± 0.06 9.79 ± 0.17 1.91 ± 0.21 0.30 ± 0.13 0.5 10 PTFE O 19.1 68.0 16.78 0.37 ± 0.06 9.16 ± 0.11 2.89 ± 0.19 0.42 ± 0.09 0.5 3 PTFE O 19.2 65.1 25.62 0.23 ± 0.00 9.20 ± 0.22 1.85 ± 0.28 0.24 ± 0.03 2 25 None 20.5 60.2 12.30 1.90 ± 0.10 8.70 ± 0.17 - - 2 18 None 19.5 88.8 14.62 1.51 ± 0.12 8.82 ± 0.06 - - 1.5 15 None 22.0 54.5 14.62 1.28 ± 0.07 7.90 ± 0.14 - - 1 15 None 19.7 78.1 15.10 0.86 ± 0.07 7.13 ± 0.16 - - 1 10 None 19.4 55.7 17.78 0.63 ± 0.07 7.39 ± 0.19 - - 0.75 15 None 18.8 83.7 15.27 0.69 ± 0.06 6.92 ± 0.07 - - 0.75 10 None 18.8 83.7 18.33 0.50 ± 0.01 7.38 ± 0.10 - - 0.75 5 None 19.4 80.6 25.20 0.31 ± 0.12 7.88 ± 0.12 - - 0.5 15 None 19.0 68.0 11.61 0.66 ± 0.11 5.82 ± 0.03 - - 0.5 10 None 19.0 68.0 14.83 0.46 ± 0.06 6.27 ± 0.16 - - 0.5 3 None 20.0 62.1 21.92 0.26 ± 0.03 7.35 ± 0.18 - - Changing obstacle shape 1.5 15 Small 30 19.2 70.0 16.34 1.03 ± 0.11 10.92 ± 0.19 3.01 ± 0.24 1.16 ± 0.13 1.5 15 Small 60 19.5 70.0 15.88 1.02 ± 0.04 12.08 ± 0.20 4.18 ± 0.24 1.15 ± 0.08 1.5 15 Small 90 20.0 82.6 16.60 1.12 ± 0.22 11.83 ± 0.17 3.92 ± 0.22 1.20 ± 0.23 1.5 15 Small 90 19.5 70.0 16.85 1.05 ± 0.09 12.43 ± 0.15 4.53 ± 0.21 1.17 ± 0.12 1.5 15 Small 120 18.5 101.4 16.99 1.02 ± 0.12 14.73 ± 0.09 6.83 ± 0.17 1.15 ± 0.14 1.5 15 Small 150 19.0 89.3 16.78 1.04 ± 0.09 16.17 ± 0.12 8.27 ± 0.19 1.16 ± 0.11 1.5 15 Small 180 19.2 70.0 17.12 1.03 ± 0.09 16.38 ± 0.20 8.47 ± 0.25 1.15 ± 0.11 1 15 Small 60 19.2 68.0 14.64 0.85 ± 0.06 8.85 ± 0.22 1.73 ± 0.28 0.85 ± 0.10 144 Table 2. continued Flux Q (mL/s) Slope α (°) Obstacle (angle θ) (°) Temperature T (°C) Viscosity µ (Pa s) Overtop?a Flow width w (± 0.20, cm) Flow velocity u (mm/s) Flow height h (mm) Wave heightb Δh (mm) Average velocityc uavg (mm/s) 1 15 Small 90 18.0 96.0 15.92 0.78 ± 0.12 10.46 ± 0.43 3.34 ± 0.46 0.82 ± 0.14 1 15 Small 120 18.3 94.3 16.06 0.79 ± 0.07 12.00 ± 0.20 4.88 ± 0.26 0.82 ± 0.10 1 15 Small 180 19.0 68.0 15.15 0.89 ± 0.07 13.92 ± 0.36 6.80 ± 0.40 0.87 ± 0.10 0.75 15 Small 30 18.0 98.2 15.79 0.64 ± 0.16 8.74 ± 0.17 1.82 ± 0.18 0.66 ± 0.17 0.75 15 Small 60 17.6 106.3 15.30 0.58 ± 0.14 9.29 ± 0.15 2.37 ± 0.17 0.64 ± 0.15 0.75 15 Small 90 18.2 101.0 14.47 0.66 ± 0.10 9.77 ± 0.50 2.85 ± 0.50 0.68 ± 0.11 0.75 15 Small 90 18.0 96.0 14.38 0.65 ± 0.15 9.24 ± 0.24 2.32 ± 0.25 0.69 ± 0.06 0.75 15 Small 120 17.6 106.3 15.61 0.65 ± 0.16 11.73 ± 0.21 4.81 ± 0.23 0.67 ± 0.17 0.75 15 Small 180 17.9 98.2 15.49 0.60 ± 0.12 13.69 ± 0.78 6.77 ± 0.79 0.64 ± 0.13 0.75 10 Small 30 18.9 82.0 18.40 0.49 ± 0.14 8.41 ± 0.11 1.03 ± 0.15 0.49 ± 0.14 0.75 10 Small 60 19.0 82.0 18.02 0.49 ± 0.05 8.59 ± 0.37 1.21 ± 0.38 0.49 ± 0.05 0.75 10 Small 90 19.1 79.0 17.25 0.53 ± 0.12 9.58 ± 0.08 2.21 ± 0.13 0.51 ± 0.12 0.75 10 Small 120 19.2 79.0 17.62 0.52 ± 0.03 10.26 ± 0.24 2.88 ± 0.26 0.51 ± 0.03 0.75 10 Small 180 19.2 79.0 18.41 0.51 ± 0.13 11.66 ± 0.18 4.29 ± 0.21 0.50 ± 0.13 0.5 15 Small 90 18.8 68.0 12.66 0.57 ± 0.07 8.62 ± 0.30 2.80 ± 0.30 0.61 ± 0.13 Changing obstacle size 2 25 Large O 20.0 74.9 13.47 1.68 ± 0.18 17.32 ± 0.10 8.62 ± 0.19 1.79 ± 0.32 2 18 Large O 20.0 84.8 15.83 1.39 ± 0.05 16.33 ± 0.33 7.51 ± 0.34 1.45 ± 0.33 1.5 15 Large O 19.0 89.3 17.25 1.04 ± 0.14 14.92 ± 0.22 7.02 ± 0.27 1.16 ± 0.34 1 15 Large O 18.1 104.0 15.18 0.78 ± 0.02 13.47 ± 0.08 6.35 ± 0.18 0.82 ± 0.35 1 10 Large O 21.0 69.8 19.61 0.60 ± 0.04 12.40 ± 0.10 5.01 ± 0.21 0.61 ± 0.36 2 25 Small O 20.0 74.9 12.95 1.86 ± 0.16 16.08 ± 0.15 7.39 ± 0.23 1.88 ± 0.37 2 18 Small O 19.5 81.6 16.09 1.36 ± 0.10 15.00 ± 0.12 6.18 ± 0.14 1.44 ± 0.38 1.5 15 Small O 20.0 84.8 16.81 1.07 ± 0.05 14.06 ± 0.09 6.16 ± 0.17 1.18 ± 0.39 1 15 Small O 20.0 84.8 16.34 0.79 ± 0.05 12.38 ± 0.10 5.25 ± 0.19 0.82 ± 0.40 1 10 Small O 19.1 90.0 20.30 0.56 ± 0.05 11.98 ± 0.12 4.59 ± 0.22 0.59 ± 0.41 145 Table 2. continued Flux Q (mL/s) Slope α (°) Obstacle (angle θ) (°) Temperature T (°C) Viscosity µ (Pa s) Overtop?a Flow width w (± 0.20, cm) Flow velocity u (mm/s) Flow height h (mm) Wave heightb Δh (mm) Average velocityc uavg (mm/s) 1.5 15 Large 30 19.0 88.2 17.26 0.98 ± 0.03 11.80 ± 0.11 3.89 ± 0.18 1.13 ± 0.42 1.5 15 Large 60 20.0 74.9 17.13 1.00 ± 0.14 12.17 ± 0.22 4.26 ± 0.26 1.14 ± 0.43 1.5 15 Large 120 18.0 101.5 17.68 0.98 ± 0.20 17.23 ± 0.14 9.33 ± 0.20 1.13 ± 0.44 1.5 15 Large 150 18.5 78.2 16.80 1.08 ± 0.17 21.42 ± 0.17 13.51 ± 0.23 1.18 ± 0.45 0.75 15 Large 60 20.0 49.0f 14.38 0.83 ± 0.08 8.60 ± 0.26 1.68 ± 0.41 0.76 ± 0.46 0.75 15 Large 120 19.0 88.2 15.26 0.67 ± 0.01 14.00 ± 0.07 7.08 ± 0.10 0.68 ± 0.47 0.75 10 Large 30 18.0 101.5 19.45 0.45 ± 0.04 10.08 ± 0.07 2.70 ± 0.12 0.47 ± 0.48 0.75 10 Large 90 21.0 58.3 17.99 0.56 ± 0.04 10.15 ± 0.32 2.77 ± 0.39 0.53 ± 0.49 0.75 10 Large 180 18.5 94.8 19.03 0.48 ± 0.04 24.37 ± 0.10 16.99 ± 0.14 0.49 ± 0.50 Changing the contact angle 1.5 15 Syrup O 20.2 77.8 16.71 1.04 ± 0.06 15.02 ± 0.26 7.12 ± 0.30 1.16 ± 0.32 1 15 Syrup O 18.9 92.7 16.41 0.77 ± 0.02 13.95 ± 0.58 6.82 ± 0.61 0.81 ± 0.33 1 10 Syrup O 20.0 84.8 20.02 0.45 ± 0.04 12.91 ± 0.45 5.51 ± 0.49 0.62 ± 0.34 1.5 15 Syrup 60 19.0 70.0 16.23 1.07 ± 0.15 11.67 ± 0.10 3.77 ± 0.17 1.17 ± 0.35 1.5 15 Syrup 90 21.0 70.2 16.54 1.10 ± 0.05 12.53 ± 0.21 4.62 ± 0.25 1.19 ± 0.36 1.5 15 Syrup 120 18.2 105.2 16.83 1.00 ± 0.05 15.19 ± 0.35 7.29 ± 0.38 1.14 ± 0.37 aY = yes, N = no, T = transitional bΔh = hexp – hcontrol cuavg = (uexp + ucontrol)/2 dAll errors are ± 2σ eTemperature was not measured, 21°C is assumed from experiments done immediately before and after. fThis viscosity is derived from the best-fitting model. The measured value of 22.4 Pa s is anomalously low. 146 Figure 4. Photographs of stationary waves on the upslope side of the obstacle. Figure 5. Flow height upslope of two 4 cm triangular obstacles (30° blue, 120° red) plotted with time. The flow height plateaus, providing a value for steady-state flow height that we use in our analysis. Flow thickness at the stationary wave shows relationships with syrup flux, slope, obstacle shape and size. The thickness relates to the flux and slope of the experiment for a constant obstacle (PTFE cylinder), reflected in the flow velocity when the syrup impacts the obstacle (Figure 6a). The flow height also increases with increasing obstacle angle, at constant flux, slope, and obstacle side length (Figure 6b). The five different flux and slope conditions (symbol colors) produce trends with obstacle angle that have different slopes, showing the combined effect of flux, slope, and obstacle angle. Obstacle size also influences the flow height, as increasing the obstacle size raises the height of the 147 wave behind the obstacle (Figure 7). The 1.7 cm increase in cylinder diameter corresponds to a ~1 mm increase in flow height over a range of experimental conditions (Figure 7a). The larger triangles yield no statistically significant change in flow height at low obstacle angles, but the difference becomes large at high angles (Figure 7b). Figure 6. Flow height measurements for (a) PTFE large cylinders plotted against flow velocity, and (b) plastic triangular obstacles plotted against obstacle angle. Error bars show two standard deviations. Figure 7. Comparison of the flow heights of stationary waves behind small and large obstacles: (a) plastic cylinders and (b) plastic triangles as in Figure 3. 148 Our results also demonstrate that surface tension plays a role in the experimental obstacle interaction. In the first experiments, flows that are thicker than the obstacle are at times unable to overtop it, suggesting the potential influence of surface tension that could hold back the flow at the scale of the thickness above the obstacle. By later using obstacles that are much taller than the flow, this problem disappears, but the contact angle between the syrup and the obstacle still shapes the intersection between the two [e.g., Baxter et al., 2009]. Comparing experimental flow heights of plastic obstacles with those that are either made of PTFE or coated in syrup highlights this (Figure 8). The PTFE obstacles yield significantly lower flow heights, on average 1.5 mm lower, than their corresponding experiments with plastic obstacles (blue symbols). However, the plastic and syrup-coated plastic experiments (red symbols) are not statistically distinguishable. The differences (or lack thereof) appear to originate from the shape of the wave against the obstacle (Figure 8 photographs), where the fully wetting syrup-coated obstacle pulls the surface of the wave upward, while the extremely hydrophobic PTFE pulls the surface downward. The wave height is not measured against the obstacle, but about 2 mm back, reducing this effect to some extent. While the contact angle impacts the measured flow height values, it is independent of the trends in the experimental results. Our experiments can therefore still provide insight into the physics of viscous flow-obstacle interactions. 5. Discussion These experimental results can help develop our understanding of the physics of obstacle interaction and its implications for lava flow behavior. We first assess the behavior of a viscous flow without an obstacle, and then analyze the factors that control 149 Figure 8. Comparison of the flow height upslope of the obstacle as measured behind obstacles with different surface properties. Flow heights in experiments with plastic obstacles are compared to those with equivalent PTFE obstacles (blue symbols) or syrup- coated obstacles (red symbols). The differences appear to correspond to the shapes of the wave as it intersects the obstacle at a different contact angle. Photographs show flat syrup surface against a plastic obstacle, a steep, wetting surface against a syrup-coated plastic obstacle, and a plunging, high contact angle against a PTFE obstacle. obstacle interaction and its influence on flow advance. Analytical models by Lister [1992] allow the calculation of flow lengths, widths, and velocities before obstacle intersection, but we must develop new theory to describe the obstacle-flow interaction and its aftermath. 5.1. Viscous Flow without an Obstacle We check our experimental data of the evolution of flows with time against the original Lister [1992] model to confirm that the equations apply. In this formulation, at early times the flow should be axisymmetric, with the down-slope (X) and cross-slope 150 (Y) coordinates growing with t1/2 as in equation (4). Later, the underlying slope beings to dominate the behavior and X ~ t7/9, while Y ~ t1/3 [equations (7) and (8)]. We compare these scaling relationships with the down-slope and cross-slope advance data from our control experiments without an obstacle (Figure 9). If plotted on a log-log plot, the slopes of the X and Y spreading represent the exponent of t in the model, and they match the theory well (Figure 9 triangles). This plot also demonstrates that once the flow reaches the obstacle (25 cm, or where it disappears underneath the caliper), the long-time sloping flow regime has been reached. We therefore use the sloping regime equations to compare to the flow behavior near the obstacle. When all of the control experiments are analyzed in this way, we find mean early-time exponents of 0.58 ± 0.07(1σ) for X (~1.17/2) and 0.56 ± 0.15 for Y (~1.11/2), and long-time exponents of 0.82 ± 0.08 for X (~7.36/9) and 0.36 ± 0.10 for Y (~1.08/3). Although the values are variable between experiments, the means match three of the expected exponents within error. The exception is the early- time exponent for X, which falls statistically above the expected value of 0.5. Therefore, the long-time scaling equations of Lister [1992] match well to our experiments and can be used to predict flow behavior. The Lister [1992] model also provides a scaling relationship for maximum flow height in the sloping regime in equation (9). While this does not predict the height of the flow at a specific point, it can be used to calculate the maximum height of the flow when the flow front has reached the obstacle, which occurs within the sloping regime. The predicted values parallel the measured steady-state flow heights (Figure 10). We use this comparison to calculate coefficients to convert equation (9) into a formula that can be 151 Figure 9. Evolution of the down-slope and cross-slope dimensions of an example experimental gravity current (1.5 mL/s at 15° slope with no obstacle). Initial spreading is proportional to t1/2 [equation (4)], while eventually down-slope flow grows as t7/9 [equation (7)] and cross-slope flow grows at t1/3 [equation (8)]. used to predict the flow height at the location of the obstacle for a given flux and slope, 5.2. Analysis of Overtopping Experiments We apply our understanding of viscous flows on a slope without an obstacle to discern how the flow parameters might influence the outcomes of obstacle interaction. We hypothesized that flows that are thicker than the obstacle height will be able to flow over it, but that run-up of the flow onto the obstacle could allow thinner flows at high € Hpredict =1.106 Q2µ3 ρg( )3 cosα sin2α $ % & & ' ( ) ) 1 9 t −1 9 + , - - . / 0 0 +1.069 . (10) 152 Figure 10. Comparison of the flow height at 25 cm down slope from the vent calculated with equation (9) [Lister, 1992] and measured as the steady-state thickness of the control experiments. The values are consistently offset. velocity to also overtop. We compare our overtopping results to flow heights and flow velocities at the location of the obstacle to test these hypotheses (Figure 11). These plots use calculated heights and velocities without the presence of an obstacle to build continuous contours, as we have shown that these represent experimental values well. In our overtopping experiments, flows were able to overtop the short, cylindrical obstacle at high flux and high slope, which corresponds best to high velocities (Figure 11b). However, it is likely that where the flow height is greater than the obstacle height (to the right of the 6.25 mm contour in Figure 11a), surface tension prevents the flow from overtopping. It is therefore clear that flow velocity exerts a strong control on obstacle overtopping, but flow height must play a role, too. 153 Figure 11. Overtopping experiment results compared to (a) flow heights calculated for all experiments with equation (10), and (b) flow velocities calculated from equation (7). Both are calculated at the location of the obstacle. 5.3. The Influence of Flow Velocity The result that flow velocity can control obstacle overtopping supports the principle of a velocity-driven run-up height. At high experimental velocities, we observe local thickening within a stationary wave upslope of the obstacle that could be produced by a run-up style interaction (Figure 6a). With experimental flow velocities ranging from 0.2 to 1.9 mm/s, equation (1) predicts run-up height (Δh) values between 10-6 and 10-4 mm. Observed Δh values (or wave heights) can be calculated from the difference between the flow height upslope of an obstacle (hexp) and the flow height of the corresponding control experiment (hcontrol) (Table 2). These values range from 1 to 17 mm, and are many orders of magnitude larger than calculated flow run-up by energy conversion. Scaling relationships between velocity (u) and wave height (Δh) also refute the application of the run-up height model. In the run-up height equation, Δh is proportional to u2 [equation (1)]. If we plot the average velocity of the two experiments 154 used to calculate the wave height, uavg, versus Δh, we instead see a power-law relationship with an exponent of 0.69 (Figure 12). Therefore, some kinetic to potential energy conversion may be producing a run-up wave feature, but it is an insignificant portion of the height of the stationary wave that we observe in these flows where viscous forces dominate over intertial forces. Figure 12. Plot of flow velocity at the location of the obstacle versus the wave height, Δh, calculated as the height of the flow behind the obstacle (hexp) minus the height of the flow with no obstacle (hcontrol, from the control experiment). Velocities are the average velocities (uavg) of the two experiments (Table 2). 5.4. Influence of Obstacle Shape and Size Obstacle shape and size, by influencing the height of the stationary wave, may also control whether an obstacle may be overtopped by a flow. While the stationary wave is strongly affected by the flow velocity, obstacle shape plays a role, as well. In results 155 from experiments with varying obstacle angles (small triangles), we observe positive linear trends between obstacle angle and flow height (Figure 6b). The slopes and intercepts of these trends vary with the experimental parameters (flux and slope, shown in the colors). This is because the flux and slope determine the flow velocity and the flow height in the absence of an obstacle [equations (7) and (9)]. To isolate the impact of the obstacle angle, we use the wave height (Δh) to remove the effect of flow height (without an obstacle) and normalize by the velocity. The observed relationship between wave height and velocity (Figure 12) suggests that dividing Δh by uavg0.69 should collapse all points onto one line with a near-zero intercept (Figure 13). It is also possible to use this result to suggest a dimensionless ratio that summarizes the behavior. The 0.69 exponent is approximately equally to two-thirds, which means the ratio has units of (mm⅓ s-⅔). Therefore, we could, for example, multiply by an expression containing gravity (g, constant in these experiments) and a length scale (l, representing obstacle size, for instance), to produce a dimensionless expression with the same distribution, This is just an example without a strong physical basis, but it hints at the possibility of summarizing the influence of obstacle angle and velocity in a way that is potentially scalable to lava flows. € Δh uavg2 3 g l2 # $ % & ' ( 1 3 . (11) 156 Figure 13. Plot of obstacle angle compared to wave height divided by average velocity0.69 (Figure 12), for the small triangle experiments. Nondimensionalization of the y-axis is discussed in the main text. Results from experiments with differently sized obstacles demonstrate that obstacle size can additionally effect stationary wave height. Large cylinders produce wave heights that are consistently offset from small cylinders at a range of fluxes and slopes (Figure 7a, hcontrol is the same for hlarge and hsmall at each point so wave height and flow height have the same offset). The difference between the wave heights upslope of large and small obstacles is the same for all represented flow velocities, suggesting that the influence of obstacle size is independent of parameters controlling the flow behavior. When we compare large and small triangles, however, we see a secondary effect from the obstacle angle (Figure 7b). In this set of experiments, the stationary wave height does not appear to depend on obstacle size at low obstacle angles, but becomes strongly dependent at high obstacle angles. This behavior, shown explicitly in Figure 14, can be 157 described by a power law relationship between obstacle angle and the difference in the flow height with the large obstacle (hlarge) and the flow height with the small obstacle (hsmall). Without an obstacle, the height difference is zero but as the obstacle angle increases, the large obstacle produces a much higher stationary wave than the small one. Moreover, the effect appears independent of flow conditions (i.e., flux and slope). Figure 14. Obstacle angle compared to the difference between the flow heights upslope of a large, 30 cm triangular obstacle and a small, 4 cm one of the same angle at the same flux and slope (hlarge - hsmall). A best-fit power law is also shown. 5.5. Interpretations from the Stationary Wave Analysis The behavior of the stationary wave offers insight into the underlying physics of the wave and therefore the physics of obstacle interaction and overtopping. We witness steady-state behavior upslope of the obstacle (e.g., Figure 5) indicating that the amount of syrup that impacts the obstacle must be balanced by the amount of syrup that flows around it. The steadiness of the stationary wave must represent a force balance where the 158 gravitational force (due to flow thickness and slope, along with gradients in syrup thickness) is balanced by the viscous resisting force. The incoming syrup is governed by the flux and slope of the experiment, along with the width of the obstacle, while the syrup making its way around the obstacle can be described as a combination of forward (down- slope) and lateral (cross-slope) flow. The important parameters we have identified – flow velocity, obstacle shape, and obstacle size – may all influence this balance and thus the height of the stationary wave. Other variables, such as viscosity and density, may cancel out, yielding the lack of viscosity effect apparent in our data. These three influential parameters govern the flux of syrup hitting the obstacle, which can be described as uHW, where u is the velocity, H is an average flow thickness, and W is the width of the syrup flow that runs into the obstacle. When the flow is wider than the obstacle, which is true of all of the experiments with small obstacles and large triangles with a 30° angle, this W is the obstacle width Wobs, where L is the triangle side length and θ is the obstacle angle. As a result, the effective flux into the obstacle (Qeff) is uHWobs. For all experiments with large triangles at obstacle angles larger than 30°, the flow is narrower than Wobs, and the entire flux of the flow (Q) intersects the obstacle. Defining Qeff helps to explain the observed positive correlations between velocity and wave height (Figure 12), obstacle angle and wave height (Figure 13), and obstacle size and wave height (Figure 7). € Wobs = 2Lsin θ 2 # $ % & ' ( , (12) 159 In the experiments where Wobs is greater than the flow width, the flux is no longer dependent on obstacle size or shape, but increasing obstacle angle still has a dramatic effect on the flow height (Figure 7b). This is likely because the obstacle angle also regulates the degree to which the flow must travel across, rather than down, the slope. The total distance Y* the flow front must move across the slope is Y* = Wobs/2 and the only driving force for cross-slope motion is gravitational collapse. Added thickening behind the obstacle would build up a thickness gradient to help propel the flow laterally. A secondary effect of forcing flow at an angle to the slope is the change in the effective slope underlying the flow. Flow along a triangular obstacle angle of θ on a surface of slope α will travel along an effective gradient (αeff), At a lower effective slope, the slope-driven component of flow decreases while the height gradient driver of flow increases. From this perspective, the two possible endmembers are: (1) flow with no obstacle (an obstacle angle of 0°) and (2) flow into a horizontal wall (the 180° obstacle). Flow advance in these two cases will be controlled by different stresses. Flow without an obstacle should be purely slope-driven, while flow along a 180° obstacle should be entirely thickness gradient-driven, and thus axisymmetric. In between, the flow is driven by both slope and thickness gradient. With the obstacle impeding downhill flow, thickening drives lateral flow, while the lower effective slope restricts sloping flow (Figure 15). The geometry also served to confine the flow against the obstacle as it flows € αeff = sin−1 sinα cos θ 2 % & ' ( ) * + , - . / 0 . (13) 160 along it. Because all of the flow along the obstacle is at an angle to the downhill direction, the syrup is confined against the obstacle on the downhill side and gravity on the uphill side. This causes both flow thickening against the obstacle and reduced flow spreading. Figure 15. Annotated picture of flow into and along a 150° large obstacle. Flow into the obstacle is driven by the experimental slope, α, while flow along the obstacle is driven by a combination of the height gradient (highlighted) and the effective slope, αeff. Gravity continues to act on the along-obstacle flow, confining it against the obstacle and limiting its ability to spread upslope. Confinement of the flow by the obstacle also affects the rate of advance of the flow along the obstacle margins. Against the obstacle, all forces driving flow front advance are projected along the obstacle, as illustrated by down-slope (X) and cross- slope (Y) flow advance data. Flow front positions in X and Y along large obstacles with angles up to 150° scale as t(7/9), as expected for X in the sloping regime [equation (7), Figure 16a]. At 180°, the behavior is different, with an exponent of ~1/2 and scales identically to the axisymmetric, gravitational collapse driven regime [equation (4)]. The 161 X and Y positions scale similarly because they follow the margins of the obstacle almost exactly, so that Y/X is approximately equal to tan(θ/2) for each experiment (Figure 16b). Figure 16. (a) Best-fit exponents of the along obstacle advance in the down-slope (X) and cross-slope (Y) directions of flows around large obstacles compared with obstacle angles. (b) Relative advance of the flow in the X and Y directions along the obstacle is dictated entirely by the obstacle angle. The ratio of advance, Y/X ≈ tan(θ/2) at every point in every experiment. The observation that X and Y scale as t(7/9) suggests that it may be possible to model flow along the obstacle with equation (7), but using half the flux and the effective slope along the obstacle. Importantly, this model assumes a fully sloping regime, and disregards the effect of thickening at the point of obstacle intersection. It also neglects the necessary condition that the flow must approach the fully axisymmetric regime, and be driven entirely by the height gradient, as the obstacle angle approaches 180˚. Furthermore, flow along the obstacle is confined, limiting the lateral spreading that is part of the formulation of equation (7). Both of these factors will increase the flow velocity along the obstacle. In fact, a comparison of the along obstacle advance (D) recorded in the experiments and modeled with equation (7) illustrates this effect (Figure 17). The experiments are distinctly faster than the model, and their offset is largely a function of 162 obstacle angle. It is possible to collapse the model results onto the one-to-one line by multiplying by a function of the obstacle angle, which projects the behavior onto the down-slope coordinate and empirically fits the experimental data. Incorporating the height gradient and confinement against the obstacle by gravity would be necessary to model this successfully, but the Lister [1992] sloping regime model serves as a good base for future work. Figure 17. Experimental data for flow advance along the obstacle (D) compared to modeled advance from equation (7) of slope-dominated flow down the effective slope (αeff) along the obstacle at half the original flux (Q/2). The model is offset, and in a manner related to the obstacle angle (θ). € −1.9 cos θ2 $ % & ' ( ) $ % & ' ( ) −1 + 3.2, (14) 163 With advance along the obstacle margin and limited flow width, the flow along the obstacle could also be considered nearly unidirectional. Indeed, some of the experiments show that X and Y may scale closer to t1, indicating that a simple model like the Jefferys equation [Jefferys, 1925], might apply. Here ualong corresponds to the depth-averaged velocity for two-dimensional, laminar flow, of constant h down a slope. Without flow height measurements along the obstacle, we use the along-obstacle velocities to calculate the flow height anticipated by equation (15). In a plot of the Jefferys height (HJefferys) against the total flow height just upslope of the obstacle, the two correlate well but are not equal (Figure 18). This implies that the increased flow velocity along the obstacle may be driven by an increase in flow thickness related to the flow height within the stationary wave. The offset between them may be because of thickening along the obstacle from confinement. The model and experimental data also diverge with increasing obstacle angle. This model assumes a fully slope-driven regime, and is therefore less and less valid at higher obstacle angles. The divergence may therefore be due to the transition to height gradient driven flow. We can additionally take a broader look at the flow advance rates over the time before the flow reaches the obstacle, along the obstacle, and past the obstacle. Example flow advance plots (Figures 19a and b) show the initial spreading of the flow, which matches the Lister [1992] equation (7) model well (Mod Q). This is followed by a period of flow along the obstacle in which the velocity is constant and the flow diverges from € ualong = ρgsinαeff h2 3µ , (15) 164 Figure 18. Plot of the flow height predicted by the Jefferys equation (15) along the obstacle margin compared to the flow height within the stationary wave upslope of the obstacle. the model; the flow then shows a reduction in velocity once the flow leaves the obstacle. We plot ratios of the flow velocities in the pre-obstacle, along obstacle, and post- obstacle segments to highlight changes along flow and deviations from the model. A comparison of the experimental along-obstacle velocity (Exp ualong) with the pre-obstacle velocity (Exp upre, measured as the velocity immediately before the obstacle) shows that the flow is faster along the obstacle at moderate obstacle angles (Figure 19c black points). A comparison of Exp ualong with the model velocity (Mod Q ualong) illustrates a similar acceleration of the flow along the obstacle at the mid-range obstacle angles, relative to the expected velocity of an identical flow that did not hit an obstacle and split (Figure 19c gray points). Faster velocities along moderate obstacle angles suggests that 165 Figure 19. (a and b) Flow advance of two examples of experimental flows before the obstacle (pre), along the obstacle, and after the obstacle (post) compared to flow modeling using equation (7). (c) Plot showing velocity ratios of the experimental velocity along the obstacle (Exp ualong) relative to the experimental velocity before the obstacle (Exp upre) and modeled velocity without an obstacle (Mod Q ualong). This highlights acceleration of the flow along the obstacle relative to the model with no obstacle and the pre-obstacle velocity at moderate obstacle angles. (d) Plot showing velocity ratios of the post-obstacle experimental velocity (Exp upost) relative to the modeled velocity without an obstacle (Mod Q upost) and the modeled velocity that divides Q by 2 at the obstacle location (Mod Q/2 upost). The experimental post-obstacle velocity is slower than the expected velocity without an obstacle, but faster than simply dividing by the flux should produce. the combined effect of slope-driven flow and height gradient-driven flow is maximized at an obstacle angle of ~90°, which is a cross-slope angle of 45°. At these angles, the flow may be thickened enough to yield greater flow velocities [e.g., equation (15)], but still has a large enough αeff to efficiently drive flow down slope [e.g., equation (13)]. 166 We use this same ratio analysis to determine the velocity difference after the flow leaves the obstacle. We model the advance of a flow with half of the original flux beyond the location of the obstacle (Mod Q/2) to demonstrate the expected impact of dividing the flow without any obstacle interaction effects like thickening or confinement. Here the experimental velocities (Exp upost) are lower than the model velocities without flow splitting (Mod Q upost), but larger than the model velocities with divided flux (Mod Q/2 upost) (Figure 19d). There is also a slight dependence on obstacle angle in these data, suggesting that while the flow may be slower once split and unconfined, factors such as flow thickening and narrowing along the obstacle may speed it up to some degree. 5.6. Qualitative Summary of Observations and Rationalizations Flow into, around, and past an obstacle is therefore a complex problem with a mixture of different flow regimes and geometric controls. We can use the parameters and models described above to separate out this range of behavior and qualitatively rationalize our observations. The increase in flow height upslope of the obstacle is the primary observation of the interaction between the flow and the obstacle, and this steady- state behavior is a function of the flow into and around the obstacle. We find that increased flow velocity causes the wave height to increase (Figure 12). An increase in flow velocity causes a greater syrup flux into the obstacle (Qeff) and consequently the flow thickness just upstream of the obstacle increases until the flux of syrup around the obstacle matches the incoming flux. The larger Δh facilitates flow around the obstacle by thickening the slope-dominated component of flow [which is proportional to h2 in 167 equation (15) and Figure 18] and by increasing the flow thickness gradient driving axisymmetric flow (e.g., Figure 15). The influence of obstacle shape and size on wave height can be similarly outlined. When the obstacle angle increases, the wave height increases (Figure 13). This can be explained by changes in both the flux into the obstacle, and movement of syrup around the obstacle. Increasing obstacle angle causes widening [equation (12)] and thus increases the amount of syrup hitting the obstacle. Moreover, with a larger obstacle angle, the effective slope along the obstacle is reduced [equation (13)] and the distance the flow must travel across the slope (Y*) increases. This requires that the flow be driven more and more by the thickness gradient, which causes the flow to thicken until lateral flow balances the incoming flux. A similar balance is caused by increasing obstacle size. Increasing either cylinder diameter or triangle side length (L) produces a wider obstacle and an equivalent increase in flux (until the total flux is captured). For a constant angle there is no change in effective slope, consistent with the similar observed behavior of 30° and 60° obstacles of different size (Figures 7b and 14). However, the height gradient needed for flow around the obstacle is expected to scale as flow height divided by obstacle width, which explains the observed increase in wave heights upslope of the high- angle, large obstacles relative to the small ones. The observations thus support a model where flux into the obstacle is balanced by flow around the obstacle and driven by a combination of slope and height gradient. We cannot, however, incorporate all three forces into a single parametric model with a simple scaling analysis. It may be possible, however, to develop end-member sloping and axisymmetric models, as in Lister [1992]. Further end-member experiments and analysis 168 to completely parameterize and nondimensionalize flow interaction with an obstacle is necessary to produce such a quantitative description. However, even our new qualitative understanding of the interaction between lava flows and underlying topography provides important insight into conditions that modulate flow advance, and effective design of lava flow barriers. 6. Conclusions These experimental results and our model framework have implications for lava flow interactions with topography, diversion barrier interventions, and velocity estimates from run-up heights. Our analysis demonstrates that the behavior of a flow intersecting an obstacle is a not simply a function of flow height or flow velocity, but also dependent on obstacle geometry. These parameters collectively control the height of a steady-state wave upslope of the obstacle, balancing the flow into the obstacle with the flow around it. These new insights inform lava flow modeling and mitigation efforts. Our results have greatly improved our understanding of lava flow response to underlying topography. The new data suggest a model for the scale and geometry of topography that may influence the flow, and how branching affects flow behavior. Prior modeling efforts have used an approximate flow height to determine the scale of topography that a flow may surmount [Favalli et al., 2005]. Our analysis indicates that flow height must play a role, but that flow velocity, flow width, and obstacle geometry are also important. Natural obstacles are not perfect cylinders or triangles, but range in shape and size considerably. The digital elevation model (DEM) used for modeling therefore needs a resolution that captures the scale of obstacle a flow may split around. 169 With further analysis to scale our results to lava flows, it may be possible to quantitatively assess this DEM resolution, given flow properties such as thickness and velocity. Our experiments also indicate the importance of modeling not simply a steepest descent path, but a flow that may branch and merge, to capture changes in flow advance rate associated with flow splitting. Future experimental work is required to similarly quantify the effects of flow merging. Collectively branching and merging could have large implications for flow lengths and advance rates that would be valuable to model [e.g., Dietterich and Cashman, in press]. The artificial topography of lava flow diversion barriers more closely matches the obstacle geometries in our experiments. The experimental results could therefore be scaled for application to predict the response of a lava flow to a V-shaped barrier. Our analysis reveals that barrier height cannot solely be determined based on flow height, but must incorporate the potential range of flow velocities and widths, along with the barrier length and internal angle. With all of these parameters, we could use a scaled model to predict the flow height upslope of a barrier with a given length and angle and calculate the required barrier height. However, this would ignore the effects of a crust that may not be negligible. The experiments also provide a test of the use of run-up height as a way to measure flow velocities in viscous flows. This method is frequently applied, but assumes a complete one-dimensional transfer of kinetic to potential energy. We find that in our syrup flows, this model yields miniscule (~10-5 mm) run-up heights, despite the mm-scale waves that we see forming upslope of the obstacle. The relationship between wave height and velocity in our data (Figure 12) also disagrees with the scaling implied by the run-up 170 height equation (1). This analysis highlights the fact that in lava and syrup, viscous forces are most often greater than intertial forces (Re < 1), and this simple intertial model fails to describe their behavior. Overall, the run-up height model ignores many parameters that we determine to be important, including obstacle shape and size, as well as parameters that we did not test a wide range of values for but may be important, including viscosity and density. Although we have delineated a number of important implications of this work, there are many avenues for building on it in the future. We determine that a balance of flow into and around the obstacle describes flow behavior during obstacle interaction, but that flow around the obstacle is governed by a combination of slope-driven and height gradient-driven flow. We require further experiments to fully capture the end-member behaviors, and a detailed mathematical analysis to encapsulate both flow regimes into a single scaled model that could be applied to lava flows interacting with both artificial and natural topography. As part of this, it is also necessary to determine the effects of cooling and crystallization, as these cause lava flows to form levees and contribute to the forces resisting flow [Kerr et al., 2006]. Other geometries present in artificial and natural obstacle interactions would also be valuable to explore experimentally, including interaction with oblique walls like those used to divert flows at Mt. Etna [Barberi et al., 2003] or those that capture flow branches at Kīlauea [Dietterich and Cashman, in press]. These topographic interactions can create flow confluences, which produce faster, longer flows that would be important to predict [Dietterich and Cashman, in press]. Our current results combined with these future directions, would develop the tools necessary for physics-based quantitative predictive models of lava flow interactions with topography. 171 CHAPTER VI CONCLUSIONS In this dissertation I have developed new techniques to collect data on lava flow emplacement in Hawai‘i and improved our understanding of how underlying topography and channel networks influence lava flow behavior. I motivated this work by outlining the goals of being able to monitor, predict, and mitigate lava flow hazards, and my research contributes to each of these. In Chapters II and III, I discuss new remote sensing methods to map active lava flows through time and collect three-dimensional views of flows that can enhance our knowledge of flow emplacement. Chapters IV and V use these data to analyze how flow interactions with underlying topography produce complex channel network geometry, and how this geometry in turn influences flow behavior. These chapters demonstrate ways to improve our ability to predict and mitigate flows by developing a fundamental understanding of the factors that govern flow emplacement. I will end by summarizing the conclusions and implications of each chapter and discussing some future directions for research. 1. Dissertation Summary In Chapter II, I illustrate a new method for mapping active lava flows with satellite remote sensing. This technique uses SAR imagery (already collected for the purposes of measuring deformation at many active volcanoes) and exploits it to build a time series of flow maps. The new technique frequently has better spatial and temporal 172 coverage than traditional field-based mapping, and provides a way to map internal flow features such as lava tubes and surface breakouts. After flow emplacement, I find that the SAR coherence data can be used to estimate flow thicknesses, providing a potential three- dimensional mapping tool. This technique could be applied to map any surface change at volcanoes or elsewhere, and can be placed in the toolbox of monitoring techniques at volcanoes around the world. In Chapter III, I describe how lidar topographic data can contribute to our record flow emplacement. With airborne laser altimetry data, we collect high resolution topographic and intensity data of the surface of lava flows. Together, these datasets can be used to map flows of different ages, and morphological structures such as flow channels and levees. The intensity of the lidar signal changes with flow age and the extent to which the flows have become vegetated after emplacement. Meter-scale surface roughness cannot distinguish flows of different ages but it does vary between ‘a‘ā and pāhoehoe morphologies and can be used to highlight internal structures within flows. In this chapter, I also use aerial photographs to create a pre-eruptive digital elevation model (DEM), which allowed me to examine the role of topography in flow emplacement and to build a three-dimensional flow model to calculate lava flow thickness and total volume. I use this to track the volume distribution between major, topographically delineated branches in the flow, as well as to discover field evidence supporting flow thinning over steep topography. I also discern flow thinning over steep topography and thickening at locations of flow branching. These results stress the importance of underlying topography as a control of lava flow behavior. 173 In Chapter IV, I use the flow data to investigate the relationship between underlying topography, the branching and merging channel network preserved in the present-day topography, and the associated flow and channel morphology and behavior. By comparing pre-eruptive slope to network complexity, I find that at higher slopes the flow thins and branches more. This suggests a relationship between flow thickness and whether the flow must split around topographic obstacles or simply overtop them. Branching because of obstacle interactions may create distributary networks. Obstacles can also cause confinement, including pre-eruptive gullies, levees, and fault scarps, which produce flow confluences, tributary networks, and long single-channel flow segments. The channel network itself can further influence the flow. Channel network connectivity and longevity control the channel morphology, such that long-lived, well- connected channels are narrower. To examine impacts of channel network geometry on flow behavior, I analyzed records from a number of Hawaiian lava flows and discovered that distributary flows are slower and shorter than confined, tributary flows. My results demonstrate the importance of considering the lava channel network in any hazard analysis. Finally, in Chapter V, I describe experiments designed to systematize the topographic controls on lava flow emplacement presented in Chapters III and IV. In Chapter III I saw flow thickening at bifurcations, which in Chapter IV I relate to field evidence for velocity-controlled lava flow run-up onto topographic obstacles dividing the flow. I also identified underlying slope and flow thinning as factors in controlling the extent of flow branching. This observation suggests flow height as a control for whether a flow splits around or overtops an obstacle. I therefore set up a series of analogue 174 experiments with a syrup lava flow that intersected a cylindrical or triangular obstacle to test these ideas and otherwise determine the important parameters that influence the interaction between a flow and an obstacle. My experimental results showed that flow height, velocity, and width, along with obstacle size and shape, all affect whether a flow overtops an obstacle. The primary control on overtopping is the size of the stationary wave that forms upslope of the obstacle. I model this behavior using the balance of syrup hitting the obstacle to the syrup flowing around it, which is a combination of flow driven from the underlying slope and flow driven by the flow thickness gradient. In agreement with the behavior of natural flows, the syrup flows also show a decrease in flow advance rate once split. These results exhibit the need to incorporate flow dynamics into predictive models of flow advance, as well as simulations used to calculate the dimensions needed for flow diversion barriers. My analyses also demonstrate the potential problems with using of lava run-up height as a means of calculating flow velocities. 2. Future Research Directions While these chapters provide much new insight into lava flow emplacement behavior, they also supply the foundation for future research directions. The results in this dissertation can inform a number of new techniques and analysis. I discuss potential research to pursue in remote sensing, channel network analysis, and experimental modeling. The remote sensing techniques described in this dissertation are new and still in development with many future potential routes. The SAR flow mapping in Chapter II 175 used data from a single SAR satellite that has since gone offline. Now, however, a pair of SAR satellites have been simultaneously imaging the same regions (TerraSAR-X and TanDEM-X), and will allow construction of DEMs through time. Combined with my initial planform flow extent mapping, we hope to use these data to produce temporal maps of flow thickness. Such a dataset would allow us to calculate lava effusion rates, and create a time series of three-dimensional maps similar to the flow thickness map described in Chapter III. Each DEM would also provide pre-eruptive topography for the next. This would provide unprecedented imaging of flow field development through time, paired with imaging of the underlying surface that would allow us to further investigate the influence on underlying topography on flow emplacement. Our use of photogrammetry to derive pre-eruptive DEMs in Chapter III can also be extended to provide very high-resolution topographic mapping of present-day flow morphology and even time-lapse DEMs. Recent advances in computer vision have enhanced traditional stereophotogrammetry to produce structure from motion, a method of assembling three-dimensional data from a large suite of two-dimensional images taken at many angles [e.g., Snavely et al., 2007]. We have begun to use this technique to take field photographs from on top of a tall pole or from a helicopter and use them to reconstruct the lava flow morphology at decimeter resolution. Paired with ground control, these data allow us to quantify small-scale flow features, such as the profile of a stationary wave located where a channel in the Mauna Loa 1984 flow branches (Figure 1). With synced time-lapse cameras placed around an area of interest, it is possible to collect sets of photographs and therefore DEMs through time. This has great potential for 176 use in the field and for laboratory experiments to capture three-dimensional flow morphology through time. Figure 1. Oblique hillshade view of a structure from motion DEM looking down flow at a 4.7 m tall stationary wave where the channel divides around an obstacle within the Mauna Loa 1984 flow. James Dietrich produced the 10 cm resolution DEM from photographs taken from a helicopter using Agisoft PhotoScan. Temporal constraints on flow emplacement would allow us to test many of the results of our channel network analysis, both about the formation of branches and confluences, and how the channel network developments and impacts flow morphology and behavior. Repeated lidar surveys during an eruption at Mt. Etna demonstrate some of what is possible, with the ability to track lava volumes down channels, flow velocities, and how channel morphology evolves with time [Favalli et al., 2010]. With pre-eruptive data, the relationship between the channel network, channel morphology, and flow behavior could be fully assessed. Repeat SAR DEMs could have a resolution as high as three meters, while time-lapse structure from motion could do even better at a small scale. Although we have focused on Hawaiian ‘a‘ā flows in our analysis, these new data 177 sources could also be used to evaluate pāhoehoe flow emplacement. Pāhoehoe flows travel quite slowly, but are common and frequently long-lived. They also can form complex lava tube networks that may significantly impact their behavior and evolution. Therefore, a similar analysis of how factors such as underlying topography influence the structure and evolution of pāhoehoe flows would be valuable. Further work could also be done to quantify topographic interactions and flow networks in lavas with different rheologies, including different compositions and crystallinities, or even extraterrestrial or subaqueous flows. An important way to extend our results to other systems is to pursue further experiments and a fully parameterized, scaled model for the interaction between a lava flow and an obstacle, and its effects in terms of flow splitting and influence on flow behavior. This requires more experiments at end-member (fully slope-driven or height gradient driven) geometries, a greater range of other potentially important parameters like flow viscosity, and quantifying the effect of cooling. With this full exploration of the problem experimentally, we can pursue a more complete model of this behavior that will apply to a whole range of flow and obstacle conditions. Utilizing the results of our syrup experiments to predict lava flow behavior is particularly dependent on incorporating the effects of cooling and crystallization. Lava flow widths, velocities, and thicknesses are all influenced by the changes in viscosity and the addition of a crust that accompany these processes [e.g., Lyman et al., 2005; Kerr et al., 2006]. However, cooling and crust formation are difficult to duplicate with syrup. We therefore have begun to perform similar experiments with molten basalt, which is heated in a furnace and then poured at a constant rate onto a slope with an obstacle. In pilot 178 experiments at the Syracuse University Lava Project, we witness stationary wave formation upslope of the obstacle along with a reduction in flow advance rate after flow splitting, validating our syrup experiments (Figure 2). We also find that the obstacle angle and flux control whether the flow is able to merge after splitting or whether both branches form levees and travel parallel to each other. Experiments at a range of fluxes, slopes, and obstacle angles are still needed for a full analysis. Figure 2. Example data from a molten basalt experiment. (a) Flow thickness map highlights flow thickening upslope of the obstacle. The image is created from time-lapse structure from motion with a DEM made by James Dietrich from cameras mounted around the experiment. (b) Flow advance plot showing a marked decrease in flow advance rate when the flow splits, with example time steps shown in (c) and (d). There are also a number of experimental geometries that would be of interest to pursue in syrup and molten basalt experiments. Most obstacles are not cylindrical or V- shaped, and many used for flow diversion are simply linear and at an angle to the flow 179 [e.g., Barberi et al., 2003]. These replicate features described in Chapters III and IV, such as the faults scarps that rerouted the Kīlauea December 1974 flow. Oblique barriers can have a significant effect on flow emplacement by confining the flow and producing flow confluences. It would therefore be particularly helpful to investigate overtopping and flow advance effects of angled linear barriers. Obstacles with non-vertical walls are also important to consider, as most artificial and natural obstacles have slopes no steeper than the angle of repose. The work described in this dissertation advances our understanding of lava flow emplacement and promotes the use of remote sensing technologies to quantify lava flow morphology and behavior. The further pursuits outlined here serve to demonstrate the applicability of the techniques and analysis presented in this dissertation to expanding our observations and knowledge of lava flows. Insights gained from this work can be employed to combat the hazards that lava flows pose to people and property through monitoring active flows, informing diversion attempts, and predicting flow paths, velocities, and lengths. Better awareness of the relationships between underlying topography, flow morphology, and flow behavior also help to interpret the features recorded in flows that have recently erupted or long been frozen. These remote sensing and analysis tools I have summarized can aid in documenting lava flow emplacement on other planets or in extracting flow histories necessary to anticipate future hazards at volcanoes around the world. With this work, we take a step toward being able to fully describe and comprehend these dynamic phenomena. 180 REFERENCES CITED CHAPTER II Briole, P., D. Massonnet, and C. Delacourt (1997), Post-eruptive deformation associated with the 1986–87 and 1989 lava flows of Etna detected by radar interferometry, Geophys. Res. Lett., 24, 37–40, doi:10.1029/96GL03705. Calabro, M. D., D. A. Schmidt, and J. J. Roering (2010), An examination of seasonal deformation at the Portuguese Bend landslide, southern California, using radar interferometry, J. Geophys. Res., 115, F02020, doi:10.1029/2009JF001314. Crisci, G. M., S. Di Gregorio, O. Pindaro, and G. Ranieri (1986), Lava flow simulation by a discrete cellular model: first implementation, Int. J. Model. Simul., 6, 137–140. Farr, T., and M. Kobrick (2000), Shuttle Radar Topography Mission produces a wealth of data, Eos Trans. AGU, 81, 583–585. Gatelli, F., A. M. Guamieri, F. Parizzi, P. Pasquali, C. Prati, and F. Rocca (1994), The wavenumber shift in SAR interferometry, IEEE Trans. Geosci. Remote Sens., 32, 855–865, doi:10.1109/36.298013. Harris, A. J. L., and S. K. Rowland (2001), FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel, Bull. Volcanol., 63, 20–44, doi:10.1007/s004450000120. Heliker, C., and T. N. Mattox (2003), The first two decades of the Pu’u ’O’o-Kūpaianaha eruption: Chronology and selected bibliography, in The Pu’u ’Ō’ō-Kūpaianaha Eruption of Kīlauea Volcano, Hawai’i: The First 20 Years, edited by C. Heliker, D. A. Swanson, and T. J. Takahashi, U.S. Geol. Surv. Prof. Pap., 1676, 1–28. Kahle, A. B., M. J. Abrams, E. A. Abbott, P. J. Mouginis-Mark, and V. J. Realmuto (1995), Remote sensing of Mauna Loa, in Mauna Loa Revealed: Structure, Composition, History, and Hazards, Geophys. Monogr. Ser., vol. 92, edited by J. M. Rhodes and J. P. Lockwood, pp. 145–170, AGU, Washington, D. C., doi:10.1029/GM092p0145. Kauahikaua, J. (2007), Lava flow hazard assessment, as of August 2007, for Kīlauea east rift zone eruptions, Hawai‘i Island, U.S. Geol. Surv. Open-File Rept. 2007-1264, 9 p. 181 Kauahikaua, J., S. Margriter, J. Lockwood, and F. Trusdell (1995), Applications of GIS to the estimation of lava flow hazards on Mauna Loa Volcano, Hawai'i, in Mauna Loa Revealed: Structure, Composition, History, and Hazards, Geophys. Monogr. Ser., vol. 92, edited by J. M. Rhodes and J. P. Lockwood, pp. 315–325, AGU, Washington, D. C., doi:10.1029/GM092p0315. Kauahikaua, J., K. V. Cashman, T. N. Mattox, C. C. Heliker, K. A. Hon, M. T. Mangan, and C.R. Thornber (1998), Observations on basaltic lava streams in tubes from Kilauea Volcano, island of Hawaii, J. Geophys. Res., 103, 27, 303–27, 323, doi:10.1029/97JB03576. Kauahikaua, J., D. R. Sherrod, K. V. Cashman, C. Heliker, K. Hon, T. N. Mattox, and J. A. Johnson (2003), Hawaiian lava–flow dynamics during the Pu’u O’o-Kupianaha eruption: A tale of two decades, in The Pu’u ’Ō’ō-Kūpaianaha Eruption of Kīlauea Volcano, Hawai’i: The First 20 Years, edited by C. Heliker, D. A. Swanson, and T. J. Takahashi, U.S. Geol. Surv. Prof. Pap., 1676, 63–87. Lu, Z., D. Mann, J. T. Freymueller, and D. J. Meyer (2000), Synthetic aperture radar interferometry of Okmok volcano, Alaska: Radar observations, J. Geophys. Res., 105, 10, 791–10, 806, doi:10.1029/2000JB900034. Lu, Z., T. Masterlark, and D. Dzurisin (2005), Interferometric synthetic aperture radar study of Okmok volcano, Alaska, 1992–2003: Magma supply dynamics and postemplacement lava flow deformation, J. Geophys. Res., 110, B02403, doi:10.1029/2004JB003148. Lu, Z., and J. T. Freymueller (1998), Synthetic aperture radar interferometry coherence analysis over Katmai volcano group, Alaska, J. Geophys. Res., 103, 29, 887–29, 894, doi:10.1029/98JB02410. Massonnet, D., and K. L. Feigl (1998), Radar interferometry and its application to changes in the Earth's surface, Rev. Geophys., 36, 441–500, doi:10.1029/97RG03139. Mattox, T. N., C. Heliker, J. Kauahikaua, and K. Hon (1993), Development of the 1990 Kalapana flow field, Kilauea volcano, Hawaii, Bull. Volcanol., 55, 407–413, doi:10.1007/BF00302000. Mazzarini, F., M. T. Pareschi, M. Favalli, I. Isola, S. Tarquini, and E. Boschi (2007),Lava flow identification and aging by means of lidar intensity: Mount Etna case, J. Geophys. Res., 112, B02201, doi:10.1029/2005JB004166. Orr, T.R. (2010), Lava tube shatter rings and their correlation with lava flux increases at Kīlauea Volcano, Hawai’i, Bull. Volcanol., 73, 335–346, doi:10.1007/s00445-010- 0414-3. 182 Patrick, M. R., and T. R. Orr (2012), Rootless shield and perched lava pond collapses at Kīlauea Volcano, Hawai'i, Bull. Volcanol., 74, 67–78, doi:10.1007/s00445- 0110505-9. Patrick M. R., T. Orr , D. Wilson, D. Dow, and R. Freeman (2011), Cyclic spattering, seismic tremor, and surface fluctuation within a perched lava channel, Kīlauea Volcano, Bull. Volcanol., 73, 639–653, doi:10.1007/s00445-010-0431-2. Poland, M., A. Miklius, T. Orr, J. Sutton, C. Thornber, and D. Wilson (2008), New episodes of volcanism at Kilauea Volcano, Hawaii, Eos Trans. AGU, 89, 37–38, doi:10.1029/2008EO050001. Rowland, S. K., A. J. L. Harris, M. J. Wooster, F. Amelung, H. Garbeil, L. Wilson, and P. J. Mouginis-Mark (2003), Volumetric characteristics of lava flows from interferometric radar and multispectral satellite data: the 1995 Fernandina and 1998 Cerro Azul eruptions in the western Galápagos, Bull. Volcanol., 65, 311–330, doi:10.1007/s00445-002-0262-x. Stevens, N. F., G. Wadge, C. A. Williams, J. G. Morley, J.-P. Muller, J. B. Murray, and M. Upton (2001), Surface movements of emplaced lava flows measured by synthetic aperture radar interferometry, J. Geophys. Res., 106, 11, 293–11, 313, doi:10.1029/2000JB900425. Wadge, G., P. A. V. Young, and I. J. McKendrick (1994), Mapping lava flow hazards using computer simulation, J. Geophys. Res., 99, 489–504, doi:10.1029/93JB01561. Wadge, G., B. Scheuchl, and N. F. Stevens (2002), Spaceborne radar measurements of the eruption of Soufrière Hills Volcano, Montserrat, in The Eruption of Soufrière Hills Volcano, Monserrat, from 1995 to 1999, edited by T. H. Druitt and B. P. Kokelaar, Mem. Geol. Soc. Lond., 21, 583–594, doi:10.1144/GSL.MEM.2002.021.01.27. Zebker, H. A., and J. Villasenor (1992), Decorrelation in interferometric radar echoes, IEEE Trans. Geosci. Remote Sens., 30, 950–959, doi:10.1109/36.175330. Zebker, H. A., P. Rosen, S. Hensley, and P. J. Mouginis-Mark (1996), Analysis of active lava flows on Kilauea volcano, Hawaii, using SIR-C radar correlation measurements, Geology, 24, 495–498, doi:10.1130/0091- 7613(1996)024<0495:AOALFO>2.3.CO;2. Zebker, H. A., P. A. Rosen, and S. Hensley (1997), Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps, J. Geophys. Res., 102, 7547–7563, doi:10.1029/96JB03804. 183 CHAPTER III Anderson-Teixeira, K. J., P. M. Vitousek, and J. H. Brown (2008), Amplified temperature dependence in ecosystems developing on the lava flows of Mauna Loa, Hawai'i, P. Natl. Acad. Sci., 105, 228–233, doi:10.1073/pnas.0710214104. Aplet, G. H. and P. M. Vitousek (1994), An Age--Altitude Matrix Analysis of Hawaiian Rain-Forest Succession. J. Ecol., 82, 137–147. Aplet, G. H., R. F. Hughes, and P. M. Vitousek (1998), Ecosystem development on Hawaiian lava flows: biomass and species composition, J. Veg. Sci., 9, 17–26, doi:10.2307/3237219. Applegarth, L. J., H. Pinkerton, M. R. James, and S. Calvari (2010), Morphological complexities and hazards during the emplacement of channel-fed ‘a‘ā lava flow fields: A study of the 2001 lower flow field on Etna, Bull. Volcanol., 72, 641–656, doi:10.1007/s00445-010-0351-1. Branca, S., E. De Beni, and C. Proietti (2013), The large and destructive 1669 AD eruption at Etna volcano: reconstruction of the lava flow field evolution and effusion rate trend, Bull. Volcanol., 75, 694, doi:10.1007/s00445-013-0694-5. Bretar, F., M. Arab-Sedze, J. Champion, M. Pierrot-Deseilligny, E. Heggy, and S. Jacquemoud (2013), An advanced photogrammetric method to measure surface roughness: Application to volcanic terrains in the Piton de la Fournaise, Reunion Island, Remote Sens. Environ., 135, 1–11, doi:10.1016/j.rse.2013.03.026. Campbell, B. A., and M. K. Shepard (1996), Lava flow surface roughness and depolarized radar scattering, J. Geophys. Res., 101(E8), 18941–18951, doi:10.1029/95JE01804. Coltelli, M., C. Proietti, S., Branca, M. Marsella, D. Andronico, and L. Lodato (2007), Analysis of the 2001 lava flow eruption of Mt. Etna from three-dimensional mapping, J. Geophys. Res., 112, F02029, doi:10.1029/2006JF000598. Crisp, J., K. V. Cashman, J. A. Bonini, S. B. Hougen, and D. C. Pieri (1994), Crystallization history of the 1984 Mauna Loa lava flow, J. Geophys. Res., 99(B4), 7177–7198, doi:10.1029/93JB02973. Deardorff, N. D. and K. V. Cashman (2012), Emplacement conditions of the c. 1,600- year BP Collier Cone lava flow, Oregon: a LiDAR investigation, Bull. Volcanol., 74, 2051–2066, doi: 10.1007/s00445-012-0650-9. 184 Deligne, N. I., K. V. Cashman, and J. J. Roering (2012), After the lava flow: the importance of external soil sources for plant colonization of recent lava flows in the central Oregon Cascades, USA, Geomorphology, 202, 15–32, doi:10.1016/j.geomorph.2012.12.009. Dietterich, H. R., M. P. Poland, D. A. Schmidt, K. V. Cashman, D. R. Sherrod, and A. T. Espinosa (2012), Tracking lava flow emplacement on the east rift zone of Kīlauea, Hawai‘i, with synthetic aperture radar coherence, Geochem. Geophys. Geosyst., 13, Q05001, doi:10.1029/2011GC004016. Favalli, M., M. T. Pareschi, A. Neri, and I. Isola (2005), Forecasting lava flow paths by a stochastic approach, Geophys. Res. Lett., 32, L03305, doi:10.1029/2004GL021718. Favalli, M., A. Fornaciai, and M. T. Pareschi (2009a) LIDAR strip adjustment: Application to volcanic areas, Geomorphology, 111, 123–135, doi:10.1016/j.geomorph.2009.04.010. Favalli, M., F. Mazzarini, M. T. Pareschi, and E. Boschi (2009b), Topographic control on lava flow paths at Mount Etna, Italy: Implications for hazard assessment, J. Geophys. Res., 114, F01019, doi:10.1029/2007JF000918. Favalli, M., A. J. L. Harris, A. Fornaciai, M. T. Pareschi, and F. Mazzarini (2010), The distal segment of Etna’s 2001 basaltic lava flow, Bull. Volcanol., 72, 119–127, doi:10.1007/s00445-009-0300-z. Favalli, M., S. Tarquini, P. Papale, A. Fornaciai, and E. Boschi (2012), Lava flow hazard and risk at Mt. Cameroon volcano, Bull. Volcanol., 74, 423–439, doi:10.1007/s00445-011-0540-6. Gaddis, L., P. Mouginis-Mark, and J. A. Hayashi (1990), Lava flow surface textures- SIR-B radar image texture, field observations, and terrain measurements, Photogramm. Eng. Rem. S., 56, 211–224. Griffiths, R. W. (2000), The dynamics of lava flows, Annu. Rev. Fluid Mech., 32, 477– 518, doi:10.1146/annurev.fluid.32.1.477. Hollingsworth, J., S. Leprince, F. Ayoub, and J.-P. Avouac (2013), New constraints on dike injection and fault slip during the 1975–1984 Krafla rift crisis, NE Iceland, J. Geophys. Res. Solid Earth, 118, 3707–3727, doi:10.1002/jgrb.50223. Horn, B. K. P. (1981), Hill shading and the reflectance map, Proceedings of the IEEE, 69(1), 14–47, doi:10.1109/PROC.1981.11918. Jackson, T. A. (1971), A study of the ecology of pioneer lichens, mosses, and algae on recent Hawaiian lava flows, Pacific Science, 25, 22–32. 185 Kauahikaua, J., D.R. Sherrod, K.V. Cashman, C. Heliker, K. Hon, T.N. Mattox, and J.A. Johnson (2003), Hawaiian lava-flow dynamics during the Pu’u Ō’ō-Kupianaha eruption: A tale of two decades, in The Pu’u ’Ō’ō-Kupaianaha Eruption of Kīlauea Volcano, Hawai’i: The First 20 Years, edited by C. Heliker et al., U.S. Geol. Surv. Prof. Pap., 1676, pp. 63–87. Lipman, P. W. and N. G. Banks (1987), Aa flow dynamics, Mauna Loa 1984, in Volcanism in Hawaii, edited by R. W. Decker et al., U.S. Geol. Surv. Prof. Pap., 1350, pp. 1527–1567. Lister, J. R. (1992), Viscous flows down an inclined plane from point and line sources, J. Fluid Mech., 242, 631–653, doi:10.1017/S0022112092002520. Lockwood, J. P., N. G. Banks, T. T. English, L. P. Greenland, D. B. Jackson, D. J. Johnson, R. Y. Koyanagi, K. A. McGee, A. T. Okamura, and J. M. Rhodes (1985). The 1984 eruption of Mauna Loa Volcano, Hawaii, Eos Transactions, 66(16), 169–171, doi:10.1029/EO066i016p00169-01. Lockwood, J. P., R. I. Tilling, R. T. Holcomb, F. Klein, A. Okamura, and D. Peterson (1999), Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawaii, U.S. Geol. Surv. Prof. Pap., 1613. Lu, Z., T. Masterlark, and D. Dzurisin (2005), Interferometric synthetic aperture radar study of Okmok volcano, Alaska, 1992–2003: Magma supply dynamics and postemplacement lava flow deformation, J. Geophys. Res., 110, B02403, doi:10.1029/2004JB003148. Marsella, M., C. Proietti, A. Sonnessa, M. Coltelli, P. Tommasi, and E. Bernardo (2009), The evolution of the Sciara del Fuoco subaerial slope during the 2007 Stromboli eruption: Relation between deformation processes and effusive activity, J. Volcanol. Geotherm. Res., 182, 201–213, doi:10.1016/j.jvolgeores.2009.02.002. Marsella, M., P. Baldi, M. Coltelli, and M. Fabris (2012), The morphological evolution of the Sciara del Fuoco since 1868: reconstructing the effusive activity at Stromboli volcano, Bull. Volcanol., 74, 231–248, doi:10.1007/s00445-011-0516-6. Mazzarini, F., M. T. Pareschi, M. Favalli, I. Isola, S. Tarquini, and E. Boschi (2005), Morphology of basaltic lava channels during the Mt. Etna September 2004 eruption from airborne laser altimeter data, Geophys. Res. Lett., 32, L04305, doi:10.1029/2004GL021815. Mazzarini, F., M. T. Pareschi, M. Favalli, I. Isola, S. Tarquini and E. Boschi (2007), Lava flow identification and aging by means of lidar intensity: Mount Etna case, J. Geophys. Res., 112, B2, doi: 10.1029/2005JB004166. 186 McKean, J., and J. Roering (2004), Objective landslide detection and surface morphology mapping using high-resolution airborne laser altimetry, Geomorphology, 57(3–4), 331–351, doi:10.1016/S0169-555X(03)00164-8. Moore, H. J., (1987), Preliminary estimates of the rheological properties of the 1984 Mauna Loa lava, in Volcanism in Hawaii, edited by R. W. Decker et al., U.S. Geol. Surv. Prof. Pap., 1350, pp. 1569–1588. Morris, A. R., F. S. Anderson, P. J. Mouginis-Mark, A. F. C. Haldemann, B. A. Brooks, and J. Foster (2008), Roughness of Hawaiian volcanic terrains, J. Geophys. Res., 113, E12007, doi:10.1029/2008JE003079. Neri, M., F. Mazzarini, S. Tarquini, M. Bisson, I. Isola, B. Behncke, and M.T. Pareschi (2008), The changing face of Mount Etna’s summit area documented with Lidar technology, Geophys. Res. Lett., 35, L09305, doi:10.1029/2008GL033740. Passalacqua, P., T. D. Trung, E. Foufoula-Georgiou, G. Sapiro, and W. E. Dietrich (2010), A geometric framework for channel network extraction from lidar: Nonlinear diffusion and geodesic paths, J. Geophys. Res., 115, F01002, doi:10.1029/2009JF001254. Pyle, D. M. and J. R. Elliot (2006), Quantitative morphology, recent evolution, and future activity of the Kameni Islands volcano, Santorini, Greece, Geosphere, 2, 253– 268, doi:10.1130/GES00028.1. Raich, J.W., A. E. Russell, and P. M. Vitousek (1997), Primary productivity and ecosystem development along an elevational gradient on Mauna Loa, Hawai‘i, Ecology, 78, 707–721. Raich, J.W., W. J. Parton, A. E. Russell, R. L. Sanford Jr, and P. M. Vitousek (2000), Analysis of factors regulating ecosystem development on Mauna Loa using the Century model, Biogeochemistry, 51, 161–191, doi:10.1023/A:1006495408992. Riker, J. M., K. V. Cashman, J. P. Kauahikaua, and C. M. Montierth (2009), The length of channelized lava flows: Insight from the 1859 eruption of Mauna Loa Volcano, Hawai‘i, J. Volcanol. Geotherm. Res., 183, 139–156, doi:10.1016/j.jvolgeores.2009.03.002. Shepard, M. K., R. A. Brackett, and R. E. Arvidson (1995), Self-affine (fractal) topography: Surface parameterization and radar scattering, J. Geophys. Res., 100(E6), 11709–11718, doi:10.1029/95JE00664. Soule, S.A., K. V. Cashman, and J. P. Kauahikaua (2004) Examining flow emplacement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i, Bull. Volcanol., 66, 1–14, doi:10.1007/s00445-003- 0291-0. 187 Stevens, N. F., G. Wadge, C. A. Williams, J. G. Morley, J.-P. Muller, J. B. Murray, and M. Upton (2001), Surface movements of emplaced lava flows measured by synthetic aperture radar interferometry, J. Geophys. Res., 106(B6), 11293–11313, doi:10.1029/2000JB900425. Tarquini, S., and M. Favalli (2010), Changes of the susceptibility to lava flow invasion induced by morphological modifications of an active volcano: the case of Mount Etna, Italy, Nat. Hazards, 54, 537–546, doi:10.1007/s11069-009-9484-y. Tarquini, S., M. Favalli, F. Mazzarini, I. Isola, and A. Fornaciai (2012), Morphometric analysis of lava flow units: Case study over LIDAR-derived topography at Mount Etna, Italy, J. Volcanol. Geotherm. Res., 235–236, 11–22, doi:10.1016/j.jvolgeores.2012.04.026. Ventura, G. and G. Vilardo (2008), Emplacement mechanism of gravity flows inferred from high resolution Lidar data: The 1944 Somma-Vesuvius lava flow (Italy), Geomorphology, 95, 223–235, doi:10.1016/j.geomorph.2007.06.005. Wolfe, E. W. and J. Morris (1996), Geologic Map of the Island of Hawaii, U.S. Geol. Surv. Misc. Invest. Ser., I-2524-A, Washington, D.C., USA. Zevenbergen, L. W. and C. R. Thorne (1987), Quantitative analysis of land surface topography, Earth Surf. Proc. Land., 12, 47–56, doi:10.1002/esp.3290120107. Zhu, A., J. E. Burt, M. Smith, R. Wang, and J. Gao (2008), The impact of neighbourhood size on terrain derivatives and digital soil mapping, in Advances in Digital Terrain Analysis, edited by Q. Zhou et al., Springer, Berlin, doi:10.1007/978-3- 540-77800-4_18. CHAPTER IV Barberi, F., and M. L. Carapezza (2004), The control of lava flows at Mt. Etna, in Mt. Etna: Volcano Laboratory, Geophys. Monogr. Ser., vol. 143, edited by A. Bonaccorso et al., pp. 357–369, AGU, Washington, D. C., doi:10.1029/143GM22. Baxter, S. J., H. Power, K. A. Cliffe, and S. Hibberd (2009), Three-dimensional thin film flow over and around an obstacle on an inclined plane, Phys. Fluids, 21, 032102, doi:10.1063/1.3082218. Bruno, B. C., G. J. Taylor, S. K. Rowland, P. G. Lucey, and S. Self (1992), Lava flows are fractals, Geophys. Res. Lett. 19, 305–308, doi:10.1029/91GL03039. Byrne, P. K., C. Klimczak, D. A. Williams, D. M. Hurwitz, S. C. Solomon, J. W. Head, F. Preusker, and J. Oberst (2013), An assemblage of lava flow features on Mercury, J. Geophys. Res. Planets, 118, 1303–1322, doi:10.1002/jgre.20052. 188 Cashman, K. V., C. Thornber, and J. P. Kauahikaua (1999), Cooling and crystallization of lava in open channels, and the transition of Pāhoehoe Lava to 'A'ā, Bull. Volcanol., 61, 306–323, doi:10.1007/s004450050299. Coltelli, M., C. Proietti, S., Branca, M. Marsella, D. Andronico, and L. Lodato (2007), Analysis of the 2001 lava flow eruption of Mt. Etna from three-dimensional mapping, J. Geophys. Res., 112, F02029, doi:10.1029/2006JF000598. Crisci, G. M., M. V. Avolio, B. Behncke, D. D’Ambrosio, S. Di Gregorio, V. Lupiano, M. Neri, R. Rongo, and W. Spataro (2010), Predicting the impact of lava flows at Mount Etna, Italy, J. Geophys. Res., 115, B04203, doi:10.1029/2009JB006431. Crisp, J., K. V. Cashman, J. A. Bonini, S. B. Hougen, and D. C. Pieri (1994), Crystallization history of the 1984 Mauna Loa lava flow, J. Geophys. Res., 99(B4), 7177–7198, doi:10.1029/93JB02973. Crown, D. A., and S. M. Baloga (1999), Pahoehoe toe dimensions, morphology, and branching relationships at Mauna Ulu, Kilauea Volcano, Hawai'i, Bull. Volcanol., 61(5), 288–305, doi:10.1007/s004450050298. Deardorff, N. D. and K. V. Cashman (2012), Emplacement conditions of the c. 1,600- year BP Collier Cone lava flow, Oregon: a LiDAR investigation, Bull. Volcanol., 74, 2051–2066, doi: 10.1007/s00445-012-0650-9. Dietterich, H. R., S. A. Soule, K. V. Cashman, and B. H. Mackey (in press), Lava Flows in 3D: Using airborne lidar and pre-eruptive topography to evaluate lava flow surface morphology and thickness in Hawai‘i, in Hawaiian Volcanoes: From Source to Surface, Geophys. Monogr. Ser., AGU, Washington, D. C. Edmonds, D. A., and R. L. Slingerland (2007), Mechanics of river mouth bar formation: Implications for the morphodynamics of delta distributary networks, J. Geophys. Res., 112, F02034, doi:10.1029/2006JF000574. Egozi R., and P. Ashmore (2008), Defining and measuring braiding intensity, Earth Surf. Proc. Land., 33, 2121–2138, doi:10.1002/esp.1658. Favalli, M., M. T. Pareschi, A. Neri, and I. Isola (2005), Forecasting lava flow paths by a stochastic approach, Geophys. Res. Lett., 32, L03305, doi:10.1029/2004GL021718. Favalli, M., A. Fornaciai, F. Mazzarini, A. Harris, M. Neri, B. Behncke, M. T. Pareschi, S. Tarquini, and E. Boschi (2010a), Evolution of an active lava flow field using a multitemporal LIDAR acquisition, J. Geophys. Res., 115, B11203, doi:10.1029/2010JB007463. 189 Favalli, M., A. J. L. Harris, A. Fornaciai, M. T. Pareschi, and F. Mazzarini (2010b), The distal segment of Etna’s 2001 basaltic lava flow, Bull. Volcanol., 72, 119–127, doi:10.1007/s00445-009-0300-z. Gleyzer, A., M. Denisyuk, A. Rimmer, and Y. Salingar (2004), A fast recursive GIS algorithm for computing Strahler stream order in braided and nonbraided networks, J. Am. Water Resour. As., 40, 937–946. doi:10.1111/j.1752- 1688.2004.tb01057.x. Guest, J. E., C. R. J. Kilburn, H. Pinkerton, and A. M. Duncan (1987), The evolution of lava flow-fields: observations of the 1981 and 1983 eruptions of Mount Etna, Sicily, Bull. Volcanol., 49, 527–540, doi:10.1007/BF01080447. Hamilton, C. W., L. S. Glaze, M. R. James, and S. M. Baloga (2013), Topographic and stochastic influences on pāhoehoe lava lobe emplacement, Bull. Volcanol., 75, 756, doi:10.1007/s00445-013-0756-8. Harris, A. J. L., and S. K. Rowland (2001), FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel, Bull. Volcanol., 63, 20–44, doi:10.1007/s004450000120. Heliker, C., G. E. Ulrich, S. C. Margriter, and J. P. Hoffmann (2001), Maps showing the development of the Puʻu ʻŌʻō - Kūpaianaha flow field, June 1984-February 1987, Kīlauea Volcano, Hawaii, Geologic Investigations Series Map I-2685, U.S. Geological Survey, Washington, D.C., USA. Heliker C., and T. N. Mattox (2003), The first two decades of the Pu‘u ‘Ō‘ō-Kupaianaha eruption: Chronology and selected bibliography, in The Pu’u ’Ō’ō-Kupaianaha Eruption of Kīlauea Volcano, Hawai’i: The First 20 Years, edited by C. Heliker, D. A. Swanson, and T. J. Takahashi, U.S. Geol. Surv. Prof. Pap., 1676, 1–28. Hidaka, M., A. Goto, S. Umino, and E. Fujita (2005), VTFS project: Development of the lava flow simulation code LavaSIM with a model for three-dimensional convection, spreading, and solidification, Geochem. Geophys. Geosyst., 6, Q07008, doi:10.1029/2004GC000869. Hon, K., J. Kauahikaua, R. Denlinger, and K. Mackay (1994), Emplacement and inflation of pahoehoe sheet flows: Observations and measurements of active lava flows on Kilauea Volcano, Hawaii, Geol. Soc. Am. Bull., 106, 351–370, doi:10.1130/0016- 7606(1994)106<0351:EAIOPS>2.3.CO;2. Huppert, H. E. (1982), Flow and instability of a viscous current down a slope, Nature, 300, 427–429, doi:10.1038/300427a0. Iverson, R. M. (1997), The physics of debris flows, Rev. Geophys., 35(3), 245–296, doi:10.1029/97RG00426. 190 Jerolmack, D. J., and D. Mohrig (2007), Conditions for branching in depositional rivers, Geology, 35, 463–466, doi:10.1130/G23308A.1. Kauahikaua, J., D.R. Sherrod, K.V. Cashman, C. Heliker, K. Hon, T.N. Mattox, and J.A. Johnson (2003), Hawaiian lava-flow dynamics during the Pu’u Ō’ō-Kupianaha eruption: A tale of two decades, in The Pu’u ’Ō’ō-Kupaianaha Eruption of Kīlauea Volcano, Hawai’i: The First 20 Years, edited by C. Heliker, D. A. Swanson, and T. J. Takahashi, U.S. Geol. Surv. Prof. Pap., 1676, 63–87. Kerr, R. C., R. W. Griffiths, and K. V. Cashman (2006), Formation of channelized lava flows on an unconfined slope, J. Geophys. Res., 111, B10206, doi:10.1029/2005JB004225. Kerr, R. C., and A. W. Lyman (2007), Importance of surface crust strength during the flow of the 1988–1990 andesite lava of Lonquimay Volcano, Chile, J. Geophys. Res., 112, B03209, doi:10.1029/2006JB004522. Lipman, P. W. and N. G. Banks (1987), Aa flow dynamics, Mauna Loa 1984, in Volcanism in Hawaii, edited by R. W. Decker, T. L. Wright, and P. H. Stauffer, U.S. Geol. Surv. Prof. Pap., 1350, 1527–1567. Lister, J. R. (1992), Viscous flows down an inclined plane from point and line sources, J. Fluid Mech., 242, 631–653, doi:10.1017/S0022112092002520. Lockwood, J. P., N. G. Banks, T. T. English, L. P. Greenland, D. B. Jackson, D. J. Johnson, R. Y. Koyanagi, K. A. McGee, A. T. Okamura, and J. M. Rhodes (1985). The 1984 eruption of Mauna Loa Volcano, Hawaii, Eos Transactions, 66(16), 169–171, doi:10.1029/EO066i016p00169-01. Lockwood, J. P., and F. A. Torgerson (1980), Diversion of lava flows by aerial bombing – lessons from Mauna Loa Volcano, Hawaii, Bull. Volcanol., 43, 727–741, doi:10.1007/BF02600367. Lockwood, J. P., R. I. Tilling, R. T. Holcomb, F. Klein, A. Okamura, and D. Peterson (1999), Magma migration and resupply during the 1974 summit eruptions of Kilauea volcano, Hawaii, U.S. Geol. Surv. Prof. Pap., 1613. Lyman, A. W., R. C. Kerr, and R. W. Griffiths (2005), Effects of internal rheology and surface cooling on the emplacement of lava flows, J. Geophys. Res., 110, B08207, doi:10.1029/2005JB003643. Macdonald, G. A. (1943), The 1942 eruption of Mauna Loa, Hawaii, Am. Jour. Sci., 241(4), 241–256, doi:10.2475/ajs.241.4.241. 191 Makaske, B., D. G. Smith, H. J. A. Berendsen, A. G. de Boer, M. F. van Nielen- Kiezebrink, and T. Locking (2009), Hydraulic and sedimentary processes causing anastomosing morphology of the upper Columbia River, British Columbia, Canada, Geomorphology, 111, 194-205, doi:10.1016/j.geomorph.2009.04.019. Malin, M. C. (1980), Lengths of Hawaiian lava flows, Geology, 8, 306–308, doi:10.1130/0091-7613(1980)8<306:LOHLF>2.0.CO;2. Marra, W. A., M. G. Kleinhans, E. A. and Addink (2013), Network concepts to describe channel importance and change in multichannel systems: test results for the Jamuna River, Bangladesh, Earth Surf. Proc. Land, 39, 766–778, doi:10.1002/esp.3482. Mazzarini, F., M. T. Pareschi, M. Favalli, I. Isola, S. Tarquini, and E. Boschi (2005), Morphology of basaltic lava channels during the Mt. Etna September 2004 eruption from airborne laser altimeter data, Geophys. Res. Lett., 32, L04305, doi:10.1029/2004GL021815. Mazzarini, F., M. T. Pareschi, M. Favalli, I. Isola, S. Tarquini and E. Boschi (2007), Lava flow identification and aging by means of lidar intensity: Mount Etna case, J. Geophys. Res.-Sol. Ea., 112, B2, doi:10.1029/2005JB004166. Moore, H. J. (1982), A geologic evaluation of proposed lava diversion barriers for the NOAA Mauna Loa Observatory Mauna Loa Volcano, Hawaii, U.S. Geol. Surv. Open-File Report 82-314. Moore, H. J. (1987), Preliminary estimates of the rheological properties of the 1984 Mauna Loa lava, in Volcanism in Hawaii, edited by R. W. Decker, T. L. Wright, and P. H. Stauffer, U.S. Geol. Surv. Prof. Pap., 1350, 1569–1588. Orr, T. R. 2010, Lava tube shatter rings and their correlation with lava flux increases at Kīlauea Volcano, Hawai‘i, Bull. Volcanol., 73, 335–346, doi:10.1007/s00445- 010-0414-3. Pelletier, J. D., and S. DeLong (2004), Oscillations in arid alluvial-channel geometry, Geology, 32, 713–716, doi:10.1130/G20512.1. Pinkerton, H., and L. Wilson (1994), Factors controlling the lengths of channel-fed lava flows, Bull. Volcanol., 56, 108–120, doi:10.1007/BF00304106. Pyle, D. M. and J. R. Elliot (2006), Quantitative morphology, recent evolution, and future activity of the Kameni Islands volcano, Santorini, Greece, Geosphere, 2, 253– 268, doi:10.1130/GES00028.1. 192 Richter, D. H., J. P. Eaton, K. J. Murata, W. U. Ault, and H. L. Krivoy (1970), Chronological Narrative of the 1959-60 Eruption of Kilauea Volcano, Hawaii, U.S. Geol. Surv. Prof. Pap., 537-E. Riker, J. M., K. V. Cashman, J. P. Kauahikaua, and C. M. Montierth (2009), The length of channelized lava flows: Insight from the 1859 eruption of Mauna Loa Volcano, Hawai‘i, J. Volcanol. Geotherm. Res., 183, 139–156, doi:10.1016/j.jvolgeores.2009.03.002. Rubinov, M. and O. Sporns (2010), Complex network measures of brain connectivity: Uses and interpretations, NeuroImage, 52, 1059–1069, doi:10.1016/j.neuroimage.2009.10.003. Shreve, R. L. (1966), Statistical law of stream numbers, J. Geol., 74, 17–37. Soule, S.A., K. V. Cashman, and J. P. Kauahikaua (2004) Examining flow emplacement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i, Bull. Volcanol., 66, 1–14, doi:10.1007/s00445-003- 0291-0. Strahler, A. N. (1957), Quantitative analysis of watershed geometry, Am. Geophys. Union Trans., 38, 913–920. Tarquini, S., M. Favalli, F. Mazzarini, I. Isola, and A. Fornaciai (2012), Morphometric analysis of lava flow units: Case study over LIDAR-derived topography at Mount Etna, Italy, J. Volcanol. Geotherm. Res., 235–236, 11–22, doi:10.1016/j.jvolgeores.2012.04.026. Tarquini, S., and M. M. Vitturi (2014), Influence of fluctuating supply on the emplacement dynamics of channelized lava flows, Bull. Volcanol., 76, 801, doi:10.1007/s00445-014-0801-2. Ventura, G. and G. Vilardo (2008), Emplacement mechanism of gravity flows inferred from high resolution Lidar data: The 1944 Somma-Vesuvius lava flow (Italy), Geomorphology, 95, 223–235, doi:10.1016/j.geomorph.2007.06.005. Vicari, A., H. Alexis, C. Del Negro, M. Coltelli, M. Marsella, and C. Proietti (2007), Modeling of the 2001 lava flow at Etna volcano by Cellular Automata approach, Environ. Modell. Softw., 22, 1465–1471, doi:10.1016/j.envsoft.2006.10.005. Walker, G. P. L. (1973), Lengths of lava flows, Phil. Trans. R. Soc. Lond. A., 274, 107– 118, doi:10.1098/rsta.1973.0030. Williams, R. S., and J. G. Moore (1983), Man against volcano: The eruption on Heimaey, Vestmannaeyjar, Iceland, U.S. Geol. Surv. General Interest Publication, 27 p. 193 Wolfe, E. W. (Ed) (1988), The Puu Oo eruption of Kilauea Volcano, Hawaii: Episodes 1 through 20, January 8, 1983, through June 8, 1984, U.S. Geol. Surv. Prof. Pap., 1463. Wolfe, E. W. and J. Morris (1996), Geologic Map of the Island of Hawaii, Miscellaneous Investigation Series, U.S. Geological Survey, Washington, D.C., USA. CHAPTER V Bagdassarov, N., and H. Pinkerton (2004), Transient phenomena in vesicular lava flows based on laboratory experiments with analogue materials, J. Volcanol. Geoth. Res., 132, 115–136, doi:10.1016/S0377-0273(03)00341-X. Barberi, F., F. Brondi, M. L. Carapezza, L. Cavarra, and C. Murgia (2003), Earthen barriers to control lava flows in the 2001 eruption of Mt. Etna, J. Volcanol. Geoth. Res., 123, 231–243, doi:10.1016/S0377-0273(03)00038-6. Barberi, F., and M. L. Carapezza (2004), The control of lava flows at Mt. Etna, in Mt. Etna: Volcano Laboratory, Geophys. Monogr. Ser., vol. 143, edited by A. Bonaccorso et al., pp. 357–369, AGU, Washington, D. C., doi:10.1029/143GM22. Baxter, S. J., H. Power, K. A. Cliffe, and S. Hibberd (2009), Three-dimensional thin film flow over and around an obstacle on an inclined plane, Phys. Fluids, 21, 032102, doi:10.1063/1.3082218. Beckett, F. M., H. M. Mader, J. C. Phillips, A. C. Rust, and F. Witham (2011), An experimental study of low-Reynolds-number exchange flow of two Newtonian fluids in a vertical pipe, J. Fluid Mech., 682, 652–670, doi:10.1017/jfm.2011.264. Castruccio, A., A. C. Rust, and R. S. J. Sparks (2010), Rheology and flow of crystal- bearing lavas: Insights from analogue gravity currents, Earth Planet. Sc. Lett., 297, 471–480, doi:10.1016/j.epsl.2010.06.051. Castruccio, A., A. C. Rust, and R. S. J. Sparks (2013), Evolution of crust- and core- dominated lava flows using scaling analysis, Bull. Volcanol., 75, 681, doi:10.1007/s00445-012-0681-2. Corti, G., A. Zeoli, and M. Bonini (2003), Ice-fow dynamics and meteorite collection in Antarctica, Earth Planet. Sc. Lett., 215, 371–378, doi:10.1016/S0012- 821X(03)00440-0. Corti, G., A. Zeoli, P. Belmaggio, and L. Folco (2008), Physical modeling of the influence of bedrock topography and ablation on ice flow and meteorite concentration in Antarctica, J. Geophys. Res., 113, F01018, doi:10.1029/2006JF000708. 194 Dietterich, H. R., S. A. Soule, K. V. Cashman, and B. H. Mackey (in press), Lava Flows in 3D: Using airborne lidar and pre-eruptive topography to evaluate lava flow surface morphology and thickness in Hawai‘i, in Hawaiian Volcanoes, From Source to Surface, Geophys. Monogr. Ser., AGU, Washington, D. C. Dietterich, H. R., and K. V. Cashman (in press), Channel networks within lava flows: Formation, evolution, and implications for flow behavior, J. Geophys. Res., doi:10.1002/2014JF003103. Favalli, M., M. T. Pareschi, A. Neri, and I. Isola (2005), Forecasting lava flow paths by a stochastic approach, Geophys. Res. Lett., 32, L03305, doi:10.1029/2004GL021718. Guest, J. E., P. D. Spudis, R. Greeley, G. J. Taylor, and S. M. Baloga (1995), Emplacement of xenolith nodules in the Kaupulehu lava flow, Hualalai Volcano, Hawaii, Bull. Volcanol., 57, 179–184, doi:10.1007/BF00265037. Harris, A. J. L., and S. K. Rowland (2001), FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel, Bull. Volcanol., 63, 20–44, doi:10.1007/s004450000120. Hayes, M., S. B. G. O’Brien, and J. H. Lammers (2000), Green’s function for steady flow over a small two-dimensional topography, Phys. Fluids, 12, 2845, doi:10.1063/1.1311970. Hon, K., J. Kauahikaua, R. Denlinger, and K. Mackay (1994), Emplacement and inflation of pahoehoe sheet flows: Observations and measurements of active lava flows on Kilauea Volcano, Hawaii, Geol. Soc. Am. Bull., 106, 351–370, doi:10.1130/0016- 7606(1994)106<0351:EAIOPS>2.3.CO;2. Huppert, H. E. (1982), The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface, J. Fluid Mech., 121, 43–58, doi:10.1017/S0022112082001797. Jeffreys, H. (1925), Flow of water in an inclined channel of rectangular section, Philos. Mag., 49, 793–807, doi:10.1080/14786442508634662. Kauahikaua, J., K. V. Cashman, D. A. Clague, D. Champion, and J. T. Hagstrum (2002), Emplacement of the most recent lava flows on Hualālai Volcano, Hawai‘i, Bull. Volcanol., 64, 229–253, doi:10.1007/s00445-001-0196-8. Kauahikaua, J., D.R. Sherrod, K.V. Cashman, C. Heliker, K. Hon, T.N. Mattox, and J.A. Johnson (2003), Hawaiian lava-flow dynamics during the Pu’u Ō’ō-Kupianaha eruption: A tale of two decades, in The Pu’u ’Ō’ō-Kupaianaha Eruption of Kīlauea Volcano, Hawai’i: The First 20 Years, edited by C. Heliker, D. A. Swanson, and T. J. Takahashi, U.S. Geol. Surv. Prof. Pap., 1676, 63–87. 195 Kerr, R. C., R. W. Griffiths, and K. V. Cashman (2006), Formation of channelized lava flows on an unconfined slope, J. Geophys. Res., 111, B10206, doi:10.1029/2005JB004225. Kerr, R. C., and A. W. Lyman (2007), Importance of surface crust strength during the flow of the 1988–1990 andesite lava of Lonquimay Volcano, Chile, J. Geophys. Res., 112, B03209, doi:10.1029/2006JB004522. Lipman, P. W. and N. G. Banks (1987), Aa flow dynamics, Mauna Loa 1984, in Volcanism in Hawaii, edited by R. W. Decker, T. L. Wright, and P. H. Stauffer, U.S. Geol. Surv. Prof. Pap., 1350, 1527–1567. Lister, J. R. (1992), Viscous flows down an inclined plane from point and line sources, J. Fluid Mech., 242, 631–653, doi:10.1017/S0022112092002520. Llewllin, E. W., H. M. Mader, and S. D. R. Wilson (2002), The rheology of a bubbly liquid, Proc. R. Soc. Lond. A, 458, 987–1016, doi:10.1098/rspa.2001.0924. Moore, H. J. (1982), A geologic evaluation of proposed lava diversion barriers for the NOAA Mauna Loa Observatory Mauna Loa Volcano, Hawaii, U.S. Geol. Surv. Open-File Report 82-314. Moore, H. J., (1987), Preliminary estimates of the rheological properties of the 1984 Mauna Loa lava, in Volcanism in Hawaii, edited by R. W. Decker et al., U.S. Geol. Surv. Prof. Pap., 1350, pp. 1569–1588. Pierson, T. C. (1985), Initiation and flow behavior of the 1980 Pine Creek and Muddy River lahars, Mount St. Helens, Washington, Geol. Soc. Am. Bull., 96, 1056– 1069, doi:10.1130/0016-7606(1985)96<1056:IAFBOT>2.0.CO;2. Richter, D. H., J. P. Eaton, K. J. Murata, W. U. Ault, and H. L. Krivoy (1970), Chronological Narrative of the 1959–60 Eruption of Kilauea Volcano, Hawaii, U.S. Geol. Surv. Prof. Pap., 537-E. Shaw, H. R. (1969), Rheology of basalt in the melting range, J. Petrol.,10, 510–535, doi:10.1093/petrology/10.3.510. Soule, S.A., K. V. Cashman, and J. P. Kauahikaua (2004) Examining flow emplacement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawai'i, Bull. Volcanol., 66, 1–14, doi:10.1007/s00445-003- 0291-0. Takagi, D., and H. E. Huppert (2010), Initial advance of long lava flows in open channels, J. Volcanol. Geoth. Res., 195, 121–126, doi:10.1016/j.jvolgeores.2010.06.011. 196 Williams, R. S., and J. G. Moore (1983), Man against volcano: The eruption on Heimaey, Vestmannaeyjar, Iceland, U.S. Geol. Surv. General Interest Publication, 27 p. Wolfe, E. W. (Ed) (1988), The Puu Oo eruption of Kilauea Volcano, Hawaii: Episodes 1 through 20, January 8, 1983, through June 8, 1984, U.S. Geol. Surv. Prof. Pap., 1463. CHAPTER VI Barberi, F., F. Brondi, M. L. Carapezza, L. Cavarra, and C. Murgia (2003), Earthen barriers to control lava flows in the 2001 eruption of Mt. Etna, J. Volcanol. Geoth. Res., 123, 231–243, doi:10.1016/S0377-0273(03)00038-6. Favalli, M., A. Fornaciai, F. Mazzarini, A. Harris, M. Neri, B. Behncke, M. T. Pareschi, S. Tarquini, and E. Boschi (2010), Evolution of an active lava flow field using a multitemporal LIDAR acquisition, J. Geophys. Res., 115, B11203, doi:10.1029/2010JB007463. Kerr, R. C., R. W. Griffiths, and K. V. Cashman (2006), Formation of channelized lava flows on an unconfined slope, J. Geophys. Res., 111, B10206, doi:10.1029/2005JB004225. Lyman, A. W., R. C. Kerr, and R. W. Griffiths (2005), Effects of internal rheology and surface cooling on the emplacement of lava flows, J. Geophys. Res., 110, B08207, doi:10.1029/2005JB003643. Snavely, N., S. M. Seitz, and R. Szeliski (2007), Modeling the World from Internet Photo Collections, Int. J. Comput. Vision, 80(2), 189–210, doi:10.1007/s11263-007- 0107-3.