Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Identifying and characterizing the spatial variability of supraglacial hydrological features on the western… King, Leonora Adele 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2018_november_king_leonora.pdf [ 111.05MB ]
JSON: 24-1.0372827.json
JSON-LD: 24-1.0372827-ld.json
RDF/XML (Pretty): 24-1.0372827-rdf.xml
RDF/JSON: 24-1.0372827-rdf.json
Turtle: 24-1.0372827-turtle.txt
N-Triples: 24-1.0372827-rdf-ntriples.txt
Original Record: 24-1.0372827-source.json
Full Text

Full Text

Identifying and characterizing the spatial variability of supraglacialhydrological features on the western Greenland Ice SheetbyLeonora Adele KingBSc, Geography, McGill University, 2008MSc, Earth and Environmental Sciences, University of British Columbia, Okanagan, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Geography)The University of British Columbia(Vancouver)October 2018c© Leonora Adele King, 2018The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance the dissertation entitled:Identifying and characterizing the spatial variability of supraglacial hydrologicalfeatures on the western Greenland Ice Sheetsubmitted by Leonora King in partial fulfillment of the requirements for the degree ofDoctor of Philosophy in GeographyExamining Committee:Marwan Hassan, Geography(Supervisor)Gwenn Flowers, SFU Earth Sciences(Supervisory Committee Member)Nicholas Coops, Forestry(Supervisory Committee Member)(University Examiner)(Departmental Examiner)(External Examiner)iiAbstractGlacier and ice sheet response to climate change is modulated in part by the temporal and spa-tial delivery of surface water to the subglacial system which can enhance ice sliding. Surface-to-bed connections are formed by near-vertical pathways through the ice called moulins. Themechanisms through which moulins form and the hydrological implications of their distribu-tion are inadequately constrained.This thesis focuses on the southwest Greenland Ice Sheet, where surface melt and relatedfeatures are abundant and where meltwater drainage through moulins has been empiricallylinked to ice speedup. In the first part, I employ flow routing over a high resolution digitalelevation model to delineate supraglacial channels and catchments. I compare my results toan independent dataset and demonstrate that flow routing is an effective tool for delineatingchannels in high detail. Whereas multispectral methods provide instantaneous discharge esti-mates and flow routing derived catchments can be used to build hourly hydrographs, the twomethods are temporally incongruous and produce discharge estimates in small catchmentsthat differ by several orders of magnitude.In the second part of the thesis, I contribute empirical insight into where moulins formand the hydrological implications of the resultant spatial organization of catchments. Usingremotely sensed data on supraglacial lake distribution, bed elevation, ice thickness, surfacevelocity and strain rate, I show that moulins form predominantly in low velocity, high strainrate ice with adverse bed slopes. Moulins formed near moulin-drained lakes are typicallylocated down-ice of the lake.The spatial distribution of moulins affects the physical characteristics of supraglacial catch-ments, including elevation, size, elongation, drainage network configuration and slope. Thisin turn impacts the temporal delivery of meltwater to the subglacial network. I use syntheticunit hydrograph theory based on catchment morphometrics to show that catchment diurnalhydrograph lag times in the study area vary by up to 11 hours. The hydrological implicationsof subglacial meltwater supply vary according to the configuration of the subglacial catch-ments, which vary temporally based on subglacial water pressure. Considering this additionallevel of complexity is a necessary next step in constraining the relationship between melt andice flow dynamics.iiiLay SummaryIt is difficult to predict how the Greenland Ice Sheet will respond to a changing climate becausewe do not understand how different processes interact within the ice sheet. For example, melt-water that reaches the base of the ice sheet can cause the ice to accelerate under certain condi-tions. This dynamic depends on several interacting processes that we have yet to fully explore.This thesis explores some of these processes by identifying 1) surface features associated withmelt; 2) controls on the spatial distribution of meltwater drainage holes (called moulins); and3) how the distribution of moulins affects the distribution of surface melt drainage in space. Ifind that moulin distribution reflects ice velocity characteristics, surface lake distribution, andthe gradient of the bed with respect to the gradient of the ice surface. This means that moulinsare unequally distributed over the ice sheet, as is water delivery to the base of the ice sheet.ivPrefaceThis thesis is original work completed by Leonora King. Guidance was given by the supervi-sory committee.This thesis includes two manuscripts, and two complementary Chapters that will be submit-ted for publication as two or more manuscripts. The published manuscripts are presented inChapter 2 and Chapter 3. Chapter 5 and Chapter 4 are the complementary chapters.A version of the work in Chapter 2 is published in Remote Sensing (King et al., 2016). Theco-authors are Marwan Hassan, Kang Yang and Gwenn Flowers. I am responsible for devel-oping the project design, performing the remote sensing analysis, and writing the text. Theco-authors provided editorial review of the manuscript prior to publication.King, L., M. A. Hassan, K. Yang and G. Flowers (2016), Flow Routing for Delineating SupraglacialMeltwater Channels, Remote Sensing, 8, 988.A version of the work in Chapter 3 is accepted for publication as a letter in Journal of Glaciology.I am the sole author.King, L. (2018), Comparing two methods of remotely estimating moulin discharge on theGreenland Ice Sheet, Journal of Glaciology.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction and motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 The link between surface melt and ice flow . . . . . . . . . . . . . . . . . . . . . 51.2 A primer on moulins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Motivation and study design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Flow routing for delineating supraglacial meltwater channels . . . . . . . . . . . . 112.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.3.1 Datasets and study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4.1 Comparison between flow routing and multispectral datasets . . . . . . 222.4.2 Comparison with Manually Digitized Datasets . . . . . . . . . . . . . . . 272.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.5.1 Moulins . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.5.2 Flow Routing in High-Elevation Catchments . . . . . . . . . . . . . . . . 342.5.3 Flow Routing in Low-Elevation Catchments . . . . . . . . . . . . . . . . 35vi2.5.4 Limitations and Considerations . . . . . . . . . . . . . . . . . . . . . . . . 362.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373 Replicating moulin discharge estimates reveals scale dependent error . . . . . . . 393.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444 The spatial variability and characteristics of moulins in southwest Greenland . . . 484.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.3 Study Site and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.3.1 Study Area and datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.3.2 Moulin classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.3.3 Lake delineation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.4.1 Overview of spatial moulin characteristics . . . . . . . . . . . . . . . . . 554.4.2 Moulin Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.4.3 Influence of lakes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.4.4 Moulin location relative to crevasse fields . . . . . . . . . . . . . . . . . . 674.4.5 Longitudinal trends in moulin distribution . . . . . . . . . . . . . . . . . 684.4.6 Comparison with modeled predictions . . . . . . . . . . . . . . . . . . . 724.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.5.1 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.5.2 Moulin Types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.5.3 Moulin Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.5.4 Influence of lakes on moulin distribution . . . . . . . . . . . . . . . . . . 774.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 775 Spatial variability of supraglacial catchment morphometrics and hydrology on thesouthwestern Greenland Ice Sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.2.1 Background: Predicting hydrograph characteristics in ungauged catch-ments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.3 Study Site and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 825.3.1 Study Area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 825.3.2 Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83vii5.3.3 Surface catchment and channel network delineation . . . . . . . . . . . . 845.3.4 Deriving catchment characteristics . . . . . . . . . . . . . . . . . . . . . . 865.3.5 Deriving catchment hydrograph characteristics . . . . . . . . . . . . . . 875.3.6 Subglacial catchment delineation . . . . . . . . . . . . . . . . . . . . . . . 895.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 905.4.1 Spatial variability of supraglacial catchment characteristics . . . . . . . . 905.4.2 Partitioning of runoff . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.4.3 Connections to subglacial catchments . . . . . . . . . . . . . . . . . . . . 1015.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1095.5.1 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1105.5.2 Catchment and channel network characteristics and variability . . . . . 1105.5.3 Partitioning of supraglacial discharge . . . . . . . . . . . . . . . . . . . . 1125.5.4 Implications of surface catchment hydrology for subglacial hydrology . 1135.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1146 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1156.1 Summary of contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166.2 Limitations, assumptions, and future work . . . . . . . . . . . . . . . . . . . . . 1196.3 Final thoughts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123A Study area characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137B Moulin Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139C Strain rate scale considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140C.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140C.2 Mean Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140C.3 Median Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143C.4 Gaussian Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146C.5 RMSE comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149D Hydraulic variables and catchment properties for derivation of the GIUH. . . . . . 150E An illustrated thesis timeline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152viiiList of TablesTable 2.1 Summary of datasets used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14Table 2.2 Summary of study catchments. . . . . . . . . . . . . . . . . . . . . . . . . . . 15Table 2.3 Moulin comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Table 5.1 Summary of catchment and network morphometrics. . . . . . . . . . . . . . 87ixList of FiguresFigure 1.1 The Greenland Ice Sheet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2Figure 1.2 NSF-funded research on Greenland surface hydrology. . . . . . . . . . . . . 3Figure 1.3 Examples of supraglacial hydrological features. . . . . . . . . . . . . . . . . 4Figure 1.4 Media representations of supraglacial channels. . . . . . . . . . . . . . . . . 5Figure 1.5 Outline of the thesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10Figure 2.1 The study area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16Figure 2.2 Workflow for channel delineation. . . . . . . . . . . . . . . . . . . . . . . . . 19Figure 2.3 Ratios of drainage density as a function of contributing area. . . . . . . . . . 24Figure 2.4 Visual comparison of channel networks. . . . . . . . . . . . . . . . . . . . . . 25Figure 2.5 The performance of flow routing relative to multispectral methods acrossinitiation thresholds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.6 Match rates between datasets. . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 2.7 A comparison of digitized, flow routed and multispectrally derived networks. 29Figure 2.8 Match percentages. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 2.9 Comparison of flow-routed and manually digitized channel networks in thelow-elevation catchment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 3.1 The study area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 3.2 Dataset comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Figure 3.3 Test catchment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 4.1 The study area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 4.2 Examples of lake drainage identification. . . . . . . . . . . . . . . . . . . . . 54Figure 4.3 Elevation distribution of moulins. . . . . . . . . . . . . . . . . . . . . . . . . 55Figure 4.4 Slope and annual melt data with respect to elevation. . . . . . . . . . . . . . 56Figure 4.5 Moulin density. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 4.6 Moulin and study area characteristics. . . . . . . . . . . . . . . . . . . . . . . 58Figure 4.7 Moulin distribution relative to ice velocity zone. . . . . . . . . . . . . . . . . 59Figure 4.8 Moulin density relative to velocity and strain bins. . . . . . . . . . . . . . . . 60Figure 4.9 Characteristics of moulin types. . . . . . . . . . . . . . . . . . . . . . . . . . . 61xFigure 4.10 Spatial distribution of moulin types relative to velocity. . . . . . . . . . . . . 62Figure 4.11 Moulin type relative to velocity regime. . . . . . . . . . . . . . . . . . . . . . 63Figure 4.12 Distance to nearest lake. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64Figure 4.13 Moulin position relative to lakes. . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure 4.14 PCA of lake drainage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure 4.15 Lake characteristics boxplots. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 4.16 Moulin position relative to crevasse fields. . . . . . . . . . . . . . . . . . . . 68Figure 4.17 Flow lines overview. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 4.18 Flow line profiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 4.19 Total annual moulin discharge relative to strain rate. . . . . . . . . . . . . . 73Figure 5.1 The study Area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Figure 5.2 Close up image of supraglacial channels. . . . . . . . . . . . . . . . . . . . . 85Figure 5.3 Maps of catchment area, lake area, and mean melt. . . . . . . . . . . . . . . 91Figure 5.4 Morphometrics of the study catchments. . . . . . . . . . . . . . . . . . . . . 92Figure 5.5 Relationship between Re and ice velocity. . . . . . . . . . . . . . . . . . . . . 93Figure 5.6 Morphometrics of channel networks. . . . . . . . . . . . . . . . . . . . . . . 94Figure 5.7 Components of the bifurcation, length, and area ratios relative to stream order. 95Figure 5.8 PCA analysis of stream ratios relative to ice characteristics . . . . . . . . . . 96Figure 5.9 Annual discharge and time to peak in study catchments. . . . . . . . . . . . 97Figure 5.10 PCA of catchment morphometrics. . . . . . . . . . . . . . . . . . . . . . . . . 98Figure 5.11 Partitioning of surface drainage. . . . . . . . . . . . . . . . . . . . . . . . . . 99Figure 5.12 The frequency of total annual catchment discharges distributed by catch-ment type. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100Figure 5.13 Partitioned annual hydrograph. . . . . . . . . . . . . . . . . . . . . . . . . . 101Figure 5.14 Subglacial catchments with different k values . . . . . . . . . . . . . . . . . . 103Figure 5.15 Changes in catchment areas for different k values. . . . . . . . . . . . . . . . 104Figure 5.16 Average catchment tp with different subglacial catchment configurations. . 105Figure 5.17 Total catchment Qyr with different subglacial catchment configurations. . . 107Figure 5.18 Catchment Geomorphic Instantaneous Unit Hydrographs. . . . . . . . . . . 108Figure 5.19 Changing catchment discharge and tp with increasing k. . . . . . . . . . . . 109Figure 6.1 Temperature and sea level rise since 1993 . . . . . . . . . . . . . . . . . . . . 115Figure A.1 Surface slope and ice thickness in the study area. . . . . . . . . . . . . . . . . 137Figure A.2 Ice surface velocity and principal strain rate in the study area . . . . . . . . 138Figure B.1 Examples of moulin identification. . . . . . . . . . . . . . . . . . . . . . . . . 139Figure C.1 Mean filtering across various kernel sizes. . . . . . . . . . . . . . . . . . . . . 143Figure C.2 Median filtering across various kernel sizes. . . . . . . . . . . . . . . . . . . 146xiFigure C.3 Gaussian Blurring across various kernel sizes. . . . . . . . . . . . . . . . . . 149Figure C.4 RMSE of different filtering methods across kernel sizes. . . . . . . . . . . . . 149Figure D.1 Hydraulic variables and catchment properties for derivation of the GIUH. . 151Figure E.1 A graphical representation of the PhD journey. . . . . . . . . . . . . . . . . . 152xiiAcknowledgmentsAs always, I would like to firstly thank my parents, Karen and Geoff, for their endless andunconditional love, their unwavering support and patience, their peculiar confidence in me,and the incredible privileges and joys they have given me in life. I’ll never be able to thankyou enough for all that you have given me.Similarly, I extend my enormous gratitude to Dave, my partner, colleague, co-author, andcomrade in academia as well as in life. Having your support, love and companionship alongthe way has made the good times better and kept the bad times in perspective.Thank you Marwan for encouraging me to pursue a PhD with you and for your kind andpatient support through my various academic meanderings and avulsions in interests. Yourencouragement and understanding have been indispensable along the way, and your pep talkshave always been well timed. Gwenn, thank you for adopting me into the world of glaciologyand giving me the guidance I needed to bring this thesis to a conclusion. It has been a privilegeto work with you. You have been a brilliant, kind and inspiring role model for how sciencecan be done with both rigor and generosity.I have an incredible network of fam-riends who are deserving of my gratitude for stickingby me through these challenging six years. Carly, Haley, Salome, Anna, Julie, Dave W, Charles,Kiely, Alexandria: thanks to all of you for your beautiful friendship - the good times we sharetogether have carried me through. To Marc in particular: thank you for enticing me into andmentoring me through the overlap between human and physical geography. You have givenme a love for Geography that has enriched this experience.Thank you to all my labmates and colleagues; Tobias, Maria, Lucy, Katie, Dave, Erik. Inparticular, thank you Shawn for your support, encouragement, collegiality and mentorshipthroughout. Thank you to the UBC Geography faculty and staff who helped and encouragedme along the way, including Vincent, Suzanne, Danny, Sandy, Sumi, Brett, Siobhan, and Eric.And thank you to Nicholas for accommodating my tight finishing timeline!I extend specific thanks to Kang Yang and Laurence Smith for their help with Chapters 2and 3, as well as to the anonymous reviewers of those chapters.xiiiDedicationTo my parents, Karen Patrice Himsl King and Geoffrey John King. I hope that this thesis andeverything else I do in this world honours you and the love you give me.xivChapter 1Introduction and motivationThe Greenland Ice Sheet (GrIS) (Sermersuaq in the Kalaallisut language of the Kalaallit Nunaatpeople of Greenland) is the second largest body of ice in the world. It is almost 2,400 kmnorth to south, and 1,100 km wide, covering a total of 1,710,000 km2, approximately a fifththe surface area of Canada (Figure 1.1). At its deepest point, the ice sheet is 3.43 km thick,and in total it contains 2,942,393 km3 of ice (Morlighem et al., 2017) (Figure 1.1), equivalent to apotential 7.4 m of global sea level rise and 10% of the Earth’s freshwater resources (Benn andEvans, 2010).There is varied scientific interest in the GrIS, ranging from ice dynamics modeling (see, forexample, Kirchner et al., 2011) to paleoclimatic studies (e.g. Appenzeller et al., 1998; Rasmussenet al., 2007; Thomas et al., 2007), to the oceanographic influence of freshwater runoff (e.g. Driess-chaert et al., 2007; Fairbanks, 1989). Increasingly, research on the GrIS revolves around a climatechange framing, particularly (but not exclusively) with regards to its potential contribution toglobal sea level rise and impact on ocean circulation. Although not without spatial and tem-poral complexity (e.g. Chylek et al., 2006), there has been a generalized warming trend overGreenland in the last three decades (Hall et al., 2013; Hanna et al., 2008). This warming has beenwidely associated with increased meltwater production (Fettweis et al., 2007; Frauenfeld et al.,2011; van den Broeke et al., 2016), culminating in 2012 when almost the entire surface of the icesheet underwent melt (Nghiem et al., 2012).Between 1991 and 2015, the GrIS has contributed an average of 0.47±0.23 mm/yr of sealevel equivalent mass loss to global sea level rise, with a peak contribution of 1.2 mm duringthe 2012 melt year (van den Broeke et al., 2016). Assuming a mid-range emissions scenario, aminimum of almost 5 cm of global sea level rise is already committed from the GrIS by the year2100, a lower estimate that will almost certainly be significantly exceeded (Price et al., 2011).There is clearly interest in projecting these contributions into the future (e.g. Church et al., 2013;Nerem et al., 2018). If increased surface temperature equated linearly to increased mass lossthrough surface melt, the relationship would be easily characterized and predictions of futuremelt would be readily made. However, like all glacial systems, the mass loss in the GrIS has acomplex and non linear relationship to climate. For example, dynamic glaciological processes1Sources: Esri, GEBCO, NOAA, National Geographic, Garmin, HERE,, and other contributors, Esri, Garmin, GEBCO, NOAANGDC, and other contributors, Esri, HERE, DeLorme, MapmyIndia, ©OpenStreetMap contributors, and the GIS user community30°0'0"W40°0'0"W50°0'0"W75°0'0"N 75°0'0"N70°0'0"N 70°0'0"N65°0'0"N65°0'0"N60°0'0"N60°0'0"NIce Thickness (m)343500 350 700175 KmFigure 1.1: The Greenland Ice Sheet and peripheral ice masses. Color gradation refers toice thickness, where data are from Morlighem et al. (2017)’s BedMachineV3 bed andice thickness model.can redistribute ice with respect to elevation in response to climate forcings (e.g. Pritchard et al.,2009). Ice redistributed to lower elevations through glacier acceleration and thinning meltsmore rapidly due to the elevational temperature gradient (Oerlemans and Hoogendoorn, 1989).Additionally, where outlet glaciers are marine terminating, increased ice discharge across thegrounding line (where glacier tongues begin to float as ice shelves) can lead to higher massloss through calving processes (e.g. Nick et al., 2009; Price et al., 2011; van den Broeke et al., 2016).On the GrIS, increased melt has been non-linearly accompanied by acceleration and dy-namic thinning on the margins and outlet glaciers of the icesheet (Krabill et al., 2004; Nick et al.,2009; Rignot and Kanagaratnam, 2006; ?, e.g.). This dynamic thinning is a poorly understoodbut crucial part of the GrIS response to climate change; ice discharge across the calving frontnow accounts for 40% of GrIS mass loss (van den Broeke et al., 2016).Nonlinear feedbacks between temperature increase and mass loss are enormously concern-ing to the scientific community. As these feedbacks have been uncovered, hyperbolic claimsand a lack of nuanced understanding of ice sheet processes have created a morass of scientific2drama, as well as the complete omission of these potential feedback effects from the 2007 IPCCFourth Assessment Report on sea level rise (O’Reilly, 2015). More refined understanding hasnow been developed, and several mechanisms have been invoked to explain the unexpectedacceleration observed in the GrIS, prominent among them lubrication of the base of the icesheet by percolating surface meltwater.     <HDU16)*UDQWV$ZDUG9DOXH0LOOLRQFigure 1.2: The grey bars indicate the number of National Science Foundation (NSF)-funded awards through time that refer to ’Greenland’ and either ’lakes’,’supraglacial’ or ’hydrology’ in their titles. The red dots indicate the value of theawards (per awards), as per the secondary y-axis. Data acquired from potential importance of surface meltwater to GrIS stability and dynamics has fosteredover a decade and a half of heightened academic and public interest in ice sheet surface hy-drology. Research interest has exploded in the last decade, with over $23 million USD awardedto National Science Foundation (NSF) researchers alone to study Greenland surface hydrology(Figure 1.2). Additionally, whereas the Antarctic ice sheet was previously thought to be toocold to be impacted by surface melt, that view is now changing (e.g. Dow et al., 2018; Langleyet al., 2016). Research on the GrIS may now be prophetic; lessons learned from the dynamicrelationships observed there might be exported to the changing surface melt conditions nowbeing observed in Antarctica.As scientists have publicized the results of their work, surface melt, epitomized by impres-sive supraglacial channels and sparkling emerald lakes (Figure 1.3), has become a visual andtangible image of glacier loss that is often injudiciously used in the media to evoke a visceralpublic reaction to climate change. Rolling Stone, for example, titled a stunningly edited arti-cle replete with images of supraglacial channels and lakes ’Greenland Melting’, as though thepublic is to understand that these systems are symptomatic of a never-before seen melt, ratherthan a seasonally repeated process (Figure 1.4).3Figure 1.3: These examples illustrate typical Greenland supraglacial features. In A) waterponds in a supraglacial lake located in a shallow depression in relatively flat ice. Thelake drains into a moulin evident to the left of the frame, and another lake is seen inthe background. In B) and C) supraglacial meltwater channels are incised into theice, fed by densely spaced, dendritic tributary channels. All images courtesy of G.Flowers.In fact, supraglacial hydrological features are widespread and globally ubiquitous wher-ever ablation zones are to be found, recurring seasonally to the great inconvenience of fieldresearchers. The extent to which these features represent or in fact contribute to ice sheet re-sponse to climate change is still the subject of investigation - an investigation for which wehave only recently acquired the methodological tools and conceptual approaches with whichto piece together the dynamics of this complex system.This thesis is situated in and cognizant of this rapidly proliferating interest in ice sheetsupraglacial hydrology. Despite a limited understanding of these systems and a paucity ofempirical data, assumption-laden numerical models have been used to make predictions onthe relative importance or unimportance of surface hydrology to the future of the GrIS (e.g.Parizek and Alley, 2004; Poinar et al., 2015). In this context, the motivation of this current workis to check methodological and empirical assumptions about these systems and to contributeempirical and descriptive work of supraglacial hydrology that has thus far been overlooked.As per the title of the thesis, the objectives are twofold. The first two chapters are method-ological, and speak to the ’identification’ of supraglacial hydrological features. With the rapid41.1. The link between surface melt and ice flowFigure 1.4: The New Yorker (top) and Rolling Stone (lower left) both ran climate-changethemed specials on supraglacial hydrology (Goodell, 2016; Kolbert, 2016). TheGuardian (lower right) uses an image of a supraglacial channel in an article aboutGreenland Ice Sheet melt, although the article does not pertain to these systems atall (Abraham, 2015).and recent release of new supraglacial catchment and river datasets (Smith et al., 2015, 2017;Yang and Smith, 2016), Chapters 2 and 3 independently replicate prior findings to verify method-ological limitations. Chapters 4 and 5 are empirical, and contribute to the ’quantification’ di-mension of the title. These chapters describe the spatial distribution of moulins and catchmenthydrology, respectively, in relation to ice conditions.To situate the reader, this introduction provides a brief but relevant overview of the historyand current understanding of the importance of surface hydrology to Greenland Ice Sheetdynamics, and then presents the specific outline and motivation of the thesis chapters.1.1 The link between surface melt and ice flowIt has been known for some time that subglacial water pressure can increase the basal slidingvelocity of a temperate glacier (Iken, 1981). However, it was only within the last 20 years thatfield observations of ice motion and ablation rate empirically linked meltwater production toseasonal and spatial variability in GrIS velocity (Zwally et al., 2002). Thereafter, a series ofobservational campaigns, primarily in southwestern GrIS, have been conducted to refine ourunderstanding of how, when, and where surface meltwater contributes to ice sheet dynamics.Consistent observations reveal that throughout a given melt season, there is a strong correla-tion between melt rate and ice velocity (Andrews et al., 2014; Palmer et al., 2011; Sundal et al.,51.1. The link between surface melt and ice flow2011; van de Wal et al., 2008), evidence that seasonal meltwater drives seasonal basal sliding.Additional evidence for the linkage between surface melt and ice velocity is inferred from thestrong correlation between rapid lake drainage and rapid ice speed-up events (Das et al., 2008;Joughin et al., 2008, 2013).There is a distinct, non-linear temporality to this relationship, such that ice velocity in-creases markedly to a maximum during melt periods, decreases to a minimum immediatelyafter peak summer melt, and gradually increases to winter velocities after melting has ceased(Andrews et al., 2014; Joughin et al., 2008; Sole et al., 2013; Zwally et al., 2002). Furthermore,although seasonal and event-specific relationships between melt and velocity are consistent,other observations suggest that at coarser, inter-annual scales, the relationship is highly com-plex and non-linear. For example, there has been a generalized decrease in ice velocities in theSW GrIS over the last several decades, despite increasing melt rates over the same time period(Tedstone et al., 2015; van de Wal et al., 2008). Similarly, Sundal et al. (2011) found that althoughpeak rates of ice acceleration were positively correlated with degree of melting (higher melt insummer led to higher velocities relative to the winter), the mean summer flow rates were not;years with lower melt were in fact associated with higher ice velocities.The proposed mechanism to reconcile these disparate observations is subglacial conduitsthat evolve seasonally to accommodate the increased summer melt. Schoof (2010) proposedthat discharge in the subglacial hydrological network below a critical value is insufficientto maintain large-scale channelization. This situation promotes distributed, high pressuredrainage through numerous small conduits and cavities; water supplied to this sub-glacial net-work will rapidly increase effective pressure and promote bed separation and sliding. Abovea critical discharge value, sub-glacial drainage evolves towards macroporous flow throughlarge channels that have evolved to convey the supplied discharge and reduce effective pres-sure, reducing ice sliding velocities. This model of the co-evolution of basal water pressureand sub-glacial conduits is consistent with the seasonal observations of meltwater supply andice velocity, and has been suggested as a causal mechanism (Andrews et al., 2014; Joughin et al.,2013; Moon et al., 2014; Sundal et al., 2011).Recent observations offer new insights into other factors affecting the relationship betweenmeltwater drainage and ice sheet velocity. Through comparisons of hydraulic head in moulinsand boreholes with ice velocity measurements, Andrews et al. (2014) observed strong diurnalcorrelations between moulin hydraulic head and ice velocity in the Sermeq Avannarleq region.However, they found that this correlation displays seasonal hysteresis; diurnal ice velocityfluctuations decreased during the melt season although moulin hydraulic head fluctuationsdid not, suggesting non-hydrologic drivers of ice motion. Similarly, several fast-flowing out-let glaciers, including the rapidly accelerating Jakobshavn Isbrae (Sermeq Kujalleq in Green-landic) glacier appear to be insensitive to meltwater supply, and in fact speed up towards themelt season’s end, after seasonal speedups on other glaciers have ended (Joughin et al., 2008;Palmer et al., 2011; Tedstone et al., 2015). Non-hydrologic factors such as bed topography and61.2. A primer on moulinstransfer of longitudinal stress have been invoked to explain some of this spatial variability(Andrews et al., 2014; Joughin et al., 2004).Thus, there are both theoretical and empirical reasons to assume a relationship betweensurface meltwater hydrology and ice dynamics. However, the relationship is both tempo-rally and spatially variable. Parts of this variability may be attributed to the transfer of stressthrough the ice in response to differences in bed topography (Andrews et al., 2014; Joughinet al., 2004). However, where subglacial hydrology has strong feedbacks on ice basal veloc-ities, the magnitude, timing and location of surface meltwater delivery to moulins may alsocritically influence sub-glacial channel evolution and thus the spatial variability of ice veloc-ity. To date, prior studies have tended to consider meltwater delivery in terms of regionalmelt rate, despite growing evidence that supraglacial catchment characteristics are spatiallyvariable (Banwell et al., 2012; Gleason et al., 2016; Lampkin and VanderBerg, 2014; Legleiter et al.,2014; Smith et al., 2015; Yang and Smith, 2013, 2016; Yang et al., 2016a) and consequently producemoulin discharges in SW GrIS that varied from 0.36 up to 17.72 m3s-1 during peak melt in 2012(Smith et al., 2015), with hydrograph peak lags that range from under 2 hours to over 8 (Smithet al., 2017). Still other evidence suggests that the distribution and style of drainage affects therapidity and quantity of water delivery in the ice, such that crevasse drainage dampens thediurnal cycle of meltwater inputs regionally (McGrath et al., 2011). This slower, steadier andmore distributed drainage contrasts with rapid and concentrated meltwater injection throughmoulins. It is possible that, given the strong relationship between rapid, high volume lakedrainage events and ice speedup, more distributed, low volume and sustained meltwater de-livery generates a different ice flow dynamics response, even for the same volume of injectedwater.1.2 A primer on moulinsThe term ’moulin’ comes from the french word for mill. These features were of early interestto geologists for the role they play in creating bedrock erosion features (e.g. Alexander, 1932).Their importance to glacier hydrology has made moulins features of extensive study in theglaciological community since the 1970s. Much of this focus has been in alpine glacier en-vironments with emphasis on the relationship between surface meltwater drainage and icevelocity (e.g. Holmlund and Hooke, 1983; Iken, 1972) as well as their internal structure and for-mative mechanisms (e.g. Dewart, 1965; Holmlund, 1988; Holmlund and Hooke, 1983; Phillips et al.,2013; Schroeder, 1998).Research in alpine environments has been pivotal for developing the relationship betweenmeltwater at the ice-bed interface and glacier acceleration (e.g. Iken, 1981), and has revealedfour different possible moulin formation mechanisms (Gulley et al., 2009a): 1) cut-and closurein which supraglacial streams incise downwards and are cut off from the surface by creepclosure of the channel (Gulley et al., 2009b), 2) exploitation of permeable structures such asmoraine debris or fractures (e.g. Fyffe et al., 2015), 3) hydrofracture in which water pressure71.3. Motivation and study designin a fracture exceeds the lithostatic pressure of the ice which acts to close the fracture (van derVeen, 2007), and 4) a combination of any of the above processes.Although all of these processes are presumably at work on the GrIS, moulin processes onthe ice sheet differ from those in alpine environments for several reasons. Firstly, ice thick-ness in the ablation zone of the GrIS frequently exceeds 1 km (Morlighem et al., 2017), whereasalpine glacier ice is typically well below 500 m in thickness (Gärtner-roer et al., 2014). Secondly,higher volumes of supraglacial melt are available for moulin formation on the ice sheet than onalpine glaciers due to large catchment sizes and the presence of extensive supraglacial lakes.Thirdly, because the ice sheet overrides underlying topography, its ice dynamics are funda-mentally different from alpine glacier and its flow direction is driven primarily by gradientsin ice thickness rather than bedrock topography(Cuffey and Paterson, 2010).For these reasons, moulin formation on the ice sheet cannot be understood by simple ex-trapolation of alpine glacier research. The complexity of conditions on the ice sheet makessimple predictions of moulin distribution difficult (e.g. Colgan and Steffen, 2009; Phillips et al.,2011). Recent research emphasizes that stochastic lake drainage events likely play an im-portant part in ice sheet moulin formation due to tensile stress propagation following lakedrainage (Christoffersen et al., 2018; Hoffman et al., 2018; Stevens et al., 2015). However, signif-icantly more empirical work is needed into the variety of conditions under which ice sheetmoulins can form.1.3 Motivation and study designThere is a tension in ice sheet research between the urgency and desire to advance a nuancedunderstanding of a complex and dynamic system, and the reality of a theoretical grasp that islimited by a dearth of empirical data and observation. This tension plays out in contemporaryice sheet hydrology research as well. Remote sensing researchers focus on quantifying thevariability of supraglacial features using multispectral imagery on scales of 2 to 15 m (Lampkinand VanderBerg, 2014; Smith et al., 2015, 2017; Yang and Smith, 2016). Models of supraglacialdrainage assume predictable patterns of moulin formation and route flow to those locationsover digital elevation models 30 m or larger in grid size (Banwell et al., 2016, 2013; Clason et al.,2012; Koziol et al., 2017). Field workers generally interest themselves with how hydrologicalprocesses play out below the ice (Andrews et al., 2014; Catania et al., 2008; Meierbachtol et al.,2013).The result of this ontological and methodological siloing is that advances in remote sens-ing of these systems are rapid, and methods are quickly becoming obsolete (Smith et al., 2015,2017). However, the scale and practical considerations of these methods are such that the mod-eling community does not integrate them into their models. Meanwhile, the modeling com-munity advances hypotheses about the projected role of surface hydrology in future ice sheetdynamics without reconciling their numerical assumptions to empirical observations that sug-gest otherwise. For example, whereas Poinar et al. (2015) predict that limited crevassing at high81.3. Motivation and study designelevations will inhibit surface-to-bed connections in the context of an expanding ablation zoneand Koziol et al. (2017) assume that lake-drainage is the dominant mode of moulin formationat high elevations, consistent empirical evidence suggests that moulins without obvious for-mation mechanisms are widespread at high elevations (Lampkin and VanderBerg, 2014; Morrisset al., 2013; Smith et al., 2015; Yang and Smith, 2016).I argue that too little empirical work has been done to test and validate the governing as-sumptions that underpin our current understanding of how supraglacial hydrological systemsactually work, and to draw together the disparate approaches to understanding their dynam-ics. Further, I perceive a gap in methods between those used for identifying supraglacial fea-tures in the modeling domain (i.e. flow routing over digital elevation models) and the currentdominant approaches used in multispectral delineation of surface features (Smith et al., 2015,2017; Yang and Smith, 2016).To pull together these disparate threads and relate numerical interests in process to empir-ical focus on variability, this thesis has four parts, laid out in Figure 1.5. In the next chapter(Chapter 2), I capitalize on the release of a new, high resolution (2 m) Digital Elevation Model(DEM) called ArcticDEM (Noh and Howat, 2015) to test whether the flow routing approachtraditionally employed by models of supraglacial meltwater transport can be substituted forthe high resolution channel delineation that is currently being done using multispectral meth-ods (Smith et al., 2015, 2017). Although multispectral methods are promising and practical inmany ways, they provide only a snapshot in time of the system and, even with high resolu-tion imagery, underestimate channel network extent. If flow routing can be successfully usedfor channel delineation work, it offers the advantage of being readily integrated into surfacehydrology models which already rely on flow routing of surface water.Another advantage of flow routing to delineate channels is that it is easily extended tocatchment delineation. It is one of the only methods currently available on the GrIS for catch-ment delineation, and is particularly promising in the typically low slopes of the GrIS giventhe high spatial resolution of ArcticDEM. In Chapter 3, I test whether catchments delineatedin this way produce moulin discharge values similar to those estimated by Smith et al. (2017),the only existing dataset of moulin discharges. I explore the ways in which this independentreplication of estimating moulin discharge provides insight into the limitations of the twomethods.The fourth chapter is the first of two results chapters. Previous attempts to identify moulinlocation using coarse resolution ice characteristic data have had limited success (Clason et al.,2012; Colgan and Steffen, 2009; Phillips et al., 2011). Without practical insight into the complexformative mechanisms of different moulin types, models therefore rely on overly simplifiedassumptions of moulin formation (Banwell et al., 2013; Clason et al., 2015; Koziol et al., 2017).In this chapter, I use flow routing to derive catchment hydrological characteristics, and I sub-sequently examine the interacting effects of hydrology and ice characteristics on the spatialdistribution of moulins. I adopt an explicitly spatial approach and look for trends in how ice91.3. Motivation and study designconditions and catchment hydrology vary with respect to one another, to moulin type, and tothe spatial distribution of each.The spatial arrangement of moulins dictates, to a large extent, the spatial variability ofsupraglacial catchment characteristics and hydrology. The fifth and final empirical chapter ex-plores the ways in which supraglacial catchment morphometrics vary in space, and expandswork by Smith et al. (2017) to infer catchment hydrograph parameters from catchment morpho-metrics. Assuming that diurnal supraglacial catchment hydrographs can significantly influ-ence the diurnal and seasonal evolution of the subglacial hydrological network, I consider thehydrological implications of this supraglacial catchment hydrological variability with respectto different configurations of subglacial catchments.Figure 1.5: A graphical overview of the thesis outline. The first section, which deals withmethodological concerns, focuses on a study area to the north of the second section,which focuses on empirical observations of moulin and catchment spatial variability.10Chapter 2Flow routing for delineatingsupraglacial meltwater channels2.1 SummaryGrowing interest in supraglacial channels, coupled with the increasing availability of high-resolution remotely sensed imagery of glacier surfaces, motivates the development and testingof new approaches to delineating surface meltwater channels. I utilized a high-resolution (2m) digital elevation model of parts of the western margin of the Greenland Ice Sheet (GrIS) andretention of visually identified sinks (i.e., moulins) to investigate the ability of a standard D8flow routing algorithm to delineate supraglacial channels. I compared these delineated chan-nels to manually digitized channels and to channels extracted from multispectral imagery. Idelineated GrIS supraglacial channel networks in six high-elevation (above 1000 m) and onelow-elevation (below 1000 m) catchments during and shortly after peak melt (July and August2012), and investigated the effect of contributing area threshold on flow routing performance.I found that, although flow routing is sensitive to data quality and moulin identification, it canidentify 75% to 99% of channels observed with multispectral analysis, as well as low-order,high-density channels (up to 15.7 km/km2 with a 0.01 km2 contributing area threshold) ingreater detail than multispectral methods. Additionally, I found that flow routing can delin-eate supraglacial channel networks on rough ice surfaces with widespread crevassing. Ourresults suggest that supraglacial channel density is sufficiently high during peak melt that lowcontributing area thresholds can be employed with little risk of overestimating the channelnetwork extent.2.2 IntroductionSupraglacial meltwater channels are ubiquitous fluvial features in the ablation zones of glaciersand ice sheets (e.g. Knighton, 1981; Rippin et al., 2015; Smith et al., 2015; Wyatt and Sharp, 2015).They are of intrinsic interest to researchers as unique fluvial systems, as well as for the impor-112.2. Introductiontant roles they play in glacial hydrology. As fluvial systems, researchers have been interestedin their morphometrics and controlling processes (e.g. Dozier, 1970; Gleason et al., 2016; Karl-strom et al., 2013; Knighton, 1981; Yang and Smith, 2016; Yang et al., 2016a), and in how theirform may encode landscape processes (Karlstrom and Yang, 2016). Hydrologically, supraglacialmeltwater channels influence the spatial and temporal distribution of surface water inputs toen- and sub-glacial hydrological systems (Banwell et al., 2013; Joughin et al., 2013; Sundal et al.,2011; Wyatt and Sharp, 2015), which in turn have been shown to influence glacier and ice sheetdynamics (e.g. Schoof , 2010; van de Wal et al., 2008; Wyatt and Sharp, 2015). This recent andgrowing interest has been concurrent with a growing availability of high-resolution remotelysensed imagery and topographic data of ice surfaces (DGGS Staff , 2013; Noh and Howat, 2015;Shean et al., 2016; USGS, 2012; Yang and Smith, 2013).Relatively little is known about the characteristics of these channel networks and their con-trols, in part because of a dearth of data on their form and behavior through space and time(Fountain and Walder, 1998; Irvine-Fynn et al., 2011). Most data pertaining to supraglacial chan-nel form and process come from spatially and temporally limited field studies carried out inthe 1970s–1980s (Dozier, 1970, 1976; Ferguson, 1973; Iken, 1981; Knighton, 1981; Marston, 1983;Parker, 1975; Thomsen, 1986). Recently, remotely sensed imagery of ice surfaces has facilitatedvarious approaches to supraglacial channel delineation. Manual delineation of meltwaterchannels and lakes has been widely used with imagery from various satellites. Landsat ETM+has been used to assist with meltwater channel digitization on the Devon Ice Cap (Wyatt andSharp, 2015) and the Greenland Ice Sheet (GrIS) (Fitzpatrick et al., 2014; Lampkin and VanderBerg,2014; Yang and Smith, 2016). Optical imagery from WorldView-1 has been used to delineatechannels on the GrIS (Das et al., 2008; Joughin et al., 2013; McGrath et al., 2011). Landsat-7/8, L-band ALOS-PALSAR, and C-band RADARSAT images were used to detect mainstem reachesof large supraglacial meltwater channels on the southwest GrIS (Poinar et al., 2015).Multispectral methods, in which channels are automatically delineated based on the spec-tral reflectance of water relative to snow and ice, have also been employed with significantsuccess. Smith et al. (2015) and Yang et al. (2015) used multispectral methodologies on high-resolution ( 2.0 m) WorldView-2 imagery of the southwest GrIS to produce the most extensivedataset to date, and Yang et al. (2016a) analyzed fluvial morphometry of these supraglacialriver networks. With topographic data, Rippin et al. (2015) used drone imagery to producehigh-resolution (0.01 m) surface Digital Elevation Models (DEMs) of an alpine glacier in Sval-bard, and identified channels using an iterative elevation threshold. DEMs have also increas-ingly been used to delineate channels through flow routing methods modified for use on icesurfaces (Andrews, 2015; Arnold et al., 2014; Banwell et al., 2012; Karlstrom and Yang, 2016; Yanget al., 2015).Each methodology has its limitations. Manual digitization is labor-intensive and sensitiveto user bias. Multispectral methods have been highly successful (Smith et al., 2015). However,at commercially available multispectral image resolution ( 2 m), they underestimate channel122.2. Introductiondensity by as much as sevenfold when compared to manually digitized channel networks(spatial resolution 0.5 m) (Smith et al., 2015). Additionally, both methods have been foundto be of limited utility at low elevations on the GrIS where supraglacial rivers are smaller,less dense, and more intersected by crevasse fields (Yang and Smith, 2016). Furthermore, thederived network does not directly contain spatial information on channel network linkages,and requires additional geometric processing to create continuous networks and derive basicconnectivity data (Yang and Smith, 2013, 2016; Yang et al., 2016a).In terrestrial landscapes, flow routing is a standard and widely accepted methodology fordelineating channels from DEMs (Tarboton et al., 1991). This method assumes an impermeablesurface with no storage and requires that sinks (low elevation pixels surrounded by higherelevation pixels) be leveled to force flow routing from drainage boundary to catchment outlet(Arnold, 2010). This processing step is problematic in ice surfaces, however, where surface wa-ter drains via moulins (sinks on the ice surface) to the englacial hydrological system (Arnold,2010; Gulley et al., 2009a; Nienow and Hubbard, 2005). Yang et al. (2015) and Rippin et al. (2015)both concluded that flow routing overestimates hydrological network extent for this reason.Studies employing flow routing to delineate supraglacial meltwater drainage pathways haveemployed various approaches to this problem. Banwell et al. (2012) defined all depressions inmoderate resolution (30 m) DEMs as supraglacial lakes, and modeled flow routing over theGrIS by employing a lake filling model. Banwell et al. (2013) and Arnold et al. (2014) furtherassumed these lakes drain through moulins, and therefore preserved all topographic depres-sions as sinks during flow routing over moderate resolution DEMs (30 and 90 m, respectively).Yang et al. (2015) tested an automated sink retention size threshold (termed a “depression areathreshold”), but advised caution when using flow routing over moderate resolution (30–40 m)DEMs to delineate supraglacial channels. To overcome the limitations of automated moulindetection, non-automated sink preservation has been used by Andrews (2015) and Karlstromand Yang (2016), in which manually identified moulin locations are preserved as sinks duringflow routing over moderate- (25 m, resampled from 10 m) and high-resolution (2 m) DEMs,respectively.Flow routing has many advantages with regard to channel delineation, even in the rela-tively complicated case of ice surfaces. Firstly, it is a widely used, relatively simple method ofchannel delineation, and a basic feature of most GIS platforms (e.g. Conrad et al., 2015; ESRI,2016; GRASS Development Team, 2015). Channel locations are derived based on a physical re-lationship, i.e., the geometric network relationship between one channel pixel and anotheris physically interpretable based on their elevation differences and positions in the network.Contrast this to, for example, manually digitized channel networks, in which it can be dif-ficult to distinguish flow direction between pixels, and spatial proximity may not be a goodindicator of network connectivity. Network connectivity is implicit in flow routing; no furtherprocessing is needed to connect channel segments and order networks. This avoids the com-plication of shape analysis (e.g. Yang and Smith, 2016) or linear feature detection (e.g. Yang and132.3. Materials and MethodsSmith, 2013), and facilitates the easy coupling of channel networks with hydrological mod-els and topographically derived morphometrics such as stream ordering, contributing area,and gradient (e.g. Maidment, 2002). Furthermore, using flow routing to delineate supraglacialchannels provides the opportunity to capitalize on the growing availability of often free, high-resolution topographic data of a variety of glacier and ice sheet surfaces (DGGS Staff , 2013;Jóhannesson et al., 2003; Noh and Howat, 2015; Shean et al., 2016).In this paper, I assess the performance of flow routing on high-resolution (2 m) ice surfaceDEMs (Noh and Howat, 2015) with targeted (non-automated) moulin identification to system-atically assess the utility and limitations of flow routing in supraglacial channel delineation. Iemploy this methodology to assess the performance of flow routing in delineating supraglacialchannels without the complicating factor of automated sink detection, i.e., under the assumedconditions of known sink locations. Both crevasses and moulin-type sinks are retained dur-ing sink filling, based on manual identification of sinks. The impacts of channel initiationthreshold on the flow routing-derived channel network extent and accuracy are investigated.This assessment is made relative to two independent datasets—one manually digitized fromWorldview-1 0.5 m resolution panchromatic imagery, and one of channels derived throughmultispectral analysis, developed by Smith et al. (2015) (see Table 2.1).Channel Extraction MethodImagery Source Acquisition Date(s)Resolution (m)ReferenceFlow Routing (FR)SETSM DEM (derived from WorldView-1)24 July 2012 (Catchments 1–3)2.0Noh and Howat (2015) 20 August 2012 (Catchments 4 and 5)20 August 2011 (Catchment 6)Multispectral (MS) WorldView-2 18–23 July 2012 2.0Smith et al. (2015)Manual Digitizing (MD)WorldView-121 July 2012 (Catchment 7) 0.5 n/a20 August 2012 (Catchment 5)Table 2.1: Summary of datasets used.2.3 Materials and MethodsSupraglacial meltwater channels were extracted from DEM datasets using a flow routing (FR)workflow described below. The performance of this workflow is evaluated by comparing theresults to two independent datasets: (1) meltwater channels derived from multispectral (MS)imagery by Smith et al. (2015); and (2) manually digitized (MD) fluvial channel networks, asdescribed below. A summary of the datasets is provided in Table 2.1. The following sectionsdescribe the data sources, channel extraction workflow, and error analysis.142.3. Materials and Methods2.3.1 Datasets and study areaThe spatial and temporal coverage of the study was constrained by the availability of data. Assuch, this study is focused on Greenland’s southwestern margin where there is a coincidenceof multispectral and topographic datasets. The study area displays extensive supraglacialchannel development (Smith et al., 2015), and has been the focus of previous studies on theGrIS’s surface meltwater systems (e.g. Fitzpatrick et al., 2014; Gleason et al., 2016; Karlstrom andYang, 2016; Smith et al., 2015). A total of seven meltwater catchments were delineated throughflow routing on the ice surface. Their characteristics are described in Table 2.2 and their spa-tial extents are shown in Figure 2.1. Six of the catchments (Catchments 1–6) were above 1000m in elevation (‘high’-elevation catchments), where the ice surface has a lower slope and thesurface topography is generally uniform (Clason et al., 2015; Yang et al., 2016a). Catchment 7(‘low’-elevation catchment) was delineated on Russell Glacier, where surface slopes are higherand numerous areas of extensional and compressive ice flow create a rough surface, impart-ing structural elements to the ice surface. To date, there has been no successful, automatedmethodology developed for delineating supraglacial channels at low elevation, often heavilycrevassed parts of the GrIS (Yang et al., 2016a). This seventh location was chosen specificallyto assess the performance of flow routing in delineating channels across a range of ice sur-face morphologies. Only one low-elevation catchment was selected as MS-derived meltwaterchannel data were not available for low-elevation locations and the comparison was thereforelimited to MD channels.CatchmentCentral CoordinatesArea (km2)Mean Slope (%)Mean Elevation (m)DatasetsFR MS MD149°12′48.853″W23.1 11.0 1172 X X67°15′3.539″N249°13′58.816″W71.5 11.9 1147 X X67°8′52.932″N349°11′56.907″W15.5 10.7 1230 X X67°13′15.191″N449°14′5.529″W78.2 11.6 1119 X X67°2′24.692″N548°52′8.931″W79.2 13.0 1244 X X X66°41′56.575″N648°18′22.375″W31.2 8.2 1528 X X66°55′17.935″N749°44′22.47″W42.2 16.1 838 X X67°6′19.874″NTable 2.2: Summary of study catchments.152.3. Materials and Methods67°N 10’66°N 40’49°WMS datasetCatchment 1Catchment 2Catchment 3Catchment 4Catchment 5Catchment 6Catchment 70 20 km1000m1500m80°N70°N60°NGreenland60°W30°WNFigure 2.1: Extent of multispectral (MS) dataset, with the six overlapping flow routing-derived catchments identified as polygons. Base map from ESRI ArcMap WorldImagery (2016) (Jóhannesson et al., 2003). Star in inset map marks approximate loca-tion.The dataset used to test the channel extraction by flow routing is a 2 m resolution DEM oflarge parts of the western GrIS. This ‘Surface Extraction with TIN-based Search-Space Mini-mization (SETSM)’ DEM is derived from overlapping DigitalGlobe Worldview-1 satellite im-ages and has reported root mean square errors of 3.8 and 2.0 m in the horizontal and vertical,respectively (Noh and Howat, 2015). Experience working with the data suggests that errors ingeo-referencing may produce local errors and offsets of several tens of meters between fea-tures in the DEM and their real-world coordinates. An offset was similarly observed by An-drews (2015). Flow routed (FR) channel networks were delineated using flow routing over the162.3. Materials and MethodsSETSM DEM in all seven study catchments, as described in Section 2.2.The six high-elevation GrIS catchments (Catchments 1–6) extracted from the SETSM DEMusing flow routing overlap in space with multi spectral-derived meltwater channel networksextracted from 2 m resolution WorldView-2 images dated to July 2012, provided by Smithet al. (2015) (Tables 1 and 2). This MS dataset represents channels containing meltwater at thetime of image acquisition, as identifiable at the 2 m resolution. These overlapping FR andMS datasets are used to compare the channels derived with multispectral and flow routingmethodologies. The FR dataset for Catchments 1–3 overlaps in time with the MS datasets, andis dated to one month after the MS dataset for Catchments 4 and 5. The FR dataset for Catch-ment 6 is from August 2011. Data availability and overlapping spatial extent for the MS andFR datasets limit the spatial extent of our study area and restrict the temporal resolution of ourdata choices. However, numerous observations suggest that supraglacial meltwater channels(particularly high-order, mainstem channels) perennially re-occupy the same spatial location(Ewing, 1970; Knighton, 1981; Smith et al., 2015), even as abandoned channels are advecteddown-glacier (Karlstrom and Yang, 2016). Furthermore, Lampkin and VanderBerg (2014) observethat the extent of the supraglacial network on the GrIS, observable at 15 m resolution, did notchange significantly between early July and early August 2007. Considering melt on the GrISin 2012 was higher than in 2007 (Nghiem et al., 2012), it follows that the fluvially incised chan-nel network was sufficiently unchanged between the July and August image acquisition datesto justify comparison between datasets from the two dates. Again, I distinguish here betweenwater-filled channels, which cannot explicitly be delineated by flow routing (the MS dataset)and whose extent is likely sensitive to timing, and the fluvially incised network, which is un-likely to change following peak melt and will persist until erased by either ablation of the icesurface or burial by seasonal snow.A third and manually digitized dataset (MD) was also developed from 0.5 m panchromaticimages from the WorldView-1 satellite, obtained for a high-elevation GrIS location that over-laps with both a MS- and FR-derived catchment (Catchment 5) and the low-elevation GrISoutlet glacier (Russell Glacier, Catchment 7). The channels in these two images were digitizedmanually for comparison with FR-delineated channels in both Catchment 5 and Catchment7, as well as with the MS dataset in Catchment 5. Digitized channels were identified as rec-tilinear features that are darker and appear recessed relative to the surrounding ice surfaces,and were systematically digitized from the terminal moulin up to the smallest channels thatcould be readily distinguished in the 0.5 m resolution imagery. They may or may not havebeen actively filled with water at the time of image acquisition—it is not objectively possibleto differentiate water from shadow, particularly in small channels, and therefore I assume thatthis dataset represents the full extent of the fluvially incised channel network, identifiable at0.5 m resolution, regardless of the presence of meltwater at the time of image acquisition. Al-though I developed this dataset to represent an objectively mapped channel network, there area number of errors inherent in manual digitization, including, for example, user bias, the scale172.3. Materials and Methodsof the image, limitations of panchromatic imagery (particularly with regards to illuminationand shading), and difficulty distinguishing flow direction and connectedness, particularly forsmall channels (e.g. Paul et al., 2013).The three datasets represent three distinct types of channel networks. The MS dataset rep-resents active meltwater channels and the MD dataset represents channels that are fluviallyincised but may or may not actively contain meltwater at the time of image acquisition. TheFR dataset represents the likely path that water takes through a landscape, where networkextent is set by a user-specified channel initiation condition. Therefore, these datasets are fun-damentally different, and comparing them in this context is not done intended to suggest thatthey are equal attempts to represent the same processes. However, the MD dataset representsfluvial incision at times of maximum melt, and therefore by definition should contain the MSdataset. For this current study, then, these datasets provide independent opportunities to as-sess flow routing for reproducing channel network structure and characteristics, rather than adirect comparison of their abilities to reproduce the same processes.Channel network extractionThe general methodological workflow is illustrated using the DEM from which Catchment 7was extracted (Figure 2.2). In the flow routing channel extraction method, flow is routed tomoulin locations in the SETSM DEM which are preserved as sinks during the sink filling stageof a flow routing workflow. Two types of potential englacial entry locations are included: sinksin crevasse fields (Figure 2.2, step 3a) and isolated moulins (Figure 2.2, step 3b). Although notall crevasses are necessarily englacial entry locations, preserving them as sinks in flow routingconservatively overestimates the presence of sinks in a DEM and limits the catchment extent(Figure 2.2, step 6). Both crevasse fields and isolated moulins were visually identified andmapped from the DEM and hill shade layers, and sinks within (crevasse fields) or at (moulins)those locations were preserved during sink filling and flow routing. Moulins were readilyidentifiable as places where channels terminated abruptly, or depressions with no outlets (seeFigure 2.2, step 3b for an example). Crevasse fields were identified as areas with abundantparallel or cross hatched and discontinuous linear depressions. A 1 by 1 km grid was overlainon the DEM, and a grid search pattern was used to assist in the visual identification of sinks(Figure 2.2, steps Steps 2, 3a and 3b ).182.3. Materials and Methods1. Raw DEM 2. 1x1 km search grid3a. Identify crevasse fields 3b. Identify moulins4. Fill sinks - preserve sinks at moulin locations and inside crevasse fields.5. Derive flow direction and flow accumulation.Crevasse FieldsCatchments (>2500m )ChannelsSearch GridHigher ElevationLower Elevation6. Generate sink catchments 7. Define stream channels210 km150 m5 kmNFigure 2.2: Workflow for channel delineation, exemplified by the extraction of Catchment7 from a DEM tile (SETSM tile 15_38_5_5). The processing steps are described inSection 2.2. Only catchments above 2500 m2 are shown in Step 6; however, manymore sinks catchments are present after flow direction processing.Flow routing was carried out using the ArcHydro toolbox in ArcGIS 10.3.1 (ESRI, 2016),which allows for sink preservation during sink filling and flow routing using a standard D8flow routing algorithm (O’Callaghan and Mark, 1984). A multi-directional flow (MDF) rout-192.3. Materials and Methodsing algorithm was also tested using the MATLAB based TopoToolbox (Schwanghart and Kuhn,2010). Similar to findings by Karlstrom and Yang (2016), channel networks delineated with MDFwere similar to those delineated with D8. I therefore employ only a D8 flow routing algorithmhere as it is coupled with sink preservation tools in ArcHydro (‘Fill Sinks’ tool using IsSinkfield, followed by ‘Flow Direction with Sinks’ tool). Flow was routed over the ice sheet DEMto visually identified sink locations, generating flow direction, flow accumulation, and sinkwatershed layers (Figure 2.2, Steps 4, 5, and 6). Watersheds of interest can then be extractedfor further analysis (e.g., Figure 2.2, Steps 6 and 7, showing the extraction of Catchment 7).Streams were defined by employing a flow accumulation area threshold (e.g., Figure 2.2, Step7). Choosing an appropriate channel initiation threshold criterion is a pervasive problem forhydrologists and geomorphologists (Montgomery and Foufoula-Georgiou, 1993). It is an unre-solved question in alluvial systems, and has not been addressed in supraglacial systems. Inthis study, I investigate the methodological effects of channel initiation area thresholds on flowrouting and meltwater channel delineation as described in the following section.Channel initiation thresholdField observations of supraglacial channel initiation suggest several different mechanisms ofchannel initiation. At the beginning of the melt season, saturated snow may transition intochannelized meltwater flow through the formation of rills (Kamintzis, 2015; Karlstrom and Yang,2016; Knighton, 1981; Stenborg, 1968), overflow from slush swamps, lakes or hollows (Brykała,1998; Dozier, 1970; Ferguson, 1973; Hodgkins, 1997), or slush avalanches (Fountain, 1996; Onesti,1987). Although these processes have been empirically described, there has been only onenumerical description of the hydrodynamic controls on rill-like supraglacial channel initiation(Mantelli et al., 2015). There is therefore no generalizable and mechanistic methodology forpredicting channel initiation points in a flow network. In terrestrial landscapes, inflections inslope/area relationships are often used (Montgomery and Dietrich, 1994, 1988; Tucker and Bras,1998); however, it is not yet clear if these relationships can be extended to supraglacial sys-tems (Karlstrom and Yang, 2016), and, furthermore, channel initiation locations are likely belowthe resolution of available imagery (Karlstrom et al., 2013; Mantelli et al., 2015). In previoussupraglacial research, Karlstrom and Yang (2016) employed a threshold of 0.02 km2 and notedthat the resultant channel networks were sensitive to this threshold. Researchers using flowrouting will likely need to make locally based judgements on the initiation threshold they em-ploy. In the present study, I avoid testing a specific initiation area and instead specify a rangeof 10 channel initiation threshold areas, bounded on the lower end by the high density MDnetworks, and on the upper end by the lower-density MS network.I assume that the MS-derived initiation points represent an upper bound for the initiationarea threshold of the active meltwater network, and the MD-derived initiation points rep-resent a lower bound on the full extent of the fluvially formed network. I calculate initiationpoint density based on the density of first-order channel endpoints in the MS and MD datasets202.3. Materials and Methods(number of initiation points divided by catchment area, giving us the average area per initia-tion point). An MS-derived initiation point density is available for each of the high-elevationcatchments, and MD-derived initiation point densities are available for Catchments 5 and 7.The MD-derived point density is applied as the lower bound on channel initiation densities inthe other high-elevation catchments, where there are no MD datasets.Our objective here is not to suggest that the MD dataset channel head densities are directlygeneralizable to other parts of the GrIS. Indeed, a generalizable methodology for identifyingchannel initiation points is a complicated issue that is highly dependent on scale and local gov-erning processes (Mantelli et al., 2015; Montgomery and Dietrich, 1988; Zhang and Montgomery,1994), and is beyond the scope of this work. In this work, I investigate the MD thresholdas a suitable and conservative channel initiation area recommendation value during or rightafter peak melt time. In the absence of process based approaches to mapping channel initia-tion locations, initiation area thresholds may be useful for the purpose of helping researchersmap general surface drainage patterns rather than to reveal the physical processes of channelinitiation.Quantification of discrepancy between datasetsSystematic and local offsets sometimes exist between the channel networks from the differentdatasets due to stitching and georectification errors, and potentially, in the case of datasets withdifferent acquisition dates, due to ice advection. This makes comparisons of the three flow ex-traction methodologies difficult. A direct spatial comparison would necessitate a large allow-able offset, which would complicate the analysis in such a dense network of channels. Becauseour interest here is an evaluation of flow routing as a means of delineating supraglacial chan-nels, and not an assessment of the datasets, the datasets were spatially adjusted to ensurespatial overlap and facilitate error quantification. This was done by linear translation usingreadily identifiable channel junctions as displacement links between datasets.Many flow routing algorithms, including D8, cannot produce realistic channel segmentsin the absence of topographic variation, including over lakes, and furthermore cannot capturechannel bifurcations. As such, channel lengths in lakes and bifurcations were removed fromthe datasets. Lakes were visually identified as areas with multiple adjacent, perfectly straightchannel segments in the FR dataset, and large braided sections were automatically identifiedby building polygons from enclosed areas in the MS dataset. Examples of these lakes andbraided areas can be seen in Figure 2.4 .Similarity between the datasets was determined as follows. (1) The lowest drainage densitydataset was identified, hereafter referred to as Dataset A (the MS dataset in Catchments 1–6and the FR dataset in Catchment 7); (2) A 15-m buffer was applied to Dataset A (buffer sizedetermination is explained below); (3) The comparison dataset (Dataset B) was clipped to thebuffer; (4) Before comparing the datasets, I removed any lengths of channel equal to or shorterthan twice the buffer width in the clipped Dataset B (i.e., any lengths of channels shorter than212.4. Results30 m). This was done to avoid artificially high match rates due to including inclusion of theends of tributary channels from Dataset B incidentally clipped to the Dataset A buffer.Similar to Yang et al. (2015), two metrics of comparison were employed: (1) the length ofthe clipped Dataset B relative to the length of Dataset A, which provides a measure of the sim-ilarity between Datasets A and B; and (2) the length of Dataset B outside of the buffer relativeto the total length of Dataset B, which provides a measure of how much additional channelnetwork is delineated in Dataset B. Yang et al. (2015) refer to these two metrics of comparisonas ‘completeness’ and ‘miscoding’; however, this terminology is not adopted here. As dis-cussed in Section 2.1, the MS dataset represents the lower-density active meltwater channelsand the MD and FR datasets target the fluvially incised network; these are therefore not directcomparisons of the same system, and I therefore avoid the implication of a direct comparisonimplicit in the terms ‘completeness’ and ‘miscoding’.Appropriate buffer size was determined by comparing match rates between the datasetsover a range of buffer sizes. A 15 m buffer size was chosen because additional increases tobuffer size produced only moderate gains in match rate and, because the data had been spa-tially translated, I wanted to employ as small a buffer size as possible while still allowing forsome unavoidable minor offset between the datasets.2.4 Results2.4.1 Comparison between flow routing and multispectral datasetsMoulin identificationBy design, the FR delineated catchments contain only one sink: large, plainly visible in theDEM, and easily identifiable moulins or crevasses drain the basins at the outlets. In all six FRdelineated catchments, the location of these moulins coincides with the location of a moulinidentified in the MS imagery. However, in three out of six cases, additional moulins were visu-ally identified in the MS dataset that were not identified from the DEMs. Table 2.3 shows thebreakdown of moulins that were identified in the MS imagery but not in the DEM inspection.The three catchments in which additional moulins were identified are all more than twice aslarge as the remaining three catchments. In all of these cases, three moulins were missed, cu-mulatively draining 1.4%, 1.9%, and 8.0% of the three FR catchments, respectively. Thus, onlyone catchment was overestimated by more than 2% as a product of failing to identify smallmoulins.222.4. ResultsCatchmentMissed Moulins MS Drainage Density (km/km2)MS Initiation Threshold (km2)Number Area Drained1 0 0.0% 2.8 0.32 3 8.0% 2.2 0.53 0 0.0% 2.4 0.34 3 1.4% 2.0 0.55 3 1.9% 1.4 0.86 0 0.0% 0.6 1.8Table 2.3: Summary of moulins identified from optical imagery in the multispectral im-agery that were not visually identified from the DEM and hill shade layers, anddrainage densities and initiation thresholds derived from the MS dataset.Comparison of Drainage NetworksDrainage density of the MS channels ranges from 1.4 to 2.8 km/km2 for Catchments 1 to 5, andis 0.6 km/km2 in Catchment 6 (Table 2.3). The channel initiation contributing area thresholdsestimated from the MS dataset (first-order channel endpoints over catchment area) vary from0.3 to 1.8 km2, although only Catchment 6 has an initiation area over 1 km2 (Table 2.3). Allof these initiation area thresholds are more than an order of magnitude higher than estimatedfrom the MD network in Catchment 5, which has an initiation contributing area of 0.01 km2. Asa result, the FR channel networks delineated across this range of initiation areas have a widerange of extents and thus a range in the match percentages with the MS dataset and the FRdrainage densities. Regardless of the initiation threshold used, flow routing produced higherchannel densities than observed with multispectral methods, with this difference increasingwith decreasing channel initiation area (Figure 2.3). This increase was particularly pronouncedin Catchment 6, in which drainage density increased by more than an order of magnitude aschannel initiation area decreased towards the measured digitized initiation area threshold.232.4. Results10–2 10–1 100 1010510152025301Catchment23456DdFR/Dd MSInitiation area (km2)Figure 2.3: The ratio of flow-routed (FR) to multispectral (MS) drainage density as a func-tion of initiation contributing area. Initiation area is bounded on the lower end bythe manually digitized (MD) initiation area for Catchment 5 and on the upper endby the measured initiation area for each MS catchment.Figure 2.4 visually compares the MS and the FR datasets in three catchments—Catchment5 with the greatest percentage overlap between the two datasets, Catchment 6 with the leastoverlap, and Catchment 2, which is intermediate. The other three catchments are includedin the supplementary material. The FR results illustrated here are those produced using thelowest initiation threshold (0.01 km2), close to the digitized initiation threshold and where thematch rates are the highest (Figure 2.3). The depicted FR channels have been clipped to theMS buffer to highlight only those lengths of channel that overlap between the two, excludingthe additional lower-order, higher-density channels in the FR dataset outside of the MS buffer.It is clear that, generally, the shape, extent, and placement of the flow routing-derived catch-ment boundaries are visually consistent with the MS dataset in all test catchments. Even inCatchment 6, with low match rates at higher channel initiation areas, the general morphologyof the network is largely similar between the two datasets.242.4. Results0 MoulinsFR channelsMS channelsBraided areasLakes0 2.5 km0NNN2.5 kmA) Catchment 2B) Catchment 5C) Catchment 61 kmFigure 2.4: A visual comparison of flow-routed (FR) and multispectral (MS) channel net-works in three of the test catchments (Catchments 2 (A); 5 (B); and 6 (C)).Figure 2.5 plots the comparisons between the FR and MS datasets. There is a range in theperformance of flow routing in the different catchments, and the match percentage dependsheavily on the initiation threshold. For example, in Catchment 5, at the (largest) MS-derivedchannel initiation area, around 72% of the MS dataset is matched by the FR dataset, and almost85% of the FR dataset is within the MS dataset buffer. In comparison, in Catchment 6, only 31%of the MS dataset is matched by the FR dataset, and only 28% of the FR dataset is within theMS dataset buffer. At the (smallest) MD-derived channel initiation area, 99% of the MS datasetis matched by the FR dataset in Catchment 5, and 75% is matched in Catchment 6. Catchment252.4. Results6 has the lowest match rates in general. It also has the largest MS channel initiation area (bymore than double), the highest mean elevation of the catchments, and a temporal offset ofa year between the MS and the FR datasets. Overall, at channel initiation areas similar tothose of the MD network in Catchment 5, all the catchments have match rates above 75%, atwhich point the majority of the FR channel network falls outside of the MS buffer (i.e., the FRnetworks have significantly higher density than the MS networks).010020406080%FR inside MS buffer10010–110–2050100%MS matched by FRInitiation area (km2)1Catchment23456Figure 2.5: The performance of flow routing (FR) relative to multispectral (MS) methodsacross a range of initiation thresholds, bounded on the low end by the digitizedinitiation threshold from Catchment 5 and on the upper end by the initiation pointdensity derived from the MS dataset for each catchment. Solid lines represent thepercent of the FR dataset within the MS buffer (left y-axis), and the dashed linesrepresent the percentage of the FR network length inside the MS buffer (right y-axis).Although flow routing captures more than 75% (and as much as 99%, in the case of Catch-ment 5) of the MS channels, these high percentages are only achieved with a small channelinitiation threshold, up to two orders of magnitude smaller than the active channel extentmeasured from the MS dataset. As a result, in all catchments, flow routing at low channel ini-tiation thresholds produces channel networks that are between five and eight times as dense(or more than 30 times as dense in the case of Catchment 6). If the full extent of these net-works were to be plotted in Figure 2.4, the result would be a much more extensive networkwith many lower-order channels (Figure 2.7 shows, for illustrative purposes, the full extent of262.4. Resultsthe FR dataset for Catchment 5 before clipping). As can be seen in Figure 2.6, the MS datasetcontains less than 10% of the FR first-order channels (using Strahler stream ordering), and lessthan 20% of the FR 2nd order channels identified with a 0.01 km2 initiation threshold.It appears that in most cases, MS channels are most similar to FR channels higher than 5thorder, based on the ordering of the FR stream network. The drop in match rates for order 7is likely a result of the fact that there are only three catchments with seventh-order channels.These sections tend to be quite short, therefore mismatched segments become proportionallymore pronounced. Manual digitization of Catchment 5 therefore serves to validate whetherthese high-density, lower-order channels are indeed real. The results of that comparison aresummarized in the following section.0204060801001204321 5 6 7Stream order% MatchFigure 2.6: The match rates between multispectral (MS) and flow-routed (FR) datasetsin the six high-elevation catchments by FR Strahler stream order, using a 0.01 km2initiation threshold.2.4.2 Comparison with Manually Digitized DatasetsComparison of the FR dataset with manually digitized datasets serves two aspects of the meth-ods evaluation. Firstly, the higher channel density and larger number of lower-order chan-nels in the FR dataset (relative to the MS dataset) can be compared to the MD dataset whichincludes low-order and potentially dry channels. Secondly, by utilizing the MD dataset forCatchment 7, the performance of flow routing on varied and rougher ice surfaces can be as-272.4. Resultssessed. In this section I first compare the results of the MD, MS and FR methods in the oneoverlapping catchment (Catchment 5), and thereafter evaluate the performance of flow rout-ing in Catchment 7.High-Elevation CatchmentAll three datasets (FR, MD, and MS) overlap in Catchment 5. A total channel length of 1590km was digitized in this catchment, and the manually digitized dataset has very high drainagedensity (20.2 km/km2), 1.7 times the drainage density identified with flow routing (12.2 km/km2based on digitized initiation point density of 0.01 km2) and 14 times the density of channels de-lineated with multispectral methods (1.4 km/km2). Both MS and FR methods reproduced thegeneral network structure (Figure 2.7), as well as channels at a smaller scale (Figure 2.7B). TheMS dataset matched around 5% more of the MD channel length than did the FR dataset, suchthat 93% of MS channels overlapped with MD channels as opposed to 88% of the FR channels.A total of 857 km of FR-derived channels overlap with MD channels (54% of the total MDchannel length), whereas the MS channels overlapped with 106 km of MD channels (7% of thetotal MD channel length), consistent with the fact that MS networks reflect the lower-densityactive meltwater channels.282.4. Results0 2 4 kmNMoulinsMS channelsFR channelsMD channels0 500 m(A)(B)NFigure 2.7: A comparison of the manually digitized (MD) network, the flow-routed (FR)network with a channel initiation condition based on the MD network, and the mul-tispectral (MS) dataset for Catchment 5. In (A), the datasets are overlain in the orderof the legend entries shown, with the lowest density network (MS) on top, and thehighest density (MD) on the bottom. This figure illustrates the relative morpholo-gies and extents of the networks. In (B), the gray square from (A) is enlarged and theorder in which the networks are overlain is reversed. The widths of the networks aredifferentiated for illustrative purposes—MS is the thickest, and MD is the thinnest.This figure illustrates the overlap between the networks.292.4. ResultsThe performance of flow routing relative to the MD dataset varied by stream order. Streamorder is scale-dependent, and I therefore use stream order here as a proxy for relative streamsize. Figure 2.8 demonstrates that for larger, higher-order channels, match rates in Catchment5 can exceed 95%; however, these match rates drop to close to 70% for first-order channels. Partof this mismatch can be attributed to the fact that low-order channels high up in the networkoverestimate the upward network extent because it is difficult to specify where channels start.This can be seen visually in the upper, westernmost reaches of Figure 2.7 where long, overpredicted low-order FR tributaries can be observed.0 1 2 3 4 5 6 75060708090100% MatchStream OrderCatchment 5Catchment 7Figure 2.8: The match percentages between manually digitized (MD) and flow-routed(FR) channels in Catchment 5 and 7, but by FR-derived Strahler stream order.Low-Elevation CatchmentThe ice surface of low-elevation Catchment 7 is heavily crevassed, yet the main stem channelhas a maximum flow length of 21.6 km and width of around 27 m at the moulin entry loca-tion, suggesting that there can be significant supraglacial channel development even in thepresence of crevasse fields. Digitizing these surface channels was complicated by the surfacemicrotopography—the flow direction in the smaller channels is not always clear to the nakedeye, highlighting the advantage of incorporating DEMs into channel network delineation. Intotal, 532 km of channel length was digitized in Catchment 7. Manual digitizing identified achannel network that was 1.6 times denser than the FR network (12.35 km/km2 as opposed to302.4. Results7.50 km/km2).Figure 2.9 illustrates the MD and FR networks in Catchment 7. The FR network is bothunderestimated relative to the MD network in one area (illustrated by the yellow square) andoverestimated in several others (as can be seen in the upper reaches of the network). Manualdelineation included an 11-km branch of channel that was attributed to an adjacent catchmentduring flow routing, highlighted by the yellow square in Figure 2.7. The reason for this dis-crepancy may be a large stitching error in the DEM, and may have resulted in flow beingrouted to an adjacent catchment. Additionally, six moulins in the upper reaches of the catch-ment were also identified during manual digitizing that were not identified during the moulingrid search for flow routing sink retention. Five of the moulins were located in a crevassed areaon the periphery of the catchment, and drained a total of 8.7% of the flow-routing delineatedcatchment. The sixth moulin was isolated from the others and drained 3% of the flow-routingdelineated catchment.312.4. Results0 3 kmMoulinsFR channelsMD channels0 500 m(A)(B)NNFigure 2.9: Comparison of flow-routed (FR) and manually digitized (MD) channel net-works in Catchment 7, the low-elevation catchment. As in Figure 2.6, in (A) the FRchannels are overlaid on MD channels to illustrate the channel networks morpholo-gies. In (B) the order is reversed, with FR on the bottom, illustrated with a thickerline, to highlight the overlap between the networks.322.5. DiscussionA total of 75% of the FR channel length matched MD channels. This is more than 10%less than the match rate in the high-elevation catchment. However, as in Catchment 5, matchrates are highly dependent on stream size/order. As seen in Figure 2.8, in the highest-orderchannels, match rates exceed 90%, and this match drops to below 65% for low-order channels.Again, as in Catchment 5, much of the mismatch in the low-order channels was located at theupper ends of channel reaches, as can be seen at the northernmost and southernmost extentsin Figure 2.9, and in the western part of the reach where a missed moulin resulted in overpredicted channel extent. This suggests that mismatch arises not from an inability to delineatechannels, but rather from the reliance on simple initiation area thresholds for specifying chan-nel initiation points, as well as moulins that were not identified during the moulin grid search.Close inspection of the more central parts of the network (for example, Figure 2.9B), show thatwhen channel initiation area is not under predicted, lower-order FR channels visually comparevery well with MD channels.2.5 DiscussionDEMs are increasingly available at high resolution for ice surfaces (DGGS Staff , 2013; Noh andHowat, 2015; Rippin et al., 2015; Shean et al., 2016), and reproducible, accurate methods of ex-tracting surface meltwater channels are necessary for quantifying their temporal and spatialcharacteristics. Flow routing has numerous advantages for delineating meltwater channels;however, previous research by Yang et al. (2015) on the use of flow routing over moderate-resolution ice surface DEMs found that (1) automated sink preservation thresholds missedover half of mapped moulins and resulted in significantly overestimated drainage network ex-tent; (2) mapped channel networks could be replicated by flow routing with match rates in ex-cess of 85% when a 600 m buffer was applied; (3) less than 40% of the FR drainage network didnot match the mapped drainage network, leading to the conclusion that these channels weremiscoded. My results using manual moulin identification and retention on high-resolutionDEMs suggest that moulin identification is indeed a persistent problem with flow routing,but that when sink location can be determined, flow routing over a high-resolution DEM canproduce channel networks that largely resemble both MS and MD channel networks, using abuffer size of only 15 m. Furthermore, our results support the suggestion of Yang and Smith(2016) that high-resolution DEMs may provide a viable solution to delineating supraglacialchannels from low-elevation GrIS catchments.2.5.1 MoulinsVisual identification of moulins is labor-intensive, and thus in some ways offsets the practicaladvantages of using flow routing over other channel delineation methodologies. However,there are no automatic moulin detection methodologies currently available, and manual iden-tification was used by Smith et al. (2015) in the production of the extensive MS dataset em-ployed here as well as by Karlstrom and Yang (2016) and Andrews (2015). I tested moulin identi-332.5. Discussionfication with DEM and hill shade layers, as optical imagery associated with the SETSM datasetwas not available to us. Digitization from topographic raster layers may be necessary in otherstudies when optical imagery is not available or accessible (e.g., LiDAR-derived DEMs), or isnot sufficiently well georeferenced to the topographic data to easily couple identified moulinswith the DEM (e.g. Doctor and Young, 2013).2.5.2 Flow Routing in High-Elevation CatchmentsIn almost all of the high-elevation catchments, flow routing did an exceptional job of delin-eating the shape, extent and placement of channels identified through multispectral methods.Match rates of up to 99% could be achieved between the FR network and the MS network. Sev-eral important limitations of the comparison are worth noting: firstly, lakes and large braidedareas were removed as they cannot be properly captured by D8 flow routing. Integrationof more sophisticated flow routing algorithms (e.g. Tarboton, 1997) with sink preservation al-gorithms (e.g. Maidment, 2002) may better handle flow dispersion and channel bifurcation.However, neither I nor Karlstrom and Yang (2016) observed significant differences in flow net-works derived with multiple flow direction or Dinf algorithms, respectively. Secondly, chan-nel density and the match rates are highly dependent on the initiation area used to define thechannels. Match rates were ubiquitously lower at large channel initiation thresholds basedon the low drainage density MS datasets, and higher at the smaller channel initiation thresh-olds generalized from the MD dataset. However, at lower channel initiation areas, most ofthe FR channel network lay outside of the MS dataset buffer because low channel initiationthresholds result in significantly higher channel densities. Therefore, although it was possibleto achieve high match rates, it was necessary to densify the network, implying that the activemeltwater extent is not solely or explicitly a function of contributing area, and a more adap-tive criterion for channel initiation may be required. Initial modeling work by Mantelli et al.(2015) suggested that slope, meltwater depth, and small scale roughness exert primary con-trols on channel spacing and initiation. Thus, accurate modeling of channel networks requiresa more flexible initiation criterion than channel initiation contributing area, and integrates theevolution of the supraglacial network over time.Our high-elevation catchment MD channel density of 20.1 km/km2 falls within the rangeof 6.0–31.7 km/km2 digitized in a catchment of similar elevation by Smith et al. (2015), whofound that digitized channel networks were at least seven times as dense as those derivedfrom multispectral imagery. Comparison of the FR network with these MD channel networkssuggests that almost 90% of the FR channels overlap with MD channels and may thus be flu-vially formed, that the full extent of the fluvial channel network is at a minimum five timesdenser than the multispectral network, and that MS channels generally represent higher-orderchannels in the full supraglcially incised networks. Therefore, for researchers interested inrepresenting channel networks, including lower-order channels, it appears that in the absenceof a more process-based and generalizable criterion for channel initiation, small channel ini-342.5. Discussiontiation areas can be prescribed without risk of overestimating the fluvially formed channelnetwork. However, the match rates between the FR and MD dataset decrease with decreasingstream order. Although match rates for first-order channels was were still above 70%, this issignificantly less confidence than can be had in higher-order channels.The MS and FR networks in Catchment 6 had the lowest match rates across the range ofchannel initiation thresholds tested. This may have been because there was a difference of ayear between the acquisition date of the DEM and MS imagery. However, the main channelsin the MS imagery matched FR channels at low channel initiation thresholds, suggesting thatthe time difference was at least not solely responsible for the mismatch. Catchment 6 differsfrom the other catchments in its higher elevation (an average of 1528 m a.s.l. versus a rangeof average elevation from 1147 to 1244 m a.s.l. in the other high-elevation test catchments)and lower slope (8.2% compared to a range of 10.7%–16.1% for the other catchments). It alsohas a much lower MS channel density than the other catchments. Very little data exist per-taining to the spatial variability of supraglacial channels or their controls; however, Yang et al.(2016a) noted a decrease in channel density with increasing elevation, supporting the notionthat at higher elevations less meltwater is available to incise channel networks. Fluvial chan-nel incision is dependent on channel slope and water discharge (e.g. Tucker and Bras, 1998); it istherefore possible that fluvial channels and networks in higher-elevation catchments, charac-terized by lower slopes and lower meltwater supply, are less developed than in lower elevationcatchments. As such, multispectral rather than flow routing methods may be better suited toidentifying meltwater channels in low-slope, high-elevation catchments.Conversely, the mismatch may arise from the presence of perennial snow and firn coverat this high-elevation catchment. The mean equilibrium altitude in this part of the GrIS isestimated as 1553 m a.s.l. (van de Wal et al., 2008). Observations suggest that extensive firn sat-uration and active surface runoff were present on snow covered ice above this average equi-librium line altitude in 2012 (Mikkelsen et al., 2016). Catchment 6 is located at a mean elevationof 1528 m, suggesting that supraglacial channel development was occurring in saturated firnand that subsurface meltwater flow may have been occurring in this catchment (Karlstromet al., 2014). This subsurface flow would not have been detectable by MS techniques, and maycontribute to some of the high mismatch in this catchment.2.5.3 Flow Routing in Low-Elevation CatchmentsAlthough the ‘low-elevation’ catchment in this study (Catchment 7) is still above 800 m, itscharacteristics differ significantly from the ‘high-elevation’ catchments. It is steeper, and widespreadadjacent zones of extensional and compressive flow create a rough surface with significantmicrotopography (as can be seen in Figure 2). Manual digitizing nonetheless reveals a well-developed and extensive supraglacial channel network. Previously, there has been limitedautomated delineation of supraglacial channels over heavily crevassed ice surfaces (e.g. Karl-strom and Yang, 2016; Smith et al., 2015), in spite of the fact that the coincidence of crevasses352.5. Discussionand supraglacial meltwater channels have been linked to high ice flow velocity and variabil-ity (Joughin et al., 2013; Wyatt and Sharp, 2015). Yang and Smith (2016) noted that multispectralmethods often fail to fully map supraglacial river networks at low elevations (below 1200 m)on the GrIS, although their methods were limited by the 15-m resolution of Landsat data.There is excellent SETSM coverage of many low-elevation sectors of the GrIS and its outletglaciers. This research therefore represents an attempt to identify a readily useable methodfor extracting surface meltwater channel characteristics from this morphologically complexice surface utilizing this new and expanding dataset. Flow routing is able to replicate the ma-jor structure and extent of a manually mapped low elevation catchment. Match rates above80% and 90% for fourth- and fifth-order channels, respectively, and above 75% for second-and third-order channels, but drop to below 65% for first-order channels. The mismatchesin the low-order channels arise primarily in the higher reaches of the catchment where smallmoulins are difficult to identify from the DEM and variable channel initiation points cannotbe adequately captured with a threshold initiation area. Although no overlapping MS datasetwas available for Catchment 7, it is reasonable to expect that multispectral methods would beof limited utility in such a rough, shadowed surface, and that flow routing provides a reason-able and physically objective method for delineating channels (Yang and Smith, 2016). Moulinidentification will likely be an ongoing challenge in rough ice, in particular, and the dangersof overestimating channel extent may be more significant here (Joughin et al., 2013).Interestingly, the digitized channel density and digitized initiation point density are lowerby 1.6 and 2.2 times, respectively, when compared to the MD network in Catchment 5. Ifboth slope and melt rate control channel initiation (Mantelli et al., 2015), higher channel den-sity would be expected in this lower elevation, higher slope catchment. Whereas Yang et al.(2016a) observed a decrease in channel density with increasing elevation, I do not observe thisbetween our low- and high-elevation catchments. It is compelling that MD initiation pointdensity in particular is lower in Catchment 7 than it is in the higher-elevation, digitized Catch-ment 5. This suggests that density is not low simply because there is more water drainingthrough the ice in this crevassed region, but rather because the points of initiation are morewidely spaced. This unexpected observation suggests that there are other controls on channelinitiation and density that need to be explored, or that the mechanisms of channel initiationdiffer between this rough surface ice and the generally crevasse-free ice surface of the high-elevation catchments.2.5.4 Limitations and ConsiderationsThere are a number of important limitations of this work that must be considered when in-terpreting the results. Firstly, manual steps incorporated into this workflow offset, to somedegree, the ease and rapidity of flow routing. These steps include manual identification ofcrevasses, moulins, lakes, and channel initiation densities. Previous work suggests approachesto automating parts of this process, including predicting moulin location (e.g. Clason et al.,362.6. Conclusions2012; Colgan and Steffen, 2009), automated crevasse detection (e.g. Jóhannesson et al., 2003), andautomated lake mapping (e.g. Fitzpatrick et al., 2014). The objective of the work presented hereis not to develop a fully automated approach to delineating the supraglacial fluvial network,but rather to employ flow routing under ideal conditions of known sink location to addressprevious reservations about flow routing over ice surfaces (Yang et al., 2015). The success offlow routing demonstrated here will support researchers in incorporating flow routing intotheir research with whatever modifications necessary.Secondly, as stated throughout, flow routing cannot directly map meltwater channel ex-tent. Although I compare it here to the MS dataset, I do so only to take advantage of a comple-mentary dataset that serves as a point of relative comparison. Given the ambiguity in channelinitiation area threshold, I suggest that researchers most interested in meltwater extent at agiven point in time are best served by multispectral methods. Although flow routing doesnot directly allow researchers to infer channel network extent, derived channels may be usefulto researchers interested in channel morphometrics that are dependent on fluvial processesacting over timescales of channel evolution. This has been done in supraglacial channels forcatchment characteristics and slope-area analysis (Karlstrom and Yang, 2016), and extensively interrestrial systems for network structure and self-similarity (see Rodriguiez-Iturbe and Rinaldo,1997), and scale independent morphometrics such as sinuosity and junction angles. Generalnetwork geometry, independent of initiation area, may also be of interest to researchers for arange of methodological purposes such as satellite image geo-registration (Yang et al., 2016b).2.6 ConclusionsIn this paper, I employ a newly available high-resolution (2 m) surface DEM of the GreenlandIce Sheet (Noh and Howat, 2015) to test flow routing as a method for delineating supraglacialchannels. I manually identify moulins and crevasse fields as sink locations and retain thesinks during flow routing. I compare the resultant channel networks to both manually dig-itized channels and a dataset of supraglacial channels delineated by Smith et al. (2015) from2 m resolution multispectral imagery. I find that flow routing is capable of replicating up to99% of the multispectrally derived channel networks when employed with a 0.01 km2 chan-nel initiation threshold, and can produce networks with significantly higher channel density.I confirm, using manually digitized channel networks, that when moulin locations are prop-erly determined, this higher channel density is real rather than an artifact of the flow routingmethodology. The accuracy of flow routing decreases with decreasing stream order, likely asa result of uncertainties regarding channel initiation area. Furthermore, I find that flow rout-ing over high-resolution DEMs provides a viable methodology to extract supraglacial channelnetworks at low elevations on the GrIS, where steep slopes and microtopography otherwisemake channel extraction challenging. Channel initiation threshold requires further investiga-tion, but researchers wishing to employ flow routing to delineate supraglacial channels maysafely employ quite low initiation thresholds (0.01 km2, in this particular study) during or after372.6. Conclusionspeak melt without risk of overestimating channel extent. Researchers will nonetheless needto pick channel initiation thresholds that reflect local surface conditions. Low-elevation catch-ments may have higher initiation thresholds because of the importance of microtopographyin determining channel heads, whereas high-elevation catchments may have higher initiationthresholds because of low melt rates.Supplementary Materials: The following are available online at, Fig-ure S1: MS and FR datasets for Catchments 1, 3 and 4.38Chapter 3Replicating moulin discharge estimatesreveals scale dependent error3.1 SummaryIn the midst of rapid proliferation of high resolution remotely sensed data, this brief con-tribution calls attention to the need for and utility of independent replication of results andthe importance of scale considerations in remote sensing analysis. In this chapter, I replicatemoulin discharge estimates published by Smith et al. (2015) using a different method and anindependent dataset. I show that the difference between the two approaches is highly scaledependent, and for small supraglacial rivers the difference can be several orders of magnitude.I suggest that the difference may be attributed to a scale incongruity between the spatial reso-lution of the two methods employed by Smith et al. (2015) – field measurements of hydraulicgeometry and remotely sensed river widths. These findings emphasize the need for scale-dependent measurements of error, as well as the methodological insight that can be gainedfrom replication studies.3.2 IntroductionSubglacial meltwater has been hypothesized to influence ice velocity in ice sheets and glaciers(e.g. Iken, 1981). Recent work has emphasized how the variability in surface meltwater supplyto basal drainage networks in the Greenland Ice Sheet (GrIS) temporally influences the devel-opment of subglacial hydrological pathways and thus the spatial and temporal relationshipbetween ice melt and velocity (e.g. Andrews et al., 2014; Banwell et al., 2016; Hoffman et al., 2011).Given the likely expansion of supraglacial meltwater networks under a warming climate (Lee-son et al., 2015; Poinar et al., 2015), there has been considerable interest in where surface melt-water reaches the bed through moulins, and how much water is supplied to these points ofentry (e.g. Koziol et al., 2017; Smith et al., 2015, 2017; Yang and Smith, 2016).Supraglacial hydrological systems span thousands of square kilometers and are difficult393.3. Methodsand expensive to study extensively in-situ (although it can and has been done for select catch-ments, see Andrews et al., 2014; Gleason et al., 2016; Smith et al., 2017). Given the importanceof moulin hydrology in determining ice response to surface melt, it is important to developmethods of accurately estimating moulin discharge remotely. However, moulin discharge es-timates derived from remotely sensed data can be problematic because, for example, largeproportions of the hydrological networks are below the spatial resolution of earth observationsatellites (e.g. Smith et al., 2015), and diurnal variability in melt, stream flow, and supraglacialchannel geometry occurs over time scales shorter that the return period of these satellites (e.g.Gleason et al., 2016).Independent replication of moulin discharge estimates can help to highlight limitationsin different contemporary approaches. This is essential for concretizing our presumed un-derstanding of these systems which, given the difficulty in accessing them, we may neversignificantly ground truth. More generally, the need for independent replication of findings isan increasingly articulated concern in the age of computational science (Peng, 2011). However,given the rapid pace of data collection, and increases in computational power and analyticalsophistication, fully independent replication is rarely undertaken (Leek and Peng, 2015).To this end, I compare two recent approaches to estimating moulin discharge on the GrIS.Each of these approaches makes use of a different high-resolution data product, and thusprovides an interesting opportunity to contrast two independent approaches to estimatingmoulin discharge. Throughout, I emphasize how insight into the limitations of a remotelysensed data product can be inferred from replication studies.3.3 MethodsThe first of the two moulin discharge estimation methods was developed by Smith et al. (2015),who combined field measurements of supraglacial hydraulic geometry with satellite remotesensing of supraglacial river widths to extrapolate a dataset of discharges at 523 moulins dur-ing peak daily melt (13:53 to 14:09 local time) between 18-23 July, 2012. River masks werederived from 2 m resolution WorldView-2 (WV2) images using multispectral methods, andchannel widths were subsequently extracted by spatially averaging measurements every 2 me-ters over 1 km reaches to reduce sensor resolution uncertainty. Empirical hydraulic geometryrelationships were developed through in-situ measurements of channel width and dischargeat 78 channel cross sections on the GrIS, and discharge was inferred from remotely sensedchannel widths using the resultant hydraulic geometry relations. This method was appliedonly to narrow, single thread channels up to 20 m wide. Moulin locations were visually deter-mined, assisted by automated river identification; abrupt channel termination was indicativeof moulin location. Instantaneous moulin discharge (m3/s) was estimated based on the up-stream contributing channel cross section. This dataset of instantaneous discharges from Smithet al. (2017) is hereafter referred to as the hydraulic geometry (HG) dataset.The second method for moulin discharge estimation makes use of the Polar Geospatial403.3. MethodsCenter’s 2-m resolution DEM which is stereo-photogrammetrically derived from 1-m reso-lution WorldView-1 images (Noh and Howat, 2015). It is therefore fully independent of WV2imagery used by Smith et al. (2015). Catchments were delineated by routing flow over the DEMto moulin locations identified by Smith et al. (2015). There was limited overlap between highquality, error free DEMs from the summer of 2012 and the Smith et al. (2015) dataset, constrain-ing the coverage of our dataset to the periphery of the Smith et al. (2015) dataset (see Figure3.1). Moulins at locations identified by Smith et al. (2015) were identified in the 2012 DEMand preserved as sinks during DEM filling and D8 flow direction calculation, all implementedthrough ArcHydro Tools (Maidment, 2002) in ArcMap 10.3.1 (ESRI, 2016). Further descriptionof the methods is provided in King et al. (2016). Catchment areas for each moulin were thus ex-tracted, and, similar to Yang and Smith (2016) and Karlstrom and Yang (2016), total daily moulindischarge was calculated by intersecting catchment area with daily runoff production for thecorresponding date from which Smith et al. (2015)’s dataset was derived (either the 18th, 21stor 23rd of July), as calculated by the RACMO2.3 climate model (Noël et al., 2015).These two datasets are not directly comparable. Whereas the RACMO2.3 surface massbalance (SMB) model used provides melt data at daily timescales, HG moulin discharge esti-mates are instantaneous and were inferred from imagery collected during the peak of the dailymelt cycle (1400h) which is not coincident with the timing of peak discharge. Additionally, ithas recently been shown that surface mass balance models such as RACMO2.3 may overesti-mate runoff production relative to moulin discharges (Smith et al., 2017). To compare the twomethods herein, the daily SMB-routed catchment discharge was converted to hourly melt hy-drographs, using a method proposed by Smith et al. (2017). Smith et al. (2017)’s contributionhighlights two main points. Firstly, SMB overestimates of daily discharge can be empiricallyscaled such that M = c ∗ M′, where M is SMB modeled melt (in dimensions of length overtime), c is an empirically derived adjustment coefficient, and M′ is effective (empirically ad-justed) melt. Secondly, they propose a synthetic unit hydrograph (SUH) curve whereby, for agiven hourly input of melt, the moulin discharge q (hr-1) at time t is equivalent to:q(t) = em[ ttp]m[e−m(ttp)]hp (3.1)Where m is an empirically derived equation shape factor, tp is the time to peak discharge(hr), and hp is the peak discharge (hr−1). Time to peak discharge is a function of the main stemlength (L, in km), distance between the point on the center flow line closest to the centroid of agiven catchment and the catchment’s moulin (Lc, in km) and an empirically defined coefficientCt such that tp = Ct(Lc · L)0.3. Similarly, hp is a function of an empirically defined coefficientCp and tp such that hp = Cp/tp. Using the above defined SUH, hourly moulin hydrographscan be estimated by convolving hourly time series of melt data with the q curve, such thatQ = M′ ∗ q where ∗ is a convolution operator (in this case, the NumPy ’convolve’ function).The resultant hourly melt rate at the moulin can be converted to volumetric discharge through413.3. Methods48°0'0"W49°0'0"W50°0'0"W67°40'0"N67°20'0"N67°20'0"N67°0'0"N67°0'0"N66°40'0"N66°40'0"NSmith et al. 2015 Rivers Overlapping Catchments^30 KmFigure 3.1: The main image provides a close up of the study area, which is located at thered star in the inset image of Greenland. Kangerlussuaq is located immediately tothe west. Grey polygons represent the supraglacial catchments that overlap betweenthe two datasets compared in this study. Note that because Smith et al. (2015) moulinestimates are based on the channel width immediately upstream of the moulin, itwas not necessary for the entirety of the Smith et al. (2015)-derived channels to over-lap with a flow routed catchment for that catchment to be relevant to this analysis.Background imagery is a 2015 Landsat-8 image dated August 25.multiplication with the catchment area.Their findings are incorporated into this study to derive the 1400h moulin discharge foreach catchment based on SMB-modeled melt. This dataset of daily SMB derived runoff dis-tributed hourly through the use of the SUH approach is henceforth referred to as the SMB/-SUH dataset. As hourly melt data is not available for 2012, daily RACMO2.3 runoff is con-verted into hourly data by fitting it to a gaussian distribution between 0600h and 0000h, peak-423.4. Resultsing at 1400h. L and Lc values and a unique SUH curve (Equation 3.1) were derived for eachflow-routing derived catchment. M′ was calculated using a c scaling factor of 0.65, based onSmith et al. (2017)’s observed difference between RACMO2.3 and in-situ measured dischargefor a specific catchment in July of 2015. This coefficient is designed to empirically adjust forwater retention processes as well as for errors inherent in the RACMO2.3 model. However,both snow and ice conditions as well as RACMO2.3 errors vary temporally, spatially, and alongan elevation gradient (Noël et al., 2015). Therefore, the use of this (and in fact any one) valueof c for the entire study study area in 2012 is unrealistic. However, it is clear that SMB runoffmust be scaled downwards to adjust for water retention processes, and given a dearth of fieldmeasurements, this value moves us towards better comparability between the datasets. Sim-ilarly, although likely not universally appropriate, the empirical values 2.1, 1.36, and 0.49 form, Ct, and Cp from Smith et al. (2017) are used.Not all Smith et al. (2015) moulins could be seen in the DEM layer, and therefore somedifferences between moulins were identified in the different datasets. In these cases, only onesink location was preserved and Smith et al. (2015) moulin discharges were summed within agiven DEM-derived catchment to compare moulin discharge values. Regardless, 82 out of 88catchments contained only one moulin as identified by either dataset, and all catchments werevisually inspected to ensure that flow-routing-identified catchment geometries and placementmirrored channel planform as identified by Smith et al. (2015) and that their rivers did not crossflow-routing-delineated catchment boundaries.3.4 ResultsHydraulic geometry-inferred discharges at 1400h ranged from 0.84 to 11.36 m3/s. In contrast,SMB/SUH 1400h moulin discharges based on RACMO2.3 spanned a similar range of values,ranging from 0.05 to 9.3 m3/s. Interestingly, the difference between the datasets was highlyscale dependent (Figure 3.2) such that the SMB/SUH estimates underestimated discharge rel-ative to the HG dataset in small catchments, and the highest agreement between the datasetswas observed in larger catchments. This relationship was significant (p-value < 0.001) basedon linear regression of log-transformed data (r2 = 0.69). Regression analysis suggests that thebest agreement between datasets was observed in catchments around 19 km2 in area.To further explore the discrepancy, SMB/SUH discharge estimates in a small catchmentwere compared with discharge estimates in the HG dataset. The catchment is 0.49 km2 andboth datasets largely agree on the geometry and extent of the hydrological system. The HGdataset estimates a 1400h discharge of 2.69 m3/s for this catchment. RACMO2.3 predicts totaldaily M’ of 14.2 mm/day (an M value of 22.0 mm/day), translating into a 1400h SMB/SUHdischarge of 0.17 m3/s (Figure 3.3), such that SMB/SUH discharge is just 6.32% of HG de-rived estimates. Peak daily discharge for this catchment from SMB/SUH measurements is0.21 m3/s, still less than 10% of the HG 1400h discharge.433.5. Discussion and Conclusion&DWFKPHQW$UHD(NPQHG/QSMB/SUH(%)5$&025XQRIIPPGD\Figure 3.2: Comparison between 1400h moulin discharges calculated by Smith et al. (2015)(HG) and as derived through SUH scaled daily RACMO2.3 data (SMB/SUH). Y-axisis the ratio of the two datasets multiplied by 100 (i.e. in %). The colors of the pointsare scaled relative to the RACMO2.3 predicted daily total runoff (mm/day).3.5 Discussion and ConclusionIn this study an instantaneous 1400h discharge for catchments on the GrIS is inferred by con-volving hourly SMB-modeled melt with catchment-specific SUH curves. These modeled dis-charges are compared to discharges inferred from hydraulic geometry relations and satelliteremote sensing of river widths on one of three different dates in July 2012. The latter methodprovides a snapshot of discharge in time, whereas the former attempts to broadly simulatedischarge throughout the day. With few tools available for the glaciological community tosimulate moulin discharge, this letter considers to what extent these methods might be com-pared and how this comparison informs us of their limitations.My findings indicate that deriving moulin discharge solely on SMB models produces sig-nificantly lower discharge estimates than those based on remote-sensing observation in smallbasins. The magnitude of the mismatch between the two datasets is substantial but is diffi-cult to interpret because of empirical assumptions. For example, the difference would changesubstantially for different values of c (the parameter that scales SMB-modeled melt, and a pa-rameter which certainly varies spatially and temporally). Additionally, the SMB model itselfdoes not account for snow and runoff processes (Irvine-Fynn et al., 2011; Smith et al., 2017) andmay produce runoff values that are lower than observed in the field. Therefore, the scale de-pendence of the difference is of more interest here than the absolute difference between the443.5. Discussion and Conclusion     7LPHKUP3 V0HOW+\HWRJUDSK60%68++\GURJUDSKP3 VFigure 3.3: Curves refer to the small 0.49 km2 test catchment. The dashed line is the dailymelt hyetograph, consisting of a gaussian curve between 0600 and 0000hr, peaking at1400hr (dashed vertical line). The solid line shows the moulin hydrograph, wherebypeak discharge is delayed relative to peak melt. Both the hyetograph and hydro-graph refer to the primary y-axis. The secondary (blue) axis is scaled according tothe difference between the SMB/SUH 1400h predictions and the HG instantaneousobservations. For example, the horizontal dashed line corresponds to the HG 1400hmeasurement of 2.69 m3/s on the right axis, and the inferred SUH/SMB value of0.17 m3/s on the left axis.datasets. There are two possibilities to explain this scale dependent mismatch. The first isthat this SMB/SUH approach underestimates instantaneous discharge in smaller rivers, andthe second is that the hydraulic geometry/multispectral approach leads to higher estimates ofdischarge in smaller rivers.Clearly, supraglacial channels were visible from satellite imagery around 1400h on thedates of acquisition, suggesting that water was flowing through channels at least several me-ters wide. The small discharges predicted by SMB/SUH approaches therefore do not seemconsistent with the observations of rivers in satellite imagery, suggesting this approach mayindeed be underestimating discharge in small rivers. SMB/SUH based approaches might in-troduce scale-dependent error in a number of ways. For example, the use of empirical coef-ficients derived from just one field study for scaling melt and building SUH curves herein isalmost certainly unrealistic, and substantially more field work will facilitate building empiri-cal datasets for hydrograph response to different catchment properties. The use of non-scaledependent empirical values as well as the gaussian-distributed melt hyetograph may fail toaccurately reproduce the daily distribution of the moulin hydrograph. Additionally and im-portantly, the discharge value is highly dependent on basin area, for which smaller errors in453.5. Discussion and Conclusiondelineation are of more relative importance in smaller basins than in larger ones.However, given the high melt values that would be needed to produce the discharges re-ported by Smith et al. (2015) in the small sample catchment above, it is reasonable to concludethat some of the discrepancy between datasets arises because a reliance on hydraulic geom-etry introduces several scale-dependent methodological limitations. Firstly, their hydraulicgeometry relation of Q = 0.10w1.84, where Q is discharge and w is width, is highly sensitiveto measurements of width. Given the 2-m resolution WV2 imagery used in their analysis,this uncertainty is higher in smaller channels, and may explain the tendency to overestimatedischarge in smaller channels, similar to conclusions by Bjerklie et al. (2005).A second potential scale-dependent limitation in using hydraulic geometry to estimatemoulin discharge arises from extrapolation of hydraulic geometry across river sizes. Obser-vations suggest that supraglacial channels are often deeply incised, and that the difference indischarge between channels of different widths is not well described by a width/dischargerelationship (Gleason et al., 2016; Marston, 1983), although Yang et al. (2016a) observed thatwidth was highly sensitive to discharge on the GrIS. Others have shown that velocity increasesrapidly in response to increasing discharge, often leading to little change in width to depth ra-tios (Gleason et al., 2016; Knighton, 1981). Clearly, there is substantial variability in the nature ofhydraulic geometry in supraglacial channels.If an insensitivity between width and discharge contributes to the discrepancy between es-timates observed here, one would expect the discrepancy to be dependent on river discharge.Indeed, Figure 3.2 suggests that on lower melt days, HG discharge estimates are higher rel-ative to RACMO2.3 estimates compared to higher melt days. It is possible that RACMO2.3underestimates melt on those days, (and) or it is possible that the w/Q hydraulic geometryrelationship is such that the width of channels does not vary substantially between high andlow discharges, and best approximates discharge on high melt days.Finally, the temporal sensitivity of these strongly diurnal systems fundamentally compli-cates a one to one comparison of modeled runoff and observed instantaneous discharge. Asteep hydrograph requires an exceptionally well parameterized SUH to accurately time thedelivery of meltwater to a moulin. This is likely particularly true in small basins with shorterpeak times. The high instantaneous discharges observed with satellite remote sensing of theriver widths suggest that the SUH parameterization used herein does not accurately capturethe rapidly rising hydrograph of these smaller systems.Accurate, scale-independent methods of moulin discharge estimation are crucial for pre-dicting the impacts of surface melt on ice sheet behaviour. Whichever method is employed,issues of temporal and spatial scale dependence are likely to be omnipresent, and researchersare advised to consider ways of developing scale dependent estimates of error and attempt-ing independent replication of their results. Further, researchers should consider the waysin which scale-dependent methodological issues become internalized into the outputs of ourresearch. For example, previous work has shown that catchment characteristics (including463.5. Discussion and Conclusionsize) are non-uniformly distributed on the ice sheet (e.g. Yang et al., 2016a), and that the spa-tial distribution of moulins influences sub-glacial channel development (Banwell et al., 2016).Smaller catchments are generally more plentiful, and particularly so at lower elevations whereseasonal melt begins the earliest (e.g. Yang and Smith, 2016; Yang et al., 2016a). Over- or under-estimating the volume and timing of their meltwater contributions to basal drainage systemsmay have significant implications earlier in the season when the link between meltwater andice velocity is strongest (Andrews et al., 2014).SMB/SUH based approaches are promising and practical for deriving moulin dischargeestimates across surface catchments. However, given the heavily empirical nature of SUHderivation and its apparent sensitivity to scale in small catchments, substantially more theo-retical and field-based work is needed before any approach can be considered to be reliableover large areas with variable catchment geometries, slopes, and ice and snow conditions.Meanwhile, as illustrated in this letter, independent replication, although seldom undertakenand institutionally disincentivized (see Nosek et al., 2015), can be an accessible tool for identi-fying a scale incongruity.47Chapter 4The spatial variability andcharacteristics of moulins in southwestGreenland4.1 SummaryMoulins convey surface water to the bed of glaciers and ice sheets, where subglacial hydrolog-ical processes can influence ice velocity. Knowing where and how moulins form is thereforea critical component of projecting future ice sheet dynamics under enhanced ice melt condi-tions. However, we have a very limited understanding of the variable conditions under whichmoulin formation occurs. Moulin formation in ice sheet hydrology models is driven by a setof over-simplifying assumptions that are not consistent with empirical observations of moulindistribution. This chapter explores the ways in which ice sheet and hydrological conditionscollectively produce the spatial arrangements of moulins in a study area in south west Green-land. I adopt an explicitly empirical and spatial approach, and observe that moulin formationis most common in low velocity, high strain rate ice regimes with adverse bed slopes. I observethat contemporary ice sheet hydrology models do not capture the range of conditions underwhich most moulins form. Additionally, I observe that moulins proximal to in-situ drainedlakes are predominantly located down-ice of the lake, and many of these moulins appear tohave no obvious formative mechanism. I interpret these findings to support the recent theorythat lake hydrofracture is a transient process related to temporary ice flexure in response torapid moulin drainage events.4.2 IntroductionMoulins are near vertical pathways through which surface meltwater can rapidly reach thebed of a glacier or ice sheet. There, if water supply overwhelms subglacial channel capacity,484.2. Introductionmeltwater can increase basal water pressure and facilitate ice sliding (e.g. Iken, 1981; van deWal et al., 2008; Zwally et al., 2002). As water enlarges subglacial channels over time (e.g. overthe course of a melt season), channel capacity increases, reducing water pressure, and icevelocity is decoupled from surface meltwater supply (e.g. Schoof , 2010). The density and dis-tribution of moulins and their hydrology throughout the melt season are important controlson the temporal evolution of the subglacial hydrological network and thus have the potentialto strongly influence seasonal ice velocity dynamics (Banwell et al., 2016; Hoffman et al., 2011).Understanding where, when and how moulins form is important for constraining the complexnon-linear relationship between ice sheet dynamics and surface melt, as well as for projectingthis relationship into a warming future (Leeson et al., 2015; Poinar et al., 2015).On the GrIS, meltwater reaches the bed through a number of mechanisms that have beenvariously investigated. Moulin formation through lake hydrofracture, in which supraglaciallakes drain rapidly (on timescales of hours) through fractures at their base, has received themost intense empirical and theoretical investigation. This focus is motivated by the well doc-umented association between rapid lake drainage and short term ice acceleration (Das et al.,2008; Hoffman et al., 2011; Joughin et al., 2008, 2013). Ice sheet hydrology models that simu-late moulin location typically do so using a combination of surface energy budget-driven meltproduction, flow routing, and lake hydrofracture initiated by a threshold volume of water(Banwell et al., 2016, 2013; Clason et al., 2015). Advancements in model sophistication haveincorporated moulin drainage through crevasse hydrofracture when routed flow interceptsthreshold surface strain rates (Poinar et al., 2015) or stress values (Clason et al., 2015; Koziol et al.,2017).These models all have several assumptions in common. Firstly, they all assume hydrofrac-ture is the dominant moulin formation mechanism. Secondly, they assume that lake drainageand crevasse hydrofracture comprise the primary mechanism of moulin formation at high ele-vations. However, mounting evidence suggests that a more nuanced understanding of moulinformation and distribution is required. Modeled and empirical studies in south west Green-land show that only a fraction of moulins are located in drained lake basins (Catania et al.,2008; Koziol et al., 2017; Smith et al., 2015; Yang and Smith, 2016). Further, only a minority ofsupraglacial rivers terminate directly in crevasse fields (Smith et al., 2015). The majority ofmoulins appear to terminate away from either of these obvious formation mechanisms (Cata-nia et al., 2008; Hoffman et al., 2018; Koziol et al., 2017; Smith et al., 2015), possibly created bylocalized and transient ice-fracture caused by ice speedup following the drainage of nearbysupraglacial lakes (Hoffman et al., 2018; Stevens et al., 2015). Additionally, there is consistent ev-idence of high elevation moulins outside of lake basins (Lampkin and VanderBerg, 2014; Poinaret al., 2015; Smith et al., 2015; Yang and Smith, 2016) despite a decreasing frequency of lakedrainage events with increasing elevation (Morriss et al., 2013; Poinar et al., 2015).These empirical studies largely focus on quantifying moulin and contributing catchmentcharacteristics in space, with a primary focus on elevation (Smith et al., 2015; Yang and Smith,494.3. Study Site and Methods2016; Yang et al., 2016a). As such, observed spatial trends are largely divorced from the un-derlying processes that control moulin characteristics and distribution. Although previouswork illuminates a complex and spatially variable supraglacial hydrological system, we stillhave little understanding of the mechanistic controls on moulin variability. Empirical predic-tive models of moulin location have had limited and mixed success (Clason et al., 2012; Colganand Steffen, 2009; Phillips et al., 2011), and describing or predicting moulin location remainsa challenge for building a comprehensive understanding of ice sheet hydrology. In light ofthe disconnect between the large variability observed in empirical studies and the embeddedsimplifying assumptions of modeling approaches, there is a need for a more holistic under-standing of where moulins occur, and why.In this study, I address this discrepancy by relating observed moulin locations to the spa-tial distribution of ice sheet characteristics and dynamics. I adopt an empirical and explicitlyspatial approach, looking for spatial trends in how these conditions vary with respect to oneanother and with respect to moulin type and location. I test whether model assumptions oflocal moulin conditions reflect observed moulin distributions. This contribution offers neededempirical insight towards building a better understanding of the variability and complexity inspatial configurations and distributions of moulins.4.3 Study Site and Methods4.3.1 Study Area and datasetsThe study area covers 2700 km2 on the southwest margin of the ice sheet, approximately 150km south east of Nuuk (Figure 4.1). It extends from 584 to 1761 m a.s.l. in elevation, with amaximum ice flow velocity of 245 m/yr (Joughin et al., 2010), and includes portions of threeoutlet glaciers that drain from two large bedrock troughs. The largest of the outlet glaciersis simply named Sermeq (Bjørk et al., 2015). Ice thickness is typically less than 1 km but is inexcess of 1.5 km in places (Morlighem et al., 2017). Maps of ice thickness, velocity, strain rate,and surface slope can be accessed in Appendix A.This study area was chosen in part because of a good coincidence in high quality, highresolution datasets in this area, but also because it differs from areas that have previouslybeen the focus of supraglacial studies. Previous work has focused extensively on Russell andIsunguata Glaciers (e.g. Clason et al., 2015; Fitzpatrick et al., 2014; Palmer et al., 2011; Smith et al.,2015, 2017; Sole et al., 2013; Sundal et al., 2009; van de Wal et al., 2008; Yang et al., 2016a) andthe Jakobshavn Isbrae\Sermeq Avannarleq\Paakitsoq region (e.g. Andrews et al., 2014; Banwellet al., 2016; Hoffman et al., 2011; Joughin et al., 2013; Koziol et al., 2017; Lampkin and VanderBerg,2014; McGrath et al., 2011; Morriss et al., 2013; Phillips et al., 2011). Very little work is available inthe southernmost part of the GrIS, despite a widespread presence of supraglacial hydrologicalfeatures and an extensive ablation zone (Nghiem et al., 2012; Noël et al., 2015). This presentwork, then, expands the geographic purview of our quantification of supraglacial features.504.3. Study Site and MethodsThe study area extent is fixed by the availability of DEM data in this area.49°W50°W51°W64°N63°40'NStudy Area^Nuuk20 Km¹Figure 4.1: Study area. The town of Nuuk is located approximately 150 km to the northwest. Background imagery is from Landsat 8, acquired on July 19, 2015. Inset datais the outline of the ice sheet, with the study area marked as a red star.Four different data products were used in this analysis and are described below: an icesurface DEM, ice velocity data, bed DEM, and runoff production data from a climate model.Ice surface topography data were available through a 2 m resolution DEM of the ice sur-face (ArcticDEM, available from This DEM isstereogrammetrically derived from overlapping WorldView-1 satellite imagery, with reportederrors of 3.8 and 2.0 m in the horizontal and vertical, respectively (Noh and Howat, 2015). De-spite some seams and blunders, it has been used to delineate supraglacial channels with highaccuracy (King et al., 2016). The DEM tiles used in this study are derived from 2015 imageryacquired on, in order of increasing elevation, September 2nd, August 25th, and July 13th.Ice surface DEMs were used to visually identify moulins and to produce moulin catchmentsthrough D8 flow routing using the ArcHydro Toolbox in ArcMap 10.3.1 (Maidment, 2002), asdescribed in Chapter 2 and King et al. (2016). Examples of moulin identification are shown inAppendix B.Ice velocity data, used to derive surface strain rates, are available in the form of 200 m res-olution gridded ice velocity rasters from 2015-2016 (MEaSUREs Greenland Ice Sheet Velocity514.3. Study Site and MethodsMap from InSAR data), available from The velocitydata refer to winter velocities, broadly defined from September 1st to May 31st. Velocity is de-rived from speckle tracking methods with Synthetic Aperture Radar data, acquired in this caselargely from Sentinel-1A (Joughin, 2002). Average x- and y-velocities errors are of 2.5 m/yr and6.3 m/yr, respectively. Following Poinar et al. (2015), I derived principal surface strain rates (e˙,in yr−1) for each pixel using the MEaSUREs surface velocity data as:e˙ =12(∂u∂x+∂v∂y)+12√(∂u∂x− ∂v∂y)2+(∂u∂y+∂v∂x)2(4.1)where u and v are velocity in the x and y (east-west and north-south) directions, respec-tively.Gridded velocity products were filtered using a 3 x 3 pixel (600 m) gaussian blur (smooth-ing) filter prior to deriving strain rates. Filter window sizes up to 3 km were tested (AppendixC). The 3 x 3 filter was chosen as it preserves the many linear high strain rate features (Ap-pendix B) and the 600 m window size approximates one average ice thickness in the studysite. Similar to (Poinar et al., 2015), velocity gradients were calculated using centered differ-ences, implemented with NumPy’s Gradient function.Bed elevation and ice thickness data were acquired from 150 m resolution BedMachineV3(Morlighem et al., 2017) available from This data productis derived from a range of radar surveys carried out between 1993 and 2015 and generalizedusing a mass conservation approach. As per the gridded BedMarchineV3 error data, the errorin bed and ice thickness in the study area is an average of 106 m.Daily meltwater runoff production for the study area is provided by the 100 m resolutiongridded daily RACMO2.3 climate model (Noël et al., 2015). Runoff production in units of mm/-day was translated into moulin discharge in m3/day by intersecting the daily modeled runoffproduction with catchment area, as described in Chapter 3. The sum of this daily catchmentdischarge for the year provides annual discharge.Crevasse fields were digitized based on visual inspection concurrently with moulin iden-tification. Crevasses were distinguishable from other linear surface features (e.g. channels) bytheir geometric straightness and the repeated pattern of parallel neighbouring crevasses. AllDEM sinks within visually identified crevasse fields were retained so that water draining tothem remained there. The crevasse dataset was not intended to be authoritative; rather its pur-pose was to 1) limit the extent of moulin catchments by capturing water draining into crevassefields, and 2) to provide context for the relationship between strain rate and crevassing in thestudy area.Much like Joughin et al. (2013), I found moulin identification more difficult at lower eleva-tions where there is an abundance microtopographic features such as crevasses. Additionally,our moulin dataset almost certainly excludes smaller moulin features that are below the reso-lution of what is distinguishable from the DEM and thus overestimates catchment area (partic-524.3. Study Site and Methodsularly higher up in the drainage network, see King et al., 2016). Therefore our dataset does notrepresent a complete inventory of moulins in the study area. However, it does reproduce thegeneral drainage patterns arising from moulins draining evident stream channels apparent inthe 2 m resolution DEM.4.3.2 Moulin classificationBased on the literature, I identified four moulin types: moulins formed in association with vis-ible crevassing, moulins formed through lake-drainage, moulins formed due to the presenceof a sediment or bedrock protrusion, and moulins for which no readily discernible formationmechanism is clear.I defined crevasse-proximal moulins as those occurring within a horizontal distance of adigitized crevasse field equivalent to one local ice thickness. To identify lake-drainage asso-ciated moulins, I visually inspected each moulin location in Landsat-8 imagery spanning thesummer of 2015 (dated June 24th, July 19th, and August 4th), and classified a moulin as lake-drainage associated if water was clearly accumulated in that spot at one point in the summer,and gone or substantially reduced at a later point in time (Figure 4.2). It is therefore difficultto specify the exact timing of the lake drainage event; the lake may have drained through sud-den hydrofracture, or it may have drained through moulin development at the end of a shortincised drainage channel. Moraine-proximal moulins refer to moulins formed at the periph-ery of bedrock or sediment (e.g. moraine) features. I have included these moulins as seperatemoulin types, as previous work has shown that moulins can form through the exploitationof permeable pathways in the ice caused by embedded sediment such as moraines (Gulleyet al., 2009a). Finally, moulins not classified into these other three categories were classified as’other’ type moulins.534.3. Study Site and Methods0.5 KmA. B. C.MoulinsFigure 4.2: These images illustrate how lake-drainage associated moulins were identified.From left to right, (A) shows moulins as visible in the DEM-derived curvature layer.Figure (B) shows a lake visible in Landsat-8 imagery from July 19th, that has mostlydrained by (C) August 4th, 2015.4.3.3 Lake delineationSupraglacial meltwater was identified from the aforementioned Landsat 8 imagery using thespectral band combination proposed by Morriss et al. (2013) for the identification of supraglaciallakes. Morriss et al. (2013) used the normalized difference lake index (NDLI), such that:NDLI =(Rr − RnirRnir + Rr)(4.2)where Rr and Rnir are the reflectance values in the red and near infrared spectra, respec-tively. Given the 15 m resolution of the imagery and the relatively low contrast between waterand ice during the melt season (Yang and Smith, 2016), this multispectral approach was notsensitive enough to delineate supraglacial channels. Water bodies that were identified by theanalysis were large and not elongate, and were therefore interpreted as surface water storage.To remove spuriously identified pixels and small water bodies (such as widenings in chan-nels), only surface water bodies exceeding 10,000 m2 were preserved for analysis. Each waterbody was inspected over the course of the summer in the three Landsat 8 images. Lakes thatwere present at one point in time and then gone or significantly reduced at a later point intime were assigned as ’draining’ lakes. Draining lakes that contained moulins were assigned544.4. Resultsas ’moulin-drained lakes’. In contrast, lakes that do not drain in-situ through moulins wereassumed to be ’pathway’ lakes (Yang and Smith, 2016), draining to downstream moulins.4.4 Results4.4.1 Overview of spatial moulin characteristicsIn total, the moulin dataset includes 327 moulins. As observed in previous studies (e.g. Smithet al., 2015; Yang and Smith, 2016), moulin density was highest at intermediate elevations, de-clining from peak density beyond 1300 m and dropping below the hypsometric distributionof the study region at around 1400 m (Figure 4.3). Only 32 moulins were present at elevationsabove 1400 m. The highest elevation moulin was observed at 1589 m, although I note that ourstudy area does not extend all the way to the upper limit of the ablation zone. Landsat imagerysuggests the presence of moulins more than 25 km inland from the upper extent of our studyarea. Moulin catchment area ranged in size from 0.32 to 58.87 km2, with a median basin sizeof 2.21 km2. Although there was a generalized trend of increasing maximum catchment sizewith elevation, both large and small catchments were observed throughout the study area.       (OHYDWLRQPDVO)UHTXHQF\,FH6XUIDFH(OHYDWLRQ 0RXOLQ(OHYDWLRQ%DVLQ$UHDNP2 Figure 4.3: Moulin distribution in different elevation ranges (blue line) relative to the hyp-sometric distribution of ice surface in the study area (grey line). The blue pointsshow the distribution of moulin basin sizes (right y-axis) relative to moulin eleva-tion.The decline in moulin density above 1300 m and then abruptly again beyond 1400 m gen-erally coincides with declining ice surface gradient and melt rates (Figure 4.4). Throughoutthe study region, median surface slopes decline with increasing elevation, but high surfaceslopes abruptly disappear above 1500 m. Annual surface melt rates decline with increasingelevation, with the decline becoming less pronounced above 1300 m where melt rates drop554.4. Resultsbelow 1 m/yr.Figure 4.4: A) Slope and B) annual melt data from the study region relative to elevation.The 95th percentile of slope and the median of melt within elevation bins are alsoplotted. The vertical dashed line marks 1300 m a.s.l.Moulins were not uniformly distributed. Join counts statistics, which give a measure ofspatial autocorrelation in binary data (e.g. moulin presence/absence), suggest statistically sig-nificant (p = 0.001) positive spatial autocorrelation between moulins in the study area. Moulinstherefore group together; where there is one moulin, there are likely to be more nearby. Thisspatial grouping does not appear to simply reflect the elevation gradient. Figure 4.5 shows aplot of the kernel density function of moulins in the study area (the probability of a moulin per305 m x 305 m pixel). There is a visually distinct spatial cluster of moulins in the north easternquadrant at low elevation, as well as in the southern half of the study area at mid-elevation(Figure 4.5). The remaining low-elevation parts of the study area have low moulin concen-trations (Figure 4.5). Collectively, the clustered and spatially autocorrelated distribution ofmoulins suggests localized controls on moulin formation.564.4. Results  'LVWDQFHNP       H Figure 4.5: Moulin kernel probability densities based on point distributions throughoutthe study area. Kernel density functions are similar to smoothed histogram, wherethe color scale shown here is the probability of finding a moulin, per unit area. Ele-vation contours are shown as black lines, and moulins are represented by blue dots.To develop preliminary insight into what might control this spatial distribution of moulins,I compared ice surface velocity, ice thickness, principal surface strain rate, and ice surfaceslope data at moulin locations to regional characteristics throughout the study area (Figure4.6). Considerable variability is present throughout the dataset, however moulins form atstatistically significantly (p < 0.05, as per Kruskal Wallis one way analysis of variance) lowervelocities than are present throughout the study area overall (medians of 37.1 vs. 57.1 m/yr).They do not appear to differ significantly in terms of the distribution of strain rates (mediansof 0.007 vs. 0.010 yr-1), ice thicknesses (medians of 493 vs. 616 m), or surface slopes (1.4 vs.574.4. Results2.1) in the study area. Moulins appear to form primarily in extensional ice environments; 86%of moulins formed under positive strain rates. The lowest surface strain rate at a moulin was-0.030 yr-1, and the highest 0.06yr-1.  HiceP)UHTXHQF\$  ViceP\U)UHTXHQF\%   \HDU 1)UHTXHQF\&   6ORSHPP)UHTXHQF\'0RXOLQV8QFUHYDVVHG,FH&UHYDVVHG,FH$OO,FHFigure 4.6: Frequency distributions of a) ice thickness, b) velocity, c) principal strain rateand d) slope, throughout uncrevassed ice, crevassed ice, all ice (including bothcrevassed and uncrevassed) and moulin locations. The bin sizes and ranges areconsistent across the series.Logistic multiple regression of binary moulin presence/absence gridded data relative to icethickness, slope, velocity and strain rate produces an R2 of 0.03 and thus offers little promise interms of understanding moulin location. Multiple linear regression using distance-to-moulinsas a dependent variable and the aforementioned independent variables produces an R2 of just0.04, and again offers little insight into the distribution of moulins.That moulins form at statistically significantly lower velocities is reflected in their spatialdistribution. Moulin clustering appears to be concentrated in areas of low velocity. Figure4.7 (A) shows moulin distribution relative to three zones of velocity (low, medium and high).These velocity zones are distinguished using Jenks Natural Breaks histogram optimization, adata clustering method designed specifically for mapping which minimizes variance within agrouping whilst maximizing variance between groupings. Whereas 47, 33, and 20% of the ice584.4. Resultssurface is in the lowest, middle, and highest velocity zones, 65, 22 and 13% of moulins are inthese zones, respectively.This pattern of higher moulin density with lower ice velocity does not appear to reflect asimple distribution of ice velocity relative to elevation (Figure 4.7 B). Although moulin densityis consistently low at all elevations in medium and high ice velocity zones, moulin density ishigher in low velocity zones throughout the elevation range, and particularly so for the middleelevation ranges, reflecting the higher moulin density at mid-elevations, as previously noted.  'LVWDQFHNP$9HORFLW\/RZ 0HGLXP +LJKP P(OHYDWLRQ%LQ0RXOLQ'HQVLW\NP2 %9HORFLW\=RQH/RZ0HGLXP+LJKFigure 4.7: A) Moulin location overlain on three zones of velocity, whereby low velocityrefers to 0-54 m/yr, medium to 42-116 m/yr, and high to >116 m/yr. B) Moulindensity within each velocity zone, as a function of elevation bins.Interestingly, moulin density appears highest in areas with high localized strain rates withinregionally low velocity zones (Figure 4.8), despite strain rate distributions being statisticallysignificantly (p < 0.001) higher in high velocity zones relative to low velocity zones.594.4. Results/RZ 0HGLXP +LJK9HORFLW\&ODVV0RXOLQ'HQVLW\NP2  UDQJH\U1!Figure 4.8: Moulin density within strain rate bins, within velocity classes. The three strainrate bins divide strain rate according to compressive, (< 0 yr−1), extensional butbelow the crevassing threshold as per (Poinar et al., 2015) (0 - 0.005 yr−1) and abovethe crevassing threshold (> 0.005 yr−1).4.4.2 Moulin TypesOf the 327 identified moulins, 36 (11%) appeared to be associated with a moraine or a nunatak(collectively termed ’moraine-proximal’), 92 (28%) were within one local ice thickness of acrevasse (crevasse-proximal moulins), 48 (15%) were associated with observed surface waterdrainage (lake-drainage), and the remaining 151 (46%) were of an ’other’ type. In general,there were few significant differences in ice or hydrological characteristics between moulintypes. According to Kruskal-Wallis tests, ice thickness, strain rates and ice velocity were sig-nificantly (p< 0.05) higher at crevasse-type moulins than at any of the other moulin types, andtotal and peak annual discharge was significantly lower. The other three moulin types did notvary significantly in relation to each other (Figure 4.9).604.4. Results&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU /DNH'UDLQDJHHiceP$&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU /DNH'UDLQDJH\U1 %&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU /DNH'UDLQDJHV iceP\U&&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU /DNH'UDLQDJH(OHYDWLRQP'&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU /DNH'UDLQDJH4aP3 \U(&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU /DNH'UDLQDJH4maxP3 /d)Figure 4.9: Characteristics of (A) ice thickness, (B) surface strain rate, (C) ice surface ve-locity, (D) ice surface elevation, (E) total annual discharge in 2015 and (F) maximumdaily discharge in 2015, for each moulin type (crevasse-proximal, n = 92; moraine-proximal, n = 36; other-type, n = 151; lake-drainage, n = 48. Upper and lower boundson boxes represent the 75th and 25th percentiles, where the upper and lower extentsof the error bars represent the 95th and 5th percentiles. Outliers are shown by points.The spatial clustering of moulins relative to velocity zones becomes particularly apparentwhen moulins are broken down by type (Figure 4.10). ’Other’-type moulins visually appearto dominate in low velocity zones, crevasse-proximal moulins appear dominant at highervelocities, lake-associated moulins appear common in low and medium velocity zones, andmoraine-proximal are, implicitly, restricted to areas adjacent to moraine or bedrock features.614.4. Results  'LVWDQFHNP&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO/DNH'UDLQDJH2WKHUFigure 4.10: The spatial distribution of moulin types relative to velocity zone. Moulinsare represented by circles, where the circle color corresponds to moulin type. Inthe background, black, grey and white coloring represents low, medium and highvelocity ice, respectively.Indeed, moraine-proximal, lake-associated and ’other’ type moulins are all quantitativelymore common in low-velocity zones, whereas only crevasse-proximal moulins appear morecommon in high-velocity zones (Figure 4.11). The density of ’other’ type moulins appears to beparticularly high in low-velocity zones, whereas lake-drainage associated moulins appear tobe more equally distributed between low and medium velocity zones. Conversely, all moulintypes appear more dense at higher positive (extensional) strain rates. No moraine-proximal624.4. Resultsmoulins are found in compressive ice regimes. Only ’other’ type moulins appear to exist atany significant spatial density in compressive ice./RZ 0HGLXP +LJK9HORFLW\5HJLPH0RXOLQ'HQVLW\NP2 &UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO/DNH$VVRFLDWHG2WKHU\U 1 \U 1 !\U 15DQJH0RXOLQ'HQVLW\NP2 Figure 4.11: Moulin density for each moulin type, relative to the three previously outlinedvelocity zones (see Figure 4.7).4.4.3 Influence of lakesLake drainage events have been suggested as important mechanisms of moulin formation,such that stochastic lake drainage generates transient high local strain rates through localizedflexure of the ice (Hoffman et al., 2018; Stevens et al., 2015). To evaluate the relationship be-tween moulins and lakes, I compiled two datasets. The first is based on visual inspectionsof Landsat-8 imagery around each moulin to identify which moulins appeared to be associ-ated with surface water accumulation and subsequent drainage (lake-drainage moulins). Thesecond dataset was built through automatic lake classification using multispectral analysesof the same Landsat images. Each of the lakes was then inspected and assigned as ’moulin-draining’, ’non-moulin draining’, and ’non-draining’. Moulin-draining lakes were assignedwhen moulins were located within a draining lake basin or immediately adjacent to the basin.634.4. ResultsMoulins not related to observed lake drainage (non lake-drainage-associated) were locateda median distance of 3.0 km from the nearest surface water body (as determined by multispec-tral water delineation), with 32% located more than 4 km from any surface water body (Figure4.12). ’Other’-type moulins tended to be located in closer proximity to lakes draining throughmoulins, with a median distance of 2.5 km, relative to median distances of 3.2 and 3.1 forcrevasse-proximal and moraine-proximal moulins, respectively.     'LVWDQFHWR1HDUHVW/DNHP&XPXODWLYH)UHTXHQF\&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHUFigure 4.12: The cumulative frequency of the distance between a non-lake-drainage-associated moulin and its closest draining lake, where ’draining lake’ refers specif-ically to lakes draining through moulins.Of the 120 water bodies identified by multispectral analysis of Landsat imagery, drainageevents were observed in 104 lakes. Of these, 44% drained in-situ through moulins locatedin the lake basin, whereas the others acted as pathway lakes, draining to terminal moulinsdownstream, presumably following overflowing or drainage channel incision (e.g. Kingslakeand Sole, 2015). Of the lakes draining through moulins in situ, 54% had moulins external tothe lake basin within a distance of 2 local ice thicknesses. In contrast, only 34% of pathwaydraining lakes and 38% of non-draining lakes had moulins within 2 local ice thicknesses. Ofthe lakes with nearby moulins (within two ice thicknesses), nearby moulins appeared to belocated predominantly downstream, with 65-68% of nearby moulins downstream of lakes,regardless of lake drainage type (Figure 4.13).Although only a limited number of lakes drained through in-basin moulins, these lakesappeared to be in closer-proximity to downstream crevasse-proximal moulins (18% of down-stream moulins), a moulin type which did not frequently appear downstream of non-terminallakes (7% of downstream moulins).644.4. Results'LVWDQFH)URP/DNHHice),FH)ORZ$,FH)ORZ%    'LVWDQFH)URP/DNHHice)'LVWDQFH)URP/DNHHice),FH)ORZ&HiceHiceHiceHice&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO/DNH'UDLQDJH2WKHUFigure 4.13: Moulin positions within 2 local ice thicknesses (in the x and y) are plottedrelative to lake location for (A) lakes that drain in-situ through moulins, (B) lakesthat drain to downstream locations, and (C) lakes that do not drain (terminal lakes).The position of moulins has been oriented relative to ice flow, rather than throughspatial location. Distances in the x and y are in units of local ice thicknesses, as il-lustrated by the concentric circles. Note that the lake-basin moulins through whichthe lakes drain are not shown.If I collectively consider all non-draining lakes as well as lakes that overtop and feed down-stream channels as ’non draining lakes’ (in contrast to in-situ draining lakes), I can questionwhy some lakes appear to hydrofracture whereas others (and indeed most) do not. Figure 4.14shows the principal components of PCA analysis on the ice thickness, lake area, local strainrate, and distance to the nearest moulin for in-situ (moulin-drained) and non moulin-drained(undrained) lakes. PCA analysis allows us to visualilze the relationship between data thatvaries in multiple dimensions (in this case, 4) in a 2-dimensional graph by plotting the datapoints along the two primary axes of variance (PC-1 and PC-2), which in this case collectively654.4. Resultsexplain 67% of the variance in the dataset. The points (scores) represent where each moulinsits in these reduced dimensions, and the magnitude and direction of the arrows (PCA coeffi-cients, i.e. loadings) represent how the original variables (e.g. lake area) contribute to each ofthe Principal Component axes. Lake area, ice thickness and distance to the nearest moulin arestrongly related to PC-1, and strain rate is associated with the secondary primary axis of vari-ability (PC-2). The distribution of the scores relative to the coefficients suggests that drainedand undrained lakes appear to be distinguished primarily with regards to distance from theclosest moulin.    3&3&HiceALakedmoulin8QGUDLQHG0RXOLQ'UDLQHGFigure 4.14: The first two principal components of a PCA analysis on local characteristicsaround draining and non-draining lakes are shown (PC-1 and PC-2) with regards toice thickness Hice, distance to the nearest moulin dmoulin, lake area ALake, and strainrate e˙. Arrows represent the direction and magnitude of the coefficients for eachindepedent variable (i.e. the PCA loads), and the points represent the transformedpoint values for each of the lakes (i.e. the PCA scores).As per Figure 4.15, non-draining lakes form, on average, at only slightly higher ice thick-nesses and at marginally smaller lake areas. Both moulin-drained and non-draining lakesoccur almost exclusively in extensional (positive) strain rates, however moulin-drained lakesoccur at higher (although not statistically significantly so) strain rates. Moulin-drained lakes664.4. Resultsoccur with statistically significantly (p < 0.001) closer neighbouring moulins.HiceP$A lakeP2 %1RQ'UDLQLQJ 0RXOLQ'UDLQHG\U1 &1RQ'UDLQLQJ 0RXOLQ'UDLQHGd moulP'Figure 4.15: The average local characteristics around draining and non-draining lakes interms of A) ice thickness, B) lake area, C) local strain rate, and D) distance to thenearest moulin (excluding moulins contained within the given lake).4.4.4 Moulin location relative to crevasse fieldsThe presence of crevasses is an important control on moulin formation by a) creating ex-ploitable weaknesses in the ice where crevasses intersect supraglacial channels or b) damagingthe ice, creating down-ice exploitable features where crevasses have closed due to a compres-sive strain regime in the ice. Figure 4.16 illustrates that most moulins of all types are locatedup-ice of crevasse fields. For lake-drainage and other-type moulins, 29 and 33% of crevassefields are down-ice of moulins, respectively, suggesting little to no difference in their formationprocesses relative to crevasse field location.674.4. Results&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU /DNH'UDLQDJHd crevP PRXOLQV,FH)ORZ0RXOLQ/RFDWLRQFigure 4.16: Moulin locations relative to crevasse fields in the study area. ’Up’ or ’down’ice is in relation to local ice flow direction, rather than cartesian coordinates. Thewidth of the violin plots is proportional to the moulin count (e.g. wider plots =more moulins). The white dot represents the median, and the thick black line rep-resents the interquartile range. The thin black line spans the 95% confidence inter-val.4.4.5 Longitudinal trends in moulin distributionMany of the processes on ice sheets are longitudinal in nature - that is, they depend primarily(but not exclusively) on up- or down- ice processes rather than on what is occuring in laterallyadjacent ice. This can be attributed to the downslope nature of gravity driven ice flow (Cuffeyand Paterson, 2010, Ch. 8.2) as well as the longitudinal nature of supra- and sub-glacial drainagenetworks. Given the difficulty in identifying conclusive relationships between moulins and icecharacteristics in 2-dimensional data, in this section I consider whether controls on moulin for-mation may result from longitudinal relationships between ice sheet and surface hydrologicalprocesses and features.To qualitatively investigate possible longitudinal trends in moulin formation, I extractedseven flow-line profiles in the study area spanning a north-south gradient (Figure 4.17). Theflow-lines are extracted along surface velocity vectors as per the MEaSUREs velocity dataset,and therefore represent the path that surface ice travels. These flow-lines represent a rangesection of conditions. Line 1 flows through the low velocity, low strain rate, moulin-densenorthern part of the study area. Lines 2 travels from a lower velocity, higher elevation areaalong the periphery of a fast flowing outlet glacier, intersecting multiple crevasse fields andareas of high strain rate towards its western end. Line 3 follows a high velocity, generallylow strain rate, moulin and crevasse density outlet glacier. Line 4 follows a more complexroute along the periphery of an outlet glacier, primarily in medium velocity ice with high684.4. Resultsstrain rates, intersecting few moulins and multiple crevasses. Line 5 follows a southern outletglacier with higher moulin density than Line 3, and higher crevasse frequency. Lines 6 and 7both terminate in the same general area, but Line 6 starts from a low velocity part of the studyarea, and generally avoids high strain rates and crevasse fields. Line 7 follows a high velocityice stream, intersecting multiple areas of high strain rate and multiple crevasse fields.  'LVWDQFHNP0RXOLQ&UHYDVVHV/DNHV     (yr 1)Figure 4.17: Flow lines overlain on strain rate field. Lines are labeled 1-7, north to south.By investigating strain rate, bedrock and ice surface topography, crevasse field presenceand lake presence along the flowlines, I looked for evidence of longitudinal relationships inmoulin distribution (Figure 4.18). The most obvious relationship in these long profiles is that694.4. Resultsareas where ice surface slope parallels bedrock slope produce few moulins. This is particularlyevident in the protracted parallel ice surface and bedrock slopes in panels F between 75 kmand 125 km and G between 100 km and 130 km. The ice surface slope at Line 7, in particular,largely mirrors the bed slope, which may explain the absence of moulins with the exception ofone crevasse-type moulin.The corollary of this relationship may explain the high occurrence of moulins in the northwestern part of the study area, where there is an adverse bedslope (Panel A). Where bed sur-face slope begins to trend upwards at around 120 km along the flow line, significant crevassingand crevasse-proximal moulins are observed. Thereafter, low strain rates predominate as ve-locity decreases, and with them high supraglacial lake density is observed. Moulins in thisarea are ’other’-type, and appear to occur down-ice of crevassing features.Line 5 (Panel E) passes through the largest number of moulins (15). The reasons for this arenot entirely clear, but a conspicuous difference in Line 5 relative to the other, lower moulin-density panels is a lack of high-relief bedrock undulations until 200 km, resulting in few signif-icant crevasse features. Although crevasse-proximal moulins were observed within a distanceof two ice thicknesses of the flow line, the lack of extensive crevassing or high strain ratesalong the flow line suggests the potential for large contributing areas to develop in this re-gion (Phillips et al., 2011). This is in stark contrast to Line 4, where extensive crevassing andhigh strain rates until 200 km appear to largely limit moulin development. This same limit onmoulin development may explain the absence of moulins along line 7, and the ends of Lines 2and 6.Lines 2 and 3 have broadly similar bed topographies, with undulations for the first 100-150km, until the flow lines intersect the bedrock trough in which the high velocity outlet glacieris located. Ice surface topography shows relatively little variability along these flow lines.Although the flow lines intersect only three lake-drainage associated moulins (Lines 2, 4,and 5), it is interesting that all three are associated with troughs in bedrock topography, withrelatively low ice thicknesses (all less than 750 m) and the nearby presence of crevasses.704.4. Results    (OHYDWLRQP$/LQH     (OHYDWLRQP%/LQH     (OHYDWLRQP&/LQH     (OHYDWLRQP'/LQH      (OHYDWLRQP(/LQH    (OHYDWLRQP)/LQH        'LVWDQFHDORQJOLQHNP(OHYDWLRQP*/LQH%HGURFN&UHYDVVHV/DNHV0RXOLQV&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO2WKHU/DNH'UDLQDJH(yr1 )v iceP\UFigure 4.18: Verticalprofiles of the sevenflow lines shown inFigure 4.17. The strainrate line represents theice surface elevation.Moulins and lakesare shown if they fallwithin two local icethicknesses of the flowline. Grayscale repre-sents surface velocity,and is not meant torepresent a velocityprofile. Crevasses arevertically offset down-wards and moulinsupwards for visualclarity.714.4. Results4.4.6 Comparison with modeled predictionsIn this final section, I ask how well assumptions of moulin formation reflect the reality of iceand hydrological conditions around moulins. In Figure 4.19, I plot different moulin types ac-cording to the relationship between the ratio of total annual water supply (Qa) to local icethickness (Hice) and local strain rate (e˙). The horizontal dashed lines represent an annual dis-charge to ice thickness ratio of 4000 and 2000 m2, equivalent to two fracture area (Fa) thresholdsfrom the literature (Arnold et al., 2014; Clason et al., 2012; Koziol et al., 2017), where Fa refers to thearea of a crevasse under a lake which could contain enough water to result in hydrofracture.If I assume a crevasse at each moulin location of area Fa, and I assume that the total annualdischarge can be accumulated at that moulin location and would therefore be available to in-duce hydrofracture, these horizontal lines represent the Qa to Hice threshold at which I wouldexpect hydrofracture to occur. If this assumption were true, any moulins plotting above thisline would be associated with hydrofracture, and all lake-drainage associated moulins wouldbe expected to plot above this line.The vertical grey line represents e˙ = 0 yr−1, and the red vertical line 0.05 yr−1, the strainthreshold at which Poinar et al. (2015) assumed crevassing to occur. If this threshold were anaccurate assumption of the representation of the conditions required to initiate moulin forma-tion, I would expect few moulins to plot to the right of this line.In Figure 4.19 there are mixed results with regards to these assumptions. Many non-lake-drainage associated moulins plot above the Fa thresholds, although this may be explained bylocal topographic limits on the storage of the total annual meltwater supply. Almost all lake-drainage associated moulins plot above the Fa = 2000 m2 line, suggesting that this simplifyingassumption of lake hydrofracture does indeed do a good job of reproducing conditions underwhich hydrofracture occurs.With regard to strain rate, crevasse type moulins in particular appear to span the widestrange of values, including many negative (compressive) values. Although most moulins doappear to be concentrated between 0 and 0.05 yr−1, a large number are found in regions ofextensional strain rates exceeding this value. Moulins that form at particularly high strainrates appear to be associated with smaller catchments, whereas larger catchments tend to plotcloser to the 0.05 yr−1 threshold, suggesting that higher water availability may initiate moulinformation at lower strain values. A single threshold value of strain rate does not seem likelyto capture the complex dynamics of moulin formation.724.4. Results4a+ice5HJLRQDO)UHT    \U 1 H 4a+ice    \U 1 H &UHYDVVH3UR[LPDO/DNH$VVRFLDWHG0RUDLQH3UR[LPDO2WKHUNP2NP2NP2NP25HJLRQDO)UHTFigure 4.19: Total annual moulin discharge relative to strain rate at a moulin. The verticalgrey line represents a strain rate of 0 yr-1, the vertical red line represents a strain rateof 0.005 yr-1 (the strain rate at which crevassing has been reported, see Poinar et al.(2015)), and the horizontal line represents the threshold at which the total annualdischarge, standardized by ice thickness, is sufficient to generate hydrofracture ifall discharge at that point was available to enlarge an existing crevasse of Fa. Fathresholds are empirical thresholds based on existing models (Arnold et al., 2014;Banwell et al., 2016) that presume an existing crevasse below lakes. Size scaling iscontinous; legend entries are examples only. Background grey histogram repre-sents regional strain rate distribution in the study area, as per the right axis label.734.5. Discussion4.5 Discussion4.5.1 LimitationsThis study encompasses a number of methodological limitations. The first of these pertainsto the visual identification of surface features. I have identified moulin and crevasse field lo-cations from visual inspection of DEMs and supraglacial lakes from Landsat-8 imagery. Theprimary limitations of this study lie in the identification of these features themselves as wellas the coarse temporal and spatial resolution of the gridded datasets. Care was taken to en-sure that supraglacial rivers do not cross the boundaries of delineated catchments; however Icannot be certain of how many moulins below the resolution of the imagery might have beenmissed. Additionally, I visually identified crevasse fields, and as such the extent of the fieldsare set by judgment of the researcher and not an objective set of quantitative criteria. This mayparticularly impact the classification of crevasse-proximal moulins.The second limitation of this study pertains to the temporal and spatial resolution and un-certainty of the input data, particularly the bed elevation and ice thickness (at a resolutionof 1 km) and the ice velocity measurements (which are derived from winter velocities only).In particular the temporal resolution of the ice velocity measurements is such that summerspeedup and associated strain rates are unknown, and transient, short lived increases in strainrate from tensile stress propagation cannot be captured in this analysis. Because of this lim-ited spatial resolution, it is not surprising that no relationship between strain rate and moulinlocation could be found, a limitation also expressed by Williamson et al. (2018). It is likely thatrelating strain rate to moulin formation can only be done with high-temporal-resolution (atleast daily) ice velocity measurements.Finally, this study is limited by its spatial extent and thus its limited number of replicates.A wider range of ice velocity, melt, and surface topography conditions as well as a largernumber of replicates would greatly enhance the findings of this study.4.5.2 Moulin TypesSmith et al. (2015) is, to my knowledge, the first to identify moulins according to their formationmechanism across a large study area. They report that lake drainage accounts for only 3% ofmoulin types, with 16% forming due to crevasse fields and the remaining moulins formingfrom undetermined mechanisms or shear fracture (36 and 44%, respectively). Conversely,depending on the melt conditions with respect to elevation, Koziol et al. (2017) found that just38% of moulin-drained surface melt drained through moulins outside of lake basins, withthe remainder draining through lake-basin moulins. However, Koziol et al. (2017)’s model isbased on Fa simulated lake hydrofracture rather than observed lake drainage, and further,their study was limited to only 45 observed moulins. My results align most closely with Smithet al. (2015)’s observations that most moulins occur outside of crevasse fields and lake basins.Whereas other studies of moulin type have reported the partitioning of surface meltwa-744.5. Discussionter between moulin types (Koziol et al., 2017) or the distribution of moulin types with eleva-tion (Smith et al., 2015), my study is the first to look explicitly at the distribution of moulintypes with regards to ice characteristics including strain rate and velocity. Although I findthat crevasse-proximal moulins are located, as expected, in areas of higher strain rate thanother moulin types, I find few other differences that provide insight into formative conditionsof the enigmatic ’other’ type moulin. This is a persistently pernicious problem with inte-grating moulins into hydrological models. Whereas lake-drainage appears to be adequatelysimulated by the use of a threshold fracture area, existing models (e.g. Banwell et al., 2012,2013; Clason et al., 2012, 2015) do not have a mechanism for integrating non-lake-drainage- orcrevasse-associated moulins. ’Other’ type moulins have only been integrated through visualidentification and integration into the routing model (Koziol et al., 2017).I did observe that ’other’ type moulins are particularly concentrated in low-velocity ice,positioned closer to lakes, and, like the other three moulin types, typically located upstreamof crevasse fields. Crevasse-type moulins were distinct from other moulin types with regardsto higher ice thickness, strain rate, elevation and velocity, and lower total annual discharges.Lower discharges may be attributable to smaller catchment areas resulting from supraglacialchannel truncation upon intersection with a crevasse. However, a single strain rate thresholddoes not appear to adequately capture the initiation of crevasse-proximal moulins. Smallercatchments and lower discharge yields relative to ice thickness may require higher strain ratesto initiate hydrofracture, for example. Additionally, crevasse-type moulins often and unex-pectedly formed in compressive ice. These moulins may be associated with splaying-typecrevasses which result from compressive ice regimes. I note, however, that our methodol-ogy of visually identifying crevasses may omit some crevasses and crevasse fields that werenot readily identifiable in the DEM, and that the attribution of moulins within one local icethickness of visually identified crevasses may erroneously attribute non-crevasse-associatedmoulins to crevasse formation. Nonetheless, the low sensitivity of moulins to strain ratesacross moulin types, coupled with the observation of moulins up to 25 km inland echoes Yanget al. (2016a)’s caution in projecting moulin formation processes into the warming future with-out proper constraints on high-elevation moulin formation (Poinar et al., 2015).4.5.3 Moulin ClusteringMoulins in our study area appear to be aggregated in clusters that are visually distinct andare statistically significantly spatially autocorrelated. This suggests that they are either a) co-producing, such that the formation of a moulin might trigger the nearby formation of othermoulins or (and) b) that they are highly dependent on local conditions and thus form in similarand spatially homogeneous environments.Indeed, although like others (Clason et al., 2012; Colgan and Steffen, 2009; Phillips et al., 2011)I did not find that regression analysis using ice characteristics was a productive approachfor predicting moulin distribution, I did observe that moulins are not randomly distributed754.5. Discussionrelative to ice conditions. Specifically, I observed that moulin distribution is highly sensitiveto ice velocity. Moulins form at statistically significantly lower ice velocities compared to thedistribution of ice velocity in the study area, and appear to be particularly concentrated inareas of localized high extensional strain rates within regionally low velocity zones. This isan important spatial distinction because velocity does not tend to vary rapidly over shortdistances, resulting in spatially homogeneous ’zones’ of relatively constant velocity. However,high strain rates, which can be of an equal magnitude between two low-velocity pixels or twohigh-velocity pixels, show significantly more spatial heterogeneity, making high strain ratesmore evenly spread throughout the study area. For example, the highest strain rate in our’low’ velocity zone was 0.11 yr−1, and was 0.10 yr−1 in our ’high’ velocity zone.This clustering appears to be dependent on moulin type, such that low-velocity zones aremost strongly associated with dense ’other’ type moulins. As ’other’ type moulins are themost common moulin type, this clustering in low-velocity zones is significant, as it is one ofthe most distinguishing characteristics of ’other’ type moulins.Why might moulins, and in particular ’other’ type moulins, form predominantly in lowvelocity zones? It is tempting to attribute higher velocities to higher strain rates, howeverhigh-velocity ice can be conspicuously free of high strain rates over large distances, just aslow-velocity ice regions can contain localized high strain rates. Such examples in our studyarea include the outlet glaciers traversed by flow lines 3 and 5, in which high surface strainrates are confined to the periphery of the ice stream where high velocity ice shears againstlower-velocity ice.The clearest example of clustering in low velocity ice is the high density cluster of moulinsin the low-velocity ice of the northwest corner of the study area, through which flowline 1runs. What stands out in this area is the inverse direction of the downward sloping ice tothe upward sloping bedrock. What stands out in areas of low moulin density is a relativelyparallel topographic profile along the ice surface and basal topography. Where the bed slopeparallels the surface slope, ice flow is governed by the topographic bed gradient and the effectof the ice thickness gradient is minimal (or absent). As can be seen in Lines 6 and 7, this createsstrain rates close to 0 yr−1, and moulin formation appears to be inhibited.Ultimately, bedrock topography fundamentally influences ice flow, and thus surface dy-namics expressed in variables such as strain rate and velocity reflect conditions at the bed.Understanding the mechanisms that control moulin distribution therefore can be informed byan understanding of how bed topography relates to moulin distribution. There are disparatetheories and observations of the relationship between moulin formation and bedrock topog-raphy. Catania et al. (2008) observed an increased frequency of moulin in regions where icethinned over bedrock ridges, and Joughin et al. (2013) proposed that stretching and thinningof ice over bedrock ridges produces high tensile strain rates and crevasses, thus providingopportunities for moulin formation through hydrofracture. These observations seem to beconsistent with our findings that low variability in bedrock slope relative to ice surface slope764.6. Conclusionproduces few moulins.I suggest that the relationship between basal topography and moulin formation is a promis-ing area of research that might provide insight into projecting future trends in surface waterdrainage that go beyond considerations of strain rate variability with elevation (e.g. Poinaret al., 2015). Work in this vein has been limited but promising (e.g. Lampkin and Vanderberg,2011), and is worthy of further researcher investment.4.5.4 Influence of lakes on moulin distributionThis is, to my knowledge, one of the most detailed empirical analyses of lake drainage, inwhich I investigate relationships between lake drainage and strain rate, lake area, ice thickness,and the distance to the nearest moulin.I found that the strongest determinant of lake drainage relative to non-drainage was thedistance to the nearest moulin. ’Other’ type moulins were typically closer to lakes, and morethan two thirds of the moulins nearest to lakes were located down-ice of the water body. Al-though, similar to findings by Fitzpatrick et al. (2014), lake size did not appear to differenti-ate between draining and non-draining lakes for larger lakes, very small lakes typically didnot drain. Additionally, lakes at higher strain rates had higher incidents of in-situ moulindrainage.Collectively, I interpret these findings to support the contemporary hypothesis that lakedrainage is not a water volume-threshold initiated process, but a dynamic process in whicheither nearby moulin formation triggers lake drainage through transient increases in stress(Stevens et al., 2015), or vice versa: lake drainage creates ice flexure, allowing moulin formationnearby, where fractures are opened up temporarily and exploited by surface water where theyintersect streams (Hoffman et al., 2018). If this is the case, it explains why annual strain ratealone is not a useful predictor of moulin location. Spikes in higher temporal resolution velocitydata may prove to be more productive in determining where and when moulins form.4.6 ConclusionThis chapter presents a detailed empirical description of the spatial distribution of moulinswith regards to the distribution of other moulins, ice sheet characteristics, and supraglacialhydrological features. Many of my observations confirm findings of previous work. For ex-ample, I observe, as have others, that moulin density peaks at intermediate elevations (Yangand Smith, 2016; Yang et al., 2015), which is likely attributable to extensive crevassing on icesheet margins and outlet glaciers as well as the difficulty that identifying moulins in crevassedice presents (Joughin et al., 2013). My findings reinforce that most moulins do not appear tobe associated with supraglacial lake hydrofracture or to crevasses, but rather are dominatedby another formative process (Catania et al., 2008; Koziol et al., 2017; Smith et al., 2015; Yang andSmith, 2016).Although, as with previous attempts to synthesize moulin distribution statistically (Clason774.6. Conclusionet al., 2012; Colgan and Steffen, 2009; Phillips et al., 2011), I did not document any conclusive con-trols on the formation of moulins, I observed that moulins are spatially clustered and autocor-related, and most densely concentrated at localized areas of high strain rate located in broaderregions of low velocity. I observed that moulin densities are lowest in areas where ice surfacetopography parallels bed slope, and I suggest that further analyses into the basal controls onmoulin distribution may prove useful for extrapolating moulin location predictions into thefuture of a changing ice sheet. Additionally, empirical evidence presented here of moulin loca-tion with respect to lake position strengthens recent hypotheses that lake drainage events arecorrelated with moulin formation. Further, the timescale of these events is sufficiently shortthat annual scale strain rate data is not a useful predictor of moulin location.Overall, this work highlights the difficulty of classifying moulins and associating theirformation with intuitive processes and readily accessible proxies. It also challenges the simpleconception of moulin formation as currently treated in models.78Chapter 5Spatial variability of supraglacialcatchment morphometrics andhydrology on the southwesternGreenland Ice Sheet5.1 SummaryThe physical characteristics of supraglacial catchments dictate where and when surface melt-water is delivered to the subglacial system on a seasonal and diurnal timescale. In this chapter,I quantify catchment morphometrics for supraglacial catchments in a study area in southwestGreenland, and use these morphometrics to derive hydrograph parameters using synthetichydrograph theory. I find that peak hydrograph lag times can range from less than a hour tomore than 11 hours. The exact impact of this range of supraglacial hydrographs on subglacialhydrology is contingent upon the spatial configuration of subglacial catchments, which de-pends on subglacial water and ice overburden pressure. Under conditions of high subglacialwater pressure, subglacial catchments towards the interior of the ice sheet become more nu-merous and elongated, and thus are fed by fewer moulins and are more sensitive to theirindividual hydrographs.5.2 IntroductionMeltwater generated on the Greenland Ice Sheet (GrIS) surface is collected in internally drainedcatchments and delivered to moulins (Koziol et al., 2017; Smith et al., 2015; Yang and Smith, 2016).There, near-vertical pathways convey the meltwater to the ice sheet-bed interface. It is welldocumented that meltwater at the bed can increase water pressure and facilitate ice slidingand acceleration (e.g. Hoffman et al., 2011; Iken, 1981; Palmer et al., 2011; van de Wal et al., 2008),795.2. Introductionand further that this effect is modulated temporally by an increase in subglacial channel con-veyance capacity as these channel networks evolve (e.g. Andrews et al., 2014; Moon et al., 2014;Schoof , 2010).Subglacial channels form through channel melt and the melt opening of channel walls, inopposition to ice overburden pressure that acts to close them through creep closure (Rothlis-berger, 1972). Once surface water supply to the subglacial network exceeds a critical dischargerelative to subglacial capacity, the drainage system evolves from distributed to channelized,and subglacial water pressure decreases, eventually decoupling ice velocity from water sup-ply (Kamb, 1987; Schoof , 2010). Where and when on the ice sheet meltwater is supplied tosubglacial networks has been shown to play an important role in subglacial channel evolu-tion; more densely spaced moulins result in an earlier seasonal transition from distributed tochannelized drainage (Banwell et al., 2016). The important role of surface water supply in in-fluencing subglacial channel development and thus ice velocity has motivated recent scientificinterest in the spatial distribution of moulins and supraglacial catchments (e.g. Banwell et al.,2016; Koziol et al., 2017; Smith et al., 2015, 2017; Yang and Smith, 2016).Empirical studies suggest considerable variability in catchment size, drainage density, net-work morphometry, drainage dynamics and moulin discharge (Koziol et al., 2017; Smith et al.,2015; Yang and Smith, 2016; Yang et al., 2016a), which collectively influence diurnal moulin hy-drographs (Smith et al., 2017). Complementary studies have shown that both ice velocity andsub-glacial channel water pressure vary substantially seasonally and diurnally (Andrews et al.,2014; Hoffman et al., 2011; McGrath et al., 2011; Meierbachtol et al., 2013). It is therefore reasonableto suspect that daily moulin hydrographs play an important role in daily ice dynamics and inthe seasonal development of the regional subglacial drainage network.However, the extent to which supraglacial catchment characteristics and hydrology vary inways that influence the regional temporal and spatial variability in meltwater supply remainsunclear. Contemporary models of supra- to sub-glacial water routing (Banwell et al., 2016, 2013;Clason et al., 2015; Koziol et al., 2017; Poinar et al., 2015) are reliant on assumptions of moulin andcrevasse formation that do not appear to reflect the complexity of moulin formation processes(see Chapter 4). What are the hydrological implications of the discrepancy between empiricalobservations of a complex and dynamic supraglacial hydrological system and the requisitesimplifying assumptions of numerical models?In this chapter, I consider whether supraglacial catchment characteristics vary in ways thatmight be important to subglacial drainage network development. I consider which catchmentcharacteristics likely constitute important controls on regional hydrograph variability, andhow this variability changes with different configurations of subglacial drainage networks. Ibuild on work by Smith et al. (2017) to employ traditional hydrological methods of hydrographprediction to investigate the empirical differences in catchment hydrology across a study areain the south western flank of the GrIS. Further, I build on work by Yang et al. (2016a) to quan-tify the variability of catchment morphometrics, but I extend this work to consider the ways805.2. Introductionin which this catchment variability might (or might not) matter in the context of the surfacemelt/ice acceleration dynamic. I derive both seasonal and diurnal mid-summer hydrographsfor moulins in our study area, and investigate the extent to which these catchment character-istics produce physically meaningful differences in hydrograph characteristics.5.2.1 Background: Predicting hydrograph characteristics in ungauged catchmentsVery few hydrographs have been constructed for supraglacial channels on the GrIS (Gleasonet al., 2016; Smith et al., 2017). Furthermore, there are thousands of catchments across the icesurface, ranging in size from only several hundreds of square meters to many tens of squarekilometers, covering an enormous range of ice conditions (Yang and Smith, 2016; Yang et al.,2016a) and evolving temporally.The problem of predicting hydrographs in ungauged basins is a familiar one in traditionalhydrology, where flood hydrographs are of significant research and practical interest. In thesesystems, many approaches have been developed through time to facilitate such predictions,ranging from simple empirical scaling with respect to a reference catchment, to statistical ap-proaches, to parameterized, distributed numerical GIS-based models (Dingman, 2015, Ch. 2).Conceptually, it is well understood that the catchment hydrograph is determined by processesoccurring in two separate spatial domains: the unchannelized hillslope dominated by diffu-sive processes, and the stream network dominated by channelized flow (Bracken and Croke,2007; Tucker and Bras, 1998). Hillslope characteristics, including length, slope, and surface/-subsurface characteristics control the amount of time required for water to travel from a givenpoint on the landscape to the stream network. Stream network characteristics including itslength, density and slope, control how rapidly water moves from where it enters the channelto the catchment outlet (Dingman, 2015).Unit Hydrograph (UH) theory is one of the oldest approaches to estimating catchmenthydrographs based on catchment characteristics (Sherman, 1932). In UH theory, a unit hy-drograph is the hydrograph resulting from one unit (inch, cm, etc...) of constant intensity anduniformly distributed rainfall in a catchment. A UH can be empirically derived from observeddischarge and used to predict catchment response to storms of different intensities by convolv-ing the UH curve with a storm hyetograph. However, when a basin is ungauged, the UH mustbe constructed (i.e. synthesized), giving rise to various approaches to creating a Synthetic UnitHydrograph (SUH).SUH theory is one of the most popular ways to predict ungauged catchment hydrographsbased on catchment-scale characteristics, and this technique has recently been employed onthe GrIS by Smith et al. (2017). There are various approaches to building a SUH (see Singhet al., 2014), with the Snyder Unit Hydrograph (Snyder, 1938) being one of the most popularand the one employed by Smith et al. (2017). The Snyder SUH relies on a number of empiricalparameters as well as three simple catchment morphometrics: catchment area, length, and ameasure of elongation.815.3. Study Site and MethodsHowever, supraglacial drainage networks are extraordinarily dense (Smith et al., 2015), andwater transport is likely dominated by the drainage network (fluvial) process domain. Refin-ing a SUH for supraglacial applications would benefit from integration of drainage networkmorphometrics. A number of such SUH approaches have been developed, collectively termedGeomorphological Instantaenous Unit Hydrographs (GIUH) (Rodriguez-Iturbe and Valdes, 1979),and make use of catchment morphometrics as outlined in Table 5.1 and described below.In the following sections, I describe the study area, the methods used to extract the catch-ments, drainage networks, and their morphometrics, and the methods used to build the GIUHcurves and hydrograph variables.5.3 Study Site and Methods5.3.1 Study AreaIn this study I focus on a 2700 km2 part of the ice surface on the south west margin of theice sheet, approximately 150 km south east of Nuuk (Figure 5.1). This is the same study areaused in Chapter 4. It extends from 584 to 1761 m in elevation, with a maximum ice velocityof 245 m/yr, and includes (portions of) three outlet glaciers that drain from two large bedrocktroughs. Ice thickness is typically less than 1 km but is in excess of 1.5 km in places. AppendixA includes images of surface slopes, ice thickness, ice velocity, and surface principal strainrates in the study area.I focus on this study area to complement the work done in Chapter 4 to continue to build aholistic understanding of how bedrock, ice, and hydrological processes collectively influencethe relationship between surface melt and ice sheet dynamics. The original motivation forchoosing this region was a coincidence of high quality data in the expansive south westernGrIS ablation zone; an area that has not previously been the focus of supraglacial hydrologystudies.825.3. Study Site and Methods49°W50°W51°W64°N63°40'NStudy Area^Nuuk20 Km¹Figure 5.1: Study area in south west Greenland. Background imagery is Landsat-8 im-agery dated August 02 and 04, 2015. Inset data is the outline of the Greenland IceSheet, where the location of the study area is shown with a red star.5.3.2 DatasetsFour different data products were used in this analysis and are described in Chapter 4 aswell as briefly below: an ice surface DEM, ice velocity data, bed elevation data, and runoffproduction data from a climate model.This study capitalizes on the recent release of 2.0 m resolution Digital Elevation Models(DEMs) for the GrIS ice surface. These DEMs are photogrammetrically derived from 0.5 mresolution WorldView-1 imagery by the Polar Geospatial Center using the Surface Extrac-tion with TIN-based Search-space Minimization (SETSM) algorithm, as described by Noh andHowat (2015), and acquired from: This DEMhas reported root mean square errors of 3.8 and 2.0 m in the horizontal and vertical, respec-tively. The DEM tiles used in this study are derived from 2015 imagery acquired on, in orderof increasing elevation, September 2nd, August 25th, and July 13th. Ice surface DEMs wereused to visually identify moulins and to produce moulin catchments through D8 flow routingusing the ArcHydro Toolbox in ArcMap 10.3.1 (Maidment, 2002), as described in Chapter 2 andKing et al. (2016).Ice velocity data, used to derive surface strain rates, were available in the form of 200 m835.3. Study Site and Methodsresolution gridded ice velocity rasters from 2015-2016 (MEaSUREs Greenland Ice Sheet Veloc-ity Map from InSAR Data, available from The velocitydata refer to winter velocities, broadly defined from September 1st to May 31st. I derived prin-cipal surface strain rates from this velocity data following Poinar et al. (2015). Gridded velocityproducts were filtered using a 3 x 3 pixel (600m) gaussian blur (smoothing) filter prior to de-riving strain rates. Velocity gradients were calculated using central differences, implementedwith Numpy’s Gradient function.Bed elevation and ice thickness data were acquired from 150 m resolution BedMachineV3data (Morlighem et al., 2017) available from This data prod-uct is derived from a range of radar surveys carried out between 1993 and 2015 and general-ized using a mass conservation approach. The error in bed and ice thickness in the study areais an average of 106 m, as per the included gridded error data fields.Daily meltwater runoff production for the study area is provided by the 200 m resolutiongridded daily RACMO2.3 climate model (Noël et al., 2015). Runoff production in dimensionsof length/time was translated into total annual discharge by intersecting the modeled runoffproduction with catchment area, as described in Chapter 3. We observed that taking the sumof all intersecting cells underestimates total runoff in catchments that are small relative tothe 100 m resolution of the runoff data. For this reason, catchment discharges were from themedian runoff of all the grid cells within a catchment.5.3.3 Surface catchment and channel network delineationCatchments and channels were extracted from the DEMs using standard flow routing algo-rithms, implemented in ArcHydro toolbox for ArcMap 10.3.1 (ESRI, 2016; Maidment, 2002). Afull description of the delineation method is provided in King et al. (2016). In summary, flowis routed over the DEM to pre-identified moulin or crevasse locations. Because there is noautomated approach to moulin identification (e.g. Yang and Smith, 2016), sinks were visuallyidentified from the DEM, facilitated with a 1 km by 1 km grid search of elevation and curva-ture layers. Moulins were identifiable as areas of abrupt channel termination, and were oftenvisible as large depressions or holes in the DEM or where channels terminated in crevasses.Additionally, crevasse fields were manually digitized as polygons. DEM sinks within thosecrevasse fields were preserved during DEM filling. Catchments that abut the DEM tile periph-ery were excluded, as they may not represent complete catchments. Limitations of the flowrouting methodology are discussed in King et al. (2016).As WorldView imagery was not available, panchromatic Landsat 8 imagery with a 15 mresolution was used to assist in moulin identification when possible. However, many sinks aretoo small to be visible in 15 m resolution imagery, and many sinks were evident in the DEMthat could not be seen in LANDSAT imagery. Much like Joughin et al. (2013), I found moulinidentification more difficult at lower elevations where abundant microtopographic featuressuch as crevasses made moulin identification difficult. Additionally, our moulin dataset al-845.3. Study Site and Methodsmost certainly excludes smaller moulin features that are below the resolution of the DEM andthus overestimates catchment area (particularly higher up in the drainage network, see Kinget al., 2016). Therefore our dataset does not represent a complete inventory of moulins in thestudy area. However, it does reproduce the general drainage patterns arising from moulinsdraining stream channels apparent in the 2 m resolution DEM.Identifying the full extent of a drainage network is not straightforward and has been thesubject of significant study through the years (e.g. Montgomery and Dietrich, 1989; Montgomeryand Foufoula-Georgiou, 1993; Montgomery and Dietrich, 1988; Tucker and Bras, 1998; Tucker et al.,2001). Several methods of stream initiation identification from DEMs have been proposed in-cluding, for example, curvature and slope scaling analysis (e.g. Tarboton et al., 1991). Inspectionof the DEM derived curvature layers suggests that curvature adequately identifies main chan-nels, but given the resolution and noise of the DEM it fails to identify smaller channels. Slopescaling analysis has similarly been shown to inadequately capture the transition from unchan-nelized to channelized domains in supraglacial environments (Karlstrom and Yang, 2016).Previous analysis of 1-m resolution WorldView-1 imagery suggests that drainage densityon the GrIS can be exceptionally high relative to terrestrial environments (Smith et al., 2015),and visual inspection indicates that under full melt conditions small tributary channels extendalmost all the way to the drainage divide (Figure 5.2). Therefore, in this study I define streamchannels using a simple accumulation area threshold of 0.04km2 (equivalent to a 200 m by 200m area) to reflect the high channel density that would be expected at peak melt. This thresholdis low enough to generate dense channel networks but sufficiently high that, based on visualinspection, it does not extend drainage channels into unchannelized areas.250 mFigure 5.2: An 800 by 1000 m closeup of a WorldView-1 panchromatic image (image ID12AUG20151700-P2AS-R2C2-055095122010-01-P001) of supraglacial channels to thenorth of our study site. Blue lines are channels specified using a 0.04 km2 contribut-ing area threshold.855.3. Study Site and MethodsThis accumulation area threshold for channel initiation introduces a methodological lim-itation into our study. It likely results in an under-estimation of the channel network extent.As can be seen in Figure 5.2, it does not capture many of the smallest tributary channels.Nonetheless, given that a universally applicable initiation area threshold does not exist, I choseto under- rather than over-estimate the network. Further, to employ a low universal initiationthreshold, I assume melt across the entire study region such that the channels are fully devel-oped and no parts of the catchments remain covered in seasonal snowpack.5.3.4 Deriving catchment characteristicsCatchment and channel network characteristics that are traditionally used in derivations of hy-drographs are shown in Table 5.1. I derived these morphometrics for supraglacial catchmentsand channels in our study area using ArcMap 10.3.1 (ESRI, 2016) functions. The catchment-scale bifurcation, length and area ratios (RB, RL and RA, respectively) were derived for eachstream order as per Table 5.1, and the catchment-scale value of each was determined as:RB =1em(5.1)RL, RA = em (5.2)where m in each case is the slope of the linear regression of the following respective rela-tionships:log Nw, log L¯w, log A¯w vs. ω (5.3)where w refers to the stream order. Unlike alluvial channels, the highest order stream of asupraglacial channel network ends abruptly in a moulin. As a result, the highest order channelis often shorter than the lower order channels. I therefore excluded the highest order channelfrom the calculation of RL.865.3. Study Site and MethodsTable 5.1: Summary of catchment and network morphometrics used in this study.The terminal moulin type for each catchment was determined as described in Chapter4. Four moulin types were identified. 1) ’Crevasse-Proximal’: when a moulin was within ahorizontal distance equivalent to one local ice thickness of a digitized crevasse field 2) ’Lake-Drainage’: a moulin associated with lake drainage based on visual inspection of Landsat-8 imagery over the 2015 melt season 3) ’Moraine-Proximal’: a moulin that appeared to beassociated with either moraine debris or nunataks, and finally 4) ’Other’: all other moulins.5.3.5 Deriving catchment hydrograph characteristicsI am particularly interested in how catchment characteristics influence daily hydrographs be-cause of the strongly diurnal velocity changes observed on the ice sheet (e.g. Andrews et al.,2014), however I am limited by a lack of hydrological data in this part of the ice sheet (indeedin most parts of the ice sheet). I make use of the Geomorphic Unit Instantaneous Hydrographbecause, given the high drainage density of these catchments, it is reasonable to expect thatwater routing is dominated by channel flow rather than hillslope processes (Montgomery andDietrich, 1989). The GIUH is explicitly based on channel network characteristics and is there-fore useful for moving beyond the SUH approach. Given the aforementioned assumption offull melt across the catchment, I consider daily hydrologic characteristics only for snow-freeconditions and full drainage network development, i.e. I assume that GIUH parameters thatare derived here apply only under snow-free conditions. This further simplifies the analysisby assuming minimal retention of runoff in seasonal snowpack, allowing us to not account forthe delay or retention of runoff in the snowpack.The GIUH relies on four morphometric parameters: the length of the main channel (L), andthe bifurcation, area, and length ratios (RB, RA and RL), all of which can be extracted from the875.3. Study Site and Methodsflow routed channel network. The GIUH requires a hydraulic variable, v, the average peakdiurnal flow velocity (characteristic velocity, in m/s). These five parameters can be used toidentify time to peak (tp, in hr−1) following Rodriguez-Iturbe and Valdes (1979) as:tp = 0.44(Lv)R0.55B R−0.55A R−0.38L (5.4)Even with field data with which to calibrate these equations, estimation and even defini-tion of v remains a challenge in GIUH methodology (Singh et al., 2014). In this study, then, Iemploy a very crude estimate of v based on the Manning’s equation:v =( 1n)R23√Sc (5.5)where R is hydraulic radius, n is the Manning’s n (empirical roughness factor) and Sc ischannel slope. Although Sc is readily derivable, hydraulic radius is not. I therefore assumethat R approximates channel depth d, and estimate d using empirical data for supraglacialstreams on the GrIS from Gleason et al. (2016) where d = 0.0001L + 0.01. Additionally, I useempirical values of n as reported by Gleason et al. (2016) where average n =0.035.Clearly this is a problematically rough estimate of v. However, previous studies havefound that supraglacial channels tend to widen rather than deepen (Yang et al., 2016a), andGleason et al. (2016) report that depth and surface velocity in supraglacial channels remain rel-atively insensitive to downstream distance, whereas channel slope appears to decrease rapidly.I therefore suggest that the crude approximation of v used herein is sufficient to facilitate a rela-tive comparison of GIUH parameters, although I do not suggest that these reported values areaccurate in absolute terms. Nonetheless, our derived v and tp values are well within reportedempirical ranges reported by Smith et al. (2017) and Gleason et al. (2016). Appendix D illus-trates relationships between estimated hydraulics and catchment properties. Field-observedlag time reported by McGrath et al. (2011) for a 1.14km2 catchment as well as the SUH esti-mates for Rio Behar by Smith et al. (2017) plot well within our calculations of tp as a functionof catchment area.To create GIUH curves q(t), I follow Kumar et al. (2007) where:q(t) =1kΓ(n)(t/k)n−1e(−t/k) (5.6)where Γ is the gamma function, and n and k are defined based on Horton’s Laws andderived semi-empirically as:n = 0.5764[RBRA]0.55· R0.05L (5.7)k =0.44Lv·[RBRA]0.55· R−0.38L ·1n− 1 (5.8)885.3. Study Site and MethodsGIUH curves require channel network morphometrics. Crevasse-drained catchments donot have channelized drainage, and therefore do not have channel network morphometrics.For these catchments, I estimated tp using the SUH approach employed by Smith et al. (2017)as outlined in Chapter 3, which only requires catchment morphometrics. In this approach, tpis a function of the main stem length (L, in km), distance between the point on the center flowline closest to the centroid of a given catchment and the catchment’s moulin (Lc, in km) andan empirically defined coefficient Ct such that tp = Ct(Lc · L) Subglacial catchment delineationI assigned supraglacial catchment water delivery through moulins to their respective sub-glacial catchments. I assumed that moulins deliver water vertically downward to the bed, i.e.that water reaches the bed immediately vertically below the moulin. Thereafter, I assumed, ashave many others (Banwell et al., 2013; Flowers and Clarke, 1999; Lewis and Smith, 2011; Lindbacket al., 2015; Pitcher et al., 2017; Shreve, 1972), that water flows subglacially at the ice-bed interfaceon impermeable bedrock, following the steepest gradient of hydraulic potential (a function ofthe combined pressure potential related to the overburden of ice and the elevation potentialproduced by subglacial topography). Subglacial hydraulic potential (φ) is calculated as:φ = kρig(h− z) + ρwgz (5.9)in which ρi is the density of ice (917 kg m−3), ρw is the density of water (1000 kg m−3), gis the acceleration due to gravity (9.81 m s−2), h and z are the elevation of the ice and bedrocksurface (m), respectively. The parameter k is the ratio of water pressure (pw) to ice overburdenpressure (pi).If there is no ice overburden pressure, k is equal to 0 (i.e. atmospheric pressure), anddrainage is simply influenced by bedrock topography. As water pressure increases and k ap-proaches 1 (i.e. pw is equal to pi), the drainage network becomes more strongly influencedby ice thickness (i.e ice surface topography). Over long timescales and large spatial scales, kcan be assumed to equal one (Shreve, 1972). However, empirical estimates of k suggest thatthe value ranges spatially and temporally (Banwell et al., 2013; Chu et al., 2016; Lindback et al.,2015). In the absence of practical constraints on values of k, prior catchment delineation worksuggests varying k to constrain catchment boundary uncertainty (Banwell et al., 2013; Chu et al.,2016; Pitcher et al., 2017). Given the large ice thicknesses found on the GrIS, I assumed a fairlyhigh value of k as ice overburden pressure is a function of ice thickness. As per Pitcher et al.(2017), I varied k between 0.8 and 1.1.Again following Pitcher et al. (2017), I delineate subglacial catchments by substituting hy-draulic potential fields for elevation fields during flow routing in ArcGIS 10.3. Although Flow-ers and Clarke (1999) note that a multiple-direction flow routing algorithm best represents ahydraulically diffuse glacier bed, in this study I employ a D8 flow routing algorithm whichproduces discrete catchment boundaries (Tarboton, 1997).895.4. Results5.4 Results5.4.1 Spatial variability of supraglacial catchment characteristicsThe derived characteristics of supraglacial catchments are shown in Figure 5.3. Catchmentsvary in size from 0.04 to 58.9 km2, with a median size of 6.5 km2 (Figure 5.3A). Catchmentsat lower elevations appeared to be smaller, with the interesting exception of catchments onthe northern outlet glacier. As described in Chapter 4, this outlet glacier is a high velocity,low-crevasse and moulin density glacier draining through a deep bedrock trough.Lakes are substantially more widespread at high elevations (Figure 5.3B). Most catchmentshave few to no lakes, however several higher elevation catchments have cumulative lake areasup to 1.56 km2. As expected, median catchment runoff production decreases almost linearlywith elevation, ranging from 4.5 m/yr to just 0.21 m/yr (Figure 5.3C).905.4. Results  'LVWDQFHNP$ &DWFKPHQWVAcNP2  'LVWDQFHNP%  ALNP2  'LVWDQFHNP& 0HOWP\UFigure 5.3: Characteristics of the supraglacial study catchments, including A) catchmentarea, B) total lake area, and C) mean catchment melt. All values are catchment-meanvalues. To accentuate differences and limit the visual effects of outliers, the colormap extent is stretched between the 25th and 75th percentiles of the data. The color-bar refers to the color gradations in the maps, and is also the x-axis of the histograms,which shows the distribution of catchments characteristics.Catchment elongation ratio (a function of the ratio of catchment area to length), surfaceslope, and the slope of the highest order stream in each catchment are shown in Figure 5.4.Smaller elongation ratios imply more elongate catchments. For a given channel slope andcatchment area, water in a catchment with a smaller Re would take longer to reach the moulinfor a given amount of melt. In our study area, elongation ratios range from a lower 5th per-centile of 0.25 to a 95th percentile of 0.60, with a median of 0.39. Catchments appear moreelongate along the outlet glaciers where ice is flowing faster (see Figure 4.10) as well as abovethe 1500m contour line (Figure 5.4 A). Catchments along the peripheries of the glaciers aswell as in the low velocity ice in the northernmost section of the study region appear morecompressed. Generally, the elongation ratio appears to mirror the distribution of ice surfacevelocity (Figure A.2).915.4. Results  'LVWDQFHNP$  &DWFKPHQWVRe  'LVWDQFHNP%  Ss  'LVWDQFHNP&  ScFigure 5.4: Morphometrics of the study catchments,including A) catchment elongationratio, B) catchment surface slope and C) the slope of the highest order channel ineach catchment. All values are catchment-mean values, and values are stretchedbetween the 25th and 75th percentiles to accentuate differences. Note that panelsB and C can only be calculated for channelized basins, and therefore the crevasse-drained catchments included in panel A are absent.The apparent spatial coincidence of low Re values and high ice velocity prompted us toexamine the relationship in more detail. Catchments with higher annual velocities do, onaverage, have lower Re values, indicating more elongated catchments (Figure 5.5). The inverserelationship is statistically significant (p < 0.001).925.4. Results  9HORFLW\P\U5HFigure 5.5: Binned median elongation ratio Re for given velocity ranges (black line). Greyshaded area represents the 25th and 75th percentiles of the binned data.The 5th and 95th percentile of catchment slopes Ss range from 1.01 to 4.8%, with medianslopes of 2.1%. Particularly high slopes were confined to the ice sheet peripheries and tocatchments draining towards moraines or nunataks (Figure 5.4 B). A band of slightly elevatedcatchment gradients is evident along the 1500m contour, consistent with a high surface slopesjust below 1500m (see Ch 4, Fig 4.4). Channel slopes SC are generally lower than catchmentsurface slopes and do not appear to vary over as wide a range of values, with a 5th percentileof 0.38 %, a 95th percentile of 3.0 %, and a median value of 1.2 % (Figure 5.4 C).The three channel network ratios of Horton’s Laws are shown in Figure 5.6. Bifucationratio RB ranges from a lower 5th percentile of 2.27 to an upper 95th percentile of 5.0, with amedian of 3.36. The length ratio RL has a lower 5th percentile of 1.0 and a 95th of 2.6, witha median of 1.6, implying that streams of a higher order are less than twice as long as theirlower-order tributaries. The area ratio, RA ranges from a lower 5th percentile of 2.1 to anupper 95th percentile of 4.8 with a median value of 3.1. No obvious spatial patterns stand outin any of these metrics (Figure 5.6).935.4. Results  'LVWDQFHNP$  &DWFKPHQWVRB  'LVWDQFHNP%  RL  'LVWDQFHNP&  RAFigure 5.6: Morphometrics of the channel networks, including A) bifurcation ratio, B)length ratio, and C) area ratio. Values are stretched between the 25th and 75th per-centiles to accentuate differences. Only moulin-drained catchments are shown, ascrevasse-drained catchments were unchannelized and therefore network morpho-metrics could not be calculated.The breakdown of the ratios and their components according to stream order is shownin Figure 5.7. N and A¯ vary exponentially with stream order, typical of terrestrial rivers. L¯varies very little with increasing stream order until exceeding stream order 4. L¯ dips at themaximum observed stream order because moulins truncate the length of the highest orderstreams, as observed by Yang and Smith (2016). The highest order stream was therefore notused in the calculation of RL. RB and RA do not vary much between stream orders for streamorders 1-3, however RA increases consistently between stream orders.945.4. ResultsN$L P&  6WUHDP2UGHUA P2 H(R B%R LP'    6WUHDP2UGHUR A)Figure 5.7: The A) number, C) average length and E) average contributing area of streamsin catchments, and the corresponding B) bifurcation, D) length and F) area ratiosbetween stream orders. The shaded area represents the range between the 25th and75th percentile of the values in all of the catchments, and the black line representsthe median.The lack of spatial patterns in bifurcation, length and area ratios begs the question whetherthese catchment morphometrics are random or reflect local conditions. To test whether, likeRe, these morphometrics might be influenced by ice sheet characteristics, I explored their rela-tionship to ice thickness (Hice), ice velocity (vice), and strain rate (e˙) using principal componentanalysis (PCA). Although data normality is not a requirement for PCA, data that were skewedwere log transformed to normalize them as much as possible. Results are shown in Figure 5.8.PCA loadings are shown by the arrows, indicating how the independent variables con-tribute to each of the two principal components shown, and how they vary with respect toeach other. The PCA scores are shown by the points, indicating where individual data pointsplot in the reduced dimensions of PC-1 and PC-2. RB, RL and RA do not appear to have anystrong relationship to ice properties. The loadings for each ratio are strongly correlated withthe second principal component, but are entirely orthogonal to any of the ice characteristic955.4. Resultscoefficients, which are very strongly correlated to the first principal component. There is alsono apparent clustering of the data points nor stratification according to catchment type.        3&3& RevHiceRBRLRA2WKHU0RUDLQHSUR[LPDO&UHYDVVHSUR[LPDO/DNHGUDLQDJHFigure 5.8: The first two components (PC-1 and PC-2) of Principal Component Analysisof RA, Re, RL, RA and ice characteristics ice thickness (Hice), ice velocity (vice) andprincipal strain rate (e˙). PCA loadings are indicated by the arrows, and demonstratethe co-variability of the variables amongst each other and their correlation with eachof the principal components (axes). The scores are shown by the points and arecolored based on the catchment type.The total volume of melt delivered to moulins or crevasses in the study area in 2015amounted to 2.87 km3, equivalent to an average runoff production of 1.33 m/yr throughoutthe study area. The distribution of discharge is spatially and temporally uneven. Becausecatchment discharge depends on both annual melt production and catchment area, large an-nual discharges are found at both low elevations (where melt is highest) and at high elevations(where melt is low but large catchments are present) (Figure 5.12 A).Time to peak discharge varies significantly. Some catchments have near instantaneous hy-965.4. Resultsdrograph peaks; the smallest tp is 0.15 hr for a crevasse-drained (unchannelized) catchmentand 0.7 hr for a moulin-drained (channelized) catchment. Other catchments have tp on theorder of several hours, with the slowest hydrograph peak taking 11.6 hrs to reach the moulin.Again, I stress that these numbers should, in the absence of suitable empirical data, be con-sidered in relative terms rather than in absolute terms, and also apply only to snow-free con-ditions. Time to peak appears to be higher at elevations above 1500 m as well as along thenorthern high velocity outlet glacier (Figure 5.12).  'LVWDQFHNP$&DWFKPHQWVQytdNP3    'LVWDQFHNP%tpKU Figure 5.9: A) total annual catchment discharge and B) the time to peak for each catch-ment. To accentuate differences and limit the visual effects of outliers, the color mapextent is stretched between the 25th and 75th percentiles of the data.I examined the interrelationship between catchment morphometrics using Principal Com-ponent Analysis (Figure 5.10). Peak time and catchment area are best explained by the firstprincipal component (and implicitly with each other), and SS and SC are best explained bythe second principal component. Bifurcation ratio is weakly correlated with tp, which appearsweakly and inversely correlated to slopes and area ratio. No grouping or clustering is apparentin the scores, even when considered by moulin type.975.4. Results      3&3&AReSsSc RBRLRAtp2WKHU0RUDLQHSUR[LPDO&UHYDVVHSUR[LPDO/DNHGUDLQDJHFigure 5.10: The first two principal components of a PCA analysis on the inter-variabilitybetween catchment area (A), time to peak (tp), bifurcation ratio (RB), average an-nual catchment melt, surface slope (SS), channel slope (SC), elongation ratio (Re),area ratio (RA) and length ratio (RL).5.4.2 Partitioning of runoffMeltwater is unequally partitioned between catchment types. Moulins drained 61% of thesurface runoff, with ’other’ type moulins dominating channelized drainage (30%), followedby crevasse-proximal moulins (14%), lake-drainage associated moulins (10%), and moraine-proximal moulins (6%). Crevasse (unchannelized) drainage accounted for 39% of surfacedrainage (Figure 5.11).985.4. Results&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO/DNH'UDLQDJH2WKHU &UHYDVVH'UDLQDJH4ytdP3 \UHFigure 5.11: Partitioning of total annual surface drainage between moulin types andcrevasse drainage. * indicates moulins, in contrast to un-channelized crevassedrainage.Although crevasses drained the most surface water of any of the catchment types, theydid so through multiple smaller catchments, each one draining less water on average thanmoulin-drained catchments do (Figure 5.12). Lake drained catchments, although they drainedrelatively little of the surface melt, did so in larger catchments than did any of the other moulintypes. Crevasse-proximal, ’other’-type, moraine-proximal, and lake-drainage moulins havemedian annual discharges of 1.58× 106, 2.06× 106e, 2.37× 106 and 2.87× 106, respectively.Conversely, more than half of crevasse catchments have total annual discharges of less than2.48e5 m3, equivalent to just 8.6% of the median lake-drainage associated moulin catchments,for scale.995.4. Results&DWFKPHQW4ytdP3\U)UHTXHQF\&UHYDVVH3UR[LPDO0RUDLQH3UR[LPDO/DNH'UDLQDJH2WKHU&UHYDVVHVFigure 5.12: The frequency of total annual catchment discharges distributed by catchmenttype. * indicates moulins, in contrast to un-channelized crevasse drainage.As a result of the spatial distribution of catchment types relative to elevation, there is aslight temporal offset in when different catchments contribute to the annual distribution ofsurface melt drainage (Figure 5.13). Drainage through crevasses, which are located at lowerelevations, begins and ends the melt season. ’Other’ type moulin-drained catchments beginto contribute to the annual hydrograph just a few days after, and have overtaken crevassesas the largest source of meltwater a week later. Contributions from lake-drainage catchments,in which drainage is initiated only when accumulated melt exceeds the local ice thicknessmultiplied by an Fa of 4000 m2, are delayed by almost 20 days relative to crevasse and ’other’moulin drained catchments.1005.4. Results     'D\QsumP2 H/DNH'UDLQDJH2WKHU0RUDLQH3UR[LPDO&UHYDVVH'UDLQDJH&UHYDVVH3UR[LPDOFigure 5.13: Partitioning of total annual catchment hydrograph between moulin typesand crevasse drainage, based on daily RACMO2.3 data. Lake-drainage dischargeis calculated on a catchment by catchment basis, such that lakes open and drain theaccumulated melt (and thereafter the daily melt) when the volume of accumulatedmelt exceeds local ice thickness multiplied by a fracture area threshold equivalentto 4000m2. * means moulins, in contrast to un-channelized crevasse drainage.5.4.3 Connections to subglacial catchmentsI employ four different flotation fraction (k) values to evaluate the spatial configuration ofsubglacial catchments, ranging from 0.8 to 1.1. The k value appears to substantially affectthe spatial configuration of subglacial catchments (Figure 5.14). At lower flotation fractionvalues, subglacial catchments are dominated by four large basins. The predominance of largecatchments decreases with increasing k. Interestingly, however, this is not simply because ofthe addition of new catchments. The number of subglacial catchments increases from 283 tojust 304 between k = 0.8 to k = 1.1. The difference in configuration appears to come from thegrowth of peripheral catchments that expand inland and erode the larger catchments, reducingtheir extent. As observed by Chu et al. (2016), the difference in spatial organization becomesparticularly pronounced above k = 1 (Figure 5.14).As a result of this reorganization of subglacial catchments, the number and distributionof moulins feeding individual subglacial catchments changes. At lower k values, subglacialcatchments are fed by an average of 60 moulins each. That number drops to 43, 30, and 16 fork values of 0.9, 1.0, and 1.1, respectively. At k = 0.8, the largest number of moulins feeding a1015.4. Resultsgiven catchment is equal to 95. This number is reduced to 37 when k = 1.1.1025.4. Results$  &DWFKPHQWV0RXOLQV%  0RXOLQV  'LVWDQFHNP&  &DWFKPHQWV  'LVWDQFHNP'  Figure 5.14: The spatialconfigurations of sub-glacial catchments underdifferent k values wherek = A) 0.8, B) 0.9, C) 1.0and D) 1.1. Histogramsrepresent the numberof moulins feeding thedifferent subglacial catch-ments.1035.4. ResultsABasinP2&DWFKPHQWVFXPXODWLYHN Q N Q N Q N Q Figure 5.15: The cumulative percentage of subglacial catchments of different areas forfour different values of k.As the subglacial catchments are reconfigured with increasing k, the distribution of thevalues of supraglacial tp within each subglacial catchment changes (Figure 5.16). The im-pact is particularly apparent spatially. At low values of k where subglacial drainage is domi-nated by several particularly large subglacial catchments, average subglacial tp from contribut-ing surface catchments is fairly uniform over the interior of the study area. As k increases,surface catchment meltwater contributions become more variable for sub-glacial catchments.Supraglacial catchments with longer hydrograph delays become concentrated, and the maxi-mum average subglacial catchment tp increases from 5.04 hr for k = 0.8 to 6.9 hr for k = 1.1 asthe number of subglacial catchments connected to supraglacial catchments through moulinsincreases from 79 to 131.1045.4. Results$   &DWFKPHQWVtp(hr)%   tp(hr)  'LVWDQFHNP&   &DWFKPHQWV  'LVWDQFHNP'   Figure 5.16: The averagevalue of contributingmoulin tp for sub-glacialcatchments under fourdifferent k values where k= A) 0.8, B) 0.9, C) 1.0 andD) 1.11055.4. ResultsSimilarly, the amount of water being supplied to subglacial catchments varies with dif-ferent values of k. As k increases, subglacial catchments become more fragmented and thesubglacial hydrology becomes more sensitive to individual moulin inputs, the median annualsubglacial catchment water supply decreases from 6.19e5 m3 to 1.04e5 m3. Spatially, dischargeto subglacial catchments becomes more distributed, with most catchments having lower totalannual discharges (Figure 5.17).1065.4. Results$ &DWFKPHQWVQyrP3\U% QyrP3\U  'LVWDQFHNP& &DWFKPHQWV  'LVWDQFHNP' Figure 5.17: The sum ofannual discharges formoulins connected tosub-glacial catchmentsunder four different kvalues where k = A) 0.8,B) 0.9, C) 1.0 and D) 1.11075.4. ResultsThe reconfiguration of subglacial catchments has a modest effect on average contributingGIUHs (Figure 5.18). The number of catchments increases with increasing k, and the averagecontributing shape becomes more variable and heterogeneous. This can be attributed to theinward growth of peripheral subglacial catchments which become connected to more surfacedrainage conduits, and the influence of individual moulins becomes more pronounced.TKU1 $ Q  % Q     WKUTKU1 & Q     WKU' Q Figure 5.18: Each line represents the average of all the GIUHs of the moulins contributingto a given subglacial catchment for four different values of k, where k = A) 0.8, B)0.9, C) 1.0 and D) 1.1.Figure 5.19 summarizes how the distribution of surface water delivery to subglacial basinschanges with changing k values. At low k values, supraglacial catchments typically have lower1085.5. Discussionaverage lag times for discharge contributed from connected supraglacial catchments (Figure5.19A). Concurrently, an increase in k increases the total amount of water supplied to indi-vidual subglacial catchments on average (Figure 5.19B). Therefore, at low k values most sub-glacial catchments are small, located at the periphery of the ice sheet, connected to few surfacemeltwater conduits, and supraglacial drainage is concentrated in just a few large basins. Assubglacial catchments expand inland under high assumed flotation fraction conditions, theybecome connected to more moulins, and individual moulin drainage becomes more importantfor individual subglacial catchments.tpKU&XPXODWLYH&DWFKPHQW)UHT$Q Q Q Q Q Q Q Qyr %Figure 5.19: Changing A) average tp and B) total subglacial annual discharge for eachsubglacial basin, as determined by the connected supraglacial basins. Each line inthe above plots pertains to a different value of k.5.5 DiscussionIn this chapter, I consider the extent to, and the ways in which supraglacial catchment mor-phometrics vary across a study area in southwest Greenland, and how the spatial distributionof these characteristics, relative to subglacial catchment configuration, might impact regionalhydrograph characteristics.1095.5. Discussion5.5.1 LimitationsThe work presented in this Chapter is predicated on several assumptions and limitations thatpropagate through the analysis. First, the simple identification of moulins and catchmentsis made based on visual observation rather than a specific set of quantitative criteria. As aresult, small catchments are very likely to have been missed. Additionally, preserving all ofthe crevasses as sinks assumes that meltwater at those locations reaches the bed, which maynot be true. However, I note that models of moulin formation use thresholds of stress andstrain to indicate crevasse formation, and assume moulin formation at those locations (Banwellet al., 2016; Koziol et al., 2017; Poinar et al., 2015). Our assumption of hydrofracture in crevassefields is therefore consistent with contemporary modeling approaches.The fixed channel initiation threshold introduces additional constraints. Although fixinga low channel initiation threshold is reasonable under the assumption of bare ice conditionsand maximum channel development (resulting in high drainage density), I recognize that thechannel initiation threshold has an effect on the derivation of the GIUH curves (Chavan andSrinivas, 2015). This consideration, in addition to the use of empirical parameters derived foralluvial catchments for deriving the GIUHs, warrants restating that the hydraulic parametersderived in this analysis (velocity, time to peak, GIUH curves) should be interpreted as relativerather than absolute. Despite these limitations the derived tp and v values are well withinreported ranges for supraglacial channels (Gleason et al., 2016; McGrath et al., 2011; Smith et al.,2017).Future work in this area must address the limitations highlighted above. Better constraintson channel network extent are necessary, and promising multispectral work has been donein this vein by Smith et al. (2015) and Yang et al. (2016a). However, multispectral approachesare limited to a snapshot in time and do not accurately depict the high drainage density ofthese systems (King et al., 2016; Smith et al., 2015). New approaches are needed that allowfor modeling the co-evolution of the channel networks and the ablating snowpack on spatialresolutions that allow discrimination of low order channels (e.g. at finer resolutions than the30m ASTER DEM used by Banwell et al., 2013).5.5.2 Catchment and channel network characteristics and variabilityCatchments in this study region vary in predictable ways that are observed elsewhere on theice sheet. Small catchments dominate at lower elevations, and larger catchments exist pre-dominantly at higher elevations (Yang and Smith, 2016; Yang et al., 2016a). Because catchmentdischarge is a function of both its contributing area and its annual melt, the emerging spatialpattern of total catchment discharge is complex: higher annual discharges are observed bothat low elevations where melt rates are higher as well as at higher elevations where melt ratesare lower but larger catchments are present. The difference in annual discharge between highand low discharge catchments is almost three orders of magnitude.The time to peak, however, which to a large extent determines the daily hydrograph, dis-1105.5. Discussionplays a somewhat more consistent spatial pattern, with longer tp above 1500m elevation, inthe north east corner of the study area, and along the northernmost outlet glacier. Our valuesof tp are well within the range of those determined by Smith et al. (2017) based on SUH theory,where tp generally ranged from 2 to 8 hours. Further, the spatial pattern of higher tp at higherelevations mirrors the pattern they observed in the Kangerlussuaq/Rusell Glacier region.The upper limit of the range of tp values reported herein (15 hr) exceeds the upper limitof Smith et al. (2017)’s SUH estimates. I suggest that the higher peak times I observe inlarge catchments may be attributed to the inclusion of channel morphometrics in hydrographderivation, an aspect of supraglacial hydrograph derivation that has not been previously ad-dressed. Larger catchments with correspondingly longer peak travel times are typically lo-cated at higher elevations, where channel slopes are also lower. Including channel morpho-metrics in the derivation of the hydrograph allows for inclusion of catchment properties thatcontribute to a wider range of catchment responses to a given hyetograph.The bifurcation ratios in our study area, which range from just over 2 to just over 5 fall wellwithin the range of alluvial values which typically range from around 2 up to 3 or 4 (Horton,1945). With a median RB of 3.4, supraglacial channels in our study area are on the higherend of what is observed in alluvial systems, and indicate highly dissected drainage basins, aswould be expected with the high drainage density of these systems. Additionally, our rangeof catchment bifurcation ratios is almost identical to the range of RB values calculated by Yangand Smith (2016) using multispectral delineation of supraglacial catchments. They report amean RB of 3.7 ± 1.9.The length RL and area RA ratios observed in our study region are on the low end of therange typically observed in alluvial systems. Although highly variable in alluvial systems,the length ratios typically range from 1.5 to about 3.5, and the area ratios range from 3 to 6(Rodriguez-Iturbe and Valdes, 1979). The median RL and RA values are 1.6 and 3.1, respectively.These low values suggest that supraglacial channels of a given order are not substantiallyshorter than channels of the next highest order, and that area does not increase as rapidlydownstream as it does in alluvial channels.Given the high drainage densities observed in supraglacial catchments (Smith et al., 2015)I attribute these low RL values to low piracy rates - channels parallel each other for a signif-icant way before joining, resulting in long but frequent low-order channels. This hypothesisis consistent with our qualitative observations of long, incised, parallel low order channel net-works, as well as field observations of regularized channel spacing (Karlstrom et al., 2014), andnumerical analysis of supraglacial channel inception (Mantelli et al., 2015). Finally, I note thatas in Figure 5.10, RL increases with channel slope, suggesting channel piracy increases withincreasing slope, as expected as per Mantelli et al. (2015). The hydrological implication of thisobservation is that flow remains in low-order channels for a larger proportion of its traveltime than would be expected in a typical alluvial basin. In low order channels, flow velocity islikely to be lower than it is higher order channels, and the hydrograph peak therefore delayed1115.5. Discussion(Gleason et al., 2016; Marston, 1983).Interestingly, despite the length of the channel segments, area does not appear to increasewith channel order to the same extent one would expect in a typical alluvial river. I attributethis to the manner in which supraglacial catchments are elongated; because they are narrowerthan they are wide, a given increase in distance downstream does not produce a rapid increasein contributing area.5.5.3 Partitioning of supraglacial dischargeCrevasse drainage is the most important source of meltwater to the subglacial hydrologicalsystem in our study region, similar to findings by Koziol et al. (2017). Whereas Koziol et al.(2017) reported that crevasse drainage amounted to 51.5% of surface drainage, I find the con-tribution to be just 39%. The discrepancy may certainly be attributed to different regionalcharacteristics of the study area. However, I also suggest that the use of a stress threshold forcrevasse initiation, as employed by Koziol et al. (2017), does not accurately capture the dynam-ics of surface-to-bed connections. Given that channels and catchments can occur in localizedstrain rates that exceed those predicted to initiate crevasse formation, it is likely that Koziolet al. (2017) overestimate the contribution of crevasse drainage.Koziol et al. (2017) observed that 85% of channelized drainage was attributed to lake-basinmoulins. Their findings differ substantially from ours, where I observe that lake drainageaccounts for just 17% of moulin drainage, and that ’other’ type moulins are the most importantcontributor of channelized meltwater. Again, this may be attributed to regional differences.However as I discuss in Chapter 4, this difference may also arise from their use of modeledrather than observed lake-drainage moulins.The partitioning of drainage between catchment types has both spatial and temporal hy-drological implications. Crevasse drainage occurs in smaller catchments, resulting in flashierhydrographs, more spatially distributed inputs, and, due to the concentrated presence ofcrevasses at lower elevations (Joughin et al., 2008; Poinar et al., 2015), earlier annual meltwa-ter supply to subglacial networks. Channelized surface drainage, in contrast, is higher involume, spatially concentrated due to larger catchment areas, and temporally delayed bothannually and diurnally relative to crevasse drainage. However, based on field measurementsof a moulin water budget, McGrath et al. (2011) suggest that crevasse drainage may have the ef-fect of dampening diurnal hydrographs due to temporary crevasse storage. It may be that thenet effect of crevasse drainage later in the season is to provide more sustained water deliveryto the subglacial network.Banwell et al. (2016) consider the implications of moulin density on subglacial channel evo-lution using a numerical model of lake drainage events. They test moulin densities rangingfrom 0.02 km−2 to 0.08 km−2, and infer that higher moulin densities lead to earlier and moreextensive subglacial channelization. They are unable to reconcile the timing of the onset oftheir subglacial channelization with dye-tracer experiments by Cowton et al. (2013) and Chan-1125.5. Discussiondler et al. (2013) who observe subglacial channelization four weeks earlier than the model pre-dicts.In our empirical study, I observe moulin densities of 0.1 km−2 and total bed to surfaceconduit densities (i.e. including crevasse drainage) of 0.38 km−2 (although I note that crevassedrainage is highly spatially concentrated in crevasse fields). I note that non-lake basin moulinsand crevasse drainage begin supplying meltwater to the bed several weeks before lake drainagebegins to be important, which may explain the earlier observed onset of channelization. ’Other’type moulins are not included in physically based models of surface-to-bed meltwater deliv-ery (Banwell et al., 2013; Clason et al., 2012), and these models therefore cannot fully capturethe temporal dynamics of meltwater drainage. Further, the large quantities of surface wateralready delivered to the bed through ’other’ type moulins and crevasses before the onset oflake drainage could create the necessary conditions for rapid lake-drainage as speculated byAndrews et al. (2014) and Hoffman et al. (2018).5.5.4 Implications of surface catchment hydrology for subglacial hydrologyThe volume of water as well as its seasonal and diurnal variability are the main ways in whichsupraglacial hydrology may impact subglacial channel development. In our study area, therange in supraglacial catchment areas creates orders of magnitude differences in catchmentdischarges and moulin spacing, creating a mosaic of low-discharge catchments interspersedwith large, high-discharge catchments. Because tp is dependent partly on flow length, largercatchments typically have more sustained and delayed hydrographs. Further, the pattern oftp broadly but not entirely mirrors the spatial pattern of discharge.Schoof (2010) numerically demonstrated that rapid inputs of large volumes of water or sub-stantial early season inputs overwhelm the capacity of established drainage networks. Thisdecreases effective pressure, leading to bed separation and ice acceleration. Conversely, thesubglacial hydrological network will evolve gradually with sustained water inputs towards ef-ficient channelization. In our study area, early season inputs are dominated by fast tp crevassedrainage, quickly being overtaken by more sustained and higher volume ’other’-type moulininputs, long before lake drainage inputs are initiated.Clearly, capturing the impact of the density of surface-to-bed connections requires con-sidering the spatial complexity introduced by a range of catchment types, and is particularlyimportant when considered relative to the spatial configuration of subglacial catchments. Theflotation fraction, k, used to delineate sub-glacial catchments using a hydraulic potentiometricsurface represents the ratio of water pressure to ice overburden pressure (Shreve, 1972). Al-though over the long term an assumption of steady state may justify an average conditionof k = 1, this assumption does not hold for the timescales over which the supra-\sub-glacialsystems co-evolve during a melt season. Existing models of supra- and sub- glacial catchmentconnections use a temporally constant and spatially uniform value of k. Banwell et al. (2016)and Banwell et al. (2013) constrained a k value based on pro-glacial discharge measurements1135.6. Conclusionfrom May to August of their study year, however current coupled supra-\sub-glacial mod-els do not allow for temporally variable subglacial catchments. On the periphery of the icesheet where moulins are dense and subglacial catchment delineation is sensitive to k valueselection, this appears to significantly affect spatial coupling between specific supra- and sub-glacial catchments.The hydrological effect of variable supra- to subglacial catchment coupling through a meltseason would be that some subglacial catchments would become isolated from high discharge,long tp supraglacial catchments with sustained hydrographs, whereas others, particularly atlow elevations, would become connected to such catchments and would transition from aflashier hydrograph to one that is more sustained. Because the spatial organization of thecatchments is particularly sensitive to high subglacial water pressure (high k), the ways inwhich surface meltwater delivery spatially and temporally affects k will introduce a dynamicfeedback between subglacial catchment evolution and supraglacial hydrology that warrantsfurther study.5.6 ConclusionThis chapter considers the extent to which supraglacial catchment morphometry varies inspace, and the ways in which this variability might influence both supra- and sub-glacialcatchment hydrology. Moulins and crevasses were visually identified, and their contributingcatchments were derived using traditional flow routing techniques. Catchment morphomet-rics were similarly extracted and used to build Geomorphic Instantaneous Unit Hydrographcurves to extract time to peak and hydrograph shape. The spatial distribution of catchmenthydrology and its partitioning between catchment types was considered throughout.Catchment hydrology and morphometrics are variable in ways that are not captured bycontemporary models of supra- to sub-glacial hydrological connection. At the same time,I note that the connections between supra- and sub-glacial basins vary significantly as sub-glacial catchments are delineated using a range of realistic k value parameters. To what extentis this variability glaciologically meaningful? Fully modeling the complexity of these systemsmay require more research and computational investment than necessary from an ice dynam-ics perspective, depending on the importance of this variability. Future work on modeling theinterplay between spatially variable catchment hydrographs and subglacial catchment con-figuration would constrain field efforts required to build a more complete picture of the co-evolution of surface and subglacial hydrological systems.114Chapter 6Concluding remarksSince beginning this PhD in September, 2012, the Greenland Ice Sheet has lost 1401 Gt of mass(Tedesco et al., 2016). Global mean temperature has risen a shocking 0.28◦ C and sea levelshave risen 20.5 mm (Figure 6.1), with around a quarter of sea level rise being contributedby Greenland’s melt (Chen et al., 2017). Beyond approximately 54 tonnes of CO2 equivalentemitted to the atmosphere, what contributions have I made to Greenland’s future during thisPhD?     7HPSHUDWXUH$QRPDO\&6HD+HLJKW9DULDWLRQPP$ %Figure 6.1: The global temperature anomaly relative to 1951-1980, and global sea heightvariation relative to the 20-year mean. The grey bar represents the 6 year spanof this PhD, and ’A’ and ’B’ reference the PhD journey, summarized in Ap-pendix E. Temperature data accessed at . Sea level data accessed at Both datasets accessed on July 19, 2018.In this final chapter, I reflect on the contributions I have made towards improving bothour understanding of the Greenland Ice Sheet’s surface hydrology as well as towards refiningthe ways in which we learn about the system. I summarize the unavoidable but important1156.1. Summary of contributionsassumptions and limitations that are integrated into the research. Finally, from my lessonslearned, I offer my view of where this research should go next.6.1 Summary of contributionsFlow routing is an effective method for extracting supraglacial channels from highresolution DEMs.When I began to explore methods to extract supraglacial catchments, there were disparate ap-proaches used in the literature, including visual identification (Das et al., 2008; Joughin et al.,2013; McGrath et al., 2011; Wyatt and Sharp, 2015), multispectral delineation (Lampkin and Van-derBerg, 2014; Smith et al., 2015; Yang and Smith, 2016), and elevation thresholds (Rippin et al.,2015). Although flow routing was routinely used in hydrological modeling, there was an ap-parent distrust of its utility for delineating channels themselves (Rippin et al., 2015; Yang et al.,2015). However, flow routing offers many advantages over multispectral methods, includingeasy integration into contemporary surface hydrology models, the ability to capitalize on thenewly released 2 m resolution ArcticDEM (as well as many other advances in photogramme-try and DEM construction), easy extraction of catchment morphometrics (as done in Chapter5), and relative temporal insensitivity compared to multispectral methods.In my second chapter, I demonstrate that traditional flow routing techniques are in factremarkably effective at delineating supraglacial channels. I verify their effectiveness by com-paring extracted flow routed channels to an independent dataset of supraglacial channelsderived using multispectral analysis. In having verified the utility of this method, I hopeto have equipped researchers with the necessary confidence to use flow routing for furtherexploration of supraglacial hydrological systems. ArcticDEM is increasingly being used inglaciological work (e.g. Barr et al., 2018; Ross et al., 2018) and commercial imagery at 0.3 mresolution is available from WorldView-3 as of August 2014, just as new open source meth-ods in photogrammetry using satellite imagery are becoming widespread (e.g. Shean et al.,2016). Catchment morphometric analysis is well poised to be better integrated into the glacialhydrologist’s toolkit.Moulin discharge estimates are limited by the spatial and temporal resolution of themethodology used.In the third chapter, I reconcile multispectral estimates of moulin discharge to estimates de-rived from the intersection of a commonly used surface mass balance model with flow-routingderived catchments. This work is motivated by a generally articulated need for better repli-cation of scientific work (Peng, 2011), as well as an interest in how different methodologicalapproaches to quantifying supraglacial hydrology compare. I observe a scale dependent dis-crepancy between approaches which, in small catchments, results in differences in moulindischarge estimates of several orders of magnitude. I attribute this difference to the combined1166.1. Summary of contributionseffect of sensor resolution limitations in the derivation of channel widths from multispectralimagery as well as the inadequacy of synthetically derived hydrographs (using available pa-rameterization) for capturing the flashiness of runoff in small supraglacial systems.Although advances in catchment hydrograph derivation are being rapidly made (Smithet al., 2017), this chapter highlights the incongruities between the spatial and temporal scalesover which the methods apply. Whereas multispectral methods capture an accurate snapshotin time, that snapshot cannot readily be extrapolated over a diurnal melt cycle. Conversely,whereas synthetic unit hydrograph theory facilitates construction of continuous hourly hydro-graphs, it does not capture the flashiness of instantaneous discharge measurements.Localized ice velocity, strain rate, lake proximity and bed surface slope all affect the spatialdistribution of moulins. Models should move beyond a focus on moulin formationprimarily driven by meltwater threshold-volume hydrofracture and consider how icedynamics contribute to the conditions necessary for moulin formation.Accurate projections of moulin formation are crucial for predicting non-linear ice dynamicsunder climate change induced expansion of ablation zones (Leeson et al., 2015; Poinar et al.,2015). However, modeling approaches have been generally unsuccessful at capturing the com-plex range of conditions that create moulins, and rely on simplified assumptions of hydrofrac-ture induced by a threshold water volume. The complex reality of where and why moulinsform remains an open question. In my fourth chapter, I empirically quantify the ways in whichhydrology and ice sheet characteristics collectively produce the spatial distribution of moulinsin a study area in south west Greenland. Whereas previous work has attempted to collapsecomplexity and to look for statistical trends, I explicitly consider the spatial nature of moulinfeatures to explore what complex dynamics might not be captured by prior work. I observethat although ice characteristics in aggregate do not explain trends in moulin distribution,localized ice conditions do exhibit some controls on moulin formation.The highest densities of moulins appear to be associated with high strain rates within lowvelocity ice. This association is particularly apparent for ’other’ type moulins, for which for-mative mechanisms have thus far been omitted by models despite their important contribu-tions to ice sheet hydrology (see Chapter 5). The importance of low velocity ice to moulinformation is not entirely clear, but may be attributed to larger incidents of crevassing in high-velocity ice which inhibit surface water channelization. We observe that ’other’ type moulinsare generally closer to moulin-drained lakes than other moulin types, and that moulins locatedclosest to lakes that drain in-situ are typically located downstream of these lakes. Collectively,we interpret those findings as supporting the hypothesis that moulin formation is associatedwith short term ice flexure during rapid surface drainage events (Andrews et al., 2014; Hoffmanet al., 2018). This temporally transient process may explain why identifying formative pro-cesses for ’other’ type moulins is difficult to do, and why temporally averaged strain rates donot appear to correlate with moulin locations.Additionally, we observe that longitudinal ice processes may be important in determining1176.1. Summary of contributionsmoulin distribution, particularly with respect to the slope of the ice surface relative to the slopeof the bed. Where these two slopes parallel one another and ice flow is governed only by thetopographic gradient of the bed rather than the ice thickness gradient, strain rate is close to 0yr−1 and moulin formation appears to be inhibited. Conversely, in areas where bed slopes areadverse to surface slopes, moulin formation appears to be enhanced. This may be due to thepresence of splaying crevasses as ice compresses, and would explain the observed occurrenceof crevasse-proximal moulins in compressive ice.Moulin formation is driven by physical processes, however identifying the signature ofthese processes from gridded data is difficult as it is non-unique: different processes can causemoulins to form under different conditions. These empirical insights advance our understand-ing of moulin formation away from a focus on water volume threshold induced hydrofracturetowards an emphasis on how ice dynamics themselves produce the spatial configuration ofmoulins.In the broader context of moulin research, this work provides limited additional insightinto the processes at work in alpine glaciers where ice thicknesses are lower (Gärtner-roer et al.,2014) and bodies of ice do not typically override underlying topography. However, the spatialscale of this analysis and the data used herein as well as its focus on processes that are commonon the GrIS makes its findings useful for understanding moulin location in the context of theAntarctic Ice Sheet. Increasing surface melt in Antarctica is leading to heightened interest insurface-to-bed connections there (e.g. Pollard et al., 2015).Catchment characteristics significantly impact the spatial and temporal delivery of waterto the subglacial system. Understanding the importance of this variability on ice dynamicsrequires better constraints on the co-evolution of the supra- and sub-glacial hydrologicalsystem.With the spatial distribution of supraglacial catchments set by the distribution of moulins,the final chapter considers the hydrological implications of variable catchment properties.Again adopting an explicitly spatial and empirical approach, we extract hydrologically rel-evant catchment morphometrics using flow-routing derived channel networks. We use thesemorphometrics to estimate hydrograph parameters for each of these catchments, and considerhow the spatial variability thereof might influence the spatial and temporal delivery of waterto subglacial catchments. We find that the combined influence of catchment length, channelslope and depth, and the bifurcation, length and area ratios create a range of hydrograph lagtimes that range from under an hour to more than 11 hours. Importantly, our findings are onlyapplicable under snow free conditions (i.e. peak melt), and the characteristics of catchmenthydrographs could be expected to vary seasonally (Yang et al., 2018).It is difficult to know how important diurnal-timescale hydrology is to ice dynamics un-less we can understand its effect on the evolution of the subglacial network. Assuming thatmoulins delivery water vertically to the bed of the ice sheet, the importance of this dailysupraglacial hydrological variability is highly dependent on subglacial catchment configura-1186.2. Limitations, assumptions, and future worktions. Contemporary models of surface-to-bed-connections model subglacial catchment con-figurations using a single and temporally fixed flotation fraction, and link moulin meltwaterinputs to these fixed catchments. However, changing the flotation fraction changes the spatialconfiguration of subglacial catchments. As peripheral subglacial catchments expand inlandunder higher k value conditions, interior subglacial catchments become connected to fewerindividual moulins, and as a result the hydrological characteristics of individual moulins be-come relatively more important for determining the hydrology of subglacial catchments. Con-versely, when the water pressure is relatively low, subglacial drainage is dominated by fewer,larger subglacial catchments with multiple surface water inputs that dampen the importanceof individual supraglacial catchment characteristics.This could mean, for example, that as the subglacial channels evolve and subglacial catch-ment extents are rapidly reconfigured due to changes in water pressure throughout the meltseason, sub-glacial catchments could switch between being fed by rapidly-delivered low-discharge supraglacial catchments to delayed, high discharge catchments, and then back again,with significant implications for the daily evolution of the subglacial channel network. Higherk values are likely to be observed under high runoff production but before subglacial channel-ization, implying moulin meltwater delivery during peak runoff is concentrated in multiplesmaller inland subglacial catchments than during periods of lower k.6.2 Limitations, assumptions, and future workInescapably, the research in this thesis internalizes many assumptions and has many limita-tions. In this section, I highlight the most important of these, and consider their implicationsto my general findings. I suggest future work that would help to overcome these limitationsand build on the findings of this thesis.Identifying where channels beginChapter 2 aims to offer new methods for extracting supraglacial channels from remotely sensedimagery. One of the main limitations that emerges from this work is the use of a universalchannel initiation threshold for delimiting the extent of the channel network. This is a limita-tion that carries forward into Chapter 5, where derivation of the Geomorphic InstantaneousUnit Hydrograph relies on defining the channel network. A way forward with this methodthat overcomes, at least to an extent, the use of an arbitrary area threshold would be to trun-cate the up-channel extent of the network based on multispectral remote sensing of snowcover. Wherever snow cover is present on the ice, it can be assumed that channelization hasnot occurred. Multi-temporal multispectral imagery could be used over the course of the meltseason to simulate seasonal expansion of the channel network.1196.2. Limitations, assumptions, and future workBuilding catchment hydrographsChapter 3, which derives catchment hydrographs using synthetic unit hydrograph theory, isseverely limited by the use of empirical values in SUH derivations. The same is true for thederivations of GIUHs in Chapter 5. Although the derived channel hydraulics appear to bein line with field observations, albeit limited ones, the values reported herein cannot be in-terpreted as absolute. However, the variability in derived hydrograph characteristics betweenthese catchments still provides insight into the relative hydrological differences between them,and is valuable for refining where more data acquisition efforts should be focused.As in all research, empirical field data are always needed to support remote sensing andmodeling research. Empirical parameters employed in estimating the GIUHs and SUHs wouldclearly benefit from more collection of field data over a wide range of snow and ice conditions.Developing individual catchment budgets for even a limited number of moulins across a rangeof conditions and throughout the melt season would go a long way towards reducing the con-siderable uncertainty in estimating supraglacial catchment hydrological parameters. Perhapseven more temporally sustainable, field efforts to build long term proglacial discharge mea-surements across a wide spatial range of proglacial streams would provide a suite of benefitsto hydrological research, particularly if these might be coupled with several dye-tracer exper-iments from supraglacial catchments.Identifying moulinsChapters 4 and 5 employ a moulin and supraglacial catchment dataset that relies on visualidentification of moulins, introducing the potential (and likelihood) for errors in the configu-ration of supraglacial catchments. It is ironic that one of the main limitations in a study aboutidentifying moulin locations is the identification of moulin locations. However, before one canfigure out why something is where it is, one must figure out where it is in the first place. In thecase of moulins, this is easier said than done. There is no automated way to identify moulinlocations, and it has always been done manually, often assisted by assuming moulin forma-tion at the truncated end of a supraglacial channel (Colgan and Steffen, 2009; Koziol et al., 2017;Phillips et al., 2011; Smith et al., 2015). In this thesis, I also identify moulins visually, but I amfurther inhibited by having no access to satellite imagery of the study area that is concurrentwith the date of acquisition of the source imagery for the DEMs. The raw elevation data as wellas the derived curvature product show remarkably clear surface features, and many moulins,particularly those with large contributing streams, and clear. However, smaller moulins arealmost certainly missed, with the net effect of over-estimating the extent of moulin-drainedcatchments.At the end of this thesis, I do not yet see a way forward to practically identify moulinsbased on ice conditions or from automated delineation. Combining multispectral remote sens-ing of channels with flow routing might be a limited way forward. The abrupt termination ofchannels often denotes the presence of a moulin. If multispectral remotely sensed channels1206.2. Limitations, assumptions, and future workcould be georegistered with flow-routed channels, flow routed channels could be truncated atthe ends of multispectral channels. This would, however, still miss small moulins and involvea fair amount of manual checking. Additionally, identifying moulins in this way does not gen-erate insight into the controls on moulin formation, and does not contribute to our ability toproject moulin distribution into the future.Limited spatial coverageAn additionally important consideration in Chapters 4 and 5 is the limited extent of the studyarea. The study area is fairly typical of the south western ablation zone, however the resultspresented here would gain context from both an inland expansion of the study area as wellas comparisons with other study areas along a latitudinal gradient. Although it is well es-tablished that the relationship between surface melt and ice velocity is spatially variable (e.g.Moon et al., 2014; Stevens et al., 2016; Tedstone et al., 2015) there is very little, if any, synthesisof the variability of supraglacial catchment characteristics across the GrIS. Beyond our lim-ited study area, how variable are supraglacial catchments? In particular, how do they varybetween different ice conditions? Are the differences we observe here between high eleva-tion (i.e. above the 1500 m drop off in gradient), peripheral ice, and outlet glacier tonguesconsistent in space, and if so, what are the implications for how subglacial processes play outbeneath these different ice regimes?There is good ArcticDEM coverage for much of the western ablation zone. A simple latitu-dinal comparison could focus on the Russell and Jakobshavn regions, where good bed topog-raphy data are available, and many comparative studies have been done with respect to themeltwater–ice acceleration relationship (e.g. Andrews et al., 2014; Banwell et al., 2016; Joughinet al., 2010; Palmer et al., 2011; Zwally et al., 2002).Future directions in modeling and theoryOver the last 15 years, ice sheet researchers have made significant progress in understandingthe magnitude and direction of the relationship between surface melt and ice sheet accelera-tion (Flowers, 2018). The challenge now lies in digging down into the enormous complexity ofthese variable and dynamic systems and developing a theoretical framework that is consistentwith the observations and therefore has predictive power. Two useful modeling directionsemerge from this work. The first is to further consider the relationship between bed slopeand ice slope in determining moulin formation. An obvious question relates to the physicallymeaningful relationship between bed and surface slope that pre-conditions the ice for moulinformation. Further work might explore the importance of basal drag as an independent vari-able. Numerical models might move beyond meltwater threshold volume-drive hydrofractureby considering ice pre-conditioning due to basal drag, localized strain rate and nearby moulindrainage.As an interim step between knowing and not-knowing the specific controls on moulin1216.3. Final thoughtsformation (a task that remains persistently elusive, e.g. Williamson et al., 2018), the insightdeveloped in this thesis could be used to develop a probabilistic approach to predicting thelocation of ’other’ type moulins. This statistical approach could be used to move us ahead inunderstanding the importance of surface hydrology in lieu of more process-based approaches.The second modeling direction is to explore the relationship between surface catchmenthydrology and subglacial catchment configuration in a more deliberate way. This could bereadily done for the existing study area by building subglacial basin hydrographs under dif-ferent k configurations. Simple lag parameters could be calculated for supraglacial inputsbased on subglacial catchment configuration and water pressure, and a more complete pictureof the diurnal impacts of different feeder supraglacial catchment characteristics could be as-sessed. Additionally, more emphasis should be placed on the use of physically-based modelsof subglacial channel evolution (Koziol and Arnold, 2018).6.3 Final thoughtsIce sheet researchers hold tremendous public power derived from the importance of the icesheets to the public and political narrative of climate change (O’Reilly, 2017). Ice sheet researchis high profile, expensive to undertake, involves an unavoidable degree of personal risk (in thefield, at least), and is conflated with competitive and masculinist approaches to science (Careyet al., 2016). None of these factors are conducive to the careful and considerate slow sciencesome have called for in the contemporary context of academic research (Lane, 2017).These researchers are lucky to be so well resourced, and indeed we as a society facingrunaway climate change are also lucky that they are. However, along with this resource in-vestment and public scrutiny comes an academic need to occupy the frontier of knowledge.My meandering academic path has taken me through many earth science subdisciplines, andnever has my bibliography been so replete with Nature articles. Clearly, there is much tobe gained by way of grants, prestige, public interest and collaborations in moving quicklytowards the next high profile project, and much to be lost if you don’t. Although I do notquestion the value or rigor of the work emerging from the field (in fact I applaud it), I do mar-vel at the range of non-scientific choices that are internalized into research priorities in thisarea (King and Tadaki, 2018). I wonder what the scientific cost is; we’ve left a lot of ontologicalcorpses behind in our bid to render the ice sheets knowable and computable, and there’s stillmeat on the bones.122BibliographyAbraham, J. (2015), Global warming is causing rain to melt the Greenland ice sheet.Alexander, H. (1932), Pothole Erosion Author, The Journal of Geology, 40(4), 305–337.Andrews, L. C. (2015), Spatial and temporal evolution of the glacial hydrologic system of thewestern Greenland ice sheet : observational and remote sensing results, Ph.D. thesis, Uni-versity of Texas at Austin.Andrews, L. C., G. A. Catania, M. Hoffman, J. D. Gulley, M. P. Luthi, C. Ryser, R. L. Hawley,and T. A. Neumann (2014), Direct observations of evolving subglacial drainage beneath theGreenland Ice Sheet, Nature, 514, 80–83, doi:10.1038/nature13796.Appenzeller, C., T. Stocker, and M. Anklin (1998), North Atlantic Oscillation DynamicsRecorded in Greenland Ice Cores, Science, 282(October), 446–450.Arnold, N. (2010), A new approach for dealing with depressions in digital elevation modelswhen calculating flow accumulation values, Progress in Physical Geography, 34(6), 781–809,doi:10.1177/0309133310384542.Arnold, N. S., A. F. Banwell, and I. C. Willis (2014), High-resolution modelling of the seasonalevolution of surface water storage on the Greenland Ice Sheet, The Cryosphere, 8(4), 1149–1160, doi:10.5194/tc-8-1149-2014.Banwell, A., I. Hewitt, I. Willis, and N. Arnold (2016), Moulin density controls drainage de-velopment beneath the Greenland ice sheet, Journal of Geophysical Research: Earth Surface,121(12), 2248–2269, doi:10.1002/2015JF003801.Banwell, A. F., N. S. Arnold, I. C. Willis, M. Tedesco, and A. P. Ahlstrøm (2012), Modelingsupraglacial water routing and lake filling on the Greenland Ice Sheet, Journal of GeophysicalResearch, 117(1), doi:10.1029/2012JF002393.Banwell, A. F., I. C. Willis, and N. S. Arnold (2013), Modeling subglacial water routing atPaakitsoq, W Greenland, Journal of Geophysical Research: Earth Surface, 118(3), 1282–1295,doi:10.1002/jgrf.20093.Barr, I. D., M. D. Dokukin, I. Kougkoulos, S. J. Livingstone, H. Lovell, J. Malecki, and A. Y. Mu-raviev (2018), Using ArcticDEM to analyse the dimensions and dynamics of debris-coveredglaciers in, geosciences, 8(216), doi:10.3390/geosciences8060216.Benn, D., and D. Evans (2010), Glaciers and Glaciations, 802 pp., Hodder Education.123BIBLIOGRAPHYBjerklie, D. M., D. Moller, L. C. Smith, and S. L. Dingman (2005), Estimating discharge inrivers using remotely sensed hydraulic information, Journal of Hydrology, 309, 191–209, doi:10.1016/j.jhydrol.2004.11.022.Bjørk, A. A., L. M. Kruse, and P. B. Michaelsen (2015), Brief communication: Getting Green-land’s glaciers right - A new data set of all official Greenlandic glacier names, Cryosphere,9(6), 2215–2218, doi:10.5194/tc-9-2215-2015.Bracken, L., and J. Croke (2007), The concept of hydrological connectivity and its contributionto understanding runoff-dominated geomorphic systems, Hydrological Processes, 21, 1749–1763, doi:10.1002/hyp.Brykała, D. (1998), Evolution of supraglacial drainage on Waldemar Glacier (Spitsbergen) inthe period 1936–1998, Polish Polar Studies: 25th International Polar Symposium, pp. 247–261.Carey, M., M. Jackson, A. Antonello, and J. Rushing (2016), Glaciers , gender , and science : Afeminist glaciology framework for global environmental change research, Progress in HumanGeography, 40(6), 1–24, doi:10.1177/0309132515623368.Catania, G. A., T. A. Neumann, and S. F. Price (2008), Characterizing englacial drainage inthe ablation zone of the Greenland ice sheet, Journal of Glaciology, 54(187), 567–578, doi:10.3189/002214308786570854.Chandler, D. M., J. L. Wadham, G. P. Lis, T. Cowton, A. Sole, I. Bartholomew, J. Telling,P. Nienow, E. B. Bagshaw, D. Mair, S. Vinen, and A. Hubbard (2013), Evolution of the sub-glacial drainage system beneath the Greenland Ice Sheet revealed by tracers, Nature Geo-science, 6(3), 195–198, doi:10.1038/ngeo1737.Chavan, S. R., and V. V. Srinivas (2015), Effect of DEM source on equivalent Horton – Strahlerratio based GIUH for catchments in two Indian river basins, Journal of Hydrology, 528, 463–489, doi:10.1016/j.jhydrol.2015.06.049.Chen, X., X. Zhang, J. A. Church, C. S. Watson, M. A. King, D. Monselesan, B. Legresy, andC. Harig (2017), The increasing rate of global mean sea-level rise during 1993 – 2014, NatureClimate Change, 7(June), 492–495, doi:10.1038/NCLIMATE3325.Christoffersen, P., R. Pettersson, A. Hubbard, S. H. Doyle, and S. Grigsby (2018), Cascadinglake drainage on the Greenland Ice Sheet triggered by tensile shock and fracture, NatureComm, 9(1064), 1–12, doi:10.1038/s41467-018-03420-8.Chu, W., T. T. Creyts, and R. E. Bell (2016), Rerouting of subglacial water flor between neigh-boring glaciers in West Greenland, Journal of Geophysical Research: Earth Surface, 1211, 925–938, doi:10.1002/2015JF003705.Received.Church, J., P. Clark, A. Cazenave, J. Gregory, S. Jevrejeva, A. Levermann, M. Merrifield,G. Milna, R. Nerem, P. Nunn, A. Payne, W. Pfeffer, D. Stammer, and A. Unnikrishnan (2013),Sea level change, in Climate change 2013: The physical science basis. Contribution of workinggroup I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change., editedby T. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. Allen, J. Boschung, A. Nauels, Y. Xia, andP. Midgley, Cambridge University Press, Cambridge, UK and New York, NY, USA.Chylek, P., M. K. Dubey, and G. Lesins (2006), Greenland warming of 1920-1930 and 1995-2005,Geophysical Research Letters, 33(11), 1–5, doi:10.1029/2006GL026510.124BIBLIOGRAPHYClason, C., D. W. F. Mair, D. O. Burgess, and P. W. Nienow (2012), Modelling the delivery ofsupraglacial meltwater to the ice/bed interface: Application to southwest Devon Ice Cap,Nunavut, Canada, Journal of Glaciology, 58(208), 361–374, doi:10.3189/2012JoG11J129.Clason, C. C., D. W. F. Mair, P. W. Nienow, I. D. Bartholomew, A. Sole, S. Palmer, andW. Schwanghart (2015), Modelling the transfer of supraglacial meltwater to the bed of Lev-erett Glacier, Southwest Greenland, Cryosphere, 9(1), 123–138, doi:10.5194/tc-9-123-2015.Colgan, W., and K. Steffen (2009), Modelling the spatial distribution of moulins near Jakob-shavn, Greenland, in IOP Conference Series: Earth and Environmental Science.Conrad, O., B. Bechtel, M. Bock, H. Dietrich, E. Fischer, L. Gerlitz, J. Wehberg, V. Wichmann,and J. Bohner (2015), System for Automated Geoscientific Analyses (SAGA).Cowton, T., P. Nienow, A. Sole, J. Wadham, G. Lis, I. Bartholomew, D. Mair, and D. Chandler(2013), Evolution of drainage system morphology at a land - terminating Greenlandic outletglacier, Journal of Geophysical Research: Earth Surface, 118, 29–41, doi:10.1029/2012JF002540.Cuffey, K. M., and W. Paterson (2010), The physics of glaciers, 4th editio ed., Butterworth-Heinemann/Elsevier, Burlington, MA.Das, S., I. Joughin, M. Behn, I. Howat, M. King, D. Lizarralde, and M. Bhatia (2008), Fracturepropagation to the base of the Greenland ice sheet during supraglacial lake drainage, Science,320(May), 778–782.Dewart, G. (1965), Moulins on Kaskawulsh Glacier, Yukon, Journal of Glaciology, pp. 2–3.DGGS Staff (2013), Elevation Datasets of Alaska: Alaska Division of Geological and Geophys-ical Surveys Digital Data Series 4, doi:10.14509/25239.Dingman, S. L. (2015), Physical Hydrology, 3rd ed., 643 pp., Waveland Press, Inc., Long Grove,IL.Doctor, D. H., and J. a. Young (2013), An evaluation of Automated Gis Tools for Delineat-ing Karst Sinkholes and Closed Depressions From 1-Meter Lidar-Derived Digital ElevationData, Thirteenth Multidisciplinary Conference on Sinkholes and the Engineering and Environmen-tal Impacts of Karst, pp. 449–458.Dow, C. F., W. S. Lee, J. S. Greenbaum, C. A. Greene, D. D. Blankenship, K. Poinar, A. L.Forrest, D. A. Young, and C. J. Zappa (2018), Basal channels drive active surface hydrologyand transverse ice shelf fracture, Science Advances, 4(June), 1–9.Dozier, J. (1970), Part II: Channel adjustments in supraglacial stream, in Studies of morphologyand stream action on ablating ice, pp. 69–117.Dozier, J. (1976), An examination of the variance minimization tendencies of a supraglacialstream, Journal of Hydrology, 31, 359–380.Driesschaert, E., T. Fichefet, H. Goosse, P. Huybrechts, I. Janssens, A. Mouchet, G. Munhoven,V. Brovkin, and S. L. Weber (2007), Modeling the influence of Greenland ice sheet meltingon the Atlantic meridional overturning circulation during the next millennia, GeophysicalResearch Letters, 34, 1–5, doi:10.1029/2007GL029516.125BIBLIOGRAPHYESRI (2016), ArcGIS Desktop: Release 10.3.1.Ewing, K. (1970), Part II: Supraglacial streams on the Kaskawulsh glacier, Yukon Territory, inStudies of morphology and stream action on ablating ice, pp. 131–167.Fairbanks, R. G. (1989), A 17,000-year glacio-eustatic sea level record: influence of glacial melt-ing rates on the Younger Dryas event and deep-ocean circulation, Nature, 342, 637–642.Ferguson, R. (1973), Sinuosity of supraglacial streams, Geological Society of America Bulletin,(January), 251–255.Fettweis, X., J. P. van Ypersele, H. Gallée, F. Lefebre, and W. Lefebvre (2007), The 1979-2005Greenland ice sheet melt extent from passive microwave data using an improved versionof the melt retrieval XPGR algorithm, Geophysical Research Letters, 34(5), 1–5, doi:10.1029/2006GL028787.Fitzpatrick, A. A. W., A. L. Hubbard, J. E. Box, D. J. Quincey, D. Van As, A. P. B. Mikkelsen, S. H.Doyle, C. F. Dow, B. Hasholt, and G. A. Jones (2014), A decade (2002-2012) of supraglaciallake volume estimates across Russell Glacier, West Greenland, Cryosphere, 8(1), 107–121,doi:10.5194/tc-8-107-2014.Flowers, G. E. (2018), Hydrology and the future of the Greenland Ice Sheet, Nature Communi-cations, 9(2729), 1–4, doi:10.1038/s41467-018-05002-0.Flowers, G. E., and G. K. C. Clarke (1999), Surface and bed topography of Trapridge Glacier,Yukon Territory, Canada: digital elevation models and derived hydraulic geometry, Journalof Glaciology, 149.Fountain, A., and J. Walder (1998), Water flow through temperate glaciers, Reviews of Geo-physics, 36(3), 299–328.Fountain, A. G. (1996), Effect of Snow and Firn Hydrology on the Physical and ChemicalCharacteristics of Glacial Runoff, Hydrological Processes, 10(August 1994), 509–521, doi:10.1002/(SICI)1099-1085(199604)10:4<509::AID-HYP389>3.0.CO;2-3.Frauenfeld, O. W., P. C. Knappenberger, and P. J. Michaels (2011), A reconstruction of annualGreenland ice melt extent , 1784 – 2009, Journal of Geophysical Research, 116, 1–7, doi:10.1029/2010JD014918.Fyffe, C. L., B. W. Brock, M. P. Kirkbride, D. W. F. Mair, N. S. Arnold, C. Smiraglia, G. Dio-lauiti, and F. Diotri (2015), An investigation of the influence of supraglacial debris on glacier-hydrology, The Cryosphere Discussions, 9, 5373–5411, doi:10.5194/tcd-9-5373-2015.Gärtner-roer, I., K. Naegeli, M. Huss, T. Knecht, H. Machguth, and M. Zemp (2014), A databaseof worldwide glacier thickness observations, Global and Planetary Change, 122, 330–344, doi:10.1016/j.gloplacha.2014.09.003.Gleason, C. J., L. C. Smith, V. W. Chu, C. J. Legleiter, L. H. Pitcher, B. T. Overstreet, A. K. Ren-nermalm, R. R. Forster, and K. Yang (2016), Characterizing supraglacial meltwater channelhydraulics on the Greenland Ice Sheet from in situ observations, Earth Surface Processes andLandforms, 41(14), 2111–2122, doi:10.1002/esp.3977.Goodell, J. (2016), Greenland Melting, Rolling Stone.126BIBLIOGRAPHYGRASS Development Team (2015), Geographic Resources Analysis Support System (GRASS)Software.Gulley, J., D. Benn, E. Screaton, and J. Martin (2009a), Mechanisms of englacial conduit forma-tion and their implications for subglacial recharge, Quaternary Science reviews, 28, 1984–1999,doi:10.1017/CBO9781107415324.004.Gulley, J. D., D. I. Benn, D. Müller, and A. Luckman (2009b), A cut-and-closure origin forenglacial conduits in uncrevassed regions of polythermal glaciers, Journal of Glaciology,55(189), 66–80, doi:10.3189/002214309788608930.Hall, D. K., J. C. Comiso, N. E. Digirolamo, C. a. Shuman, J. E. Box, and L. S. Koenig (2013),Variability in the surface temperature and melt extent of the Greenland ice sheet fromMODIS, Geophysical Research Letters, 40(10), 2114–2120, doi:10.1002/grl.50240.Hanna, E., P. Huybrechts, K. Steffen, J. Cappelen, R. Huff, C. Shuman, T. Irvine-Fynn, S. Wise,and M. Griffiths (2008), Increased runoff from melt from the Greeland Ice Sheet: a responseto global warming, Journal of Climate, 21, 331–341, doi:10.1175/2007JCLI1964.1.Hodgkins, R. (1997), Glacier hydrology in Svalbard, Norwegian High Arctic, Quaternary Sci-ence Reviews, 16(9), 957–973, doi:10.1016/S0277-3791(97)00032-2.Hoffman, M. J., G. A. Catania, T. A. Neumann, L. C. Andrews, and J. A. Rumrill (2011), Linksbetween acceleration, melting, and supraglacial lake drainage of the western Greenland IceSheet, Journal of Geophysical Research: Earth Surface, 116(4), 1–16, doi:10.1029/2010JF001934.Hoffman, M. J., M. Perego, L. C. Andrews, S. F. Price, T. A. Neumann, J. V. Johnson, G. Catania,and M. P. Lüthi (2018), Widespread moulin formation during supraglacial lake drainages inGreenland, Geophysical Research Letters, 45.Holmlund, P. (1988), Internal geometry and evolution of moulins, Storglaciaren, Sweden, Jour-nal of Glaciology, 34(117), 242–248.Holmlund, P., and R. L. Hooke (1983), High Water-Pressure Events in Moulins, Storglaciären,Sweden, Geografiska Annaler. Series A. Physical Geography, 1(2), 19–25.Horton, R. (1945), Erosional development of streams and their drainage basins; hydrophysicalapproach to quantitative morphology, Geological Society of America Bulletin, 56, 275–370.Iken, A. (1972), Measurements of water pressure in moulins as part of a movement study ofthe White Glacier, Axel Heiberg Island, Northwest Territories, Canada, Journal of Glaciology,61, 53–58.Iken, A. (1981), The effect of the subglacial water pressure on the sliding velocity of a glacierin an idealized numerical model., Journal of Glaciology, 27(97), 407–421.Irvine-Fynn, T. D. L., A. J. Hodson, B. J. Moorman, G. Vatne, and A. L. Hubbard (2011),Polythermal glacier hydrology: A review, Reviews of Geophysics, 49(4), 1–37, doi:10.1029/2010RG000350.Jóhannesson, T., H. Björnsson, F. Pálsson, O. Sigurðsson, and o. Þorsteinsson (2003), LiDARmapping of the Snæfellsjökull ice cap , western Iceland, Jökull, pp. 19–32.127BIBLIOGRAPHYJoughin, I. (2002), Ice-sheet velocity mapping: A combined interferometric and speckle-tracking approach, Annals of Glaciology, 34, 195–201, doi:10.3189/172756402781817978.Joughin, I., W. Abdalati, and M. Fahnestock (2004), Large fluctuations in speed on Greenland’sJakobshavn Isbræ glacier, Nature, 432(7017), 608–610, doi:10.1038/nature03130.Joughin, I., S. B. Das, M. A. King, B. E. Smith, I. M. Howat, and T. Moon (2008), SeasonalSpeedup Along the Western Flank of the Greenland Ice Sheet, Science, 320(5877), 781–783.Joughin, I., B. Smith, M. Howat, T. Scambos, and T. Moon (2010), Greenland Flow Variabilityfrom Ice-Sheet-Wide Velocity Mapping, Journal of Glaciology, 56(197), 415–430, doi:doi:10.3189/002214310792447734.Joughin, I., S. B. Das, G. E. Flowers, M. D. Behn, R. B. Alley, M. A. King, B. E. Smith, J. L.Bamber, M. R. Van Den Broeke, and J. H. Van Angelen (2013), Influence of ice-sheet geometryand supraglacial lakes on seasonal ice-flow variability, Cryosphere, 7(4), 1185–1192, doi:10.5194/tc-7-1185-2013.Kamb, B. (1987), Glacier surge mechanism based on linked cavity configuration of the basalwater conduit system, Journal of Geophysical Research, 92(B9), 9083–9100.Kamintzis, J. (2015), The spatial dynamics of an annual supraglacial meltwater channel in theablation zone of Haut Glacier d ’ Arolla , Switzerland, Msc. thesis, Aberyswyth University,doi:10.13140/RG.2.1.2142.2484.Karlstrom, L., and K. Yang (2016), Fluvial supraglacial landscape evolution on the GreenlandIce Sheet, Geophysical Research Letters, 43, 2683–2692.Karlstrom, L., P. Gajjar, and M. Manga (2013), Meander formation in supraglacial streams,Journal of Geophysical Research: Earth Surface, 118(3), 1897–1907, doi:10.1002/jgrf.20135.Karlstrom, L., A. Zok, and M. Manga (2014), Near-surface permeability in a supraglacialdrainage basin on the Llewellyn Glacier, Juneau Icefield, British Columbia, The Cryosphere,8(2), 537–546, doi:10.5194/tc-8-537-2014.King, L., and M. Tadaki (2018), A Framework for Understanding the Politics of Science (CoreTenet #2), in The Palgrave Handbook of Critical Physical Geography, edited by R. Lave, C. Bier-mann, and S. Lane, pp. 67–88, Palgrave Macmillan, Cham.King, L., M. A. Hassan, K. Yang, and G. E. Flowers (2016), Flow Routing for Delin-eating Supraglacial Meltwater Channel Networks, Remote Sensing, 8(988), doi:10.3390/rs8120988.Kingslake, J., and a. Sole (2015), Modelling channelized surface drainage of supraglacial lakes,Journal of Glaciology, 61(225), 185–199, doi:10.3189/2015JoG14J158.Kirchner, N., K. Hutter, M. Jakobsson, and R. Gyllencreutz (2011), Capabilities and limitationsof numerical ice sheet models : a discussion for Earth-scientists and modelers, QuaternaryScience Reviews, 30, 3691–3704, doi:10.1016/j.quascirev.2011.09.012.Knighton, A. (1981), Channel form and flow characteristics of supraglacial streams, AustreOkstindbreen, Norway, Arctic and Alpine Research, 13(3), 295–306.128BIBLIOGRAPHYKolbert, E. (2016), Greenland Is Melting, The New Yorker.Koziol, C., and N. Arnold (2018), Modelling seasonal meltwater forcing of the velocity of land-terminating margins of the Greenland Ice Sheet, The Cryosphere, 12, 971–991, doi:10.5194/tc-12-971-2018.Koziol, C., N. Arnold, A. Pope, and W. Colgan (2017), Quantifying supraglacial meltwaterpathways in the Paakitsoq region, West Greenland, Journal of Glaciology, 63(239), 464–46,doi:10.1017/jog.2017.5.Krabill, W., E. Hanna, P. Huybrechts, W. Abdalati, J. Cappelen, B. Csatho, E. Frederick, S. Man-izade, C. Martin, J. Sonntag, R. Swift, R. Thomas, and J. Yungel (2004), Greenland Ice Sheet :Increased coastal thinning, Geophysical Research Letters, 31, 1–5, doi:10.1029/2004GL021533.Kumar, R., C. Chatterjee, R. D. Singh, A. K. Lohani, and S. Kumar (2007), Runoff estimation foran ungauged catchment using geomorphological instantaneous unit hydrograph ( GIUH )models, Hydrological Processes, 1840(May), 1829–1840, doi:10.1002/hyp.Lampkin, D. J., and J. Vanderberg (2011), A preliminary investigation of the influence of basaland surface topography on supraglacial lake distribution near Jakobshavn Isbrae, westernGreenland, Hydrological Processes, 25(21), 3347–3355, doi:10.1002/hyp.8170.Lampkin, D. J., and J. VanderBerg (2014), Supraglacial melt channel networks in the Jakob-shavn Isbrae region during the 2007 melt season, Hydrological Processes, 28(25), 6038–6053,doi:10.1002/hyp.10085.Lane, S. N. (2017), Slow science, the geographical expedition, and Critical Physical Geography,Canadian Geographer, 61(1), 84–101, doi:10.1111/cag.12329.Langley, E. S., A. A. Leeson, C. R. Stokes, and S. S. R. Jamieson (2016), Seasonal evolution ofsupraglacial lakes on an East Antarctic outlet glacier, Geophysical Research Letters, 43, 8563–8571, doi:10.1002/2016GL069511.Received.Leek, J. T., and R. D. Peng (2015), Opinion: Reproducible research can still be wrong: Adoptinga prevention approach, Proceedings of the National Academy of Sciences, 112(6), 1645–1646,doi:10.1073/pnas.1421412111.Leeson, A., A. Shepherd, K. Briggs, I. Howat, X. Fettweis, M. Morlighem, and E. Rignot (2015),Supraglacial lakes on the Greenland ice sheet advance inland under warming climate, Na-ture Climate Change, 5(January), 51–55, doi:10.1038/NCLIMATE2463.Legleiter, C. J., M. Tedesco, L. C. Smith, a. E. Behar, and B. T. Overstreet (2014), Mappingthe bathymetry of supraglacial lakes and streams on the Greenland ice sheet using fieldmeasurements and high-resolution satellite images, The Cryosphere, 8(1), 215–228, doi:10.5194/tc-8-215-2014.Lewis, S. M., and L. C. Smith (2011), Hydrologic drainage of the Greenland Ice Sheet, Hydro-logical Processes, 23, 2004–2011, doi:10.1002/hyp.Lindback, K., R. Pettersson, A. L. Hubbard, S. H. Doyle, D. Van As, A. B. Mikkelsen, and A. A.Fitzpatrick (2015), Subglacial water drainage, storage, and piracy beneath the Greenland icesheet, Geophysical Research Letters, 42, 7606–7614, doi:10.1002/2015GL065393.129BIBLIOGRAPHYMaidment, D. (2002), Arc hydro: GIS for water resources.Mantelli, E., C. Camporeale, and L. Ridolfi (2015), Supraglacial channel inception: Modelingand processes, Water Resources Research, 51, 7044–7063, doi:10.1002/2014WR016259.Marston, R. A. (1983), Supraglacial Stream Dynamics on the Juneau Icefield, Annals of the As-sociation of American Geographers, 73(4), 597–608.McGrath, D., W. Colgan, K. Steffen, P. Lauffenburger, and J. Balog (2011), Assessing the sum-mer water budget of a moulin basin in the sermeq avannarleq ablation region, Greenlandice sheet, Journal of Glaciology, 57(205), 954–964, doi:10.3189/002214311798043735.Meierbachtol, T., J. Harper, and N. Humphrey (2013), Basal drainage system response toincreasing surface melt on the Greenland ice sheet., Science, 341, 777–779, doi:10.1126/science.1235905.Mikkelsen, A. B., A. Hubbard, M. Macferrin, J. Eric Box, S. H. Doyle, A. Fitzpatrick, B. Hasholt,H. L. Bailey, K. Lindbäck, and R. Pettersson (2016), Extraordinary runoff from the Greenlandice sheet in 2012 amplified by hypsometry and depleted firn retention, Cryosphere, 10(3),1147–1159, doi:10.5194/tc-10-1147-2016.Montgomery, D., and W. Dietrich (1989), Source areas, drainage density, and channel initiation,Water Resources Research, 25(8), 1907–1918.Montgomery, D., and W. Dietrich (1994), Landscape dissection and drainage area-slope thresh-olds, in Process Models and Theoretical Geomorphology, edited by M. Kirkby, p. 417, John Wiley& Sons Ltd.Montgomery, D., and E. Foufoula-Georgiou (1993), Channel network source representationusing digital elevation models, Water Resources Research, 29(12), 3925–3934.Montgomery, D. R., and W. E. Dietrich (1988), Where do channels begin?, Nature, 336, 232–234,doi:10.1038/336232a0.Moon, T., I. Joughin, B. Smith, M. R. Van Den Broeke, W. J. Van De Berg, B. Noël, and M. Usher(2014), Distinct patterns of seasonal Greenland glacier velocity, Geophysical Research Letters,41, 7209–7216, doi:10.1002/2014GL061836.Morlighem, M., C. N. Williams, E. Rignot, L. An, J. E. Arndt, J. L. Bamber, G. Catania,N. Chauché, J. A. Dowdeswell, B. Dorschel, I. Fenty, K. Hogan, I. Howat, A. Hubbard,M. Jakobsson, T. M. Jordan, K. K. Kjeldsen, R. Millan, L. Mayer, J. Mouginot, B. P. Y. Noël,C. Ó. Cofaigh, S. Palmer, S. Rysgaard, H. Seroussi, M. J. Siegert, P. Slabon, F. Straneo,M. R. van den Broeke, W. Weinrebe, M. Wood, and K. B. Zinglersen (2017), BedMachine v3:Complete bed topography and ocean bathymetry mapping of Greenland from multi-beamecho sounding combined with mass conservation, Geophysical Research Letters, pp. 51–61,doi:10.1002/2017GL074954.Morriss, B. F., R. L. Hawley, J. W. Chipman, L. C. Andrews, G. A. Catania, M. J. Hoffman, M. P.Lüthi, and T. A. Neumann (2013), A ten-year record of supraglacial lake evolution and rapiddrainage in West Greenland using an automated processing algorithm for multispectral im-agery, Cryosphere, 7(6), 1869–1877, doi:10.5194/tc-7-1869-2013.130BIBLIOGRAPHYNerem, R. S., B. D. Beckley, J. T. Fasullo, B. D. Hamlington, D. Masters, and G. T. Mitchum(2018), Climate-change-driven accelerated sea-level rise detected in the altimeter era, Pro-ceedings of the National Academy of Sciences, pp. 1–4, doi:10.1073/pnas.1717312115.Nghiem, S. V., D. K. Hall, T. L. Mote, M. Tedesco, M. R. Albert, K. Keegan, C. A. Shuman, N. E.DiGirolamo, and G. Neumann (2012), The extreme melt across the Greenland ice sheet in2012, Geophysical Research Letters, 39(20), 6–11, doi:10.1029/2012GL053611.Nick, F. M., A. Vieli, I. M. Howat, and I. Joughin (2009), Large-scale changes in Greenlandoutlet glacier dynamics triggered at the terminus, Nature Geoscience, 2(2), 110–114, doi:10.1038/ngeo394.Nienow, P., and B. Hubbard (2005), Surface and Englacial Drainage of Glaciers and Ice Sheets,Encyclopedia of Hydrological Sciences, pp. 1–12, doi:10.1002/0470848944.Noël, B., W. J. Van De Berg, E. Van Meijgaard, P. Kuipers Munneke, R. S. Van De Wal, andM. R. Van Den Broeke (2015), Evaluation of the updated regional climate model RACMO2.3:Summer snowfall impact on the Greenland Ice Sheet, Cryosphere, 9(5), 1831–1844, doi:10.5194/tc-9-1831-2015.Noh, M.-J., and I. M. Howat (2015), Automated stereo-photogrammetric DEM generation athigh latitudes: Surface Extraction with TIN-based Search-space Minimization (SETSM) val-idation and demonstration over glaciated regions, GIScience & Remote Sensing, (August),1–20, doi:10.1080/15481603.2015.1008621.Nosek, B., G. Alter, G. Banks, D. Borsboom, S. Bowman, S. Breckler, S. Buck, C. Chambers,G. Chin, G. Christensen, M. Contestabile, A. Dafoe, E. Eich, J. Freese, R. Glennerster, D. Go-roff, D. Green, B. Hesse, M. Humphreys, J. Ishlyama, D. Karlan, A. Kraut, A. Lupia, P. Mabry,T. Madon, N. Malhotra, E. Mayo-Wilson, M. McNutt, E. Miguel, E. Levy Paluck, U. Si-monsohn, C. Soderberg, B. Spellman, J. Turitto, G. VandenBos, S. Vazire, E. Wagenmakers,R. Wilson, and T. Yarkoni (2015), Promoting an open research culture, Science, 349(6242),1422–1425, doi:10.1126/science.aab3847.O’Callaghan, J. F., and D. M. Mark (1984), The extraction of drainage networks from digitalelevation data, Computer Vision, Graphics, and Image Processing, 28(3), 323–344, doi:10.1016/S0734-189X(84)80011-0.Oerlemans, B. J., and N. C. Hoogendoorn (1989), Mass-balance gradients and climatic change,Journal of Glaciology, 35(121), 399–405.Onesti, L. J. (1987), Slushflow release mechanism: A first approximation, Avalanche formation,movement and effects, IAHS, 162(162), 331–336.O’Reilly, J. (2015), Glacial dramas: Typos, projections and peer review in the fourth assessmentof the Intergovernmental Panel on Climate Change, in Climate Cultures: Anthropological Per-spectives on Climate Change, edited by J. Barnes and M. Dove, University Press ScholarshipOnline, doi:10.12987/yale/9780300198812.001.0001.O’Reilly, J. (2017), The technocratic Antarctic: An ethnography of scientific expertise and environmen-tal governance, Cornell University Press, Ithaca and London.131BIBLIOGRAPHYPalmer, S., A. Shepherd, P. Nienow, and I. Joughin (2011), Seasonal speedup of the GreenlandIce Sheet linked to routing of surface water, Earth and Planetary Science Letters, 302(3-4), 423–428, doi:10.1016/j.epsl.2010.12.037.Parizek, B. R., and R. B. Alley (2004), Implications of increased Greenland surface melt underglobal-warming scenarios : ice-sheet simulations, Quaternary Science Reviews, 23, 1013–1027,doi:10.1016/j.quascirev.2003.12.024.Parker, G. (1975), Meadering of supraglacial melt streams, Water resources research, 11(4), 551–552.Paul, F., N. E. Barrand, S. Baumann, E. Berthier, T. Bolch, K. Casey, G. Nosenko, H. Frey, S. P.Joshi, V. Konovalov, R. L. E. Bris, N. Mo, S. Steffen, and S. Winsvold (2013), On the accuracyof glacier outlines derived from remote-sensing data, Annals of Glaciology, 54(63), 171–182,doi:10.3189/2013AoG63A296.Peng, R. D. (2011), Reproducible Research in Computational Science, Science, 334, 1226–1228,doi:10.1126/science.1213443.Phillips, E., A. Finlayson, and L. Jones (2013), Fracturing, block faulting, and moulin devel-opment associated with progressive collapse and retreat of a maritime glacier: Falljokull,SE Iceland, Journal of Geophysical Research: Earth Surface, 118, 1545–1561, doi:10.1002/jgrf.20116.Phillips, T., S. Leyk, H. Rajaram, W. Colgan, W. Abdalati, D. McGrath, and K. Steffen (2011),Modeling moulin distribution on Sermeq Avannarleq glacier using ASTER and WorldViewimagery and fuzzy set theory, Remote Sensing of Environment, 115(9), 2292–2301, doi:10.1016/j.rse.2011.04.029.Pitcher, L. H., L. C. Smith, C. J. Gleason, and K. Yang (2017), CryoSheds : a GIS modelingframework for delineating land-ice watersheds for the Greenland Ice Sheet, GIScience &Remote Sensing, 53(6), 707–722, doi:10.1080/15481603.2016.1230084.Poinar, K., I. Joughin, S. B. Das, M. D. Behn, J. T. M. Lenaerts, and M. R. Van Den Broeke (2015),Limits to future expansion of surface-melt-enhanced ice flow into the interior of westernGreenland, Geophysical Research Letters, 42(6), 1800–1807, doi:10.1002/2015GL063192.Pollard, D., R. M. Deconto, and R. B. Alley (2015), Potential Antarctic Ice Sheet retreat drivenby hydrofracturing and ice cliff failure, Earth and Planetary Science Letters, 412, 112–121, doi:10.1016/j.epsl.2014.12.035.Price, S. F., A. J. Payne, I. M. Howat, and B. E. Smith (2011), Committed sea-level rise for thenext century from Greenland ice sheet dynamics during the past decade, Proceedings of theNational Academy of Sciences, 108(22), doi:10.1073/pnas.1017313108.Pritchard, H. D., R. J. Arthern, D. G. Vaughan, and L. A. Edwards (2009), Extensive dynamicthinning on the margins of the Greenland and Antarctic ice sheets, Nature, 461, 971–975,doi:10.1038/nature08471.Rasmussen, S. O. Ã., B. M. Vinther, H. B. Clausen, and K. K. Andersen (2007), Early Holoceneclimate oscillations recorded in three Greenland ice cores, Quaternary Science Reviews, 26,1907–1914, doi:10.1016/j.quascirev.2007.06.015.132BIBLIOGRAPHYRignot, E., and P. Kanagaratnam (2006), Changes in the velocity structure of the Greenland IceSheet, Science, 311(5763), 986–990.Rippin, D. M., A. Pomfret, and N. King (2015), High resolution mapping of supra-glacialdrainage pathways reveals link between micro-channel drainage density, surface roughnessand surface reflectance, Earth Surface Processes and Landforms, doi:10.1002/esp.3719.Rodriguez-Iturbe, I., and J. B. Valdes (1979), The geomorphologic structure of hydrologic re-sponse, Water Resources Research, 15(6), 1409–1420.Rodriguiez-Iturbe, I., and A. Rinaldo (1997), Fractal River Basins: Chance and Self-Organization,Cambridge University Press, Cambridge.Ross, N., A. J. Sole, S. J. Livingstone, Á. Igneczi, and M. Morlighem (2018), Near-margin icethickness and subglacial water routing, Leverett Glacier, Greenland, Arctic, Antarctic, andAlpine Research, 50(1), doi:10.1080/15230430.2017.1420949.Rothlisberger, H. (1972), Water pressure in intra- and subglacial channels, Journal of Glaciology,11(62), 177–203.Schoof, C. (2010), Ice-sheet acceleration driven by melt supply variability., Nature, 468(7325),803–6, doi:10.1038/nature09618.Schroeder, J. (1998), Hans glacier moulins observed from 1988 to 1992, Svalbard, NorwegianJournal of Geography, 52(2), 79–88, doi:10.1080/00291959808552387.Schwanghart, W., and N. Kuhn (2010), TopoToolbox: a set of Matlab functions for topographicanalysis., Environmental Modelling & Software, 25, 770–781, doi:OI:10.1016/j.envsoft.2009.12.002.Shean, D. E., O. Alexandrov, Z. M. Moratto, B. E. Smith, I. R. Joughin, C. Porter, and P. Morin(2016), An automated, open-source pipeline for mass production of digital elevation mod-els (DEMs) from very-high-resolution commercial stereo satellite imagery, ISPRS Journal ofPhotogrammetry and Remote Sensing, 116, 101–117, doi:10.1016/j.isprsjprs.2016.03.012.Sherman, L. K. (1932), The relation of hydrographs of runoff to size and character of drainage-basins, Transactions, American Geophysical Union, 13(1), 332–339.Shreve, L. (1972), Movement of water in glaciers, Journal of Glaciology, 11(62), 205–214.Singh, P. K., S. K. Mishra, and M. K. Jain (2014), A review of the synthetic unit hydrograph:from the empirical UH to advanced geomorphological methods, Hydrological Sciences Jour-nal, 59(2), 239–261, doi:10.1080/02626667.2013.870664.Smith, L. C., V. W. Chu, K. Yang, C. J. Gleason, L. H. Pitcher, A. K. Rennermalm, C. J. Legleiter,A. E. Behar, B. T. Overstreet, S. E. Moustafa, M. Tedesco, R. R. Forster, A. L. LeWinter, D. C.Finnegan, Y. Sheng, and J. Balog (2015), Efficient meltwater drainage through supraglacialstreams and rivers on the southwest Greenland ice sheet., Proceedings of the National Academyof Sciences of the United States of America, 112(4), 1001–6, doi:10.1073/pnas.1413024112.Smith, L. C., K. Yang, L. H. Pitcher, B. T. Overstreet, V. W. Chu, K. Rennermalm, J. Ryan, M. G.Cooper, C. J. Gleason, M. Tedesco, J. Jeyaratnam, D. van As, M. van den Broeke, W. J. van deBerg, B. Noel, P. L. Langen, R. I. Cullather, B. Zhao, M. J. Willis, A. Hubbard, J. E. Box, B. A.133BIBLIOGRAPHYJenner, and A. E. Behar (2017), Direct measurements of meltwater runoff and retention onthe Greenland ice sheet surface, Proceedings of the National Academy of Sciences, 114(50), 1–40,doi:10.1073/pnas.1707743114.Snyder, F. F. (1938), Synthetic Unit-Graphs, Transactions of the Institute of British Geographers,19(1), 447–454.Sole, A., P. Nienow, I. Bartholomew, D. Mair, T. Cowton, A. Tedstone, and M. A. King (2013),Winter motion mediates dynamic response of the Greenland Ice Sheet to warmer summers,Geophysical Research Letters, 40, 3940–3944, doi:10.1002/grl.50764.Stenborg, T. (1968), Glacier drainage connected with ice structures, Geografiska Annaer. SeriesA, Physical Geography, 50(1), 25–53.Stevens, L. A., M. D. Behn, J. J. McGuire, S. B. Das, I. Joughin, T. Herring, D. E. Shean, andM. A. King (2015), Greenland supraglacial lake drainages triggered by hydrologically in-duced basal slip, Nature, 522(7554), 73–76, doi:10.1038/nature14480.Stevens, L. A., M. D. Behn, S. B. Das, I. R. Joughin, B. Noel, M. R. van den Broeke, and T. Her-ring (2016), Greenland Ice Sheet fl ow response to runoff variability, Geophysical ResearchLetters, 43(11), 1–9, doi:10.1002/2016GL070414.Received.Sundal, A. V., A. Shepherd, P. Nienow, E. Hanna, S. Palmer, and P. Huybrechts (2009), Evo-lution of supra-glacial lakes across the Greenland Ice Sheet, Remote Sensing of Environment,113(10), 2164–2171, doi:10.1016/j.rse.2009.05.018.Sundal, A. V., A. Shepherd, P. Nienow, E. Hanna, S. Palmer, and P. Huybrechts (2011), Melt-induced speed-up of Greenland ice sheet offset by efficient subglacial drainage., Nature,469(7331), 521–524, doi:10.1038/nature09740.Tarboton, D. (1997), A new method for the determination of flow directions and upslope areasin grid digital elevation models, Water Resources Research, 33(2), 309–319.Tarboton, D. G., R. L. Bras, and I. Rodriguiez-Iturbe (1991), On the extraction of channel net-works from digital elevation data, Hydrological processes, 5, 81–100.Tedesco, M., J. Box, J. Cappelen, R. Fausto, X. Fettweis, T. Mote, C. Smeets, D. van As,I. Velicogna, R. S. W. van de Wal, and J. Wahr (2016), Greenland Ice Sheet, in NOAA Arc-tic Report Card 2016.Tedstone, A. J., P. W. Nienow, N. Gourmelen, A. Dehecq, D. Goldberg, and E. Hanna (2015),Decadal slowdown of a land-terminating sector of the Greenland Ice Sheet despite warming,Nature, 526(7575), 692–695, doi:10.1038/nature15722.Thomas, E. R., E. W. Wolff, R. Mulvaney, J. P. Steffensen, S. Johnsen, C. Arrowsmith, J. W. C.White, B. Vaughn, and T. Popp (2007), The 8.2 ka event from Greenland ice cores, QuaternaryScience Reviews, 26, 70–81, doi:10.1016/j.quascirev.2006.07.017.Thomsen, H. H. (1986), Photogrammetric and satellite mapping of the margin of the inlandice, West Greenland, Annals Of Glaciology, 8, 164–167.Tucker, G., and R. Bras (1998), Hillslope processes, drainage density, and landscape morphol-ogy, Water Resources Research, 34(10), 2751–2764.134BIBLIOGRAPHYTucker, G. E., F. Catani, A. Rinaldo, and R. L. Bras (2001), Statistical analysis of drainage densityfrom digital terrain data, Geomorphology, 36(3-4), 187–202, doi:10.1016/S0169-555X(00)00056-8.USGS (2012), High Resolution Lidar Digital Elevation Models and Low Resolution ShadedRelief Maps of Antarctica from USGS.van de Wal, R., W. Boot, M. van den Broeke, C. Smeets, C. Reijmer, J. Donker, and J. Oerlemans(2008), Large and rapid melt-induced velocity changes in the ablation zone of the Greenlandice sheet, Science, 321(July), 111–114.van den Broeke, M. R., E. M. Enderlin, I. M. Howat, P. Kuipers Munneke, B. P. Noël,W. Jan Van De Berg, E. Van Meijgaard, and B. Wouters (2016), On the recent contribu-tion of the Greenland ice sheet to sea level change, Cryosphere, 10(5), 1933–1946, doi:10.5194/tc-10-1933-2016.van der Veen, C. J. (2007), Fracture propagation as means of rapidly transferring surfacemeltwater to the base of glaciers, Geophysical Research Letters, 34(1), L01,501, doi:10.1029/2006GL028385.Williamson, A. G., I. A. N. C. Willis, N. S. Arnold, and A. F. Banwell (2018), Controls on rapidsupraglacial lake drainage in West Greenland : an Exploratory Data Analysis approach,Journal of Glaciology, doi:10.1017/jog.2018.8.Wyatt, F. R., and M. J. Sharp (2015), Linking surface hydrology to flow regimes and patternsof velocity variability on Devon Ice Cap, Nunavut, Journal of Glaciology, 61(226), 387–399,doi:10.3189/2015JoG14J109.Yang, K., and L. C. Smith (2013), Supraglacial Streams on the Greenland Ice Sheet DelineatedFrom Combined Spectral – Shape Information in High-Resolution Satellite Imagery, IEEEGeoscience and remote sensing letters, 10(4), 801–805.Yang, K., and L. C. Smith (2016), Internally drained catchments dominate supraglacial hy-drology of the southwest Greenland Ice Sheet, Journal of Geophysical Research: Earth Surface,121(10), doi:10.1002/2016JF003927.Yang, K., L. C. Smith, V. W. Chu, C. J. Gleason, and M. Li (2015), A Caution on the Use ofSurface Digital Elevation Models to Simulate Supraglacial Hydrology of the Greenland IceSheet, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 8(11),5212–5224, doi:10.1109/JSTARS.2015.2483483.Yang, K., L. C. Smith, V. W. Chu, L. H. Pitcher, C. J. Gleason, A. K. Rennermalm, and M. Li(2016a), Fluvial morphometry of supraglacial river networks on the southwest GreenlandIce Sheet, GIScience & Remote Sensing, 53(4), 459–482, doi:10.1080/15481603.2016.1162345.Yang, K., L. Karlstrom, L. C. Smith, and M. Li (2016b), Automated High-Resolution SatelliteImage Registration Using Supraglacial Rivers on the Greenland Ice Sheet, pp. 1–12.Yang, K., L. C. Smith, L. Karlstrom, M. G. Cooper, M. Tedesco, D. V. As, X. Cheng, Z. Chen,and M. Li (2018), Supraglacial meltwater routing through internally drained catchments onthe Greenland Ice Sheet surface, The Cryosphere Discussions, in review, 1–32.135BIBLIOGRAPHYZhang, W., and D. Montgomery (1994), Digital elevation model grid size, landscape represen-tation, and hydrologic simulations, Water resources research, 30(4), 1019–1028.Zwally, H. J., W. Abdalati, T. Herring, K. Larson, J. Saba, and K. Steffen (2002), Surface Melt –Induced Acceleration of Greenland Ice-Sheet Flow, Science, 297, 218–222.136Appendix AStudy area characteristics  'LVWDQFHP$ Ssurface  'LVWDQFHP%  HicePFigure A.1: A) surface slope (Noh and Howat, 2015) and B) ice thickness (Morlighem et al.,2017) in the study area used in Chapters 4 and 5.137  'LVWDQFHP$  viceP\U  'LVWDQFHP% \U 1)Figure A.2: A) MEaSUREs ice surface velocity (Joughin et al., 2008) and B) derived princi-pal strain rates for the study area used in Chapters 4 and 5.please get rid of this extra page.138Appendix BMoulin Identification1 KmB1 Km1 Km1 KmACDFigure B.1: The above four plotsshow examples of moulins that areevident in the curvature layer of theArcticDEM. Moulins are denoted witha blue dot. Channels are generallyreadily distinguishable from crevassesand other linear features.139Appendix CStrain rate scale considerationsC.1 SummaryError in velocity measurements over the study area is low, with average x- and y-velocitieserrors of 2.5 m/yr and 6.3 m/yr, respectively. Nonetheless, cell to cell differences can be ex-pected and carry through into strain rate in significant ways. Thus, it is pragmatically usefulto smooth the velocity layers to reduce noise in the data. Further, cell specific strain rates maynot be as physically meangingful as strain rates averaged over areas that better capture localice dynamics. Smoothing of velocity data is therefore necessary for this analysis, howeverchoosing an appropriate scale over which to smooth the data is problematic. In this appendix,I include visual and RMSE comparisons of the three most common filtering algorithms overmultiple scales of smoothing.C.2 Mean FilterThe simplest filter is a simple sliding window of a given dimension, which replaces the centervalue in the window with the mean value of all the pixels in the window. This is repeated forevery pixel in the array. The only adjustable parameter in this filter is the size of the kernel.Outputs for kernel sizes ranging from 200 m (unfiltered) to 3000 m are shown in Figure C.1.140    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1C.3. Median FilterFigure C.1: Mean filtering across various kernel sizes.C.3 Median FilterSimilar to the Mean Filter, this filter is a simple sliding window of a given dimension, whichreplaces the center value in the window with the median value of all the pixels in the window.This is repeated for every pixel in the array. The only adjustable parameter in this filter isthe size of the kernel. Outputs for kernel sizes ranging from 200 m (unfiltered) to 3000 m areshown in Figure C.2.143    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1C.4. Gaussian FilterFigure C.2: Median filtering across various kernel sizes.C.4 Gaussian FilterWhereas the mean filter averages all of the values within the kernel, the kernel in a gaussianfilter applies weights to pixel values within the window based on a gaussian distribution,whereby pixels closer to the middle have higher weights than pixels towards the edges of thewindow. This filter has two adjustable parameters: the kernel size, and the standard deviationσ. Although a gaussian distribution has, theoretically, an infinite width, the values on eitherside of the distribution become essentially zero some distance away from the mean. The dis-tance at which this happens depends on the standard deviation of the distribution. For smallerσ, this distance will be closer to the mean, given the tighter distribution of values. The oppo-site is true for larger values of σ. Thus, the appropriate kernel size varies with σ. The tails ofgaussian functions become effectively zero beyond three standard deviations from the mean,meaning that three times the standard deviation is an appropriate rule of thumb for specifyingthe kernel size. Outputs for kernel sizes ranging from 200 m (unfiltered) to 3000 m are shownin Figure C.3.146    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1    :LQGRZVL]H P    :LQGRZVL]H P6WUDLQ5DWH\U1C.5. RMSE comparisonsFigure C.3: Gaussian Blurring filtering across various kernel sizes.C.5 RMSE comparisonsBy visually comparing the outputs of the different filtration mechanisms across various filtersizes (Figures C.1,C.2, and C.3) , it is clear that gaussian filtering does the best job at visuallypreserving high values and the spatial characteristics of the strain rate field, particularly withregards to edges and locally elevated strain rates. Further insight can be gained by comparingthe root mean square error (RMSE) between filtered layers and unfiltered layers (Figure C.4).It is clear that mean filtering produces the highest RMSE across all kernel sizes, and medianfilter the lowest.     .HUQHO6L]HP506(\U1 0HDQ)LOWHUHice2 *Hice(1/2) *Hice(1/2) *Hice     .HUQHO6L]HP506(\U1 0HGLDQ)LOWHUHice2 *Hice(1/2) *Hice(1/2) *Hice     .HUQHO6L]HP506(\U1 *DXVVLDQ%OXUHice2 *Hice(1/2) *Hice(1/2) *HiceFigure C.4: RMSE of different filtering methods across kernel sizes.149Appendix DHydraulic variables and catchmentproperties for derivation of the GIUH.For visual clarity, the figures for this appendix are shown on the following page.150&DWFKPHQW$UHDNP29HORFLW\PV$&DWFKPHQW$UHDNP2WSKU%Sc9HORFLW\PV&&DWFKPHQW$UHDNP2d'Figure D.1: The correlation between hydraulic variables and catchment properties usedin the derivation of the GIUH in Chapter 5, including A) stream velocity, B) time topeak, C) stream velocity with respect to channel slope, and D) channel depth withrespect to catchment area. In panel B, the grey line represents the catchment areaand tp reported by McGrath et al. (2011), and the black line represents the same asreported by Smith et al. (2017) for Rio Behar.151Appendix EAn illustrated thesis timeline.Figure E.1: A graphical summary of the PhD journey. In A) the author celebrates thebeginning of the PhD process in the fall of 2012. In B) a recent selfie shows theauthor celebrating submission of her thesis in the summer of 2018.152


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items