Pattern–forming instabilities in the coupling of ice sheets andsubglacial drainage systemsbyArran WhitefordBSc. Hons. Mathematics, Meteorology Victoria University of Wellington, 2013A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THEREQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Geophysics)The University of British Columbia(Vancouver)September 2018c© Arran Whiteford, 2018The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled:Pattern–forming instabilities in the coupling of ice sheets and subglacial drainage systemssubmitted by Arran Whiteford in partial fulfillment of the requirements for the degree of Masterof Science in Geophysics.Examining Committee:Christian SchoofSupervisorGwenn FlowersSupervisory Committee MemberSusan AllenAdditional ExaminerAdditional Supervisory Committee Members:Brian WettonSupervisory Committee MemberiiAbstractSharp spatial changes discovered in the basal conditions of an ice sheet do not always have anobvious source. By modelling instabilities in the coupling of an ice sheet and subglacial drainagesystem, we describe physical feedback mechanisms that force the formation of sharp spatialstructures in basal conditions and ice flow. This model predicts the spontaneous formation ofperiodic subglacial ‘sticky spot’–lake pairs, that correspond in shape to previous empirical andmodelled descriptions of similar structures. The instability that forms this structure is driven bya feedback whereby periodic humps in ice thickness redirect subglacial water to slippery spotsthat lie immediately downstream of the ice humps: the slippery regions increase ice flux into theice humps, making them grow. Scaling a one–dimensional model ice sheet coupled to a basaldrainage system, we find conditions for the instability with linear stability analysis. Solutions inthe full nonlinear model are simulated numerically, using operator splitting and finite differencemethods. The instability requires a bed permeability weakly dependent on water pressurechanges, negligible bed slopes, and a water velocity much greater than ice velocity. The ‘stickyspot’–lake pairs are predicted to form with periodic spacing and migrate upstream.iiiLay summaryUnderneath the Antarctic ice sheet—between ice and the ground—water can be found. Thiswater affects the flow of ice by lubricating the base of the ice, allowing it to slide faster. Fasterflowing ice will become thinner, whereas if ice slows down it will pile up and thicken. Thesechanges shift the flow of water: squashed water flows away from the heavy weight of thick ice towhere there is less pressure beneath thinner ice. So there is a feedback loop, where the ice speedchanges the ice thickness which redistributes the water underneath which affects the speed ofice. A simple one dimensional model ice sheet describing this feedback revealed the spontaneousformation of patterned structures—made of repeating regions of sticky and slippery bed. Thesticky regions show a dry, high friction bed and correspond to thicker ice. The slippery areasshow water pooling, little bed friction and thinner ice.ivPrefaceThis dissertation represents original work of the author, Arran Whiteford, and their supervisor,Christian Schoof.vContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vContents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xThesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Feedbacks between ice sheets and subglacial drainage systems . . . . . . . . 11.2 Description of instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Subglacial drainage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Subglacial lakes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.5 Patterns in paleo–ice sheet beds . . . . . . . . . . . . . . . . . . . . . . . . 81.6 Thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 1D coupled ice sheet–drainage model . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Introducing the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Scaling the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 Instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.1 Linear stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.2 Reduced model with unstable steady state . . . . . . . . . . . . . . . . . . 18vi3.3 A stability boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.4 Discussion of the mechanism for instability . . . . . . . . . . . . . . . . . . 254 Numerical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.1 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.2 Discussion of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.3 Limitations and research outlook . . . . . . . . . . . . . . . . . . . . . . . . 385.4 Comparison of results with existing evidence . . . . . . . . . . . . . . . . . 406 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48viiList of Figures1 Ice flows from the centre of an ice sheet to the sea. . . . . . . . . . . . . . . . . . 12 Schematic depiction of ice–stream flow over ‘sticky spot’–lake– regular bed. . . . . 33 Basal shear stress derived for MacIS . . . . . . . . . . . . . . . . . . . . . . . . . 44 Schematic of the model domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Grow rates from full and reduced dispersion relations . . . . . . . . . . . . . . . . 236 Growth rates for a range of κN . . . . . . . . . . . . . . . . . . . . . . . . . . . . 247 Properties of the linear solution show different phase shifts. . . . . . . . . . . . . 278 Schematic showing the phase shift between variables that forces a positive feedback. 289 Schematic showing a phase shift between variables forcing a negative feedback. . . 2910 Comparison of linear and non linear solutions . . . . . . . . . . . . . . . . . . . . 3111 Nonlinear evolution of an unstable solution, first parameter regime . . . . . . . . 3212 Developed structure in ice sheet and subglacial drainage . . . . . . . . . . . . . . 3313 Nonlinear evolution of an unstable solution, second parameter regime . . . . . . . 4514 A stability boundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4615 Nonlinear evolution of an unstable solution, third parameter regime . . . . . . . . 4716 Schematic describing deepening lake model extension . . . . . . . . . . . . . . . . 47viiiAcknowledgementsA warm, huge thank you to Christian Schoof, my supervisor, for his guidance throughout mytime at UBC. His great explanations, interesting tangents, hours editing, useful recommendations—I have learnt immense amounts from his wealth of knowledge and understanding of the field.Thank you to my committee: Gwenn Flowers, for Kaskawulsh boogies and great glaciologydiscussion. Brian Wetton, for kind help with the numerics, patience and teaching. Thank youHenry Dengate Thrush for help drawing diagrams. To the glaciology group: Noel, Camilo,Gabriela, Mekdes, Maryam, Sam, Natalita, Jeff. For laughs, research help, discussions, ears ofsupport, and camaraderie.My gratitude extends to the Universities of NZ, Edward and Isabel Kidsen Trust for financialassistance, and flexibility; to Huw Horgan for inspiration and an introduction to glaciology, aswell as help with finding a position and applications.Very thankful for ongoing, limitless support from Mum and Dad, Andrew, Emma and Lily. Awholehearted cheers to the extended UBC: to the Jungle House, the VOC, fondues, mountains,and all the friends in between, constant discussions with friends aboard. To the Salish Sea,for unconditional refreshing of the soul, from the beginning to the end. Joane, a terrific teammate—morning coffees with hummingbirds, cold dips, delicious dinners, and alpine adventuresthroughout the MSc.ixTo the crows in Skwxwu´7meshx1. INTRODUCTION1 IntroductionThe Antarctic ice sheet contains 27 million km3 of ice [Fretwell et al., 2013]. Snow fall augmentsice mass, and is balanced by the discharge of ice, which flows from the interior of the continentto the sea (Fig. 1) [Rignot and Thomas, 2002, Paterson, 2016]. Understanding the dynamics ofice flow is a crucial part of modelling the evolution of the Antarctic ice sheet, which itself is anintegral part of predicting future climate conditions [DeConto and Pollard, 2016, Allen et al.,2014]. In current global sea level models, ice dynamics remains one of the greatest uncertaintiesin predicting sea level rise [Church et al., 2013].The flow of ice is strongly forced by its coupling with the underlying bedrock [e.g. Rose, 1979,Engelhardt et al., 1990]. Changes in friction at this interface can drastically affect ice flux ratesby allowing ice to slide [Budd et al., 1979]. A faster flow discharges more ice from the centre ofan ice sheet to the sea. In Antarctica, regions of fast flowing ‘ice streams’ are characterised bylow basal friction and velocities of 100s of metres per year. Slower areas where ice is stuck tothe bed flow via viscous deformation at velocities of 10s of metres per year [Rignot et al., 2011,Morlighem et al., 2013].Figure 1: Ice flows from the centre of an ice sheet to the sea. Not to scale. [From [Bindschadler, 2014]]1.1 Feedbacks between ice sheets and subglacial drainage systemsWater between ice and the bed is thought to be a fundamental factor governing changes inbasal friction [Weertman, 1957, Iken and Bindschadler, 1986, Alley, 1989]: more water under11. INTRODUCTIONice at higher pressure decouples the ice from the bed, allowing it to slide more easily. Due to itsinfluence on the flow of ice, subglacial hydrology has been recognized as an important aspect ofglaciology [Flowers, 2015]. While most alpine glaciers show an evolution in subglacial drainagesystems through seasonal change in water input, Antarctic drainage systems are completelydisconnected from seasonal surface conditions [Ashmore and Bingham, 2014]. As a result, ob-served changes in basal friction or subglacial drainage are thought to be caused by varying basalconditions or by feedbacks within the ice sheet. There have been attempts to describe changesin ice dynamics through a variety of such feedbacks, for instance, the feedback driving fastflowing ice streams [Fowler and Johnson, 1996, Bennett, 2003, Kyrke-Smith et al., 2014], thespatial patterning of ice stream formation [Fowler and Johnson, 1996], ice stream stagnation[Alley et al., 1994], or reoccurring subglacial floods [Wingham et al., 2006].There is indirect evidence that bed conditions of areas of the Antarctic ice sheet may havea striped pattern, with stripes being oriented at right angles to ice flow. The clearest evidencefor such patterning comes from Sergienko and Hindmarsh [2013], and Sergienko et al. [2014].Through inverse modelling Sergienko et al. [2014] discovered distinctive laterally symmetric pat-terning in the driving and basal stresses in Antarctica and Greenland. They found the patternsto be co–located with highs in the gradient of hydraulic potential, and subsequently suggestthat subglacial water may control the evolution of these high shear stress ribs. It is feasiblethat such patterning could be formed by changes in basal slipperiness, potentially associatedwith variations in basal water pressure. Sergienko and Hindmarsh [2013] have considered howan instability in coupled ice sheet, basal sediment and subglacial drainage system flow couldcreate such patterning, but without an in–depth analysis.In a separate, theoretical study Sergienko and Hulbe [2011] explore the feedbacks between theice and subglacial drainage. They present a two–dimensional model that provides an explanationfor the common close coexistence of subglacial lakes and topographical features or sharp changesin basal conditions (Fig. 2). Sergienko and Hulbe prescribe a sticky spot and model numericallythe resulting changes to ice flow. The sticky spot causes a thickening of ice followed immediatelydownstream by ice thinning. Sergienko and Hulbe models the consequential subglacial water21. INTRODUCTIONflow and thermodynamic effects to find increased basal melt over the sticky spot and waterpooling immediately downstream. The pooled water could be interpreted as a lake the modeldoes not define a difference. Sergienko and Hulbe [2011] note that a coupled sticky and slipperyspot form the same feedback irrespective of their order in the flow direction. This suggests thatsuch a feedback could be repeated periodically, with banded slippery and sticky regions thatforce ice or subglacial water into the neighbouring region, respectively.Figure 2: Schematic depiction of ice–stream flow over ‘sticky spot’–lake– regular bed. If variations in basal to-pography are small, the gradient of hydraulic potential Φx is primarily determined by the ice–thickness gradient.The ice–thickness gradient is due to variations in the ice–flow velocity caused by the variations in basal traction.Note that the depicted situation is independent of a direction of ice flow. [From Sergienko and Hulbe [2011]]In an empirical study based on multiple remote–sensing techniques Fricker et al. [2010] foundthree lakes that support the results of Sergienko and Hulbe [2011]. The lakes are co–located withsticky spots, which are manifested by the slowing and thickening of the ice followed downstreamby thinning and speed-up of ice over the lakes. Fricker et al. [2010] hypothesise that all threelakes are likely formed as the result of feedbacks with basal shear stress, hydraulic gradientsand ice flux and as described by Sergienko et al. [2007], in a mechanism that corresponds tothat modelled in Sergienko and Hulbe [2011].1.2 Description of instabilityIn this thesis, we expand on the work of Sergienko and Hulbe [2011] and Sergienko et al. [2014]in an attempt to show that feedbacks between ice flow and subglacial drainage can cause the31. INTRODUCTIONFigure 3: Basal shear stress derived for MacIS using the method of Joughin and others (2004) with image-enhanced DEM; outlines of MacIS subglacial lakes are overlaid. [From Fricker et al. [2010]]spontaneous formation of patterned spatial structure and even lake formation. By using acoupled one–dimensional ice sheet–subglacial drainage model we show that such an instabilitycould exist, and deduce a mechanism for the instability that in essence corresponds to thatdescribed by Sergienko and Hulbe [2011], but with a key difference. Our model shows thespontaneous growth of the structure described, which relies on a slight phase shift betweenstriped basal conditions and oscillating ice thickness: regions of faster flowing ice must extendinto ice humps to make them grow.The instability forms as a result of feedbacks in the coupling of ice sheet and subglacialdrainage. Under an ice sheet, water flows both downhill and away from the weight of the ice.A wet region or lake can form where water flows downhill to fill up a topographical low, butcan also form where water flows away from heavier (thicker) ice to fill up a low in hydraulicpressure. Clarke [2005] and Iken and Bindschadler [1986] show that changes in ice thicknesshave a much greater effect on subglacial hydraulic gradients than changes in bed slope. Whilethe weight of a changing ice sheet affects the distribution and routing of subglacial water, thewater changes the flow of ice due the effect of the changes of water pressures on the friction of41. INTRODUCTIONsliding of the ice. Lastly, the speed of the ice and ice thickness are linked through conservationof mass: in a region where ice slows, ice will thicken. Where ice speeds up, it will thin. Throughthese interactions, the subglacial drainage system is intrinsically coupled to the ice sheet.Few previous studies have looked at feedbacks in this coupling in detail. This may be anissue of scales: It is common to assume that drainage is much faster than ice flow, and so whilevelocity and drainage are tightly coupled, it is easy to ignore the effect of an evolving surfacetopography on drainage in detail, assuming instead that ice geometry can be treated as fixedat the time scale on which the drainage system evolves [e.g. Hindmarsh, 1998b, Fowler, 2000,Schoof et al., 2014].1.3 Subglacial drainageMost established knowledge of subglacial hydrology comes from theoretical models which arebased on existing physical ideas of groundwater flow, that have been applied and modified tosuit subglacial hydrology in alpine glaciers, and later adjusted to ice sheets [Flowers, 2015]. Herewe list a few of the most studied, idealized components of the subglacial environment, describedin more detail by Flowers [2015]: Cavities are spaces that open in the lee of bedrock protrusions,formed as the forward motion of ice opens the space faster than it creeps closed. Sheets areidealised layers of water flowing under the ice that could describe a separation between the iceand bed but can also be an averaged description of linked cavities. Porous flow is the flowof water through porous subglacial till. R-Channels and canals are enclosed rivers that havemelted or eroded into the ice or underlying till, respectively, and evolve as water supply changes.Full drainage network models integrate these individual components, attempting to de-scribe the evolution of subglacial drainage and its effect on basal friction and ice flow [e.g.Flowers and Clarke, 2002, Kessler and Anderson, 2004, Schoof, 2010, Hewitt, 2011, Hewitt et al.,2012, Schoof et al., 2012, Hewitt, 2013, Schoof et al., 2014]. In particular, such models attemptto describe the rapid spatial and temporal changes in basal sliding by capturing changes inthe efficiency of drainage networks. An inefficient drainage system is modelled as a spatiallydistributed network of interconnected cavities, sheet flow or flow through porous bed substrate,51. INTRODUCTIONand is typically associated with low water flux and low basal friction. Higher water flux leadsto efficient drainage, modelled as an arterial network of channels.Changes in the flow of the Antarctic ice sheet have been attributed to changes in subglacialwater distribution and drainage systems [Engelhardt et al., 1990]. The majority of research onsubglacial drainage in Antarctica is focused on ice streams: dynamic, fast flowing regions of ice.Although ice streams only make up a small fraction of the area of Antarctica, they contribute90 % of the total dynamical ice loss from the continent [Bamber et al., 2000, Rignot et al.,2011]. Consequently, modelling their dynamics is crucial to calculating Antarctic’s ice massbudget. The high velocity of ice streams is associated with the presence of basal meltwater, anddeformable, saturated sediment slurries, which allow the ice to slide with little friction. Thishas been confirmed with research detailing cores, radar and seismic surveys [Blankenship et al.,1986, Alley et al., 1986, 1987, Robin et al., 1970, Engelhardt and Kamb, 1997, Kamb, 2001,Hodson et al., 2016].Drainage systems across the rest of Antarctica are thought to be variable. Difficult to ac-cess beneath kilometres of ice, empirical evidence of Antarctic drainage systems come fromobserving paleo–ice sheet beds, remote sensing, inverse modelling and sparse borehole drilling[e.g. Ha¨ttestrand, 1997, Engelhardt and Kamb, 1997, Lewis et al., 2006, Schroeder et al., 2013,Morlighem et al., 2013, Simkins et al., 2017]. Bed conditions such as intermittent pools of wa-ter or saturated sediments, lakes, channelisation of water, and areas with frozen beds withno drainage system have all been found to be widespread [Na¨slund et al., 2000, Carter et al.,2009, Langley et al., 2011, Schroeder et al., 2013, Young et al., 2016]. To encapsulate variedand often unknown drainage conditions, all large scale models of subglacial drainage in Antarc-tica have assumed that areas of sliding have a sheet or porous flow [Flowers, 2015]. This isa broad assumption that aims to approximate bed conditions and is normally based on twofindings: Firstly, Alley [1989] predicted that if there was no channelised flow, sheet flow woulddominate drainage systems under ice sheets. Later, Alley [1996] showed that in the absence ofinformation about basal conditions, modelling drainage as sheet flow is a reasonable approxi-mation for a wide variety of bed conditions. A pure sheet flow model can reasonably accurately61. INTRODUCTIONdescribe relationships between basal sheer stress and water pressure, but fails to show morecomplex evolution in the subglacial drainage, for instance changes in drainage efficiency or bedsubstrate.Canals are thought to be the most representative model for ‘discrete’ drainage elementsin Antarctica [e.g. Walder and Fowler, 1994, Simkins et al., 2017], though compared to otherdrainage models, their theory is perhaps the least developed. As opposed to R–channels, canalsare thought to exist at low effective pressure [Walder and Fowler, 1994] defined as overburdenminus basal water pressure. Model canals support a relationship with effective pressure de-creasing as discharge increases, whereas R–channels support the opposite [Ro¨thlisberger, 1972,Shreve, 1972, Walder and Fowler, 1994, Ng, 1998]. It is believed that this negative relationshipis a requirement for the feedbacks that drive ice streams, whereby low effective pressures andfast sliding lead to more frictional melt [Fowler and Johnson, 1996], supporting the existenceof canal–like or distributed drainage systems over R–channels under ice streams. R–channels,because of their tendency to increase effective pressure as discharge is increased, favour slowice flow. Additionally, it is not clear that the low hydraulic gradients and low melt rates underAntarctica are sufficient to support persistent R–channels; though R–channels have been arguedto exist at the outlets of draining lakes [Evatt et al., 2006, Pattyn, 2008], as well as under theshear margins of ice streams, where there are high melt rates [Perol and Rice, 2011].1.4 Subglacial lakesMost empirical observations of subglacial water in Antarctica are related to documenting the ex-istence of subglacial lakes [e.g. Carter et al., 2007, Wright and Siegert, 2012, Siegfried and Fricker,2018]. Airborne radar surveys first revealed subglacial lakes in the 70s [Robin et al., 1970]. Sincethen at least 379 subglacial lakes have been found [Wright and Siegert, 2012] with radar–echosounding reflections [e.g. Carter et al., 2007] or through repeated altimetry that shows local-ized drops or uplifting in the ice surface that are the result of lakes draining or filling [e.g.Siegert et al., 2005, Gray et al., 2005, Wingham et al., 2006, Fricker et al., 2007]. These obser-vations have led to a new understanding of subglacial lakes and their role in ice sheet dynamics.71. INTRODUCTIONIn particular, lakes can fluctuate in volume through quick or slow flooding and refilling events.As upstream lakes drain we can see lower lakes fill [Gray et al., 2005, Wingham et al., 2006],and rapid acceleration of ice in the drainage area [Stearns et al., 2008]. These results parallelfindings from the heavily documented jo¨kulhlaups or ‘Glacial Lake Outburst Floods’ seen inalpine glaciers [Bjo¨rnsson, 2003] and in the paleo-ice sheet record [Shaw and Gilbert, 1990].Most of the documented, identified lakes are large permanent features formed by topograph-ical basins [Wright and Siegert, 2012]. In addition, more detailed radar surveys continue to findmore water–storage features that are not topographical lakes. In radar surveys Carter et al.[2007] and Young et al. [2016] find numerous reflectors they suspect may indicate saturatedsediments or “swamps”. Similarly, Lake Whillans is an poorly defined lake that fills and drainsdynamically, and is formed in a very shallow basin ≈ 6m deep [Horgan et al., 2012].1.5 Patterns in paleo–ice sheet bedsTo date, there is no evidence of water storage structures forming quasi–periodic patterns. Theonly confirmed examples of spatial patterning in ice sheets are found in the bed of paleo–icesheets, which present a variety of well documented sediment patterns, for instance, drum-lins or bed ribs [e.g. Kinahan and Close, 1872, Ha¨ttestrand, 1997]. In order for bed topog-raphy to be transmitted to the surface it must have a wavelength greater than ice thickness[Gudmundsson, 2003]. It is therefore plausible that these sediment patterns could give rise tosignificant surface expressions on the ice sheet such as those on which Sergienko et al. [2014]build their work. Hindmarsh [1998b,a], Fowler [2000], Schoof [2007], Dunlop et al. [2008], andFowler and Chapwanya [2014] use similar models to describe pattern formation through sedi-ment transport. This research finds instabilities in a coupled ice–till environment, modelling thetwo as deformable viscous fluids divided by a shearing layer, where ice thickness and velocityare prescribed. Using a deformable till sliding law independent of effective pressure, they pro-duce patterns that could relate to documented sediment patterns, but results are inconclusive[Schoof, 2007]. To date, few studies have modelled instabilities formed from feedbacks with theice sheet itself.81. INTRODUCTION1.6 Thesis structureThe objective of this thesis is to describe a pattern forming instability in the coupling of ice sheetand subglacial drainage, and to show that such an instability could cause lakes to spontaneouslyform. In doing so, we show that such an instability could feasibly exist, describe the structurethat it forms and identify a mechanism that describes the feedback: A patch of high waterpressure, or low effective pressure, at the bed can lead to faster ice flow that causes ice to pileup downstream of the patch, leading to a reduction in hydraulic gradient and hence a reductionin water flux that, in turn, causes water to accumulate where effective pressure is already low.We define a subglacial lake as a subglacial region of water with an effective pressure of zero thatcontinues to accumulate water.In this thesis, we first introduce the coupled ice sheet, basal drainage model and describeits parts in detail. To find parameter regimes in which the model’s steady state is unstable,we carry out a linear stability analysis on the full model and find a dispersion relation whichdescribes the growth rate of different modes of the associated Fourier series solution. Next, byreducing the model to describe certain balances of terms, we force unstable solutions for certainparameter regimes. We find that the dependence of permeability on effective pressure createsa stability boundary through which the instability can be ‘switched off’. Observing unstablesolution modes, we deduce the mechanism for the instability and describe the physical feedbacksthat force the growth of solutions. Lastly, we present a finite difference scheme that simulatesthe full nonlinear model. Using operator splitting, the simulation confirms that the behaviourof the nonlinear solution follows the linear where the solution is small, then shows that theevolution of the unstable nonlinear solution continues to move away from the fixed point untileffective pressure reaches zero, when the model breaks down. This result describes the evolutionof spatially patterned structures with localised lake formation.92. 1D COUPLED ICE SHEET–DRAINAGE MODEL2 1D coupled ice sheet–drainage model2.1 Introducing the modelWe model a simple one–dimensional ice sheet coupled to a basal drainage system using threeequations. x is the position along the flowline in the downstream direction, oriented parallelto the bed, while t denotes the time passed (Fig. 4). The model’s three dependent variablesare ice velocity u(x, t), ice thickness h(x, t) and effective pressure N(x, t) — defined as thedifference between the overburden ice pressure and the pressure of subglacial water. Our spatialdomain covers a subsection of the flowline, and is small relative to the size of an ice sheet, butlarge relative to ice thickness. This scale allows us to make certain simplifying assumptions:by modelling a relatively short subsection of the ice sheet, we can apply periodic boundaryconditions, assuming that the next similarly–sized subsection of the ice sheet effectively behavesthe same as the subsection under consideration. The large length–to–depth ratio allows us toassume vertical velocities are negligible compared to the horizontal. Negligible vertical velocitiesallow us to use a depth–integrated ice flux divergence.Figure 4: The domain is a one dimensional section of ice along the flowline, where x is distance along the bed inthe downstream direction. Three variables ice thickness h, ice velocity u and effective pressure N are modelledalong the domain. θ is the angle of the bed from the horizontal.We also assume that sliding dominates internal deformation as a mechanism for ice transport.This justifies our assumption that the velocity u depends only on horizontal position x and timet, but not on distance above the bed. The appropriate depth–integrated form of the Stokes102. 1D COUPLED ICE SHEET–DRAINAGE MODELequations in this case is [MacAyeal, 1989, Schoof and Hindmarsh, 2010]4(ηhux)x − τb(u,N)− ρgh(hx − sin θ) = 0. (1)4(ηhux)x is the gradient of depth–integrated extensional stress in the ice, where η is viscosityand subscripts x and t denote partial derivatives. τb(u,N) is friction at the base of the ice sheet,which we assume depends on sliding velocity u and effective pressure N . The justificationfor such a friction law can be found either in the classical theory of sliding over hard beds[Lliboutry, 1968, Budd et al., 1979, Bindschadler, 1983, Fowler, 1986, 1987] or in models ofbasal motion facilitated by a layer of water–saturated, deforming till [Blankenship et al., 1987,Boulton and Hindmarsh, 1987, Kamb, 1991, Iverson et al., 1998]. Common to specific frictionlaws describing these two modes of sliding is that friction is a non-decreasing function of slidingvelocity (∂τb/∂u ≥ 0), and that friction increases with effective pressure (∂τb/∂N > 0). Acommonly used (though by no means unique) model that satisfies these constraints is the frictionlawτb = CuaN b,with b and C positive and a non–negative. This is sometimes referred to as a ‘GeneralizedWeertman Sliding Law’. The dependence on effective pressure will be key to our model. Forhard beds, a lower effective pressure implies water is more easily able to force its way between iceand bed, causing larger water–filled cavities to form. These reduce ice–bed contact and thereforereduce the amount of friction at the bed. For deformable beds, a lower effective pressure reducescontact forces between sediment grains and allows them to move past each other more easily.The last term in equation (1) is the ‘driving stress’, equal to the depth–integrated sum of thedownslope component of gravity and bed–parallel cryostatic pressure gradient. Here θ is thebed slope, ρ is the density of ice and g is acceleration due to gravity.We assume that accumulation and melting of ice are negligible at the length scales of interest,so that depth–integrated mass conservation can be written asht + (hu)x = 0. (2)112. 1D COUPLED ICE SHEET–DRAINAGE MODELLastly, we couple the ice sheet model with a subglacial drainage system by describing theconservation of water as it is flows under the ice:hwt(N) + qwx = 0, (3)qw = κ(N)ψ, ψ = ρwg sin θ +Nx − ρghx, (4)where qw and ρw are the flux and density of water respectively, and ψ is hydraulic gradient.Equation (3) assumes there is no water supply within the domain from the melting of ice orsurface runoff reaching the bed. Water flux is given by a form of Darcy’s law qw = κ(N)ψ, whichis used to describe fluid flow through a porous medium. Darcy’s law describes flux as beingforced by hydraulic gradient (ρwg sin θ +Nx − ρghx) and restricted by the depth–integratedpermeability κ(N) of the water–filled till column [e.g. Bear, 2013, Hillel, 2013]. When a water–filled till column is compressed, its porosity decreases. For flow through a granular medium,permeability depends on porosity so we expect κ to be a decreasing function of N (∂κ/∂N ≤ 0),for instance κ(N) = κ0e−mN . Changes in porosity are similarly related to hw, water storage perunit area of the bed: larger pores relate to more water storage. It follows that hw is assumed tobe a function of effective pressure hw = hw(N) which has a negative derivative ∂hw(N)/∂N < 0,for example hw(N) = hw0e−dN .The drainage model (3) is however not limited to describing flow through subglacial till, itcan also serve as a model for flow along the ice–till interface through a ‘sheet’ as proposed bySchoof et al. [2012] and Hewitt et al. [2012]. This is approximately true if the sheet thicknessevolves much faster than the water storage, where the sheet thickness (h in the notation ofthese papers) is a function of effective pressure N and sliding velocity u. If we use a simplifiedmodel in which sliding velocity is a function of effective pressure, then sheet thickness andthe resulting permeability are also decreasing functions of effective pressure. In any case, theminimal assumptions that our model rests on are that water storage at the bed or in basal icejust above the bed [Werder et al., 2013, Schoof et al., 2014] decreases with effective pressure (orequivalently, increases with basal water pressure), and that water flow is driven by the hydraulicgradient ψ, modulated by a permeability that can be either constant or decrease with effective122. 1D COUPLED ICE SHEET–DRAINAGE MODELpressure. Alley [1996] make equivalent assumptions about the prescribed relationships betweenflux, water storage, and the ability of the system to transmit water in a ‘sheet flow’ distributeddrainage model that is commonly used in large scale ice sheet models [Flowers, 2015].We have a minimal model in which ice flow, changes in ice thickness and the flow of water atthe bed are all coupled. Our primary aim is to determine whether these interactions can leadto the spontaneous formation of spatial structure. Take a spatial perturbation to an initiallyuniform effective pressure distribution. Consider a location where N is reduced. On its own, thiswould lead to reduced basal friction τb(u,N), but in practice, we expect that u will increase inthe same locations, either to counterbalance the effect of N on τb, or to provide sufficiently largeextensional stress gradients. A localized increase in u will lead to mass transport away from thesame location, presumably leading to a positive perturbation in ice thickness downstream of thelocation of reduced effective pressure. This surface perturbation in turn affects the hydraulicgradient that drives water flow. In fact, from this sketch, it is conceivable that a local depressionin N leads to a positive ice thickness gradient downstream, which corresponds to a reduction inhydraulic gradient ψ and therefore to a reduction in water flow out of the depression in N . Thisshould then cause water to accumulate in that depression, further lowering N . Our primaryaim in what follows is to investigate mathematically whether there are conditions under which apositive feedback of this type can indeed operate and lead to the emergence of spatial structurein basal water pressure, ice velocity, and surface elevation. We will do so by investigating thestability of the trivial, spatially uniform solution to equations (1)–(3) with respect to smallperturbations by means of a linear stability analysis. Subsequently, we will address whetherfinite perturbations can grow sufficiently for ice to begin to float, and separate from the bed.The remainder of this thesis is organised as follows: Firstly, we scale the model to a dimen-sionless form. Secondly, we analyze the stability of a steady state solution and find parame-ter regimes and wavelengths that support instability by performing a linear stability analysis.Thirdly, we examine features of a unstable linear solution and deduce the mechanism drivingthe instability. Lastly, we introduce a numerical solver and show the growing evolution of thefull nonlinear solution in a regime that forces instability.132. 1D COUPLED ICE SHEET–DRAINAGE MODEL2.2 Scaling the modelFor each variable in our model, we define a scale (or natural unit), denoted by the symbol forthe variable in square brackets, so [x] is the scale for x, [t] the scale for t etc. We then define adimensionless version of each variable by scaling asx∗ =x[x], t∗ =t[t], u∗ =u[u], h∗ =h[h], N∗ =N[N ], (5)τ ∗b (u∗, N∗) =τb(u,N)τb([u], [N ]), κ∗(N∗) =κ(N)κ([N ]), h∗w(N∗) =hw(N)hw([N ]). (6)[h] is prescribed from initial conditions as the mean ice thickness. Assume that the meanflux of water Q¯ through the domain is given, at least initially. This allows us to define theremaining scales throughτb([u], [N ]) = ρg[h] sin(θ),[h][t]=[h][u][x], κ([N ])ρwg sin θ = Q¯,[N ][x]= ρwg sin θ. (7)With [h] prescribed, these four relations define the remaining scales [N ], [u], [x] and [t](though to compute them explicitly, the functional forms of κ and τb need to be prescribed).The scheme for computing the remaining scales in principle is easy to describe: Given a meanwater flux Q¯, the drainage system permeability κ must accommodate that flux, correspondingto some value [N ]. A length scale [x] can then be computed if we require that the contributionof effective pressure gradients [N ]/[x] to hydraulic gradient be comparable to the direct effectof gravity. Given [N ], a velocity scale [u] can also be computed by balancing basal friction withthe driving stress, while [x] and [u] together define an advective time scale [t].We obtain the four dimensionless parametersδ =[h][x] sin(θ), =η[h][u][x]2ρg[h] sin(θ)γ =hw([N ])[x][t]Q¯, r =ρwρ. (8)It is straightforward to interpret the dimensionless parameters in the following way: is a142. 1D COUPLED ICE SHEET–DRAINAGE MODELdimensionless extensional viscosity, δ is the ratio of ice thickness to bed elevation change overthe length scale [x], and γ is the ratio of a water residence time scale to the advective time scale[t]. r is a material constant, with a fixed value close to unity. In dimensionless form, omittingthe asterisks on the dimensionless variables, the model takes the form:4(hux)x − τb(u,N) + h(1− δhx) = 0, (9)ht + (uh)x = 0, (10)γhwt(N) + (κ (N) (1 +Nx − rδhx))x = 0. (11)We denote the spatially uniform steady state solution to equations (9)–(11) by u(x, t) ≡ u¯,h(x, t) ≡ h¯, N(x, t) ≡ N¯ , where the only non-trivial equation in (9) – (11) that must be satisfiedby the constants u¯, h¯ and N¯ isτb(u¯, N¯) = h¯.Through our choice of scales, we have made sure that h¯ = 1, the mean ice thickness is scaledto unity. Additionally, we have set N¯ = 1 by choosing scales such that dimensionless N = 1corresponds to dimensionless κ(N) = 1, and by requiring a dimensionless mean water flux ofunity, where water flux is simply κ(N) in the absence of gradients in h or N . Lastly, we havealso ensured through our definition of the dimensionless friction law τb that τb(u¯, 1) = 1 is solvedby u¯ = 1, so the steady state is (u, h,N) ≡ (1, 1, 1), and in the steady state, κ = τb = 1.In the next two subsections, we will analyze the stability of this trivial steady state: ifa small perturbation is added to (u¯, h¯, N¯), under what circumstances will that perturbationgrow? What parameter regimes and wavelengths give rise to such growth, or instability? Toanswer that question, we will perform a linear stability analysis.153. INSTABILITIES3 Instabilities3.1 Linear stability analysisHere we analyze the stability of the steady state solution (u, h,N) ≡ (u¯, h¯, N¯) = (1, 1, 1) of thedimensionless coupled ice–flow–drainage model (9)–(11). We perturb the solution by writingu = u¯+ u′, h = h¯+ h′, N = N¯ +N ′, (12)where the primed variables are perturbations, assumed to be small compared with unity.We substitute these into (1) –(3), expanding functions τ(u,N) and κ(N) as first order Taylorseries approximations. Terms of O(u′2, h′2, N ′2) or smaller are omitted, resulting in the first order(linear) approximation of our model:0 = 4u′xx − τbuu′ − τbNN ′ + h′ − δh′x, (13)h′t = −u′x − h′x, (14)γhwNN′t = −N ′xx + rδh′xx − κNN ′x, (15)whereκN =∂κ∂N∣∣∣∣N=N¯, τbu =∂τb∂u∣∣∣∣u=u¯,N=N¯, τbN =∂τb∂N∣∣∣∣u=u¯,N=N¯.Note that by construction of the model, we assume that κN ≤ 0, τu ≥ 0, and τN > 0. Equations(13)–(15) are a set of linear partial differential equations with constant coefficients posed ona periodic domain, and can be solved by means of a Fourier series [Fourier, 1808]. Thus, werepresent the solution in the formu′(x, t) =∑kuˆ(t) exp(ikx), h′(x, t) =∑khˆ(t) exp(ikx), N ′(x, t) =∑kNˆ(t) exp(ikx),(16)where the sum is over k = 2pin/L and n is an integer in the range (−∞,∞).Substituting for u′, h′ and N ′ and subsequently using the orthogonality of distinct Fourier163. INSTABILITIESmodes, we obtain0 =(−4k2 − τu) uˆ+ (1− i) hˆ− τNNˆ , (17)dhˆdt= −ikhˆ− ikuˆ, (18)γhwNdNˆdt=(k2 − iκNk)Nˆ − rδk2hˆ. (19)This reduces the problem to solving a coupled set of linear, ordinary differential equations withconstant coefficients, of the formBv˙ = Av,where B and A are matrices, and v = (uˆ, hˆ, Nˆ)T is the solution vector. Except in the case of adegenerate spectrum, we can write the solution in the formv =∑ivi exp(σit),where each σi solves the generalized eigenvalue problemdet(A− σiB) = 0, (20)and vi is the corresponding eigenvector, satisfying (A− σiB)vi = 0.Because A is dependent on wave number k, we can find a relationship between growth rateσ and k, σ = σ(k), known as a dispersion relation. If the real part of σ is positive, then thecorresponding Fourier mode grows unstably, while a negative real part implies that it will shrinkto zero over time. The imaginary part of −σ/k gives the wave speed of the associated Fouriermode.A single unstable mode k (i.e. a single wavenumber k for which it is possible to find apositive growth rate Re(σ) for a given set of model parameters) signals instability, whereasall wavenumbers must correspond to at least non–positive growth rates to assure stability ofa steady state solution. In practice, the characteristic equation det(σB − A) = 0 defines aquadratic in σ, and is solveable in closed form. We solve explicitly for the dispersion relation173. INSTABILITIESand obtain two roots σ1(k) and σ2(k), given byσ1,2(k) =−b±√b2 − 4ac2a, (21)wherea = −4γhwNk2 − γhwN τu, (22)b = k2τu − kκNτui− γhwNki− k3κN4i+ 4k4 − δγh2wNk2 − γhwNkτui− γhwNk34i, (23)c = k2κN + k3i+ k2κNτu + k3τui− δ2k3κN i+ δ2k4 + 4k4κN + k54i− δk3rτN i. (24)Given a point in the space of parameters (, γ, δ, hwN , κN , τbu , τbN ) and wave numbers (k) wecan easily calculate the growth rates σ1 and σ2. However, the inverse of this problem, findingthe subset of parameters and corresponding wavenumbers that generate an instability (at leastone root with Re(σ) > 0) is not straightforward. The argument of the square roots in (21)are complex, and although a closed form expression for Re(σ) can be derived, this involves thesquare root of a combination of terms involving up to the fifth power of k.Consequently, we further simplify the model to focus on a specific balance of terms in themodel, or equally, balance of processes, that produces an unstable solution. In finding anunstable parameter regime, we can identify conditions suitable for generating an instability,and are able to deduce the physical mechanism driving that instability. We will then returnto the full dispersion relation (21) to test the predictions of the simplified model we are aboutderive.3.2 Reduced model with unstable steady stateOur original model is based on the horizontal length scale of interest [x] being short comparedwith the length L of the ice sheet, so that the tendency of the ice sheet to thin as the marginis approached does not manifest itself at the scale under consideration, and we can treat thedomain as being approximately periodic. We have allowed for a finite bed slope θ in orderto drive the flow at this scale. In the unperturbed steady state, θ is then also the surfaceslope, and for a realistic ice sheet, the surface slope is comparable with [h]/L. One of our scale183. INSTABILITIESrelations, equation (7), requires [N ]/[x] = ρwg sin θ. Setting θ = [h]/L, we have, with small θ,[N ] ≈ ρwg[h][x]/[L], and the requirement that [x] be much less than L translates, roughly, into[N ] being much less than the ice overburden ρig[h], which is typically required for any kind ofappreciable sliding [Schoof, 2006].Given [x] L is then equivalent to [x] sin(θ) [h], and consequently δ 1, all things beingequal, an O(1) variation in h∗ over an O(1) displacement x∗ would lead to a thickness gradientthat overwhelms all other terms in the force balance (11). This apparent inconsistency is themotivation for a rescaling of (9)–(11) that we perform below. We make our assumption morespecific than simply putting δ 1. Recall that γ is the ratio of a water residence timescale toan ice advection time scale. In many drainage systems with higher discharge than in Antarctica,water flows much faster than ice (e.g. Chandler et al. [2013]), implying that γ 1. We assumethat the same still holds true for Antarctica, and focus on the specific parameter regimeδ−2 γ δ−1 1.We also assume that the extensional stress parameter is relatively large (meaning that slidingvelocities are large enough as to make extensional stresses large), and specifically assume thatγ2δ2 = O(1). The choice of this parameter regime is of course not accidental: as we will seeshortly, it leads to a simpler version of the original model in which an instability mechanism isrelatively straightforward to identify.In the parameter regime we have identified, the following re–scaling is appropriate:x∗∗ = γδx∗, t∗∗ = γδ2t∗, u∗∗ = u∗, h∗∗ = h∗, N∗∗ = N∗. (25)Substituting the rescaled variables and dropping the double asterisks, we obtain:193. INSTABILITIES4γ2δ2 (hux)x − τb(u,N) + h(1− δ2γhx)= 0, (26)δht + (hu)x = 0, (27)hwt(N) +(κ (N)(1γδ+Nx − rδhx))x= 0. (28)To arrive at this model directly by scaling equations (1)–(4), without rescaling through equations(9)–(11) we would use the following scales:x∗∗ =x[x]′, t∗∗ =t[t]′, h1 =h1[h1]′, (29)[x]′ =κ([N ])[N ]2hw([N ])[h][u]ρwg, [t]′ =κ([N ])[N ]3hw([N ])[h]2[u]2ρ2wg2, [h1]′ =[N ]ρwg. (30)In equations (26)–(28), as δ2γ 1 and δ 1, we can see that to leading order h(x, t) isconstant: (26) shows to leading order that hx = 0 and (27) shows that ht = 0. As a result, weexpand h(x, t) as h(x, t) = h0 + δ−1h1(x, t), where h0 is then a constant. It is then immediatelyobvious that h0 is also the steady state ice thickness, and by our choice of original scales, wecan set h0 = h¯ = 1.Applying the same approach to equation (28), we obtain at leading order that∂κ∂N∂N∂x= 0.This allows for two conclusions: either N is spatially uniform, or ∂κ/∂N vanishes at leadingorder. The former is clearly not interesting, as a uniform N would prevent a positive feedbackbetween flow, thickness evolution and drainage to arise, in the sense that there cannot be anydrainage–driven spatial variations in friction τb(u,N) in that case. Consequently, we assumeinstead that permeability κ is only weakly dependent on effective pressure, and specifically, thatwe can put κ = κ0 + γδκ1(N) with κ0 constant. By our choice of scales, it then follows thatκ0 = 1.203. INSTABILITIESAt leading order in u and N and first order in h, the model then reduces to4γ2δ2uxx − τb(u,N) + 1 = 0, (31)h1t + ux = 0, (32)hwt(N) + κ1(N)x +Nxx − rh1xx = 0. (33)Note that the only nonlinear part of these equations are the constitutive laws that we prescribe:τb(u,N), hw(N), and κ1(N). If we were to choose, somewhat unrealistically, linear functions,the resulting model would be directly solvable by Fourier series without further simplification.We carry out a linear stability analysis again as described in Section 3.1, treating (u, h1, N)as the dependent variables of the model and perturbing around the steady state (u, h1, N) ≡(1, 0, 1). Following the same steps as in section 3.1 we get the Fourier domain representationof the model. For a single Fourier mode with wavenumber k and assuming exponential timedependence with coefficient σ, we write (u, h1, N) = (uˆ, hˆ, Nˆ) exp(ikx+ σt) and obtain−(4γ2δ2k2 + τu)uˆ− τNNˆ = 0, (34)σhˆ+ ikuˆ = 0, (35)hwNσNˆ − k2Nˆ + rk2hˆ+ ikκNNˆ = 0, (36)whereκ1N =∂κ1∂N∣∣∣∣N=N¯.The corresponding dispersion relation isσ(k) =−b±√b2 − 4ac2a, (37)a = hwN , b = iκ1Nk − k2, c =irτNk34γ2δ2k2 + τu, (38)4ac = 4hwNirτNk34γ2δ2k2 + τu. (39)Next, we determine analytically the behaviour of <(σ) near the origin k = 0, and for largek. The former allows us to show that an instability can indeed occur, while the latter ensuresthat the model does not predict any pathological behaviour, with short–wavelength oscillations213. INSTABILITIESbeing damped rather than growing.For small k, the solution displays positive growth rates. Assuming κ1N = 0:<(σ(k)) ∼ ±12√−2rτNk3hwNτu, k → 0. (40)hwN < 0 and so the argument of the square root is positive. <(σ) has a positive and negativesolution, denoting an instability (this limit is plotted next to <(σ) in Figure 5). Next we showthat the solution has no positive growth rate for large k. We can find a Taylor series solutionfor the growth rate by expanding the square root in (37). From (37)–(39) it follows thatσ =−b2a± b2a√1− 4acb2, (41)and from (27) we see that when k → ∞ it is also true that 4ac/b2 → 0. As a result we canexpand the square root from (41) in Taylor series about one. Expressed in full:σ(k) =−k22hwN[−1∓(1− 124hwN irτNγ2δ2k3− 14(4hwN irτN4γ2δ2k3)2− 116(4hwN irτN4γ2δ2k3)3− ...)], (42)where the ± sign is flipped as the negative b (for k → ∞) is factored out of the square root.We have two cases: picking the upper sign of ∓ in (42) gives σ ∼ k2/hwN at leading order,which is negative as hwN < 0. Picking the lower sign in ∓ in (42) for the second root, we getσ ∼ irτN4γ2δ2k− 2hwNr2τ 2N/(4γ2δ2k)2 which has a positive real part, but rapidly approaches zeroas k →∞ indicating a stable solution for large wavenumbers.Figure 5 confirms these results are reflected in the full dispersion relation, equations (21)–(24), by selecting numerical values for parameters that satisfy the conditions for instability andcalculating the growth rates <(σ1) and <(σ2) over a range of wavenumbers k. Alongside, wecompute growth rates from the reduced model. For the parameters regimes described, we cansee instability (σ > 0) for small k. As k grows, the positive growth rate approaches zero.3.3 A stability boundaryIn Section 3.2 we set κN 1. We next plot the growth rate for a range of κN ≤ 0, showingthat a small κN maintains an instability while a more negative κN will switch off the instability.223. INSTABILITIESFigure 5: With parameter values = 105, γ = 10−7, δ = 105, hwN = −e−1, κN = 0, τbu = 1/2, τbN = 1/2, wecalculate growth rates σ1 and σ2 for a range of wave numbers for both the full and reduced dispersion relationsequations (21)–(24) and (37)–(39) respectively. The full dispersion relation is adjusted to have the same scales asthe reduced. Both dispersion relations show near identical dispersion relations. An approximation of <(σ)→ 0from equation (40) is shown to accurately depict σ for small k.This is only true for the full dispersion relation (21)–(24).We can approach the case of finite κ1N in a similar way by looking at the limiting behaviourof the real part of the reduced dispersion relation for small k.<(σ1) ∼ −rτNκN 1τuk2, <(σ2) ∼(rτNκN 1τu+1hwN)k2, k → 0. (43)Although the leading order behaviour near k = 0 is quite different now, we see that theinstability invariable persists: one root (σ1) always possesses a positive real part. For each root,we also see that there is an imaginary part=(σ1) ∼ k3(rτNκN 21τu+hwNr2τ 2NκN 31τ2u), =(σ2) ∼ κN 1khwN, k → 0. (44)corresponding to waves travelling at velocity −=(σ/k) along the x-axis.In equation (43) the growth rates shrink as |κ1N | grows, but this equation only approximatesgrowth rates for |κ1N | ∼ 1. In order to see stabilization, we need to go to much more negativevalues of κN . In Section 3.2 to reduce the model, we made the assumption that permeability κ is233. INSTABILITIESonly weakly dependent on effective pressure N , which allowed us to make significant analyticalprogress. As we have seen, there is no stablization for the weak dependence of κ onN representedby κ1N (note that this is not the same as κN , in fact, κ1N = κN/γδ), and so the reduced modelis always unstable.To test whether we obtain stabilization for even more significant dependence of permeabilityon effective pressure, in Figure 6 we plot dispersion relations obtained from the full model withthe same parameter values as in Figure 5 but using a range of values of κN . We find that as|κN | grows, the σ(k) becomes negative for all k, and the instability is ‘switched off’.Figure 6: With parameter values = 1, γ = 10−3, δ = 102, hwN = −e−1, τbu = 1/2, τbN = 1/2, we calculategrowth rate for a range of wave numbers k ∈ (0, 0.3) with a range of values for κN ∈ (−1.26,−0.01) using thefull dispersion relation, equation (21). We can see that unless κN is small, we do not obtain positive growthrates. Specifically, with this parameter regime if κN < −0.45 the growth rate indicates a stable steady state andno instability.We have shown that the steady state is unstable if and only if κN is small. As discussedin section (2.1), small κN represents a model with a stiff–pored bed where changes in effectivepressure have little effect on the permeability. This result was found through a linear stabilityanalysis: we next take a closer look at the analytic linear solution presented to show that a phaseshift between dependent variables provides a feasible explanation for the physical mechanism243. INSTABILITIESdriving the instability. Firstly, we describe the analytic solution and secondly we examine itsdetails.3.4 Discussion of the mechanism for instabilityThe analytic solution we find is a linear approximation only valid when close to the steadystate (u¯, h¯, N¯). Assuming a spatially periodic domain with period L, we have represented theperturbation away from the steady state (u¯, h¯, N¯) as a Fourier series in the spatial variable x.We find that each Fourier mode is the superposition of two eigenfunctions of the eigenvalueproblem (13)–(15). Specifically:u(x, t)h(x, t)N(x, t) =u¯h¯N¯+∞∑n=−∞(c1nv1neσ1(kn)t + c2nv2neσ2(kn)t)eiknx, (45)where kn = 2pin/L, v1 and v2 are the eigenvectors that solve equation (20), and σ1 and σ2 arethe corresponding eigenvalues.We can impose initial conditions for h(x, 0) and N(x, 0), wheras u(x, 0) must satisfy theelliptic equation (33). To find c1nv1n and c2nv2n, we project h(x, 0) and N(x, 0) onto theeigenbasis for each Fourier mode n. To do this we first Fourier transform the initial conditionsh0n =1L∫ L0h(x, 0)e−iknxdx, N0n =1L∫ L0N(x, 0)e−iknxdx,and next solve the following algebraic equation for c1 and c2 for each Fourier mode n, using v1and v2 found in the eigenvalue problem.h0n = c1nv1n,2 + c2nv2n,2,N0n = c1nv1n,3 + c2nv2n,3,where v1n = (v1n,1, v1n,2, v1n,3) and v2n = (v2n,1, v2n,2, v2n,3). With c1 and c2 we can solve forall variables for all x and t. A particular Fourier mode of the solution has a wavelength of253. INSTABILITIES2pi/k, growth rate of <(σ), and a wavespeed of −=(σ/k). The solution to equation (45) is thesum of two eigenfunctions v1neiknx and v2neiknx. For a particular eigenfunction and Fouriermode, the phase shift θ of u(x, 0), h(x, 0) and N(x, 0) relative to eiknx is equal to the respectivecomponents of arg(v1n). Using unstable parameter–wavenumber regimes with the full linearmodel (13)–(15), we see different phase shifts between dependent variables for the unstableeigenfunction compared to the stable eigenfunction, shown in Figure 7.For an unstable solution, velocity and effective pressure are found to be roughly antiphasewhile a maximum in ice thickness is slightly upstream from the velocity minimum (see Figure7). From this we can provide a mechanism for the feedback that creates unbounded growth ofthe solution’s amplitude over time: a hump of ice slows water flow in the positive x–direction,with that slowing being most pronounced where the slope perturbation h′x is large and negative.This causes water to accumulate upstream of locations of negative slope perturbations, that is,near minima of h. If these correspond to minima in N , the water accumulation causes a furtherdrop in N , causing the pertubation in N to grow. The low effective pressure reduces ice–bedfriction increasing sliding and creating a high in velocities in this area. Note that although Nand h are approximately in phase, N lags h somewhat. This then causes u to lead h in phase:velocity maxima occur upstream of maxima in h, causing more ice to flow into the hump of icethan out and so the hump also grows. Figure 8 describes this in a schematic.The second, stable eigenfunction reverses the phase shifts found in its unstable counterpart:h and N are in phase, so water flows away from minima of effective pressure, and u lagsh, causing more ice to flow out of the hump than in, flattening it out. The existence of thisunstable eigenfunction is however irrelevant. When the trivial steady state (u¯, h¯, N¯) is randomlyperturbed, the perturbation will contain non–zero contributions from the unstable eigenfunction,and the solution (u, h,N) will evolve away from the steady state. Figure 9 describes this in aschematic.263. INSTABILITIESFigure 7: Eigenfunctions and subglacial drainage properties of the linear solution show different phase shifts.Calculated from Equation (45) (full linear model) for a single dominant wavenumber k = 0.0553 with parametervalues = 1, γ = 10−3, δ = 102, hwN = −e−1, κN = 0, τbu = 1/2, τbN = 1/2. Blue ‘x’s and ‘o’s respectivelydenote maximums and minimums of the unstable solutions.274. NUMERICAL ANALYSISFigure 8: Schematic showing the phase shift between variables that forces a positive feedback. A low in effectivepressure (shown by larger pores) and high in ice velocity (at the left hand dotted line) lie upstream from a humpin the ice thickness (right hand dotted line). A decrease in velocity as ice flows through the ice hump will causethe hump to grow.4 Numerical analysisThe linear stability analysis we have just performed confirms that a trivial steady state solutioncan be unstable, in which case we expect spatial structure to evolve. We are not able todetermine what the final result of that instability is from a linearized theory alone. In particular,we want to determine whether the instability could lead to the spontaneous formation of largersubglacial water bodies: will the instability simply lead to more water being stored in till insome locations than others, or will it actually cause a macroscopic, water-filled subglacial cavityor lake to form? In order to answer that question, the full model must be solved, and this isonly possible numerically.Our numerical scheme uses operator splitting as seen in MacNamara and Strang [2016]. Letui, hi and N i denote the solutions for u, h and N at the ith time ti. By ‘splitting’ vi where iis the time step, we solve for each variable ui, hi, N i separately and sequentially. Suppose thatui−1, hi−1 and N i−1 are known, given hi−1 and N i−1, ui−1 can be computed from the solutionof the elliptic problem (9). At the ith time step, we first semi-discretise equation (10) using284. NUMERICAL ANALYSISFigure 9: Schematic showing a phase shift between variables forcing a negative feedback. A low in effectivepressure and high in ice velocity (right hand dotted line) lie downstream from a hump in the ice thickness (lefthand dotted line). An increase in velocity as ice flows through the ice hump will cause the hump to shrink.a backward Euler step to update h from hi−1 to hi, treating u = ui−1 as fixed. Subsequently,we semi–discretise equation (11) using a backward Euler step in N , from N i−1 to N i, treatingh = hi as fixed. Given the updated hi and N i, ui can then again be computed by solving theelliptic problem (9).We discretise in space using a finite difference scheme, defining a regular grid of points xispaced a distance ∆x apart. Subscripts then indicate the value of a variable at the correspondinggrid point, so for instance uij = u(xj, ti). Under the operator splitting scheme described above,equation (9) then becomes4hij+1/2(uij+1 − uij)− hij−1/2(uij − uij−1)∆x2− τb(uij, N ij) + hij(1− δhij+1/2 − hij−1/2∆x)= 0, (46)where we have defined h at half grid points by an algebraic average,hij+1/2 = (hij+1 − hij)/2.Our domain is periodic, and if there are n grid points, then values at a point with spatialsubscript j = n+ 1 are identified with their values at j = 1. Equation (46) is solved for the uijvalues. Note that (46) is nonlinear, and we use Newton’s method to solve for the uij, using the294. NUMERICAL ANALYSISvalues ui−1j from the previous time step as initial guesses where available. In practice, we use apower–law friction law of the form τb(u,N) = CuaN b.Equation (10) is semi–discretised in space using a centre difference scheme, and with theoperator splitting used here, the backward Euler step for h becomeshi+1j = hij −∆t2∆x(hi+1j+1ui+1j+1 − hi+1j−1ui+1j−1), (47)which is solved for the hi+1j , using the uij computed previously for given hij and Nij .The operator splitting then goes on to solve for effective pressure N . Equation (11) isdiscretised in space using a centred scheme, givingγhw(N i+1j)− hw (N ij)∆t+1∆xκ(N i+1j+1/2)(1 +N i+1j+1 −N i+1j∆x+ rδhi+1j+1 − hi+1j∆x)− 1∆xκ(N i+1j−1/2)(1 +N i+1j −N i+1j−1∆x− rδhi+1j − hi+1j−1∆x)= 0. (48)Here N i+1j+1/2 is again given by an algebraic mean asN i+1j+1/2 = (Ni+1j +Ni+1j+1)/2.We treat the hi+1j as given by the previous step in the operator splitting, and solve for the Ni+1j .As with equation (46), equation (47) is nonlinear in the N i+1j , and we use Newton’s method tosolve it, using the N ij as initial guesses. In practice, we use constitutive relations of the formhw(N) = hw0e−dN and κ(N) = κ0e−mN .304. NUMERICAL ANALYSISFigure 10: Both linear analytic (Eq. (45)) and nonlinear numeric solutions (Sec. 4) for effective pressure plottedfor comparison. Model simulated with = 4, γ = 10−5, δ = 104 and functions τb(u,N) = u1/2N1/2, hw(N) =e−N , κ(N) = e−0.1N . While perturbations from the steady state (1) are small, the simulated nonlinear andlinear solutions behave the same, confirming the reliability of the nonlinear numerical solver.4.1 Numerical ResultsWe next use the numerical model described above to compute solutions up to finite amplitude.We present the scaled solution simulated with two sets of parameter values, the first with = 4,γ = 10−5, δ = 104, hw(N) = e−N , κ(N) = e−0.1N , τbu = 1/2, and τbN = 1/2 (Fig. 11), and thesecond with κ(N) = e−0.16N (Fig. 13), a permeability relationship closer to the shut–down ofthe instability. The reliability of the solver is confirmed in Figure 10, where the linear solutionsis shown to describes the nonlinear system well while perturbations are small compared withthe mean.314. NUMERICAL ANALYSISFigure 11: Nonlinear evolution of an unstable solution for effective pressure shows the development of spatialstructure. The solution grows until ice reaches the point of flotation N = 0 where the model breaks down.Parameters are described in Section 4.1 with κ(N) = e−0.1N , solution simulated with numerical model describedin Section 4 .The permeability law used falls within bounds of empirical studies. Empirical studies findrelationships close to κ(N) = κ0e−mN , κ0 = 2.9 × 10−13m2 and m = 5.6 × 10−6Pa−1 forpermeability on a soft bed [Boulton and Dent, 1974, McKinlay et al., 1975]. When scaledfor the typically low values of effective pressure found in Antarctic ice sheets, < 100kPa[Engelhardt and Kamb, 1997, Alley, 2001], the relationship shows a permeability weakly de-pendent on effective pressure similar to the relationship we used. The power sliding lawτb = uaN b with a = b = 1/2 falls within commonly predicted values of 0 < a, b < 1 [Lliboutry,1968, Budd et al., 1979, Bindschadler, 1983, Fowler, 1986, 1987, Boulton and Hindmarsh, 1987,Kamb, 1991, Iverson et al., 1998]. Water storage (hw) has been found to vary between at least324. NUMERICAL ANALYSISFigure 12: Velocity, ice thickness, effective pressure changing water storage, water flux and water storage cor-responding to Figure 11 at the point in time where N = 0. We can attribute the growth in the solution to thephase shift between velocity and effective pressure as described in Section 3.4. Ice flux over the maximum in icethickness (x = 35), has negative gradient and so the hump continues to grow.334. NUMERICAL ANALYSIS2 mm [Engelhardt and Kamb, 1997] and 2 m [Tulaczyk et al., 2000, Engelhardt et al., 1990,Kamb, 1991]. Scaled for small [N ] and resulting hw([N ]), the magnitudes for water storage usedour simulation lie within this range.We next demonstrate that the parameters used reflect realistic variable scales. Empiricalstudies show Antarctic scales of [u] of order 1 to 100 m/a [Rignot et al., 2011], [h] of order1000 m [Fretwell et al., 2013], and [N ] in the order of 100 kPa [Engelhardt and Kamb, 1997,Alley, 2001]. [x] = 10 km is reasonable, and corresponds to a length scale of interest to theinstability [x]′ = 100 km (see equation (29)). The computational domain was set as 10[x]′. Usingthese scales, the size of δ = 104 used corresponds to a very flat bedslope θ of order 10−5. Thisis realistic for large swaths of Antarctica including ice streams [Fretwell et al., 2013]. A typicalice viscosity of η = 1014 Pa s [Greve and Blatter, 2009], and g = 10 [e.g. Reilly and Whiteford,1965] gives an of order 1 as used in the simulation. Recall that γ is the ratio of a waterresidence timescale to an ice advection time scale. This is equivalent to the ratio of ice to watervelocities u/uwater. With γ = 10−5, ice velocities of 1 to 100 m/a correspond to water velocitiesof 0.001 to 0.1 m/s. Chandler et al. [2013] found velocities of an order of magnitude higher (0.1to 1 m/s) from subglacial drainage fed by surface melt in Greenland. As we have assumed thatthere is no contribution to subglacial water from seasonal surface melt, and very low bedslopes,scaling for water velocities of 0.001 to 0.1 m/s is reasonable.The evolution of the solution sees the development of a wave of dominant wavelength thatcontinues to grow in amplitude reaching an effective pressure of zero. Most importantly, we seethat the instability is not bounded, and as expected, a pattern with a characteristic wavelengthemerges early in simulation. As the perturbations in N and u reach finite amplitude, there isstrong evidence of nonlinear effects. Rather than evolving as smooth sinusoids, the solutionevolves varying asymmetric spatial structure. We see localized ‘sticky spot’ regions with highN , slow u and thick h, followed immediately downstream by ‘lake’ regions with low N , large uand thin h, followed downstream by more steady flowing ice.The simulation is continued until N = 0 is reached, at which point the model breaks down.Mathematically, we cannot continue the computation because the friction law does not in gen-345. DISCUSSIONeral compute a real–valued τb. Physically, the model breaks down because N = 0 not onlycorresponds to bed friction vanishing, but also to ice separating from the bed: water then fullysupports the overburden, and a cavity between ice and bed can form — the instability can leadto the formation of a lake. The simulated solutions confirm our finding from subsection 3.2,that scaled hx is small compared to Nx or ux: the full amplitude of h at t = 0.14 is less than0.0002 whereas the amplitude of u and N are much larger, of order 1.The nonlinear solution shows a bifurcation predicted in the linear model (see Section 3.3).For κN < −0.29 the spatially uniform steady state solution is stable, the instability occursfor less negative values. For −0.29 < κN < −0.15, the solution evolves into an apparentlystable travelling wave solution without N = 0 being reached, shown in Figure 15. Only for−0.15 < κN does growth occur up to the point at which ice goes afloat. Figure 14 shows themaximum and minimum values of N exhibited by the travelling wave solution as a functionof κN . The difference between the maximum and minimum values of N could be interpretedas the ‘amplitude’ of the travelling wave, and that amplitude appears to grow smoothly as κNincreases past the critical value of −0.29 at which the bifurcation to instability occurs. Notethat in the numerical simulations used to find this bifurcation point, stability is dependent ondomain size. Unstable wavelengths must be a unit fraction of the domain. Because of this,we cannot compare exact bifurcation values from nonlinear simulations to those from lineardispersion relations described in Section 3.3.5 Discussion5.1 SummaryWe have looked at a coupled ice flow–melt drainage model in one horizontal dimension, withthe aim of identifying whether any pattern–forming instabilities can appear, and what thoseinstabilities might lead to. We identified one such instability that leads to the formation ofpatterned ‘sticky spot’–lake pairs, and described the mechanisms that force the instability: Apatch of high water pressure, or low effective pressure, at the bed can lead to faster ice flow355. DISCUSSIONthat causes ice to pile up downstream of the patch, leading to a reduction in hydraulic gradientand hence a reduction in water flux that, in turn, causes water to accumulate where effectivepressure is already low. This feedback is related to the relationship between subglacial lakesand ice surface expression described in Sergienko and Hulbe [2011], where an induced slipperyspot causes a pile up of ice, and downstream pooling of water. The structure predicted bySergienko and Hulbe [2011] matches that which is simulated in our nonlinear model.We have been able to confirm the existence of this feedback by using a systematically reducedmodel as described in Section 3.2. The feedback leads to an instability when the residence timefor water in the system is much shorter than that for ice, and when bed conditions support apermeability only weakly dependent on effective pressure. In addition, the bed gradient wasvery small, and changes in ice thickness were assumed to be very small compared to steady stateice thickness. Additional assumptions on ice sheet and bed conditions went into the formulationof the original ice sheet–drainage model, described in Section 2.1.While the reduced model we describe is always unstable, its formulation is restrictive due tothese assumptions. The complexity of the full model of Section 3.1 and corresponding dispersionrelation (21)–(24) makes it difficult to analytically map out regions in parameter space that cor-respond to instability. This can, in principle, be done numerically by solving simultaneously forthe wavenumber k at which the maximum of <(σ) is attained, and requiring the corresponding<(σ) to equal zero. These two equations in principle allow k to be eliminated and give a singlerelationship between all the model parameters (, γ, δ, hwN , κN , τbu, τbN) at neutral stability,that is, at the bifurcation point where the trivial solution becomes unstable.In Section 3.3 we discussed a stability boundary where the sensitivity of permeability κ toeffective pressure N is a mechanism through which the instability can be suppressed. Thereduced model with unstable steady state made the assumption that permeability varied onlyweakly with effective pressure. When solving for growth rates with the full model, we foundthat with large enough values of ∂κ/∂N , the instability eventually disappeared. The physicalmechanism for the suppression of the instability can be described in relation to phase shifts ofthe solution as in Section 3.4. With a more negative ∂κ/∂N , h continues to spatially lag N .365. DISCUSSIONHowever, in comparison to the unstable solution, u is in phase with N , and so leads h. As aresult, flux is greater on the downstream side of the ice hump and so the perturbation producesa negative feedback.5.2 Discussion of resultsThe results of the nonlinear simulation present a number of features of interest. Firstly, forspecific parameter regimes the nonlinear solution is unstable, and grows in amplitude until itattains a zero effective pressure at some point in space. This describes the ice sheet being atthe point of flotation: the weight of the ice is supported by the drainage system pressure, notby the bed. At this point the bed and ice sheet could begin to separate to form a deepeningsubglacial lake. In previous work modelling basal pattern formation by Hindmarsh [1998a] andFowler [2000] flotation also typically occurs, though through a different process. Interestingly,regions of low effective pressure or water accumulation lie just downstream of an ice hump, andso the instability appears to be closest to a sticky–spot feedback as described by Sergienko et al.[2007], as opposed to an ice dam feedback described in jo¨kulhlaups models [e.g. Bjo¨rnsson, 2003]which would cause water to accumulate upstream of thick ice.Secondly, the regions of the nonlinear solution where N becomes small and goes to zeroare quite localized, and are separated by much larger regions in which N has become elevatedrelative to the effective pressure in the spatially uniform steady state. The effect of the instabilityis that water is taken from larger areas and concentrated in smaller pockets. In the nonlinearsolution at the point in time where the model ends, illustrated by Figure 12, mean effectivepressure across the domain is 1.1696, above the steady state average. This does not, however,lead to a net slow–down in ice. Pockets of low effective pressure lead to very fast flowing ice(maximum u = 2.4), whereas drained areas of the bed do not slow ice to the same extent(minimum u = 0.6). The large–scale effect of the instability is a net speed–up of ice: averagevelocity is 1.1037 when the model reaches flotation. Changes in ice thickness are very small, sothis corresponds to a net increase in the ice flux over the domain. Note that we are describinga transitional solution, we cannot describe the long term trend of the distribution of water or375. DISCUSSIONvelocity, without allowing for the continuation of the model in time. This requires a model thatallows the formation of deepening lakes under floating ice, as described in a suggested modelextension in Section 5.3.Thirdly, in accordance with the linearized theory, the pattern continues to move upstreamas these nonlinear effects manifest themselves. The wavespeed of the linearized pattern is−=(σ/k), only accurate while the pattern is small. In the linear analysis we found that thepattern velocity is small in comparison to the growth rate for a set of parameters that foster aninstability: this implies that the instability will grow to an appreciable size before it moves along distance relative to its own wavelength. It follows that unless the instability is very weak,we can expect the pattern to develop to a significant size before reaching the end of the idealisedfinite domain.5.3 Limitations and research outlookThe model used in this thesis contains some poorly constrained assumptions, primarily regardingwater flow at the bed. While we have assumed that relationships governing water storage andpermeability are dependent entirely on effective pressure, there are more sophisticated drainagemodels with additional degrees of freedom that future work in this area could consider. Inparticular, a more representative model drainage system should allow for the changing the natureof drainage systems from relatively inefficient distributed drainage networks to form efficientdrainage channels, such as that proposed by Hewitt [2013]. This model describes a distributedwater system with a permeability dκ/dt = f(κ, u,N), which can switch to a channelised drainagesystem with permeability dependent on channel size and channel roughness. Water storage hwis dependent on cavity size, derived from changes in ice velocity and water pressure. With thisextension to the drainage model we would expect to see an increase in drainage efficiency inregions of high hydraulic potential gradient. This could potentially exacerbate the instabilityby creating ‘R–channels’ that could slow ice downstream of areas of low effective pressure, orcould slow the feedback by smoothing perturbations in effective pressure more efficiently. Wenote that this drainage model still has flaws: in particular, the threshold for channelisation is385. DISCUSSIONnot well resolved, and the model fails to account for sediment erosion which could also lead tochannelisation.Model extension: Two dimensionsEven with the functional form of permeability and water storage in the model as posed, there aretwo obvious extensions to consider. Firstly, the model is posed in only one horizontal dimension.Extending the model to two horizontal dimensions would make it clear whether the instabilityis intrinsically one–dimensional, producing patterns at right angles to the ice and water flowdirection, or whether it breaks symmetry somehow. The three model equations (1)–(3) can begeneralized to two horizontal dimensions as follows:(4ηhux + 2ηhvy)x + (ηh (vx + uy))y − τbx = ρgh (h+ b)x , (49)(4ηhvy + 2ηhux)y + (ηh (uy + vx))x − τby = ρgh (h+ b)y , (50)ht + (hu)x + (hv)y = 0, (51)hw(N)t + (κ(N)χ)x + (κ(N)ψ)y = 0, (52)Where u(x, y, t) and v(x, y, t) are the x and y components of ice velocity respectively, h(x, y, t)is the ice thickness, N(x, y, t) is the effective pressure, b(x, y, t) is the bed surface, and χ =ρwgbx+Nx−ρghx and ψ = ρwgby+Ny−ρghy are the x and y components of hydraulic gradientrespectively.Model extension: Water–pocketsAs described in Section 4 the model breaks down when N = 0. A separate extension allowsthe one–dimensional model solution to continue to evolve by accounting for water storage inpockets between ice and water at the bed, where N = 0. This extension incorporates furtherconditions into the model that describe a region of deepening subglacial water column with afree boundary.Let hL be the depth of the subglacial water column, the distance between the ice bottomand the bed. Where N > 0, the ice is grounded, hL = 0 and the model as previously formulated395. DISCUSSION(equations (1)–(4)) holds. We expect that extended regions in which N = 0 can form, but donot allow N to become negative. In the regions where N = 0 we allow hL ≥ 0, describing theformation of a water–pocket between the ice and bed. Here the following amended model holds.4(ηhux)x + ρgh ((sin θ − (hx + hLx) = 0, (53)ht + (uh)x = 0, (54)ψ = 0 = ρwg sin θ + ρwghLx − ρghx. (55)N , hL, u, ux and N should be continuous across the free boundaries of regions where N = 0.In these regions, equations (53) and (54) solve for u and h while equation (55) solves for hLx interms of hx.5.4 Comparison of results with existing evidenceThe laterally symmetric ‘tiger stripes’ found by Sergienko and Hindmarsh [2013] and Sergienko et al.[2014] through inverse modelling are similar in structure to the patterns predicted in Section4.1 when viewed in one dimension: both describe striping or oscillating areas of high and lowbasal shear stress in the downstream direction. While this thesis does not attempt to describethe specific the structures found by Sergienko et al. [2014], it models a more general mechanismfor the spontaneous formation of spatial structure with periodic striping. The dominant wave-length of oscillations predicted is much larger (500 km) than those observed by Sergienko et al.[2014] (20 km), and unlike existing observations structures are predicted to migrate upstream.It is possible that the stripes identified by Sergienko et al. are formed by an instability similarto that described in this thesis, with other feedbacks contributing to restricting the wavelengthand wavespeed: for instance sediment transport as in Fowler and Chapwanya [2014], or ther-modynamic feedbacks as in Sergienko and Hulbe [2011]. It is also feasible that the predictedmigration is restricted or ‘pegged’ by fixed spatial changes in bed properties such as variabletopography, or geothermal heatflux, both commonly varying phenomena [e.g. Simkins et al.,2017, Whiteford and Graham, 1994].405. DISCUSSIONSergienko and Hulbe [2011] model feedbacks similar to those in this thesis as an explana-tion for the common close coexistence of subglacial lakes and topographical features or sharpchanges in basal conditions. They conclude that a sticky spot will induce an ice hump andwater pooling immediately downstream, a result that corresponds directly to the structureshown in figure 12, but found through different means. It follows that the instability wemodel may present a mechanism for the initiation, formation and growth of the sticky–spotfeedback that Sergienko and Hulbe [2011] describe. Additionally, Sergienko and Hulbe [2011]found that adding thermodynamic processes to their model further enhanced the instabilitythey modelled. Frictional melting of ice at sticky spots provided water that flowed to slipperyspots, enhancing lubrication (Fig. 2). With a spatial structure equivalent to that modelled bySergienko and Hulbe [2011], we would expect similar feedbacks from thermodynamic processes.In an empirical study based on multiple remote–sensing techniques Fricker et al. [2010] foundlakes co–located with sticky spots. They hypothesise the lakes formed as the result of thesame feedback with basal shear stress, hydraulic gradients and ice flux and as described bySergienko et al. [2007] and Sergienko and Hulbe [2011]. The ‘sticky spot’–lake pairs in thesmall region surveyed by Fricker et al. [2010] appear to be quasi–periodic, repeating four timeswith a wavelength approximately 20km long (Fig. 3). These structures correspond in shape tothe structures we predict, shown in Figure 12: an ice hump piles up due to sticky spot, justupstream of a low in ice thickness with an associated subglacial lake, followed downstream bymore average flowing ice, with these three distinct regions repeating periodically.While this thesis predicts the inception of the formation of a subglacial lake, it does notmodel the growth of a deepening lake. Our results suggest that feedbacks could spontaneouslydistribute water into patterned structures, forming of moving regions of low effective pressureor lakes, on a flat bed. This is in contrast to most of the 379 identified lakes which are large per-manent features formed by topographical basins [Wright and Siegert, 2012]. While our theorydoes not account for the typically large, individual lakes that are generally observed in practice,it does not attempt to: we have deliberately not introduced basal topography but sought out apattern–forming instability. Additionally, lakes of the kind we propose therefore would not show415. DISCUSSIONup in altimetry observations looking for drainage events (rapid, localized surface lowering dueto emptying of a lake [e.g. Siegfried and Fricker, 2018]): the lakes we describe are unlikely tofill and drain through mechanisms described in jo¨kulhlaups models [e.g. Nye, 1976, Evatt et al.,2006], because our drainage model does not permit channelization.Due to the fact that the identification of less distinct subglacial lakes and water systems withradar sounding is complicated by variable radio–wave attenuation in the overlying ice, there arefew clear examples of the kind of moving lakes or shallow areas of water that our model predicts[Young et al., 2016]. However, some identified lakes show very little variation in bed topography,suggesting a large permanent topographic basin is not a prerequisite for lake formation. Whilethere have been no explicit discoveries of lakes formed from feedbacks in the ice sheet or ofmigrating lakes, there are signs that suggest a wide variety of subglacial lakes exist. Numerous“fuzzy” or indefinite lakes suspected to be saturated sediments or “swamps” have been found byCarter et al. [2007] and Young et al. [2016]. Similarly, Horgan et al. [2012] theorises that LakeWhillans could be formed and maintained by temporary sediment erosion and deposition basedon evidence that the lake has disconnected or transient active and inactive portions. Thesekinds of indeterminate lakes may be more representative of the flat bedded water accumulationpredicted in this thesis. Even if lakes like Whillans or those described in Fricker et al. [2007]are not formed by instabilities similar to those modelled in this thesis, is is feasible that suchinstabilities catalyse the growth of existing features to form these lakes.The well documented sediment patterns in paleo–ice sheet remains [e.g. Kinahan and Close,1872, Ha¨ttestrand, 1997] such as drumlins, bed ribs or mega scale glacial lineations are notexplicity related to the ice sheet–drainage instability we have shown. Attempts to model theformation of these features [e.g. Hindmarsh, 1998b,a, Fowler, 2000, Schoof, 2007, Dunlop et al.,2008, Fowler and Chapwanya, 2014] differ from the modelling in this thesis. We have attemptedto consider a mechanism though which patterns can emerge that is not dependent on modelsfor sediment transport. As a result our model cannot simulate the formation of bed–forms.However, it is conceivable that the instability we describe leads to the formation of longitudinallypatternend features like drumlins or bed ribs by first making a pattern in basal conditions,426. CONCLUSIONcreating water–filled cavities at the ice–bed interface spontaneously, and subsequently fillingthem in with sediments to create the bed–forms.With existing data–sets [e.g. Mouginot et al., 2017] it may be possible to corroborate findingsof patterns in this thesis by conducting a high resolution Fourier analysis of Antarctic surfacevelocities, similar to that conducted on basal data by Spagnolo et al. [2017], but with spatialpatterning of velocities in the downstream direction.6 ConclusionThis thesis describes potential instabilities in the coupling of an ice sheet and its subglacialdrainage. My analysis shows that instabilities that spontaneously form patterning in subglacialdrainage are physically feasible, describes the spatial structures produced by the instability, givesus insight into the mechanisms driving the instabilities, and suggests some of the conditions thatmay foster an instability.Starting with a coupled 1D ice sheet–drainage model, we used linear stability analysis tofind parameter regimes for which the trivial fixed point is unstable. Solving the full, nonlinearmodel numerically, we found three regimes: one in which the spatially uniform steady state isstable, one in which there is growth to form a travelling wave of finite amplitude, and one inwhich growth is unbounded until flotation occurs and the model breaks down. These regimesare differentiated by the role played by effective pressure in setting the permeability of thedrainage system: the more sensitive the drainage system is to effective pressure, the less likelyan instability is to occur. Where instability occurs, a higher sensitivity to effective pressureis more likely to lead to bounded growth and a travelling wave being formed before flotationoccurs.By rescaling the coupled 1D ice sheet–drainage model we reduced the model to one whichhas an unstable steady state solution. The steps made to construct the reduced model show usconditions necessary to force the instability: changes in ice thickness are small, permeability isweakly dependent on effective pressure, the bed gradient is very shallow, and water velocity is436. CONCLUSIONmuch faster than ice velocity.Under these conditions, perturbations in the model solution grow in amplitude, developinga periodic spatial structure with localized pockets of low effective pressure that correspond tofast flowing ice, just downstream from humps in ice thickness formed by a ‘sticky spot’. Thiscorresponds to previous empirical and modelled sticky–slippery spot findings, and suggests thatsuch structures could form spontaneously in the right conditions. The instability is driven bya growing ice hump that forces more water to localized areas, reducing basal friction, which inturn increases ice flux into the hump. Although drawing specific predictions from the modelis not useful due to the large number of free parameters and lack of empirical justification forparameter values, the fundamental mechanisms, forces and structure dominating the identifiedinstability are likely key parts of our understanding of feedbacks in the coupling between icesheets and subglacial drainage systems.446. CONCLUSIONFigure 13: Evolution of a solution with parameter regime as in Figure 10 except for a change in permeabilitylaw to κ(N) = e−0.16N , closer to the shutdown of the instability described in Section 3.3.456. CONCLUSIONFigure 14: Maximum and minimum N(x, tf ) plotted against variations in κN , where tf is the point at whichthe amplitude of N has stopped growing. Parameter regime as in Figure 11 but with = 1. For κN < −0.29,(u¯, h¯, N¯) is stable. For −0.29 < κN < −0.15 the fixed point (u¯, h¯, N¯) is unstable and the solution trends towardsa second stable fixed point, with an amplitude that grows smoothly as |κN | grows.466. CONCLUSIONFigure 15: Evolution of a solution with parameter regime as in Figure 10 except for a change in permeabilitylaw to κ(N) = e−0.26N , the solution evolves into an apparently stable travelling wave solution without N = 0being reached.Figure 16: In a model extension, regions where N = 0 have an evolvable lake of depth hL, shown by the areabetween and including the black dots. This region is governed by equations (53)–(55). The black dots representfree boundaries where jump conditions must hold.47BibliographyM. R. Allen, V. R. Barros, J. Broome, W. Cramer, R. Christ, J. A. Church, L. Clarke, Q. Dahe,P. Dasgupta, N. K. Dubash, et al. IPCC fifth assessment synthesis report-climate change 2014synthesis report. 2014.R. Alley. Water-pressure coupling of sliding and bed deformation: I. water system. Journal ofGlaciology, 35(119):108–118, 1989.R. Alley, D. Blankenship, C. R. Bentley, and S. Rooney. Till beneath ice stream B: 3. tilldeformation: evidence and implications. Journal of Geophysical Research: Solid Earth, 92(B9):8921–8929, 1987.R. Alley, S. Anandakrishnan, C. Bentley, and N. Lord. A water-piracy hypothesis for thestagnation of Ice Stream C, Antarctica. Annals of Glaciology, 20:187–194, 1994.R. B. Alley. Towards a hydrological model for computerized ice-sheet simulations. HydrologicalProcesses, 10(4):649–660, 1996.R. B. Alley. The West Antarctic ice sheet: behavior and environment, volume 77. AmericanGeophysical Union, 2001.R. B. Alley, D. D. Blankenship, C. R. Bentley, and S. Rooney. Deformation of till beneath IceStream B, West Antarctica. Nature, 322(6074):57–59, 1986.D. W. Ashmore and R. G. Bingham. Antarctic subglacial hydrology: current knowledge andfuture challenges. Antarctic Science, 26(6):758–773, 2014.48BIBLIOGRAPHYJ. L. Bamber, D. G. Vaughan, and I. Joughin. Widespread complex flow in the interior of theAntarctic ice sheet. Science, 287(5456):1248–1250, 2000.J. Bear. Dynamics of fluids in porous media. Courier Corporation, 2013.M. R. Bennett. Ice streams as the arteries of an ice sheet: their mechanics, stability andsignificance. Earth-Science Reviews, 61(3-4):309–339, 2003.R. Bindschadler. The importance of pressurized subglacial water in separation and sliding atthe glacier bed. Journal of Glaciology, 29(101):3–19, 1983.R. A. Bindschadler. Meet Antarctica, Jun 2014. URL https://lima.nasa.gov/antarctica/.H. Bjo¨rnsson. Subglacial lakes and jo¨kulhlaups in Iceland. Global and Planetary Change, 35(3-4):255–271, 2003.D. D. Blankenship, C. R. Bentley, S. Rooney, and R. B. Alley. Seismic measurements reveala saturated porous layer beneath an active Antarctic ice stream. Nature, 322(6074):54–57,1986.D. D. Blankenship, C. R. Bentley, S. Rooney, and R. B. Alley. Till beneath Ice Stream B: 1.properties derived from seismic travel times. Journal of Geophysical Research: Solid Earth,92(B9):8903–8911, 1987.G. Boulton and R. Hindmarsh. Sediment deformation beneath glaciers: rheology and geologicalconsequences. Journal of Geophysical Research: Solid Earth, 92(B9):9059–9082, 1987.G. S. Boulton and D. Dent. The nature and rates of post-depositional changes in recentlydeposited till from south-east iceland. Geografiska Annaler. Series A. Physical Geography,pages 121–134, 1974.W. Budd, P. Keage, and N. Blundy. Empirical studies of ice sliding. Journal of glaciology, 23(89):157–170, 1979.49BIBLIOGRAPHYS. P. Carter, D. D. Blankenship, M. E. Peters, D. A. Young, J. W. Holt, and D. L. Morse. Radar-based subglacial lake classification in Antarctica. Geochemistry, Geophysics, Geosystems, 8(3), 2007.S. P. Carter, D. D. Blankenship, D. A. Young, and J. W. Holt. Using radar-sounding datato identify the distribution and sources of subglacial water: application to Dome C, EastAntarctica. Journal of Glaciology, 55(194):1025–1040, 2009.D. Chandler, J. Wadham, G. Lis, T. Cowton, A. Sole, I. Bartholomew, J. Telling, P. Nienow,E. Bagshaw, D. Mair, et al. Evolution of the subglacial drainage system beneath the Greenlandice sheet revealed by tracers. Nature Geoscience, 6(3):195, 2013.J. Church, P. Clark, A. Cazenave, J. Gregory, S. Jevrejeva, A. Levermann, M. Merrifield,G. Milne, R. Nerem, P. Nunn, A. Payne, W. Pfeffer, D. Stammer, and A. Unnikrishnan. SeaLevel Change, book section 13, page 1137–1216. Cambridge University Press, Cambridge,United Kingdom and New York, NY, USA, 2013. ISBN ISBN 978-1-107-66182-0. doi: 10.1017/CBO9781107415324.026. URL www.climatechange2013.org.G. K. Clarke. Subglacial processes. Annu. Rev. Earth Planet. Sci., 33:247–276, 2005.R. M. DeConto and D. Pollard. Contribution of Antarctica to past and future sea-level rise.Nature, 531(7596):591, 2016.P. Dunlop, C. D. Clark, and R. C. Hindmarsh. Bed ribbing instability explanation: testing anumerical model of ribbed moraine formation arising from coupled flow of ice and subglacialsediment. Journal of Geophysical Research: Earth Surface, 113(F3), 2008.H. Engelhardt and B. Kamb. Basal hydraulic system of a West Antarctic ice stream: constraintsfrom borehole observations. Journal of Glaciology, 43(144):207–230, 1997.H. Engelhardt, N. Humphrey, B. Kamb, and M. Fahnestock. Physical conditions at the base ofa fast moving Antarctic ice stream. Science, 248(4951):57–60, 1990.50BIBLIOGRAPHYG. Evatt, A. Fowler, C. Clark, and N. Hulton. Subglacial floods beneath ice sheets. Philosoph-ical Transactions of the Royal Society of London A: Mathematical, Physical and EngineeringSciences, 364(1844):1769–1794, 2006.G. E. Flowers. Modelling water flow under glaciers and ice sheets. In Proc. R. Soc. A, volume471, page 20140907. The Royal Society, 2015.G. E. Flowers and G. K. Clarke. A multicomponent coupled model of glacier hydrology 1. theoryand synthetic examples. Journal of Geophysical Research: Solid Earth, 107(B11), 2002.J. Fourier. Me´moire sur la propagation de la chaleur dans les corps solides. Nouveau Bulletindes Sciences de la Socie´te´ Philomathique de Paris, 6:112–116, 1808.A. Fowler. A sliding law for glaciers of constant viscosity in the presence of subglacial cavitation.In Proceedings of the Royal Society of London A: Mathematical, Physical and EngineeringSciences, volume 407, pages 147–170. The Royal Society, 1986.A. Fowler. Sliding with cavity formation. Journal of Glaciology, 33(115):255–267, 1987.A. Fowler. An instability mechanism for drumlin formation. Geological Society, London, SpecialPublications, 176(1):307–319, 2000.A. Fowler and C. Johnson. Ice-sheet surging and ice-stream formation. Annals of Glaciology,23:68–73, 1996.A. C. Fowler and M. Chapwanya. An instability theory for the formation of ribbed moraine,drumlins and mega-scale glacial lineations. In Proceedings of the Royal Society of LondonA: Mathematical, Physical and Engineering Sciences, volume 470, page 20140185. The RoyalSociety, 2014.P. Fretwell, H. D. Pritchard, D. G. Vaughan, J. Bamber, N. Barrand, R. Bell, C. Bianchi,R. Bingham, D. Blankenship, G. Casassa, et al. Bedmap2: improved ice bed, surface andthickness datasets for Antarctica. The Cryosphere, 7(1), 2013.51BIBLIOGRAPHYH. A. Fricker, T. Scambos, R. Bindschadler, and L. Padman. An active subglacial water systemin West Antarctica mapped from space. Science, 315(5818):1544–1548, 2007.H. A. Fricker, T. Scambos, S. Carter, C. Davis, T. Haran, and I. Joughin. Synthesizing multipleremote-sensing techniques for subglacial hydrologic mapping: application to a lake systembeneath MacAyeal Ice Stream, West Antarctica. Journal of Glaciology, 56(196):187–199,2010.L. Gray, I. Joughin, S. Tulaczyk, V. B. Spikes, R. Bindschadler, and K. Jezek. Evidence forsubglacial water transport in the West Antarctic ice sheet through three-dimensional satelliteradar interferometry. Geophysical Research Letters, 32(3), 2005.R. Greve and H. Blatter. Dynamics of ice sheets and glaciers. Springer Science & BusinessMedia, 2009.G. H. Gudmundsson. Transmission of basal variability to a glacier surface. Journal of Geophys-ical Research: Solid Earth, 108(B5), 2003.C. Ha¨ttestrand. Ribbed moraines in Sweden—distribution pattern and palaeoglaciological im-plications. Sedimentary Geology, 111(1-4):41–56, 1997.I. Hewitt. Seasonal changes in ice sheet motion due to melt water lubrication. Earth andPlanetary Science Letters, 371:16–25, 2013.I. Hewitt, C. Schoof, and M. Werder. Flotation and free surface flow in a model for subglacialdrainage. part 2. channel flow. Journal of Fluid Mechanics, 702:157–187, 2012.I. J. Hewitt. Modelling distributed and channelized subglacial drainage: the spacing of channels.Journal of Glaciology, 57(202):302–314, 2011.D. Hillel. Fundamentals of soil physics. Academic press, 2013.R. C. Hindmarsh. Ice-stream surface texture, sticky spots, waves and breathers: the coupledflow of ice, till and water. Journal of Glaciology, 44(148):589–614, 1998a.52BIBLIOGRAPHYR. C. Hindmarsh. The stability of a viscous till sheet coupled with ice flow, considered atwavelengths less than the ice thickness. Journal of Glaciology, 44(147):285–292, 1998b.T. Hodson, R. Powell, S. Brachfeld, S. Tulaczyk, R. Scherer, et al. Physical processes in Sub-glacial Lake Whillans, West Antarctica: inferences from sediment cores. Earth and PlanetaryScience Letters, 444:56–63, 2016.H. J. Horgan, S. Anandakrishnan, R. W. Jacobel, K. Christianson, R. B. Alley, D. S. Heeszel,S. Picotti, and J. I. Walter. Subglacial lake whillans—seismic observations of a shallow activereservoir beneath a West Antarctic ice stream. Earth and Planetary Science Letters, 331:201–209, 2012.A. Iken and R. A. Bindschadler. Combined measurements of subglacial water pressure andsurface velocity of Findelengletscher, Switzerland: conclusions about drainage system andsliding mechanism. Journal of Glaciology, 32(110):101–119, 1986.N. R. Iverson, T. S. Hooyer, and R. W. Baker. Ring-shear studies of till deformation: Coulomb-plastic behavior and distributed strain in glacier beds. Journal of Glaciology, 44(148):634–642,1998.B. Kamb. Rheological nonlinearity and flow instability in the deforming bed mechanism of icestream motion. Journal of Geophysical Research: Solid Earth, 96(B10):16585–16595, 1991.B. Kamb. Basal zone of the West Antarctic ice streams and its role in lubrication of their rapidmotion. The West Antarctic ice sheet: behavior and environment, pages 157–199, 2001.M. A. Kessler and R. S. Anderson. Testing a numerical glacial hydrological model using springspeed-up events and outburst floods. Geophysical Research Letters, 31(18), 2004.G. H. Kinahan and M. H. Close. The general glaciation of Iar-Connaught and its neighborhood,in the counties of Galway and Mayo. Hodges, Foster, 1872.T. Kyrke-Smith, R. Katz, and A. Fowler. Subglacial hydrology and the formation of ice streams.Proc. R. Soc. A, 470(2161):20130494, 2014.53BIBLIOGRAPHYK. Langley, J. Kohler, K. Matsuoka, A. Sinisalo, T. Scambos, T. Neumann, A. Muto, J.-G.Winther, and M. Albert. Recovery lakes, East Antarctica: Radar assessment of sub-glacialwater extent. Geophysical Research Letters, 38(5), 2011.A. R. Lewis, D. R. Marchant, D. E. Kowalewski, S. L. Baldwin, and L. E. Webb. The ageand origin of the labyrinth, western Dry Valleys, Antarctica: Evidence for extensive middlemiocene subglacial floods and freshwater discharge to the southern ocean. Geology, 34(7):513–516, 2006.L. Lliboutry. General theory of subglacial cavitation and sliding of temperate glaciers. Journalof Glaciology, 7(49):21–58, 1968.D. R. MacAyeal. Large-scale ice flow over a viscous basal sediment: Theory and application toice stream B, Antarctica. Journal of Geophysical Research: Solid Earth, 94(B4):4071–4087,1989.S. MacNamara and G. Strang. Operator splitting. In Splitting Methods in Communication,Imaging, Science, and Engineering, pages 95–114. Springer, 2016.D. McKinlay, A. McGown, A. Radwan, and D. Hossain. Representative sampling and testingin fissured lodgement tills. In Proceedings of the Midland Soil Mechanics and FoundationEngineering Society Symposium on The Engineering Behaviour of Glacial Materials, pages143–55, 1975.M. Morlighem, H. Seroussi, E. Larour, and E. Rignot. Inversion of basal friction in Antarcticausing exact and incomplete adjoints of a higher-order model. Journal of Geophysical Research:Earth Surface, 118(3):1746–1753, 2013.J. Mouginot, E. Rignot, B. Scheuchl, and R. Millan. Comprehensive annual ice sheet velocitymapping using Landsat-8, Sentinel-1, and RADARSAT-2 data. Remote Sensing, 9(4):364,2017.J.-O. Na¨slund, J. Fastook, and P. Holmlund. Numerical modelling of the ice sheet in western54BIBLIOGRAPHYdronning maud land, East Antarctica: impacts of present, past and future climates. Journalof Glaciology, 46(152):54–66, 2000.F. S. Ng. Mathematical modelling of subglacial drainage and erosion. PhD thesis, University ofOxford, 1998.J. Nye. Water flow in glaciers: jo¨kulhlaups, tunnels and veins. Journal of Glaciology, 17(76):181–207, 1976.W. S. B. Paterson. The physics of glaciers. Elsevier, 2016.F. Pattyn. Investigating the stability of subglacial lakes with a full stokes ice-sheet model.Journal of Glaciology, 54(185):353–361, 2008.T. Perol and J. R. Rice. Control of the width of West Antarctic ice streams by internal meltingin the ice sheet near the margins. In AGU Fall Meeting Abstracts, 2011.W. Reilly and C. Whiteford. Gravity map of New Zealand, 1: 4,000,000. Bouguer anomalies.Wellington, New Zealand, Department of Scientific and Industrial Research, 1965.E. Rignot and R. H. Thomas. Mass balance of polar ice sheets. Science, 297(5586):1502–1506,2002.E. Rignot, J. Mouginot, and B. Scheuchl. Ice flow of the Antarctic ice sheet. Science, 333(6048):1427–1430, 2011.G. d. Q. Robin, C. Swithinbank, and B. Smith. Radio echo exploration of the Antarctic icesheet. International Association of Scientific Hydrology Publication, 86:97–115, 1970.K. Rose. Characteristics of ice flow in Marie Byrd Land, Antarctica. Journal of Glaciology, 24(90):63–75, 1979.H. Ro¨thlisberger. Water pressure in intra-and subglacial channels. Journal of Glaciology, 11(62):177–203, 1972.55BIBLIOGRAPHYC. Schoof. Variational methods for glacier flow over plastic till. Journal of Fluid Mechanics,555:299–320, 2006.C. Schoof. Pressure-dependent viscosity and interfacial instability in coupled ice–sediment flow.Journal of Fluid Mechanics, 570:227–252, 2007.C. Schoof. Ice-sheet acceleration driven by melt supply variability. Nature, 468(7325):803, 2010.C. Schoof and R. C. Hindmarsh. Thin-film flows with wall slip: an asymptotic analysis of higherorder glacier flow models. The Quarterly Journal of Mechanics & Applied Mathematics, 63(1):73–114, 2010.C. Schoof, I. J. Hewitt, and M. A. Werder. Flotation and free surface flow in a model forsubglacial drainage. part 1. distributed drainage. Journal of Fluid Mechanics, 702:126–156,2012.C. Schoof, C. Rada, N. Wilson, G. Flowers, and M. Haseloff. Oscillatory subglacial drainage inthe absence of surface melt. The Cryosphere, 8(3):959–976, 2014.D. M. Schroeder, D. D. Blankenship, and D. A. Young. Evidence for a water system transitionbeneath Thwaites Glacier, West Antarctica. Proceedings of the National Academy of Sciences,110(30):12225–12228, 2013.O. Sergienko, D. MacAyeal, and R. Bindschadler. Causes of sudden, short-term changes inice-stream surface elevation. Geophysical Research Letters, 34(22), 2007.O. Sergienko, T. T. Creyts, and R. Hindmarsh. Similarity of organized patterns in driving andbasal stresses of Antarctic and Greenland ice sheets beneath extensive areas of basal sliding.Geophysical Research Letters, 41(11):3925–3932, 2014.O. V. Sergienko and R. C. Hindmarsh. Regular patterns in frictional resistance of ice-streambeds seen by surface data inversion. Science, 342(6162):1086–1089, 2013.O. V. Sergienko and C. L. Hulbe. ‘Sticky spots’ and subglacial lakes under ice streams of theSiple Coast, Antarctica. Annals of Glaciology, 52(58):18–22, 2011.56BIBLIOGRAPHYJ. Shaw and R. Gilbert. Evidence for large-scale subglacial meltwater flood events in southernOntario and northern New York State. Geology, 18(12):1169–1172, 1990.R. Shreve. Movement of water in glaciers. Journal of Glaciology, 11(62):205–214, 1972.M. J. Siegert, S. Carter, I. Tabacco, S. Popov, and D. D. Blankenship. A revised inventory ofAntarctic subglacial lakes. Antarctic Science, 17(3):453–460, 2005.M. R. Siegfried and H. A. Fricker. Thirteen years of subglacial lake activity in Antarctica frommulti-mission satellite altimetry. Annals of Glaciology, pages 1–14, 2018.L. M. Simkins, J. B. Anderson, S. L. Greenwood, H. M. Gonnermann, L. O. Prothro, A. R. W.Halberstadt, L. A. Stearns, D. Pollard, and R. M. DeConto. Anatomy of a meltwater drainagesystem beneath the ancestral East Antarctic ice sheet. Nature Geoscience, 10(9):691, 2017.M. Spagnolo, T. C. Bartholomaus, C. D. Clark, C. R. Stokes, N. Atkinson, J. A. Dowdeswell,J. C. Ely, A. G. Graham, K. A. Hogan, E. C. King, et al. The periodic topography of icestream beds: insights from the fourier spectra of mega-scale glacial lineations. Journal ofgeophysical research: earth surface, 2017.L. A. Stearns, B. E. Smith, and G. S. Hamilton. Increased flow speed on a large East Antarcticoutlet glacier caused by subglacial floods. Nature Geoscience, 1(12):827–831, 2008.S. Tulaczyk, W. B. Kamb, and H. F. Engelhardt. Basal mechanics of Ice Stream B, WestAntarctica: 2. Undrained plastic bed model. Journal of Geophysical Research: Solid Earth,105(B1):483–494, 2000.J. S. Walder and A. Fowler. Channelized subglacial drainage over a deformable bed. Journalof Glaciology, 40(134):3–15, 1994.J. Weertman. On the sliding of glaciers. Journal of glaciology, 3(21):33–38, 1957.M. A. Werder, I. J. Hewitt, C. G. Schoof, and G. E. Flowers. Modeling channelized anddistributed subglacial drainage in two dimensions. Journal of Geophysical Research: EarthSurface, 118(4):2140–2158, 2013.57BIBLIOGRAPHYP. Whiteford and D. Graham. Conductive heat flow through the sediments in Lake Rotomahana,New Zealand. Geothermics, 23(5-6):527–538, 1994.D. J. Wingham, M. J. Siegert, A. Shepherd, and A. S. Muir. Rapid discharge connects Antarcticsubglacial lakes. Nature, 440(7087):1033–1036, 2006.A. Wright and M. Siegert. A fourth inventory of Antarctic subglacial lakes. Antarctic Science,24(06):659–664, 2012.D. A. Young, D. Schroeder, D. Blankenship, S. D. Kempf, and E. Quartini. The distribution ofbasal water between Antarctic subglacial lakes from radar sounding. Phil. Trans. R. Soc. A,374(2059):20140297, 2016.58
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Pattern-forming instabilities in the coupling of ice...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Pattern-forming instabilities in the coupling of ice sheets and subglacial drainage systems Whiteford, Arran 2018
pdf
Page Metadata
Item Metadata
Title | Pattern-forming instabilities in the coupling of ice sheets and subglacial drainage systems |
Creator |
Whiteford, Arran |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | Sharp spatial changes discovered in the basal conditions of an ice sheet do not always have an obvious source. By modelling instabilities in the coupling of an ice sheet and subglacial drainage system, we describe physical feedback mechanisms that force the formation of sharp spatial structures in basal conditions and ice flow. This model predicts the spontaneous formation of periodic subglacial `sticky spot'-lake pairs, that correspond in shape to previous empirical and modelled descriptions of similar structures. The instability that forms this structure is driven by a feedback whereby periodic humps in ice thickness redirect subglacial water to slippery spots that lie immediately downstream of the ice humps: the slippery regions increase ice flux into the ice humps, making them grow. Scaling a one-dimensional model ice sheet coupled to a basal drainage system, we find conditions for the instability with linear stability analysis. Solutions in the full nonlinear model are simulated numerically, using operator splitting and finite difference methods. The instability requires a bed permeability weakly dependent on water pressure changes, negligible bed slopes, and a water velocity much greater than ice velocity. The `sticky spot'-lake pairs are predicted to form with periodic spacing and migrate upstream. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-05-14 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-ShareAlike 4.0 International |
DOI | 10.14288/1.0366280 |
URI | http://hdl.handle.net/2429/65843 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-sa/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_september_whiteford_arran.pdf [ 1.38MB ]
- Metadata
- JSON: 24-1.0366280.json
- JSON-LD: 24-1.0366280-ld.json
- RDF/XML (Pretty): 24-1.0366280-rdf.xml
- RDF/JSON: 24-1.0366280-rdf.json
- Turtle: 24-1.0366280-turtle.txt
- N-Triples: 24-1.0366280-rdf-ntriples.txt
- Original Record: 24-1.0366280-source.json
- Full Text
- 24-1.0366280-fulltext.txt
- Citation
- 24-1.0366280.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
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"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0366280/manifest