UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modelling the migration of ice stream margins Haseloff, Marianne 2015

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

Item Metadata


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

Full Text

Modelling the migration ofice stream marginsbyMarianne HaseloffDipl. Phys., Humboldt-University Berlin, 2009A Thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of PhilosophyinThe Faculty of Graduate and Postdoctoral Studies(Geophysics)The University of British Columbia(Vancouver)July 2015c©Marianne Haseloff, 2015AbstractThe Siple Coast ice streams are long, narrow bands of ice within the Antarctic ice sheet. Theymove significantly faster than the surrounding ice ridges, and therefore discharge significantlymore ice. Observations suggest that their fast flow is due to sliding along a water-saturatedbed, while the bed of the neighbouring ridges generally appears to be frozen. The ice streamvelocities and widths vary on decadal to centennial time scales, and these variations includethe migration of the ice stream margins, where the fast flow slows down to the speed of thesurrounding ice.In this thesis I show that conventional thin film models, which are used to calculate theevolution of ice sheets on continental scales, are only able to reproduce the inwards migrationof ice stream margins and the subsequent shutdown of an ice stream. These processes are theresult of an insufficient heat dissipation and freezing at the bed. Conversely, I find that thewidening of ice streams into regions where the bed is frozen can only be modelled by takingsmall-scale heat transfer processes in the ice stream margin into account. Previous research hasshown that ice stream widening results from an interplay of heating through lateral shearingin the ice stream margin and inflow of cold ice from the adjacent ridges. However, the relativeimportance of the different effects on the migration speed has not yet been quantified. Toaccount for these processes, I derive a new boundary layer model for ice stream margins. Thenumerical solution of this model provides the margin migration speed as a function of large-scaleice stream properties such as ice stream width, ice thickness, and geothermal heat flux.The influence of different basal boundary conditions and temperate ice properties on themargin migration velocity is also investigated. To derive a parameterization of ice streamwidening that can be used in continental-scale models, I consider asymptotic solutions withhigh heat production rates and high advection velocities, a limit that likely applies in real icestream margins.iiPrefaceThis thesis consists of two complementary studies that form the basis of three individual pub-lications.A version of chapter 2 is ready for submission to a geophysical research journal. The co-authors are Marianne Haseloff (first author) and Christian Schoof. I was responsible for con-ducting the research, interpreting the results, and composing the manuscript. Christian Schoofprovided continuous support during all of these stages.A version of chapter 3 with excerpts of sections 4.2–4.4 has been submitted to a peer-reviewedresearch journal. The co-authors are Marianne Haseloff (first author), Christian Schoof andOlivier Gagliardini. I was responsible for conducting the research, interpreting the results, andcomposing the manuscript. Christian Schoof provided continuous support during all of thesestages. Olivier Gagliardini provided the first setup of the numerical model in Elmer/Ice, advisedon the numerical solution, and provided editorial comments on the manuscript.A version of sections 4.5–4.7 will be submitted to a geophysical research journal. The co-authors are Marianne Haseloff (first author), Christian Schoof and Olivier Gagliardini. I wasresponsible for conducting the research, interpreting the results, and composing the manuscript.Christian Schoof provided continuous support during all of these stages. Olivier Gagliardiniprovided the first setup of the numerical model in Elmer/Ice and advised on the numericalsolution.Section C.4 was kindly provided by Christian Schoof, who performed the calculation andwrote the section; it is included for completeness.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 General properties of the Siple Coast ice streams . . . . . . . . . . . . . 11.1.2 Temporal variations of the Siple Coast ice streams . . . . . . . . . . . . . 41.2 State of knowledge and open research questions . . . . . . . . . . . . . . . . . . 51.3 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Ice stream flow and basal freeze-on . . . . . . . . . . . . . . . . . . . . . . . . . 92.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2.1 Ice flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2.2 Sediment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.3 Energy balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3 Ice stream evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.1 Outwards migration and evolution to steady states . . . . . . . . . . . . 182.3.2 Spurious margin migration . . . . . . . . . . . . . . . . . . . . . . . . . . 20iv2.3.3 Inwards migration and ice stream shutdown . . . . . . . . . . . . . . . . 222.4 Ice stream steady states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.4.1 General steady state behaviour and stability . . . . . . . . . . . . . . . . 242.4.2 The effect of lateral shear heating . . . . . . . . . . . . . . . . . . . . . . 262.4.3 The effect of subglacial drainage . . . . . . . . . . . . . . . . . . . . . . . 282.4.4 The effect of the basal friction law . . . . . . . . . . . . . . . . . . . . . . 302.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313 Boundary layer model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.1 Description of the geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.2 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.3 Ice stream . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.3.1 Leading order ice stream model . . . . . . . . . . . . . . . . . . . . . . . 443.3.2 Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.4 Ice ridge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.4.1 Leading order ice ridge model . . . . . . . . . . . . . . . . . . . . . . . . 493.4.2 Depth integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.5 Boundary layer scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.5.1 Boundary layer model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.5.2 Series expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.6 Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.6.1 Matching the stream with the boundary layer . . . . . . . . . . . . . . . 593.6.2 Matching the ridge with the boundary layer . . . . . . . . . . . . . . . . 613.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634 Solution of the boundary layer model . . . . . . . . . . . . . . . . . . . . . . . . 674.1 The finite-element solver Elmer/Ice . . . . . . . . . . . . . . . . . . . . . . . . . 704.2 Ice flow and heat production . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.3 Temperature field and migration velocity . . . . . . . . . . . . . . . . . . . . . . 734.4 Migration velocity as a function of forcing parameters . . . . . . . . . . . . . . . 764.5 Asymptotic solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78v4.5.1 The velocity field close to the transition point . . . . . . . . . . . . . . . 794.5.2 Temperature solutions and migration velocity for α 1 . . . . . . . . . 814.5.3 Temperature solutions and migration velocity for α 1 and Pe 1 . . . 834.6 Sub-temperate sliding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.7 Temperate ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 904.7.1 Local solution and numerical solution . . . . . . . . . . . . . . . . . . . . 934.7.2 Heat production at large distances from the origin . . . . . . . . . . . . . 984.7.3 Asymptotic behaviour for α 1 . . . . . . . . . . . . . . . . . . . . . . 1004.7.4 Asymptotic behaviour for α 1 and Pe 1 . . . . . . . . . . . . . . . . 1034.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1065 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1125.1 Outwards migration of ice stream margins . . . . . . . . . . . . . . . . . . . . . 1125.2 Inwards migration of ice stream margins . . . . . . . . . . . . . . . . . . . . . . 1155.3 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118A Appendix to chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131A.1 Integrated energy balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131A.2 The upper bound on ice stream velocity . . . . . . . . . . . . . . . . . . . . . . 132A.3 Ice stream flow on a fully permeable bed . . . . . . . . . . . . . . . . . . . . . . 132A.4 Solutions with small lateral shear stresses . . . . . . . . . . . . . . . . . . . . . 133A.5 Solutions with large lateral shear stresses and small drainage . . . . . . . . . . . 136A.6 The role of subglacial drainage . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139A.6.1 Ice streams without subglacial drainage and σ > 0 . . . . . . . . . . . . . 139A.6.2 Bifurcation diagrams for subglacial drainage parameter K . . . . . . . . 141A.6.3 The effect of the drainage parameter a . . . . . . . . . . . . . . . . . . . 143B Appendix to chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144viC Appendix to chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149C.1 Local behaviour for the transverse velocity . . . . . . . . . . . . . . . . . . . . . 149C.2 Boundary layer model for α 1 and Pe 1 . . . . . . . . . . . . . . . . . . . 151C.2.1 Temperature solutions and migration velocity for α 1 . . . . . . . . . 153C.2.2 Temperature solutions for α 1 and Pe 1 . . . . . . . . . . . . . . . 154C.3 Change of model for sub-temperate sliding . . . . . . . . . . . . . . . . . . . . . 155C.4 Local solution for free cold-temperate boundary (provided by Christian Schoof) 156viiList of Tables2.1 Table of variables. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Table of parameters, their values (where applicable) and their units. . . . . . . . 153.1 Parameter values and characteristic scales. . . . . . . . . . . . . . . . . . . . . . 373.2 Aspect ratios of the boundary layer model. . . . . . . . . . . . . . . . . . . . . . 423.3 Dimensionless groups. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45viiiList of Figures1.1 Velocity field of Antarctica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Velocity field of the Siple Coast region in West Antarctica. . . . . . . . . . . . . 31.3 Sketch of the ice stream geometry considered in this thesis. . . . . . . . . . . . . 72.1 Sketch of the ice stream/sediment model. . . . . . . . . . . . . . . . . . . . . . . 132.2 Temporal evolution of ice stream flow that is initially localized by a trough ineffective pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.3 Example of grid-size dependent margin migration rate. . . . . . . . . . . . . . . 202.4 Shutdown of an ice stream. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.5 Steady state solutions and stability. . . . . . . . . . . . . . . . . . . . . . . . . . 252.6 Effect of shear heating on ice stream steady states. . . . . . . . . . . . . . . . . 272.7 Dependence of steady state solutions on permeability K. . . . . . . . . . . . . . 292.8 Influence of basal friction law on the effective pressure and steady states. . . . . 303.1 Setup of boundary layer model and definition of geometric scales. . . . . . . . . 364.1 Example of a mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.2 Velocities and heat production in the boundary layer. . . . . . . . . . . . . . . . 724.3 Numerical scheme to determine vm. . . . . . . . . . . . . . . . . . . . . . . . . . 754.4 Migration velocity as function of forcing parameters. . . . . . . . . . . . . . . . 774.5 The effect of shear heating and advection on the temperature field. . . . . . . . 784.6 Comparison of numerical near-field velocity solutions with benchmark solutions. 804.7 Boundary layer structure for α 1 and Pe 1. . . . . . . . . . . . . . . . . . 824.8 Asymptotic behaviour of Vm = vm/α for α 1 and Pe 1. . . . . . . . . . . . 854.9 Relation between α, Pe and vm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.10 Influence of temperature field in the ridge on the migration rate. . . . . . . . . . 874.11 The effect of sub-temperate sliding on heat production rates and temperature. . 884.12 Migration speed against Tc for different values of Pe. . . . . . . . . . . . . . . . 89ix4.13 Geometry for the boundary conditions for temperate ice. . . . . . . . . . . . . . 924.14 Geometry of the local solution for temperate ice. . . . . . . . . . . . . . . . . . . 944.15 Sample solution with additional boundary condition for temperate ice. . . . . . 964.16 Numerical mesh and temperature along the bed for the solution shown in figure4.15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 974.17 Margin migration rates vm with additional boundary condition (4.36). . . . . . . 984.18 Sample profiles of temperate ice solutions for α and Pe of O(1). . . . . . . . . . 994.19 Far field heat production rates. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1004.20 Boundary layer geometry for temperate ice with small advection. . . . . . . . . . 1014.21 Asymptotic results for temperate ice with small advection. . . . . . . . . . . . . 1034.22 Boundary layer geometry for temperate ice with large advection. . . . . . . . . 1054.23 Dimensional migration velocities against ice stream width and ice stream thickness.1084.24 Temperature field in the boundary layer and the adjacent ice stream. . . . . . . 110A.1 Locations of asymptotic solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . 134A.2 Sample plots of semianalytical solutions in the limit of low and high lateral shear. 135A.3 Solutions for K = 0 and σ > 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . 141A.4 Unphysical solutions at small permeabilities and high velocities. . . . . . . . . . 142A.5 Unphysical solutions at small permeabilities and low velocities. . . . . . . . . . . 142A.6 Influence of drainage parameter a on ice stream steady states. . . . . . . . . . . 143C.1 Comparison of numerical near-field velocity solutions with benchmark solutions. 150xAcknowledgementsNo equation would be scaled and no sentence written here without the persistent and compre-hensive support of my adviser Christian Schoof, who not only conceived the project that formsthe basis of this thesis, but who also patiently and enthusiastically taught me how to executeit. Working with him was incredibly inspiring and enjoyable, and I have always felt extremelylucky about having such a knowledgeable and involved adviser. His support has by far exceededanything that I could have hoped for, and I cannot give enough credit and thanks to him.I was also very lucky with an all-star PhD committee consisting of Mark Jellinek, NeilBalmforth, and Gwenn Flowers, that has supported me throughout the entire course of myPhD.The enthusiasm and knowledge Mark Jellinek radiates for what seems to be every possibleresearch field has always amazed me. He has opened my eyes to many exciting research areasoutside of glaciology, taught me how to comfortably discuss the things I know and don’t know,that learning begins with admitting that there is something still to be learned, and that abalance between work and life is possible.The second chapter of this thesis would probably not exist without what seemed to be arelatively harmless question asked by Neil Balmforth in my second committee meeting. Thisdiscussion and the work resulting from it has helped me to understand much better why I wasactually doing what I was doing, and I am very grateful for the advice and helpful feedback, aswell as an introduction to perturbation methods.My committee was not complete before Gwenn Flowers joined it, bringing her glaciologicalexpertise, encouragement, constructive criticism, and a great sense of humour. She provided allof these things generously whenever they were needed.The fourth chapter of this thesis owes a great deal of its existence to Olivier Gagliardini,who introduced me to Elmer/Ice, set up the original geometry of the boundary layer model,and advised me whenever I had trouble with the numerical solution.Without Garry Clarke’s support I would have never made it to UBC, and I am also verygrateful for his insightful advice.xiAlthough not a member of my committee, Valentina Radic has always given me helpfuladvice on how to master the little pitfalls of grad school and life, and without her friendship Iwould have been a little lost sometimes (and would have slept in a sleeping bag on the floor forprobably a very long time).Thanks also to Jenn, Jean-Francois, Christophe, Kathrin, Ian, Mauro & Ursula, Aram, andthe current and former members of the UBC glaciology group for making my time here soenjoyable!Laurens: Thank you for everything – patience, edits, encouragement, and all the other thingsthat bring so much happiness into my life.Finally, many thanks to my parents who have always supported me in every respect, andwho remain close at any distance.Financial support was provided through NSERC (grant number 357193-13) and by theDepartment of Earth, Ocean and Atmospheric Sciences in form of a Four-Year-Fellowship, anEgil H. Lorntzsen scholarship, and a W. H. Matthews scholarship. Most of the numericalsolutions in chapter 4 were computed on ComputeCanada resources.xiiFu¨r Reiner & Renate.xiii1 | IntroductionThe Antarctic continent is covered by an ice sheet that is thousands of kilometres wide and about3 km thick in the interior. It gains mass through the accumulation of snow, which densifies andturns into ice. Driven by gravity, the ice moves as a viscous fluid towards the margins of the icesheet, where it is predominantly lost through the calving of icebergs into the ocean. Figure 1.1shows a velocity map of Antarctica derived from satellite measurements, with fast regions inpurple/red (Rignot et al., 2011).Ice streams are narrow bands of very fast flowing ice embedded in slowly moving ice. Thebulk of an ice stream flows at surface velocities of hundreds of meters per year, while theneighbouring ice moves at surface velocities of only a few meters per year. Despite the fact thatthey make up less than 10% of the Antarctic ice sheet, ice streams are the main conduits forice transport from the interior of the ice sheet to the coast (Bamber et al., 2000).Several of the streaming features in figure 1.1 are fast-flowing due to topographic or geologicalcontrols. However, this is true only to a limited extent for the part of West Antarctica enlargedin figure 1.2. These are the five Siple Coast ice streams which provide the motivation formy research project. As summarised below, they are characterised by very similar physicalproperties.1.1 Observations1.1.1 General properties of the Siple Coast ice streamsAll the Siple Coast ice streams have similar dimensions and speeds. They are approximately 1km thick, 50 km wide, and several hundred kilometres long (see figure 1.2). They are separatedby slowly moving inter-ice-stream regions called ice ridges. Surface elevation maps show thatthese ridges have higher surface elevations than the ice streams themselves (Shabtaie and Bent-ley, 1988; Retzlaff et al., 1993). At first glance, the slow (light-blue) region marked as Kamb IceStream (C) in figure 1.2 appears to be an exception to that pattern, but apart from its velocity1Figure 1.1: Velocity field of Antarctica. From Rignot et al. (2011). Reprinted with permissionfrom AAAS.profile it shares the geometric properties of the adjacent ice streams: it is an elongated regionwith a relatively low-angled surface, bordered by steeper inwards–sloping ridges.The Siple Coast ice streams have high basal velocities, and exhibit little vertical variationin their velocity profiles (Engelhardt and Kamb, 1998). This fast flow seems to be enabledby water at the ice stream bed, which allows rapid basal slip (Blankenship et al., 1986, 1987;Engelhardt and Kamb, 1997; Gray et al., 2005; Kamb, 2001; Fricker et al., 2007; Winberry et al.,2009). Conversely, the bed outside of the ice stream generally appears to be frozen (Bentleyet al., 1998; Catania et al., 2003), which suggests that the transition from fast to slow velocitiescan be linked to a change in thermal conditions at the ice base.The transition from fast ice stream flow to the slow flow of the ice ridges occurs in the icestream margin, which is a narrow, heavily crevassed zone at the sides of an ice stream (Bind-schadler and Scambos, 1991; Echelmeyer et al., 1994; Bindschadler et al., 1996; Jacobson and2Figure 1.2: Velocity field of the Siple Coast region in West Antarctica. Ice streams from top tobottom (South to North): Mercer (A), Whillans (B), Kamb (C), Bindschadler (D), MacAyeal(E). Figure from Le Brocq et al. (2009).Raymond, 1998; Joughin et al., 2002). The margin is only a few kilometers wide, and obser-vations suggest high localized strain rates and shear stresses in this area (Bindschadler et al.,1987; Whillans et al., 1993; Echelmeyer et al., 1994; Jackson and Kamb, 1997). Temperaturemeasurements in the margin of Whillans ice stream have shown that this leads to significantshear heating (Harrison et al., 1998).The location of the ice stream margins appears to be only weakly controlled by geological andtopographic conditions. Soft sediments have been found underneath the onset of some ice streamtributaries, while some tributaries appear to be topographically constrained (Anandakrishnanet al., 1998; Joughin et al., 2002). However, radar profiles have shown that the fully evolveddownstream regions of ice streams are not confined to troughs in the bed topography and thatsubglacial soft sediments extend beyond the ice stream margins (Bindschadler et al., 1996;Nereson and Raymond, 2007). Moreover, the direct and indirect observations summarizedbelow indicate that ice stream extent as well as ice stream velocities can change over time. Thissuggests that the Siple Coast ice streams are only weakly controlled by geological or topographicconditions, and that dynamic processes play an important role in their formation and evolution.31.1.2 Temporal variations of the Siple Coast ice streamsObservations suggest temporal variability in the flow of the Siple Coast ice streams. Recentvelocity changes have been observed directly through survey of stakes on the ice surface andthrough satellite measurements. These measurements show that ice stream margins can migrate(Whillans et al., 1987; Stephenson and Bindschadler, 1988; Bindschadler and Vornberger, 1998;Hulbe and Whillans, 1997; Echelmeyer and Harrison, 1999; Joughin et al., 2002). An exampleof such temporal variability is the recent widening of the tributaries of Whillans ice stream.From its geometrical properties, Kamb ice stream (figure 1.2) is widely interpreted to be an‘inactive’ ice stream: the depth of burial of its marginal crevasses suggests that the ice streamwas last active approximately 160 years ago (Retzlaff and Bentley, 1993). Ice-penetrating radarimaging has suggested that its margins migrated inwards as the ice stream was shutting down.Moreover, changes in flow patterns during the past 1000 years have been deduced fromtracking patterns that characterize the ice stream margins. The marginal crevasse patterns ofthe Siple Coast ice streams are advected downstream into the Ross ice shelf, and the distortionof these features demonstrates that discharge of the ice streams from which they originatedmust have changed over the last 1000 years. Radar measurements of distorted internal icelayers have been used to identify another former ice stream (the ‘Siple ice stream’), that wasflowing at an angle to the directions of the current ice streams approximately 420 years ago (seee.g. Bindschadler et al., 2000; Fahnestock et al., 2000; Conway et al., 2002; Stearns et al., 2005;Catania et al., 2006; Hulbe and Fahnestock, 2007; Catania et al., 2012).On much longer timescales, marine sediment records imply that episodes of ice streaming inthe Laurentide ice sheet have led to the discharge of large amounts of ice into the North Atlanticduring the last glacial period between 60 and 10 thousand years ago (these are called Heinrichevents, Heinrich, 1988; Bond et al., 1992; Broecker et al., 1992). It has been hypothesized thatthese ice discharges were sufficiently large to impact the ocean circulation system through theaddition of cold, fresh water, and to have altered atmospheric circulation through changes inice sheet elevation (Stokes and Clark, 2001).41.2 State of knowledge and open research questionsBased on the observations described above, several key research areas can be identified:1. Including ice streams into continental-scale simulationsIce streams discharge the majority of the Antarctic ice, and therefore affect the evolutionof the ice sheet. The last report of the International Panel on Climate Change (IPCC)states that the contribution of the ice sheets remains a significant source of uncertaintyin predicting future sea level rise (Church et al., 2013). Quantifying such effects requiresincluding ice stream dynamics in continental-scale simulations, which model the evolutionof ice sheets over hundreds to thousands of years (e.g. Bueler and Brown, 2009). At presentdifferent ice sheet models predict very different ice stream behaviour under prescribedforcing (Calov et al., 2010).2. Spatial patterningAll Siple Coast ice streams have similar physical properties and they form a pattern sug-gestive of a fingering instability. The formation of this kind of pattern may be the resultof a positive feedback that amplifies small perturbations in flow speed in an otherwise uni-form flow field. Models implementing feedbacks of this type, based on thermomechanicalcoupling, have been able to reproduce ice streams as flow instabilities (more details aboutthe relevant physics are provided in the introduction of the next chapter, section 2.1; seealso Fowler and Johnson, 1996; Payne and Dongelmans, 1997).However, this feedback alone is not sufficient to account for the apparent spatial wave-length selection in the width of the Siple Coast ice streams. Recent work has explaineda stabilization (i.e. damping) of short wavelength perturbations through the prevailingstress conditions in ice streams (Sayag and Tziperman, 2008; Hindmarsh, 2009). It hasbeen hypothesized that dynamic thinning due to increased flow speed and the limitedavailability of mass (generally, a wider ice stream would transport more mass) could sup-press the growth of large wavelengths in ice stream width (Fowler and Johnson, 1996;Hindmarsh, 2009; Kyrke-Smith et al., 2014).5However, the details of the positive feedbacks at work and some aspects of the non-linearevolution of these features are still poorly understood.3. Temporal variabilityThe origin of the temporal variations outlined above might be closely related to the for-mation of patterns: the same positive feedback that can lead to patterning can cause alaterally confined ice stream flow to oscillate in time (MacAyeal, 1993; Sayag and Tziper-man, 2011; Robel et al., 2013).The work in this thesis addresses one aspect of the non-linear evolution of a fully-formed icestream: What controls the migration of ice stream margins? In particular, the transition froma slowly flowing ice ridge to a fast flowing ice stream may correspond to a transition in thermalconditions at the ice stream bottom. How does the location at which the bed goes from thawedto frozen evolve over time? Modelling efforts during the 1990s and more recent analytical workare our main source for understanding the relevant processes, but these models are incompleteand cannot be upscaled to quantify temporal changes at the scale of an ice stream or ice sheetas a whole (Raymond, 1996; Jacobson and Raymond, 1998; Schoof, 2004, 2012; Suckale et al.,2014).The purpose of this research project is therefore(i) to better understand the physical processes that lead to the migration of ice stream mar-gins, and(ii) to derive a parametrization of these processes that can be used to better quantify icestream dynamics in continental-scale ice sheet models.1.3 Thesis outlineSince water plays a key role in facilitating the fast flow of ice streams, in chapter 2 we start byinvestigating the role of subglacial drainage in ice stream dynamics. In this chapter we couple athin film model for the ice stream flow to a simple model of the ice stream bed in which watercan move laterally and where freezing and melting can occur.6Figure 1.3: Sketch of the ice stream geometry considered in this thesis. A fast flowing ice streamis bordered by slowly moving ice, the ice ridge.We find that subglacial drainage cannot fully explain the widening of ice streams, as freezingin the subglacial sediment likely creates a thermal barrier for subglacial water flow. Where iceflows quickly (in the ice stream center) or velocity gradients are high (in the ice stream margin)heat is dissipated and melting at the ice stream bed is high, while little heat is dissipated inregions of almost stagnant flow. This leads to freezing of the ice stream bed in slowly movingregions and eventually to a completely frozen bed in the regions of slow flow.We also show that widening due to heat dissipation in the ice stream margin cannot bemodelled with thin film models. To model the outwards migration of ice stream margins,thermal and mechanical processes at the ice thickness scale must be included, but thin-filmmodels do not resolve these processes. To address this problem, we derive a boundary layermodel for ice stream margins in chapter 3.To illustrate the purpose of the boundary layer model, we consider the ice stream geometryin figure 1.3. An ice stream flowing with a velocity u in the x-direction is bordered by slowlymoving ice at its sides. This ice moves towards the ice stream, providing an inflow qin, seefigure 1.3. Modelling the evolution of the ice stream requires knowledge of the ice streamvelocity field u, the inflow qin and how the width of the ice stream changes over time, i.e.,knowledge of ∂W/∂t.In chapter 2 we analyse the dependence of u on different physical processes, and how W mayevolve due to subglacial drainage while the bed under the ice ridge is not completely frozen.7Once the bed under the ice ridge is completely frozen, the widening rate ∂W/∂t cannot bemodelled in this way. It is then controlled by the processes in the margin, for which a detailedmodel is necessary. Such a model is derived in chapter 3 by considering the margin as a boundarylayer between ice stream and ice ridge.We identify two processes that determine the temperature field in the ice stream margin:heat dissipation through lateral shear, which is controlled by the far field properties of theice stream, and inflow of cold ice from the ridge. In the boundary layer model, each of theseprocesses is parameterized through a single parameter. The balance between these processesalso determines the margin migration velocity.By solving the coupled mechanical and thermal problem posed in the boundary layer inchapter 4, we can determine that migration velocity. In order to derive a simple parameteri-zation that can be used in continental-scale models, we then consider the limit in which heatdissipation through lateral shear and inflow of cold ice are large, a limit that is likely to applyin real ice stream margins. Finally, we also consider the effect of different mechanical basalboundary conditions at the transition between ice stream and ice ridge, and how the formationof temperate ice (ice that is at the melting point) affects these migration velocities. We presentconclusions and avenues for future work in chapter 5.82 | Ice stream flow and basal freeze-on2.1 IntroductionTwo qualitatively different mechanisms for generating fast ice flow have been proposed in theliterature. The first is through topographic control. Thicker ice flows faster because the samevertical strain rate will lead to greater velocities in a thicker ice column. Thicker ice is alsoassociated with larger driving stresses and, by providing additional thermal insulation, withhigher temperatures in basal ice. Both of these effects lead to larger strain rates and furtherincrease ice flow velocities (see e.g., Bennett, 2003).The second category of streaming mechanisms relies on positive thermomechanical feedbacksrather than on pre-existing topography. In general, faster ice flow will lead to more heatdissipation, which can weaken either the ice itself or the ice-bed contact, permitting even fasterflow. Observations suggest that this type of mechanism is responsible for the existence of theSiple Coast ice streams in West Antarctica, which appear to be only weakly topographicallycontrolled (e.g., Shabtaie and Bentley, 1987; Bindschadler et al., 1996).The nature of the positive thermomechanical feedback depends on the assumed ice flowdynamics and the properties of the bed. For instance, in colder ice where there is no basalsliding, an increase in ice temperature leads to a lowering of the ice viscosity. This allows iceto flow faster and dissipate more energy through strain heating, warming the ice further andpotentially leading to streaming features (e.g., Clarke et al., 1977).Once the melting point has been reached, rapid sliding becomes possible and heat dissipationcan be focused at the bed. The dominant positive feedback mechanism now involves melting atthe base of the ice: if faster sliding leads to more heat dissipation at the bed, then this producesadditional melt water which ‘lubricates’ the bed and permits even faster sliding.There are two basic concepts of how this ‘lubrication’ works. Common to both is that theproduction of melt water raises water pressure at the bed, or equivalently, reduces effectivepressure, defined as the difference between overburden and basal water pressure. A reduction9in effective pressure is known to reduce friction at the bed for both, hard beds (e.g., Iken andBindschadler, 1986) as well as the deformable sediments found under the Siple Coast ice streams(e.g., Tulaczyk et al., 2000a).The difference between the two concepts lies in the link between melt water production andthe drop in effective pressure. In one version of the positive feedback, the bed is assumed tobe effectively ‘undrained’. Any melt generated is assumed to be stored locally at the bed, forinstance in subglacial till. If additional water storage (e.g., through a higher porosity of the bed)requires a decrease in effective pressure, then we have the necessary ingredients for a positivefeedback between faster sliding and lower effective pressure (e.g., Tulaczyk et al., 2000b; Sayagand Tziperman, 2008).A second explanation assumes that any generated melt is evacuated via an active drainagesystem. The production of more water requires more water to be drained, therefore the capacityof the drainage system must increase with decreasing effective pressure in order for the positivefeedback to operate. If an increase in drainage capacity requires a lower effective pressure, thenan increase in water production will lead to a lower effective pressure and therefore faster sliding(e.g., Fowler and Johnson, 1996; Kyrke-Smith et al., 2014).In the context of these positive thermomechanical feedbacks, the formation of ice streamscan be understood as a ‘fingering instability’ (Hindmarsh, 2009). The observed patterning ofice streams is a result of some regions flowing faster due to the positive feedback, while otherregions experience the reverse: slower ice motion leads to less heat dissipation and hence toeven slower ice motion. Successful attempts to model this include Hindmarsh (2009), who usesthe warmer ice/lower viscosity variant, and Kyrke-Smith et al. (2014), who employ the moremelting/higher drainage capacity variant of the positive thermomechanical feedback.Ice stream flow speeds vary on time scales of hundreds to thousands of years, including shut-down and reactivation cycles of fast ice stream flow (Hulbe and Fahnestock, 2007; Catania et al.,2012). Several studies have shown that the same positive feedback that explains the patternformation of ice streams can also explain this temporal behaviour. Faster ice flow, which elevatesheat dissipation rates, results in dynamic thinning (MacAyeal, 1993; Fowler and Johnson, 1996;Payne and Dongelmans, 1997; Fowler and Schiavi, 1998; Bougamont et al., 2003a; Christoffersenand Tulaczyk, 2003). When the ice thickness becomes too small, the amount of conductive10cooling at the bed will exceed the heat dissipation due to fast flow and the ice stream will shutdown. Robel et al. (2013) and Robel et al. (2014) showed that it is possible to reproduce bothsteady streaming and oscillatory (‘binge-purge’) behaviour with a zero-dimensional box modeland a similar flow line model that couple the evolution of the ice stream thickness to a simplifiedrepresentation of the energy balance at the bed.The feedback between increased heat dissipation and lubrication of the ice stream bed alsooperates in the reverse. In the ‘undrained’ view described above, insufficient heat dissipationleads to freezing at the bed. If freezing occurs by removing water from the sediment andfreezing it cleanly on at the base of the ice, then a reduction in porosity, and hence an increasein effective pressure, suppresses sliding further and leads to an even further reduction in heatdissipation (Christoffersen and Tulaczyk, 2003; Bougamont et al., 2003b). There is howeveranother possibility, namely that freezing occurs in the sediment pore space itself (Rempel,2009), with subtly different consequences. Rather than gradually strengthening basal resistanceto flow, freezing can cause a thermal barrier to form in the shape of a completely frozen sedimentlayer, which may block any residual lateral water transport, and which may abruptly suppresssliding.Ice streams are typically bordered by slowly flowing ice, which suggests that the feedbackabove operates in both directions in close spatial proximity. This is likely to result in largegradients in effective pressure between the ice stream bed and the surrounding areas. Even inthe undrained view this suggests that one may not be able to ignore drainage altogether, atleast in the lateral direction. Here we are interested in how lateral drainage driven by theseeffective pressure gradients can cause established ice streams to widen, and how this competeswith different freezing mechanisms.Although our model is built on the positive feedback between water storage and dissipation,we investigate the role of lateral drainage in widening ice streams by using a ‘weakly drained’model. This model allows for lateral drainage of melt between the ice stream margins and theice stream interior, but neglects drainage in the along-flow direction by essentially assuming aninfinitely long domain. Our geometry of an ice stream cross-section is similar to that in Sayagand Tziperman (2008).A disadvantage of this approach is that there cannot be net ice transport out of the domain11Description Variable Unitsvoid ratio eviscosity η Pa sthickness of wet sediment at e = 0 Hs mice stream thickness H mbasal melt rate m m/seffective pressure N Pawater flux in sediment layer qw m2/sbasal shear stress τb Pavelocity u m/sice stream width W mlateral coordinate y mTable 2.1: Table of variables.and we are unable to build ice ridges whose adverse surface slopes could prevent pressure-drivendrainage (Kyrke-Smith et al., 2014). The advantage is that we are left with a minimal modelof an ice stream cross-section in which dissipation, drainage, and different freezing processescompete to shape ice stream width and ice stream velocity field, and we can study the effect ofdifferent physical processes in the model.This chapter is laid out as follows: we describe the various components of our thin filmmodel in section 2.2. Subsequently, we study the evolution of ice stream width and velocity insection 2.3. One key outcome is that the model predicts the formation of quasi-steady velocityprofiles. We identify some of the physical and modelling subtleties involved and a possiblesource of numerical error leading to spurious continued ice stream widening in section 2.3.2. Insection 2.4 we study how different physical processes affect the steady states in more detail, andconclude with a discussion of our results in section Model descriptionWe consider an ice stream flowing with velocity u(y, t) in the x-direction, see figure 2.1. Ouranalysis is restricted to an ice stream cross section along the y–axis from y = −W0/2 to y =W0/2, perpendicular to the ice flow direction. Periodic boundary conditions are assumed onboth lateral sides at y±W0/2, and we assume a uniform ice thickness across the entire domainwith all time and space dependent variables depending only on the transverse coordinate y and12Figure 2.1: Sketch of the ice stream/sediment model. Left: A uniform slab of ice moves ontop of a sediment layer of thickness H0(e + 1) and width W0. Periodic boundary conditionsare assumed on both sides of the geometry at y = ±W0/2. Yellow/black: frozen sediment,blue/black: wet sediment. Panels A and B: To initialize freezing into the sediment, the effectivepressure N has to exceed a critical value Nc. While N < Nc, water is removed from the sedimentand frozen onto the base of the ice stream (marked as ‘basal ice’, see panel A). Once N ≥ Nc andfreezing into the sediment has begun, ice will intrude into the void space between the sedimentgrains and a frozen fringe will form, see panel B.time t, but not the downstream coordinate x. This assumption implies that no ice is advectedsideways.The ice stream sits on a sediment layer of thickness H0(e+ 1), where e(y, t) is the sedimentvoid ratio. H0 is the thickness the sediment would occupy if all pore space was removed (i.e.,at e = 0). H0 is hereafter called the ‘unexpanded sediment thickness’. The sediment layercan be completely unfrozen, partially frozen or completely frozen. In the former two cases weassume the unfrozen part of the sediment to be always fully saturated with water. We trackthe thickness of the unfrozen portion of the sediment layer, which is denoted by Hs(y, t) with0 ≤ Hs ≤ H0. Consequently, the volume of water per unit area bed at point y equals eHs.2.2.1 Ice flowIn the presence of a water saturated sediment layer, i.e. for Hs > 0, the motion of ice is modelledas a plug flow or ‘shallow stream’ (MacAyeal, 1989). The velocity u is controlled by a balance13of the down slope component of gravity, f , basal shear stress τb and lateral shear stress∂∂y(ηH∂u∂y)+ f − τb(u,N) = 0 if Hs > 0, (2.1)with the viscosity given by a regularized form of Glen’s law (Paterson, 1994)η =B¯21/n(∣∣∣∣∂u∂y∣∣∣∣2+D20)(1−n)/(2n), (2.2)where B¯ is a temperature-dependent viscosity factor, averaged through the ice thickness. Tokeep the model simple, we assume that B¯ does not depend on the lateral coordinate y. Thiseffectively corresponds to the assumption that the depth-averaged temperature is constant alongthe width of the ice stream. D0 is a regularization constant that ensures a finite viscosity in thelimit where lateral shear stresses are small. H is the thickness of the ice stream, τb is the basalshear stress, and n = 3 is the exponent in Glen’s flow law (Paterson, 1994). If Hs = 0, the iceis frozen to the bed, and we assume that this implies zero sliding velocity. Hence u = 0 is alsothe boundary condition for regions with Hs > 0 at the transition to Hs = 0. The driving stressis f = ρgH sin θ, with ρ the density of ice, g the gravitational constant, and sin θ = ∂H/∂x thesurface slope in downstream direction.2.2.2 SedimentThe feedback between fast flow and heat dissipation described in the introduction to this chap-ter requires that the basal shear stress decreases with decreasing effective pressure (N). Theeffective pressure is defined as the difference between ice overburden pi and water pressure pwin the sediment, i.e.∂τb∂N> 0 with N = pi − pw. (2.3)We consider two possible friction laws (constitutive relations) for τb. One is a Coulomb-plasticrelation as suggested by Kamb (1991)τb = µNu√u2 + u20, (2.4)14Description Symbol Value Unitspermeability parameter a ≥ 0temperature dependent viscosity factor (value for n = 3) B¯ 4.9× 10−25 s−1 Pa−3compressibility of sediment C 0.12constant in the power-law friction law c 1regularization of viscosity D0 5× 10−9 s−1regularization of void ratio δ 0.1reference value of void ratio e0 0.69total thickness of sediment at e = 0 H0 1 mpermeability parameter K > 0 m3/s Paa−1thermal conductivity k 2.3 W m−1 K−1along-flow length scale L 1000 kmlatent heat Lh 334× 103 J/kgfriction coefficient µ 1rheology parameter n 3reference value of effective pressure N0 103 Pacritical value for freezing Nc 104 Pageothermal heat flux qgeo 40× 10−3 W/m2exponent in basal shear stress r ∈ (0 1)density of water ρw 1000 kg/m3temperature difference ∆T 42 Kregularization of velocity u0 3× 10−7 m/sTable 2.2: Table of parameters, their values (where applicable) and their units.in a regularized form, where µ is the friction coefficient and u0 is a small regularization param-eter. Additionally, we consider a power-law friction law,τb = cNur (2.5)where c is a constant and 0 < r < 1.The effective pressure is related to the water balance in the sediment through the void ratioe. The positive feedback between dissipation and fast ice flow discussed in the introductionrelies not only on τb increasing with N but also on e decreasing with N , ∂e/∂N < 0. Tulaczyket al. (2000a) empirically found such a relationship between the void ratio and the effectivepressure, of the form e = e0 − C log(N/N0). We adopt their expression, but regularize it to15prevent the void ratio from becoming negative for large values of N > N0 exp(e0/C), i.e.,e = e0 − C log(NN0 + δN). (2.6)δ is a regularization parameter.The specific bed water volume eHs can evolve over time due to a lateral water flux qw in thesediment and at the ice-sediment interface, and due to melting or freezing:∂(eHs)∂t+∂qw∂y= m. (2.7)We assume that the water flows in the direction of water pressure gradients (and equivalentlyeffective pressure gradients)qw =κ(N)∂N∂y if Hs > 0,0 if Hs = 0,(2.8)where the effective permeability κ can be a function of the effective pressure, reflecting changesin the sediment void ratio or the opening and closure of gaps along the ice-sediment interface(Hewitt and Fowler, 2008; Creyts and Schoof, 2009; Schoof et al., 2012). The case for Hs = 0in (2.8) accounts for the fact that water cannot move through a solidly frozen sediment matrix.Consequently, we also use qw = 0 as the relevant boundary condition on (2.7) at the edgeof regions where Hs = 0. We assume a simple power-law for the dependence of the effectivepermeability on effective pressure Nκ = KN−a, (2.9)with a ≥ 0. For a = 3 this model is identical to a model proposed by Walder and Fowler (1994)for channelized flow over a deformable bed.The thickness Hs can change over time due to melting and freezing processes. Assume thatthe sediment layer is initially completely unfrozen, so that ice and water are in direct contactin the pore throats at the upper boundary of the sediment layer, see figure 2.1A. To extend theice surface into the pore space between the grains, surface tension has to be overcome, whichrequires the effective pressure to exceed a critical value Nc (see e.g., Rempel, 2009). If freezing16occurs while N ≤ Nc, then water is removed from the sediment and frozen onto the base of theice stream. Freezing is initiated when N = Nc and a frozen fringe will form, i.e. ice will intrudeinto the void space between the sediments, see figure 2.1B.Once a frozen fringe has formed, any further freezing will occur inside the sediment and nofurther water will be removed from it. Likewise, if there is a frozen fringe, any melting thatoccurs will release the frozen-on sediment at the void ratio e(Nc) at the critical effective pressurefor fringe initialization Nc. The evolution of the unexpanded thickness Hs of unfrozen sedimentis therefore modelled ase(N)∂Hs∂t=0 if N < Nc, and Hs = H0, and m < 0,m if (N ≥ Nc or 0 < Hs ≤ H0), and Hs > 0, and m < 0,e(N)e(Nc)1+e(Nc)1+e(N) m if 0 ≤ Hs < H0, and m > 0,0 otherwise.(2.10)2.2.3 Energy balanceAssuming that vertical heat conduction is the dominant form of energy transport in the ice, themelt rate m can be calculated as (Tulaczyk et al., 2000b):m =1ρwLh[qgeo − k∆TH+ τbu+ σηH(∂u∂y)2], (2.11)with ρw the density of water and Lh the specific latent heat. qgeo is the geothermal heatflux, k∆T/H describes the conductive cooling through vertical temperature gradients, k is thethermal conductivity of ice, and ∆T is the temperature difference between ice stream surfaceand base. τbu is the heat dissipation through basal sliding with τb the basal shear stress.Additionally, we have introduced a term σηH(∂u/∂y)2, which does not appear in Tulaczyket al. (2000b). Our results in chapter 3 show that we expect heating through lateral shearstresses and the advection of cold ice from the ridge to affect the temperature field in the icestream margin. ηH(∂u/∂y)2 is the vertically integrated heat dissipation rate through lateralshear, see appendix A.1 for details. If conduction alone was responsible for heat transport in17the slab, then half of this heating would reach the bed, so σ = 1/2. We also allow fractions ofσ/2 in order to parameterize crudely the effect of advection removing heat from the domain.2.3 Ice stream evolutionWe first consider the temporal evolution of the model above. A natural first step would beto find spatially uniform steady states and consider their stability. This has already beendone in Schoof and Hewitt (2013), using a version of the model that excludes the elaboratefreezing dynamics of (2.10). The authors confirm that the positive feedback here is one betweendissipation and water storage: although we have a drainage mechanism, instability occurs if andonly if e is a decreasing function of N . Models that mediate the instability through drainagemust allow for net drainage in the downstream (x–) direction (Kyrke-Smith et al., 2014), whichwe do not. A further complication of an effectively infinitely long domain in the x–directionis that it also does not allow for the suppression of the instability on long wavelengths as hasbeen previously found by Hindmarsh (2006) and Sayag and Tziperman (2008). Our aim hereis however different: we are interested in the evolution of a fully formed ice stream, and wewant to investigate whether drainage can widen such an ice stream and how this interacts withfreezing of the bed.Starting from a prescribed set of initial conditions, we solve the sediment model consistingof (2.4)–(2.11) numerically using finite volumes with an implicit Euler method, and (2.1)-(2.2)together with (2.4) and (2.5) using finite elements. The sediment model and ice flow modelare coupled through the basal shear stress, which depends on both the effective pressure andthe velocity. To deal with this coupling, we use an operator splitting. The model is integratedforward in time until temporal variations are below machine accuracy, at which point we assumethat a steady state has been reached.2.3.1 Outwards migration and evolution to steady statesWe start our simulation from an initially completely unfrozen sediment layer, see panel a infigure 2.2. The initial value for the effective pressure is chosen to be low (varying between 518Figure 2.2: Temporal evolution of ice stream flow that is initially localized by a trough ineffective pressure. Each panel from top to bottom shows velocity, effective pressure, melt rate,and a visualization of the sediment layer, respectively. Panel a: initial state, panel b: after 250years, panel c: 500 years, panel c: steady state configuration. σ = 0.01, a = 2.kPa and 5.5 kPa) in the centre of the domain. Outside of the centre, the effective pressure waschosen to be at a large value of about 100 kPa. These initial conditions lead to a localizedflow in the centre, as the effective pressure determines the basal resistance to sliding throughequations (2.4) and (2.5).As time advances, water starts to seep in the direction of effective pressure gradients, leadingto lower effective pressures at the sides of the ice stream, see panel b in figure 2.2. The rowsof figure 2.2 are snapshots of the velocity u (top row), effective pressure N (second row), meltrate m (third row), and a 1 m thick sediment layer (bottom row). Panel a shows these statesat initialization, panel b after t = 250 years, panel c after t = 500 years, and panel d showsthe final steady state. The positive range of the melt rate is displayed on a logarithmic scale.The sediment layer is shown in blue if it is unfrozen and water saturated, while yellow indicatesfrozen sediment.Outward seepage of water leads to a speed up of the flow in the regions immediately adjacentand consequently a widening of the ice stream. However, where the velocities are slow the meltrate is negative, causing the sediment to freeze. Once the sediment is completely frozen (Hs = 0,19Figure 2.3: Example of the grid-size dependent margin migration rate with a centred differencescheme for heat production. Panel A: Position of ice stream margins against time for differentgrid sizes. Panel B: velocity shown at individual grid points at some instant during the cal-culation. With a centred difference scheme, there is a finite strain rate computed at the firstgrid point where u = 0 (vertical line). Panel C: melt rate, panel D: frozen sediment thicknessaround the corresponding point.panel c of figure 2.2) the ice is stuck to the bed in these regions and no further drainage drivenmargin migration is possible. The stream evolves into a steady state.In the steady state (figure 2.2, panel d), the regions of high heat dissipation, in this casethe margins, supply water to regions with higher effective pressure where freezing occurs. Notethat this pattern is the reverse of that shown in panel c. The reason for this is that in steadystate the high concentration of shear in the margins leads to high heat production. Even withσ = 0.01 this dominates. The resulting steady state is characterized by a balanced overallenergy budget, i.e., the integrated melt rate over the whole range of flow is zero. Thereforethe high heat dissipation in the margin needs to be offset by freezing in the ice stream center,facilitated by near-zero basal friction τb and negligible dissipation there. Outside the streamthe ice is stagnant and no heat is dissipated.2.3.2 Spurious margin migrationDespite the proximity of regions with high heat dissipation to regions with a frozen bed, ourthin film model cannot capture the necessary lateral heat transport to simulate any outwardmargin migration once the sediment is frozen. We study the latter problem in chapters 3 and4.20Figure 2.4: Temporal evolution of a simulation similar to figure 2.2 but with an initially narrowertrough in effective pressure, leading to the shutdown of the ice stream. Same plotting schemeas in figure 2.2. σ = 0.05, a = 2.A ‘false’ margin migration into the frozen sediment layer can only arise due to discretizationartifacts. This is illustrated in figure 2.3, which shows close-ups of the margin in panels B –D. The dashed vertical line in these panels marks the point at which the sediment column isfrozen throughout. At this point, no sliding is possible and the velocity is zero. Consequentlyno heat can be dissipated here and the melt rate should be zero as well. If however the velocityderivative in (2.11) is calculated through some kind of averaging scheme, for example a centreddifference scheme∂u∂y∣∣∣∣i=ui+1 − ui−12∆y, (2.12)where i denotes the finite volume cell center at which the relevant quantity is computed, thenthis will result in a non-zero melt rate where the velocity is zero. This is a purely numericalartifact, which would become apparent through a grid-size dependent migration rate computedfrom the centred difference scheme (2.12) as shown in figure 2.3A. We avoid this effect (whichmight conceivably be present in some more elaborate ice sheet simulation codes) through aninward scheme at the boundaries between temperate and frozen sediment.212.3.3 Inwards migration and ice stream shutdownThe final steady state shown in figure 2.2 depends on the assumed initial conditions: the bed farfrom the regions of initially fast flow froze on a time-scale of a few hundred years. This preventeda further widening of the ice stream into these regions via subglacial drainage. Subsequently,the confined ice stream evolved into a steady state.To investigate how different initial conditions affect the final steady state, we choose a differ-ent set of initial conditions in figure 2.4. As before, we consider an effective pressure perturbationon an unfrozen, water-saturated sediment, but the lateral extent of this perturbation is now onlyabout 30 km wide. Again, the sediment starts to freeze in the regions of low velocity/high ef-fective pressure. In contrast to the previous case, the melt rate does not reach positive valuesin the regions of faster flow. This leads to a continuous increase in effective pressure as wateris removed from the sediment for the formation of basal ice (see panel b in figure 2.4). Oncethe effective pressure reaches the critical effective pressure for freeze-on, Nc (marked with thedashed red line), the sediment starts to freeze in the centre of the domain (panel c). This leadsto a completely frozen sediment and shutdown of the ice stream.We can understand this behaviour in terms of heat dissipation: the melt rate (2.11) dependson the amount of heat dissipated at the bed (τbu) and in the margin (σηH(∂u/∂y)2), and onthe background heat flux (q0 = qgeo − k∆T/H). If this background heat flux is less than zero,q0 < 0, freezing of the bed can only be prevented by heat dissipation. Where heat dissipationis consistently less than the background heat flux, as in the example in figure 2.4, the entiresediment experiences freezing and the ice stream will shut down.It is straightforward to show that at constant effective pressure both, the ice stream velocityin the center and the heat dissipation in the margin increase with ice stream width to the powerof n + 1, i.e. u ∼ W n+1 and η(∂u/∂y)2 ∼ W n+1. We therefore intuitively expect wider icestreams to be able to dissipate more heat. Conversely, narrower ice streams are expected todissipate less heat. This explains qualitatively why the example in the last section, figure 2.2,developed into a steady state of constant streaming, while the example in this section, figure2.4, led to an ice stream shutdown. In the next section, we explore the dependence of steadystates on the model parameters in more detail.222.4 Ice stream steady statesSo called ‘flow line models’ of ice streams assume that the mechanics of lateral shear that we havedescribed above can be integrated straightforwardly to produce a rate of ice discharge along theflow line as a function of ice stream width, thickness, surface slope, and possibly an additionaldynamic variable such as effective pressure N at the center line of the ice stream (e.g., Vander Veen and Whillans, 1996; Raymond, 2000; Bougamont et al., 2003b; Robel et al., 2013). Inparticular, Robel et al. (2013) use this approach to compute a quasi-steady relationship betweendischarge and ice stream geometry alone, which they then use to explain temporal oscillations inice stream flow and thickness. We reinvestigate this quasi-steady relationship between dischargeand geometry, but pay closer attention to variations in effective pressure, ice velocity and heatdissipation across the width of the ice stream.The steady states computed above correspond to ice streams flowing over a relatively weakbed that is bordered by a fully frozen bed. We therefore define the width W of an ice streamas the lateral extent of unfrozen sediment over which the ice stream flows. At the lateral icestream margins, the velocity u and water flux qw are zero. Centering steady state solutionsaround the origin, the ice stream margins are then at y = ±W/2 and the boundary conditionsfor steady state ice streams areqw =∂N∂y= u = 0 at y = ±W/2. (2.13a)These equations require symmetry around the centre of the ice stream which is at y = 0, sothat∂u∂y=∂N∂y= 0 at y = 0. (2.13b)In line with Robel et al. (2013), we also parameterize the slope sin θ = ∂H/∂x which appearsin our model above as sin θ = H/L, where L is ice stream length.23For ease of reference, the remainder of the steady state model is:B¯H21/n∂∂y[∣∣∣∣∂u∂y∣∣∣∣1/n−1 ∂u∂y]− τb + ρgH2L= 0, (2.14a)K∂∂y(N−a∂N∂y)−1ρwLh[qgeo − k∆TH+ τbu+ σB¯H21/n∣∣∣∣∂u∂y∣∣∣∣1/n+1]= 0. (2.14b)For K =∞ and σ = 0, this model is identical to the model for ice velocity described by Robelet al. (2013). Here we can go further and investigate the dependence of ice stream velocity andhence ice stream discharge not only on geothermal heat flux and on a temperature difference∆T as was done in Robel et al. (2013), but also on K and σ. This helps us identify the rolesplayed by different processes in ice stream flow. In the next sections, we will calculate thesesteady states both semi-analytically and numerically by solving (2.13a)–(2.14b) as a boundaryvalue problem using a shooting method with the Matlab ODE solver suite.2.4.1 General steady state behaviour and stabilityBefore we solve (2.13a)–(2.14b) below it is worth noting two bounds on steady state solutions.Physically, it is not possible for the basal drag to be less than zero, which imposes an upperlimit on the ice stream velocity, marked as a red line in figure 2.5 (see appendix A.2 for details).Moreover, for H > k∆T/qgeo, the melt rate is always positive. In this case, the magnitudeof the geothermal heat flux exceeds the magnitude of the conductive cooling term and morewater is produced than can be frozen on. Therefore, no steady state solutions exist beyondthis point. Physically, we expect in this case that the ice and bed become disconnected. Thiswould lead to the opening of drainage pathways which can discharge the excess melt waterdownstream. However, the model used here does not resolve these processes, and in its originalversion without any lateral confinement fixing W would lead to melting everywhere and anindefinite widening of the ice stream. This limit is marked by a vertical grey dashed line infigure 2.5.In this section we assume that K →∞ and σ = 0. For u0 = 0 in the Coulomb-plastic basalshear stress, the model can then be solved analytically (see appendix A.3). These solutions areshown as blue lines in figure 2.5. With an infinite permeability, the effective pressure is constant240 500 1000 1500 2000 2500100102104H in mu c in m/year  maximum velocityunstablestableFigure 2.5: Steady state solutions and stability. Plotted is velocity in the ice stream center,uc = u(y = 0), against ice stream thickness H. The red line is the upper bound on the velocityobtained by setting τb to zero. The vertical grey dashed line corresponds to the critical thicknessat which geothermal heat flux can just be balanced by conductive heat loss. The circles andarrows indicate the results of dynamical computations. The open circle is the initial condition,the filled dot the final steady state. a = 0, σ = 0, and K → ∞ for the analytical solution,K = 106 Pa−1 m3 s−1 for the numerical calculation.and we obtain the model of Robel et al. (2013). As in the latter paper, we find that there is nonon-zero velocity steady state below a critical ice thickness Hc: conductive cooling is then toostrong and frictional dissipation too weak to obtain an overall heat balance.One high and one low velocity branch emerge from the saddle-node bifurcation at Hc. Theexistence of these two steady solutions is due to the dual role of the effective pressure in bothlubricating the flow (low effective pressure N leads to high velocities u) and producing melt (thebasal shear stress τb increases with N). The dissipation rate τbu required to offset the conductivecooling at the bed can either be generated by a low effective pressure but high velocities (upperbranch in figure 2.5), or a high effective pressure and slow velocities (lower branch).The lower branch is unstable in the same way as in Robel et al. (2013) through the basicfeedback between dissipation, effective pressure and velocity as indicated by arrows in figure2.5. On this branch, a decrease in effective pressure leads to a sufficiently large increase invelocity so that net dissipation also increases, causing runaway acceleration of the ice streamuntil the upper branch is reached. This positive feedback is possible because lateral shearinggenerally plays a minor on the lower branch and velocities are primarily set by basal friction25and the driving stress (see appendix A.4). On the upper branch the feedback is suppressedbecause velocities are limited by lateral shearing. Further decreases have negligible effect on icevelocity and therefore actually reduce dissipation rates τbu.We see that in our model lateral shearing plays an essential role in stabilizing fast ice streamflow, that is in allowing the upper branch to exist. This can be contrasted with other modelsfor fast ice stream flow in which the positive feedback that makes our lower branch unstableis mediated by changes in the permeability of the subglacial drainage system, rather than thebed water storage capacity (e.g., Fowler and Johnson, 1996). In models of that kind, increasedwater production and lower effective pressures need an increased capacity to drain water outof the system (which we have suppressed by having an effectively infinitely long domain), andthis eventually stabilizes the velocity even in the absence of lateral shear, see also Kyrke-Smithet al. (2014).The multi-valued structure of the velocity against thickness relationship with an upper,stable and a lower, unstable branch can be used to explain ice stream surges: such surges occurif the ice discharge needed to account for net accumulation over the ice stream corresponds toa value of velocity on the lower branch (Robel et al., 2013).Our model includes two effects not considered by Robel et al. (2013): the effect of lateralshear heating, parameterized through σ, and a finite lateral hydrological permeability of thebed, parameterized through δ. In the next two sections, we investigate how inclusion of theseeffects changes the steady state solutions considered thus far.2.4.2 The effect of lateral shear heatingHeat production due to lateral shearing has a significant impact on water budget of an icestream. A steady state requires net energy (or water) balance at the bed with melting in someparts of the bed balanced by freezing in others. In the absence of melt supplied due to shearheating (as in the previous section), we know that this balance fails when geothermal heat fluxoverwhelms the conductive heat loss to the surface when H = k∆T/qgeo.The presence of the additional heat production term can cause the net balance to fail (sothere is no steady state solution) for values of H < k∆T/qgeo. In that case the stable upper260 1000 2000100102104H in mu c in m/yearaupper limitցa1900 1000500100020003000H in mu c in m/year a1←σ=0σ=1/4րցσ=1/8 −50050q w in m2 /yearσ=0b1−20 0 20−10010y in kmm in mm/year b2−50050q w in m2 /yearσ=1/8c1−20 0 20−10010y in kmm in mm/year c2Figure 2.6: Effect of shear heating on ice stream steady states. Panel a: center line velocitiesuc against ice thickness H for σ = 0, σ = 1/8, and σ = 1/4. Solid lines indicate stable steadystates, dashed lines unstable ones. Panels b1–b2: subglacial water flux (panel b1) and basalmelt rate (panel b2) for σ = 0 at the point marked with a black star in panel a1. Panels c1–c2:same as panels b1–b2 but for σ = 1/8 at the point marked with a magenta circle in panel a1.Calculations were done with K = 10−3 Pa−1m3s−1 and a = 0.solution branch of figure 2.5 has a shorter extent, with the saddle-node bifurcation at a lowervalue of Hc because the addition of shear heating permits steady states to exist even if thicknessis smaller and conductive cooling is stronger, see figure 2.6a1. When there is significant heatproduction due to lateral shearing, the distribution of the melt rate m(y) can also become thereverse of the pattern observed when there is no contribution from lateral shear (σ = 0). Forσ = 0, the only source of heat dissipation is frictional sliding, which is strongest at the centerof the ice stream where velocities are large (figure 2.6b2). Melting then occurs in the centerof the ice stream with freezing in the margins, and net water redistribution from the center tothe sides, see panel b1 in figure 2.6. The opposite is the case for steady states with moderateσ: very little friction and sufficiently strong cooling lead to net freezing in the center of the icestream (panel c2), so that water is produced in the margins and redistributed to the center, seepanel c1 in figure 2.6.At larger values of σ (but still less than 0.5, see the σ = 1/4 case in figure 2.6a), the stablebranch can disappear. This is however different from the failure at H = k∆T/qgeo. In the caseof H ≥ k∆T/qgeo the bed experiences melting everywhere, so that the ice stream is no longerthermally confined. Conversely, for σ > 0 we can fail to have net melt balance in the bed under27the ice stream and still maintain a frozen bed outside. In dynamic solutions of the latter casethe effective pressure simply reaches zero, so ice and bed decouple but water is still producedby lateral shear heating. In reality we expect that net downstream drainage starts to play arole before this happens, which we have excluded in our model.2.4.3 The effect of subglacial drainageThe previous section has already suggested that drainage is an important component in theenergy balance of an ice stream, even though a number of ice stream models ignore or onlypartially incorporate it, because the main feedback in these models (as in ours) is between iceflow, heat dissipation and local water storage (Tulaczyk et al., 2000a; Bougamont et al., 2003b;Robel et al., 2013).Our simple model cannot give a complete understanding of the role of drainage, because anydownstream flux of water in our model would have zero divergence by construction. However,the results in figure 2.6 suggest at least that lateral water transport is important in balancingexcess water production in some parts of the ice stream with net freezing elsewhere. We can infact show that no steady states exist if K = 0 and there is no drainage.In the absence of shear heating (σ = 0), we haveB¯H21/n∂∂y(∣∣∣∣∂u∂y∣∣∣∣1/n−1 ∂u∂y)− τb + ρgH2L= 0 (2.15a)qgeo −k∆TH+ τbu = 0. (2.15b)The second of these equations allows us to eliminate τb, independently of the friction law.Introducing τ = |∂u/∂y|1/n−1∂u/∂y and using the chain rule, we getB¯H21/nτ |τ |n−1dτdu−(k∆TH− qgeo)1u+ ρgH2L= 0. (2.16)Integrating givesB¯H21/n|τ |n+1 =(n+ 1)H[ρgH2Lu+(k∆TH− qgeo)log u]. (2.17)2810−6 10−4 10−2920940960980100010201040K in Pa−1 m3 s−1u c in m/year  a−20020q w in m/year  b−20 −10 0 10 20234y in kmN in kPacstableunstableKcKcblue circles in aFigure 2.7: Dependence of steady state solutions on permeability K. Panel a: Center line icestream velocity against K. Panels b and c: steady state water flux (panel b) and effectivepressure (panel c) profiles for different values of K ≥ Kc (marked with circles in panel a). Forlarge K, N becomes approximately spatially uniform as only very small pressure gradients areneeded to drive O(1) fluxes. Calculations were done with H = 1024 m, W = 50 km, a = 0,σ = 0 and the values listed in table 2.2.With H < k∆T/qgeo the right hand side of (2.17) is negative for sufficiently small values of u,contradicting the fact that the left hand side must be non-negative. It is therefore impossibleto meet a u = 0 boundary condition when K = 0. There is insufficient heat production in themargins as u goes to zero, and a certain amount of drainage is needed to supply water to theseregions. Figure 2.7a confirms this: steady solutions (one stable and one unstable branch asbefore) are possible only for K above a critical value Kc. For sufficiently large K in figure 2.7,velocity becomes insensitive to changes in permeability. Note that figure 2.7 suggests that thestable and unstable states connect in a saddle-node bifurcation. This is not always the case asthe condition N ≥ 0 is implicit in the model (otherwise dissipation could be locally negative andtherefore violate the Clausius-Duhem inequality), which can prevent the stable and unstablesolution to collapse into a single solution at the bifurcation point (see appendix A.6.2 for anexample of this behaviour).For non-zero heat production in the margins (σ > 0), we again find that steady states arenot possible unless there is some drainage or W and H satisfy a specific functional relationshipthat corresponds to the marginal existence of the ice stream (meaning a narrower ice stream290 1000 2000100102104H in mu c in m/year  aupper limit←plastic, K=∞plastic, K=5×10−4 Pa−1m3s−1r=1/4, K=5×10−4 Pa−1m3s−1r=1/3, K=5×10−4 Pa−1m3s−1r=1/2, K=5×10−4 Pa−1m3s−11.61.611.62  c1plastic−20 0 2001000y in kmu in m/year b0.520.530.54  c2r=1/40.360.370.38                 N in kPa  c3r=1/3−20 0 in km  c4r=1/2Figure 2.8: Influence of the basal friction law on the effective pressure and steady states. Panela: Center line velocity against ice stream thickness H for different choices of the basal frictionlaw. The red line again shows the upper bound on the velocity, set by τb = 0. Panel b: velocityprofiles at the point marked with a red cross in panel a. They are indistinguishable. Panelsc1-c4: the corresponding effective pressure profiles. σ = 0 and a = 0 for all calculations.would shut down, see appendix A.6.1). Note that the drainage parameter a appears not to havea significant influence on the steady states (see appendix A.6.3).2.4.4 The effect of the basal friction lawFinally, we also consider the influence of the basal friction law τb(u,N). All of the computationspresented previously were done using a Coulomb friction law (2.4). We compare the stablebranch of the u–H relationship for a Coulomb friction law with analogous calculations for thepower law (2.5). The results are virtually identical when we consider the velocity u, whileeffective pressures differ between the results for the different friction laws. This echoes ourfailed attempt to find a steady state solution without drainage, equations (2.15a)–(2.17). Thereτb was simply determined by the need to maintain the energy balance of the bed, regardlessof the basal friction law. Friction simply adjusted itself in a steady state to produce sufficientdissipation. Even though we were ultimately unable to find a solution because the requiredlocal heat balance could not be maintained in the margins, the basic insights seem to hold forreal steady states with drainage: the precise form of the basal friction law is not particularlyrelevant.302.5 DiscussionIn this chapter, we have considered the effects of basal freeze-on and subglacial drainage on icestream evolution. We used a simple model of an ice stream cross-section flowing on a weaklydrained bed, allowing for storage of melt water and lateral but not along-flow drainage of meltwater. Melt water can be produced by dissipation due to basal friction and lateral shearing,though we parameterize the relevant heat fluxes in the ice somewhat crudely. A further keycomponent of the model is that we allow freezing to cause either the removal of water from theopen pores or in situ freezing inside the pore space. Thus freezing at the bed can both changethe effective pressure and the portion of the bed through which melt water can move.We find that subglacial drainage can drive the outwards migration of ice stream marginsas long as the bed outside the ice stream is unfrozen. This migration is the result of waterseeping from the ice stream bed, where effective pressure is low, to the regions outside the icestream, where effective pressure is high (or equivalently water pressure is lower than below theice stream). Outside the ice stream, velocities are slow and little heat is dissipated, so that thebed experiences freezing.In our model, we assume that once all the pore space in the subglacial sediment layer isfilled with ice, no further drainage is possible and sliding is suppressed. Our results show thatfreezing the entire sediment layer (assuming the latter is finite) eventually happens outside theice stream and this stops any further widening. In reality there may be other processes that stopwidening, for instance the formation of thickened ice ridges, whose surface slopes also suppressdrainage-driven widening (Kyrke-Smith et al., 2014).The erosion of the thermal barrier presented by a fully frozen sediment layer may be possiblethrough lateral heat distribution from the ice stream. However, this cannot be captured by athin film model of the kind we have used, or by models (including higher order models) that donot faithfully account for the mechanical and thermal physics at the ice thickness scale in the icestream margins (Jacobson and Raymond, 1998; Schoof, 2012). As our results in section 2.3.2illustrate, a thin film model can produce an apparently plausible ongoing margin migration,but in any thin film model this is likely an artifact of how velocity gradients are interpolatedacross the ice stream margin and is consequently grid dependent. Care is therefore advisable31when interpreting model results or implementing schemes for margin migration. In the nextchapter, we develop a boundary layer model for ice stream margins that can capture the physicalprocesses that lead to ice stream widening.By carefully handling the transition between frozen and unfrozen bed in our model, wegenerally obtain steady ice stream widths and velocity profiles. We refer to these as quasi-steady in the introduction of this chapter, because these states assume a constant ice thickness,which is hard-wired into our model, but in realistic ice streams the thickness H will change overtime.In steady state, the bed is in overall thermal equilibrium. If we were to persist with afully undrained model (there is no lateral water redistribution at all), this global balance hasto become a local balance. This is however not possible in general, as shown in section 2.4.3.Lateral drainage must provide melt to the margins if we ignore the role of lateral shearing inbasal energy balance altogether. Conversely, if we allow shear heating in the ice to contributeto the basal energy balance, excess melt production in the margins must be redistributed to theice stream center. Some amount of drainage is therefore essential to maintain a global energybalance.There is however a significant caveat here, because our representation of heat fluxes at thebed is oversimplified as we do not account for advection in the ice in any realistic way. Alsoour conclusion comes out of a thin film model for an ice stream, which ceases to be valid in themargins, where the energy balance problem inevitably occurs. It is therefore not impossible thatsteady states exist in a more sophisticated model that resolves heat transport in the horizontaland vertical direction.Considerable controversy exists about the appropriate friction law for ice stream flow, withboth, Coulomb-plastic and power-law friction laws having appeared in the literature (e.g., Boul-ton and Hindmarsh, 1987; Alley et al., 1987; Hindmarsh, 1997; Hooke et al., 1997; Iverson et al.,1998; Tulaczyk et al., 2000a). Our results indicate that if a positive feedback between ice flowand effective pressure is at least in part mediated by water storage in the bed, then the icestream velocity is primarily controlled by lateral shear. The basal shear stress then adjustsitself through the energy balance at the bed, and the particular form of the basal friction lawdoes not appear to affect the ice stream velocities significantly.32In addition to our simple handling of heat fluxes in the ice, the biggest simplification inour model has been the assumption of an infinitely long domain, so there is not net ice orwater transport out of the domain. At least net ice transport out of the domain appears tobe essential in order to generate ice streams through a pattern forming instability (Hindmarsh,2006, 2009; Sayag and Tziperman, 2008; Schoof and Hewitt, 2013; Kyrke-Smith et al., 2014).Net melt drainage out of the domain may equally be essential (Fowler and Johnson, 1996;Anandakrishnan and Alley, 1997; Le Brocq et al., 2009; van der Wel et al., 2013; Christoffersenet al., 2014; Kyrke-Smith et al., 2014). In fact, observations suggest that complex drainagenetworks are likely to exist at ice stream beds (Blankenship et al., 1986, 1987; Engelhardt andKamb, 1997; Kamb, 2001; Gray et al., 2005; Fricker et al., 2007). An extension of our modelto these additional physics is therefore desirable as is the incorporation of the detailed physicsthat allow margins to migrate even after the bed outside the stream becomes fully frozen.333 | Boundary layer modelIn this thesis, we have so far considered margin migration as the result of two different processes.One possible cause of margin migration is subglacial drainage. Assuming that the high ice streamvelocities can be attributed to the presence of water, seeping of water into less lubricated regionscan speed up the ice flow there and result in ice stream widening. This process is limited byfreezing at the ice base, and in particular it is not possible for water to flow into frozen sediment.Observations in the Siple Coast (Bentley et al., 1998; Catania et al., 2003) and our model resultsin section 2.3 show that the base of the ice ridges are likely to be frozen. Drainage driven marginmigration may therefore be the exception rather than the rule in causing the widening of theSiple Coast ice streams.The other condition under which margin migration is possible is a negative energy balance.In ice streams, heat can be dissipated along the bed due to frictional heating and due to lateralshear in the ice stream margins. Together with the geothermal heat flux, this heat counteractsthe cooling of the ice surface due to low atmospheric temperatures. If the ice stream bedexperiences net cooling, existing water at the bed will freeze. If additionally water supply fromother regions of the ice stream is insufficient to compensate this freezing, a loss of lubrication willresult, which in turn slows down the ice stream and can lead to its narrowing and its eventualshutdown.In this chapter and the next, we will concentrate on the second feedback, but operating inthe reverse direction. Heat dissipation in the margins may exceed the cooling effect of the icesurface, warming the previously frozen bed outside of the ice stream and allowing the ice streamto widen. In the previous chapter, we have seen that heat dissipation in a depth-integrated thinfilm model occurs on the ice stream side of the margin, but models of this kind typicallylack a lateral heat transport mechanism that would allow the dissipated heat to cause marginmigration.Conduction of heat at the scale of the ice thickness or below along with concentrated dissi-pation of heat at the same scale (Raymond, 1996; Jacobson and Raymond, 1998; Schoof, 2004,2012; Suckale et al., 2014) can lead to margin migration, but a model that resolves the relevant34physics cannot be depth-integrated. The migration speed then needs to balance the effect ofdissipation thawing the bed at low migration speeds, and the translation of the margin coolingthe bed at high migration velocities. Schoof (2012) uses this to compute a unique migrationspeed as a function of the heat dissipation and far field thermal conditions.A complication, first identified by Jacobson and Raymond (1998), is that lateral advection ofcold ice from the adjacent ice ridges may significantly affect the temperature field in the marginand therefore also margin migration. Schoof’s model ignores the inflow altogether, while othermodels have parameterized the inflow by prescribing a simple form of the transverse velocityfield in the margin (Jacobson and Raymond, 1998; Suckale et al., 2014).In this chapter, we derive a model for the transverse velocity field as well as the velocitycomponent parallel to the ice stream in the margin from first principles. We do so by developinga boundary layer description for an ice stream shear margin that separates a rapidly sliding plugflow from a slow vertically shearing flow. In the process, we also derive the appropriate heatequation for the shear margin, incorporating the effects of margin migration on temperature. Inchapter 4 we solve this model with a combination of numerical and asymptotic methods, leadingto a parameterization of the margin migration velocity in terms of characteristic variables ofthe ice stream and the ice ridge.3.1 Description of the geometryWe consider an ice stream flowing in the positive x-direction with velocity u, bordered at itssides at y = ±ym(x, t) by slowly moving ice, the ice ridges. The geometry we have in mindis illustrated in figure 3.1. The domain has length L in the x–direction and side boundariesat y = ±W , where we assume periodic boundary conditions to hold, in addition to symmetryabout y = 0. We also treat W as the characteristic width of the ice ridges. In addition to thegeometrical parameters L and W we also assume that the characteristic half width Ws of the icestream is known. z is the third coordinate, pointing vertically upwards. The basis of the leadingorder model we develop below will be that ice streams and ice ridges are much longer than theyare wide (W  L and Ws  L) and that their thicknesses are much less than their widths.This is physically plausible: for the Siple Coast ice streams in Antarctica we have typical values35x yzshearing flow v plug flow us(x,y,t)marginstream ridgeridge Ws Wmarginym-ymLFigure 3.1: Geometry of the problem. We consider an ice stream whose principal direction offlow is along the x–axis. The y–direction is transverse to the flow, z points vertically upwards.At its sides, the ice stream is bordered by slowly moving ice or ‘ice ridges’. The transitionbetween ice stream and ice ridge is referred to as the ice stream margin. It is located aty = ±ym. The ice ridges are symmetrical about y = ±W .Ws = 25 km, W = 50 km, L = 1000 km, and thicknesses between 1000 m and 1200 m.As the ice stream moves much faster than the surrounding ridges, it discharges more ice.This lowers the ice surface of the ice stream relative to the ice ridge, see e.g. the surface elevationmaps by Shabtaie and Bentley (1988). This difference in elevation drives inflow of ice from theridges into the ice stream.Mechanically, we assume that the ice flow of the ice stream occurs mostly by rapid slip at thewell-lubricated base of the ice stream (Blankenship et al., 1986, 1987). In contrast, we assumethat the underlying bed at the base of the ice ridge is frozen, permitting no sliding and ice flowsby vertical shear. In this chapter, we will will derive leading order models that describe both icestreams and ice ridges. The flow regimes in the two regions are mechanically distinct and thetransition between these regimes occurs in a boundary layer around the ice stream margin. Aswe will show, the mechanical and thermodynamic behaviour of the boundary layer determinesthe migration rate of ice stream margins, i.e., how ym(x, t) evolves in time.In the following section, we present the governing equations and boundary conditions forice flow and energy balance. We then proceed by constructing leading order models for theice stream (section 3.3), the ice ridge (section 3.4) and the boundary layer (section 3.5). Inthe process the model is non-dimensionalised in different ways in all three regions. The scale36Description Symbol Value Unitsfluidity factor of ice A 1 · 10−16 kPa−3 s−1specific heat capacity c 2 kJ kg−1 K−1acceleration due to gravity g 9.81 m s−2thermal conductivity k 2.3 W m−1 K−1geothermal heat flux qgeo 40× 10−3 W m−2density of ice ρ 920 kg m−3accumulation [a] 0.3 m year−1ice stream length L 1000 kmsurface temperature T0 240 Kmelting point temperature Tm 270 Kice stream half width Ws 25 kmdomain half width W 50 kmTable 3.1: Parameter values and characteristic scales.applied to variables is indicated by a squared bracket and subscripts denote the region underconsideration. For example, [u]r is the scale for the velocity component u in the ice ridge, whichwe generally non-dimensionalise as u = [u]rur unless otherwise indicated. ur is then the scaledvelocity in the ice ridge, [u]s is the scale for the same velocity component in the ice stream and[u]BL is the scale in the boundary layer. In the ice stream and ice ridge, scaled variables areindicated by a subscripted s and r, respectively. In the boundary layer we use capital lettersfor scaled quantities. In other words, the x–component of the velocity in the ice stream is us,ur in the ice ridge, and U in the boundary layer. Where we deviate from this convention, weexplicitly specify the definition of the relevant dimensionless variable in the respective section.Typical values of parameters and the characteristic scales which we use to guide our analysisare given in table 3.1.373.2 The modelWe assume that ice is moving as a three-dimensional Stokes flow∂τxx∂x+∂τxy∂y+∂τxz∂z−∂p∂x= 0, (3.1a)∂τxy∂x+∂τyy∂y+∂τyz∂z−∂p∂y= 0, (3.1b)∂τxz∂x+∂τyz∂y+∂τzz∂z−∂p∂z= ρg, (3.1c)where τij are the components of the deviatoric stress tensor, p is the pressure, ρ the density ofice and g is the gravitational constant. The stresses are linked to strain rates Dij through anon-Newtonian constitutive relationDij = Aτn−1τij. (3.2a)n = 3 for typical glaciological ice-flow problems (Paterson, 1994). τ is the effective stressτ 2 = τijτij/2, (3.2b)where the summation convention has been used. A is a temperature dependent parameter.Strain rates are defined asDij =12(∂ui∂xj+∂uj∂xi)(3.2c)with ui the components of the velocity vector u. We use u = (u1, u2, u3) and u = (u, v, w)interchangeably. Mass conservation requires∇ · u =∂u∂x+∂v∂y+∂w∂z= 0. (3.3)Now turning to the boundary conditions for the ice flow, the upper surface at z = s(x, y, t)is assumed to be stress freeτijnj − pni = 0. (3.4)38The outward pointing normal vector n = (nx, ny, nz) is defined asn =(− ∂s∂x ,−∂s∂y , 1)√(∂s∂x)2+(∂s∂y)2+ 1. (3.5)The surface moves according to the kinematic boundary conditionw + a = u∂s∂x+ v∂s∂y+∂s∂t. (3.6)a is the rate of mass gain (through snowfall, if positive) or loss (through melting or sublimation,if negative) at the surface.The ice rests on sediment or bedrock and we assume that bedrock undulations and variationsare small relative to the spatial scales under consideration, so that the base of the ice is at z = 0.This is in line with our goal to model the evolution of topographically unconfined ice streams.Generalization to an arbitrary bed is straightforward.We assume that different boundary conditions apply under the ice stream, where the bed isunfrozen and sliding is possible, and under the ice ridge, where the bed is frozen and no slidingis possible. At the base of the ice stream, |y| < ym, z = 0, we assume that the direction of thetangential component of traction is given by the sliding velocityτijnj − τklnknlni = τbui|u|, (3.7)where we assume |u| > 0. τb is the magnitude of basal friction which we assume to be giventhrough a friction law as a function of the sliding speed |u| and possibly of other variables thatwe do not model here, such as effective pressure N , as described in chapter 2 (see also Iken andBindschadler, 1986; Fowler, 1987; Engelhardt and Kamb, 1997; Iverson et al., 1999; Tulaczyket al., 2000b; Kamb, 2001). We do not specify a particular friction law at this point, but willlater simply assume that τb remains low enough everywhere in the ice stream to facilitate rapidsliding. n is again the upward pointing unit vector, for a flat bedn = (0, 0, 1). (3.8)39We furthermore assume that no penetration of ice into the bed is possible,u · n = 0 (3.9)and there is no basal melting.At the base of the ice ridges we assume that there is no slip,u = 0 (3.10)for |y| > ym, z = 0. The abrupt switch from slip to no slip is known to potentially engendersingularities that may be problematic (see, for example, Hutter and Olunloyo, 1980; Barcilonand MacAyeal, 1993; Moore et al., 2010), and there may be residual sliding even where the bedis frozen (Fowler, 1986; Echelmeyer and Zhongxiang, 1987; Cuffey et al., 1999). We discuss thisfurther in section 4.6.We associate the switch from no slip to slip with a switch in basal temperature, whichwe consequently need to solve for. Within the ice, 0 < z < s(x, y, t), heat is both advectedand conducted, and produced through strain heating, while in the bed, z < 0, there is onlyconduction (assuming negligible shear in the sediment):ρ c(∂T∂t+ u · ∇T)−∇ · [k∇T ] = Dijτij for 0 < z < s(x, y, t), (3.11a)ρbed cbed∂T∂t−∇ · [kbed∇T ] = 0 for z < 0. (3.11b)c and cbed are the specific heat capacities of ice and bed, respectively, ρbed is the density of thesediment, and k and kbed are the thermal conductivities of ice and sediment.We assume a constant temperature T0 at the surface z = s(x, y, t) and a constant incominggeothermal heat flux qgeo:T = T0 at z = s(x, y, t), (3.12a)−kbed∂T∂z→ qgeo as z → −∞. (3.12b)40At the ice-bed interface we again have to distinguish between ice stream and ridge. At the baseof the ice stream, the bed is temperate, while at the base of the ice ridge, temperature is belowthe melting point and the normal component of the heat flux is continuous:T = Tm for |y| < ym, z = 0, (3.13a)T < Tm and − k∂T∂z∣∣∣∣++ kbed∂T∂z∣∣∣∣−= 0 for |y| > ym, z = 0. (3.13b)In the next sections we will non-dimensionalise the model in the three regions previously iden-tified and use the assumption of low aspect ratios to derive simplified leading order approxima-tions.3.3 Ice streamObservational data suggest that there is little vertical variation in the velocity of an ice stream,which resembles the plug flow of an ice shelf (e.g. Engelhardt and Kamb, 1998). This is due torapid sliding at the base of the ice. Models for the plug flow of ice streams and ice shelves havebeen derived by e.g. Muszynski and Birchfield (1987) and MacAyeal (1989). These are basedon the ‘shallow’ nature of ice sheet and ice shelf flow: the aspect ratio between vertical andhorizontal scales is small. We follow these authors, but additionally take the lateral confinementof the ice stream flow into account, which leads to further simplifications.We assume that the ice stream is much longer ([x]s = L ≈ 1000 km) than wide ([y]s = Ws ≈25 km), leading to a plan view aspect ratios =WsL(3.14)of about 0.025. In addition to scales for length and width, we assume that an accumulationrate scale [a] is given, as are the material parameters A, n, ρ and acceleration due to gravity g.As we will see, it is possible to derive from these the appropriate scales for velocities, stressesand ice thickness. Key to developing a shallow ice stream model is then that the scale [z]s forice thickness in the stream is much less than Ws.41Definitions = WsLr = WLδs =[z]sWs=([a]A(ρg)nLn+1r2n+3s) 1n+1Definitionδr =[z]rW = 2n+32n+2s − n+22n+2r δ12sλ ={ [τyz ]BL[τxy ]BL= 2sδ−1s forsδs. 1[τxy ]BL[τyz ]BL= 1+ 1ns δ− 1ns for sδs  1γ = [z]s[z]r = − 12n+2s − n2n+2r δ12sTable 3.2: Aspect ratios used in the matched asymptotic expansions. The three ratios on theleft, s, r and δs can be chosen independently, in which case the remaining three ratios δr, λand γ can be expressed as indicated. Our asymptotic expansions assume not only that s  1,r  1, δs  1 but also that δr  1, λ 1 as well as γ . 1. This imposes additional constraintson the ratios s, r and δs in the form δs  − 2n+3n+1s n+2n+1r , δs . 1n+1s nn+1r , and either s . δs or1+1/ns  δs  s. It is easy to see that a parameter regime like this can at least in principleexist.We assume that pressure scales with its hydrostatic magnitude [p]s = ρg[z]s. A balancebetween the longitudinal pressure gradient and lateral as well as vertical shear stresses leads to[τxy]s/Ws = [τxz]s/[z]s = [p]s/L, while balancing terms in the velocity divergence gives [u]s/L =[v]s/Ws = [w]s/[z]s. Assuming flow dominated by lateral shearing, we put [u]s/Ws = A[τxy]ns ,while the normal stress components scale as A[τxy]n−1[τxx]s = [u]s/[x]s, A[τxy]n−1s [τyy]s = [v]s/[y]sand A[τxy]n−1s [τzz]s = [w]s/[z]s. The vertical shear stress associated with transverse flow balancesextensional stresses in that region, leading to [τyz] = [z]s[τyy]s/Ws. The basal shear stress τb isassumed to scale as the vertical shear stress, [τb]s = [τxz]s. We still need an additional scale foreither velocity or ice thickness in the stream. Taken together with the width of the ice stream,these give us an ice discharge rate [u]s[z]sWs, which we assume to balance net accumulation overthe domain, [u]s[z]sWs = [a]WL. As a temperature scale we naturally choose [T ] = Tm − T0.We can now compute the scale of the ice thickness in the ice stream as[z]s =LWs([a]A(ρg)nWWs) 1n+1. (3.15)This relation reflects mass conservation: For constant temperatures and accumulation rateswider ice streams have to be thinner as they discharge more mass. Wider ice ridges, on theother hand, provide more inflow of mass, therefore they allow for thicker ice streams.42We have already defined the aspect ratio s from the known geometry of the domain and inaddition a second plan aspect ratio can be defined asr =WL. (3.16)From (3.15) it is clear that [z]s is not determined purely by the length scales L, Ws, and W andwe therefore obtain a third, vertical aspect ratioδs =[z]sWs(3.17)that is independent of s and r. With the typical parameter values given in table 3.1, wefind δs = 0.04 and this allows us to develop a thin flow model for the ice stream. In terms ofthe given geometrical parameters Ws, W , L, and the scale [a], the remaining scales for the icestream flow problem are given by[x]s = L, [y]s = Ws, [z]s =1s([a]A(ρg)nrs) 1n+1, (3.18a)[u]s =(rs[a]A1/nρg) nn+1Ws, [v]s =(1/ns r[a]A1/nρg) nn+1 Ws = s[u]s, (3.18b)[w]s =(1/ns r[a]A1/nρg) nn+1 Ws = sδs[u]s, (3.18c)[p]s =1s(ρg[a]Ars) 1n+1= ρg[z]s, [τ ]s = [τxy]s =(ρg[a]Ars) 1n+1= s[p]s, (3.18d)[τxx]s = [τyy]s = [τzz]s = s(ρg[a]Ars) 1n+1= s[τxy]s, (3.18e)[τxz]s = [τb]s =1s([a]2A2(ρg)n−1W n+1s2r2s) 1n+1= δs[τxy]s, (3.18f)[τyz]s =([a]2A2(ρg)n−1W n+1s2r2s) 1n+1= sδs[τxy]s, (3.18g)[t]s =1r(1A(ρg)n[a]nrs) 1n+1, [T ]s = Tm − T0. (3.18h)433.3.1 Leading order ice stream modelWe non-dimensionalise as u = [u]sus, v = [v]svs, etc., and put ym = [y]syM , T = Tm + [T ]sTs.At leading order, omitting terms of O(δ2s) and O(2s), the stress balance is∂τxy,s∂ys+∂τxz,s∂zs−∂ps∂xs= 0, (3.19a)∂ps∂ys= 0, (3.19b)∂ps∂zs= −1. (3.19c)The effective stress only depends on the lateral shear stress τxy,sτs = |τxy,s| , (3.20)where we have again omitted terms of O(δ2s) and O(2s). The constitutive relations also simplifysignificantly:12∂us∂ys= |τxy,s|n−1 τxy,s,∂us∂zs= 0,∂vs∂zs= 0,∂us∂xs= |τxy,s|n−1 τxx,s,∂vs∂ys= |τxy,s|n−1 τyy,s,∂ws∂zs= |τxy,s|n−1 τzz,s.(3.21)As in the ice stream model derived by MacAyeal (1989), the vertical gradients of both hori-zontal velocities are negligible. However, our model deviates from MacAyeal’s original ‘shallowshelf approximation’ in that we also neglect the gradients ∂v/∂x in downstream direction, aconsequence of the small plan aspect ratio s = Ws/L (the traditional approach is to assumethat [z] [y] = [x]). The equation of mass conservation keeps its form in the scaled variables∂us∂xs+∂vs∂ys+∂ws∂zs= 0. (3.22)The leading order surface boundary conditions at zs = ss(x, y, t) are (omitting terms ofO(2s))∂ss∂ys= 0, ps = 0, τxz,s = 0. (3.23a)44Definitionαs = 2A[τxy ]n+1s [z]2sk[T ] = 2ρg[a]L2k[T ] srδ2sνs =qgeo[z]skbed[T ]= qgeoLkbed[T ]sδsPes =[a]ρcLk rδsΓ = ρbedρcbedcκ = kbedkDefinitionαr = 2A[τyz ]n+1r [z]2rk[T ] = n+2n+1s −1n+1r δ−1s αsαBL = 2A[τ ]n+1BL [z]2BLk[T ]BL={αs for sδs . 11+1/ns δ−1−1/ns αs for sδs  1νr =qgeo[z]rkbed[T ]= 12n+2s n2n+2r δ− 12s νsPer =ρc[a][z]rk = δrδ−1s Pes = 2n+32n+2s − n+22n+2r δ−1/2s PesPeBL =ρc[a]Wk = δ−1s PesTable 3.3: Dimensionless groups. αs, Pes and νs are independent dimensionless groups in themodel along with s, r and δs. The remaining model parameters can be expressed throughthese, as indicated on the right.The first of these equalities shows that the surface is flat in the across-stream direction. Theevolution of the surface is determined byws +sra = us∂ss∂xs+ vs∂ss∂ys+∂ss∂ts= us∂ss∂xs+∂ss∂ts(3.23b)on account of (3.23a). The scaled boundary conditions at the base z = 0 are:τxz,s = τb,sus|us|, τyz,s = τb,svs|us|, ws = 0, (3.24)where we have restricted ourselves to |ys| < yM .At leading order (again omitting terms of O(δ2s)) the heat equation becomesPes(∂Ts∂ts+ us∂Ts∂xs+ vs∂Ts∂ys+ ws∂Ts∂zs)−∂2Ts∂z2s= αs |τxy,s|n+1 for 0 < zs ≤ ss(xs), (3.25a)PesΓ∂Ts∂ts− κ∂2Ts∂z2s= 0 for zs ≤ 0, (3.25b)where αs = 2ρg[a]r[z]2s/(sAk[T ]s) is the ratio between dissipative heat production and con-ductive cooling. The Peclet number measures the relative magnitude of heat advection toconduction and is defined as Pes = ρc[a]Wδs/k. With the parameter values in table 3.1, we getαs ≈ 6 and Pes ≈ 14. Note that a summary of dimensionless groups may be found in table 3.245and table 3.3.The scaled boundary conditions for the heat equation areTs = −1 at zs = ss(xs), (3.26a)Ts = 0 at zs = 0, (3.26b)∂Ts∂zs→ −νs as zs → −∞, (3.26c)where νs = [z]sqgeo/kbed[T ]s is a scaled geothermal heat flux, with νs ≈ 0.8 using the valuesfrom table 3.1. In writing down (3.25) we have again restricted ourselves to |y| < yM . Note thatPes, αs, νs are independent of each other and of s, r and δs, as they depend on the additionalconstants k, c and qgeo.3.3.2 IntegrationWe can immediately integrate (3.19c) using (3.23a)2 to show that ps = (ss− z). From (3.21) wesee that us depends on xs, ys and t only. Integrating (3.19a) from zs = 0 to zs = ss and using(3.23a)3 and (3.24)1 yields121/n∂∂ys[ss∣∣∣∣∂us∂ys∣∣∣∣1−nn ∂us∂ys]− τb,sus|us|= ss∂ss∂x, (3.27)where we have also used (3.23a)1 which shows that ss depends on xs and ts only. Because weassume a flat bed at zs = 0, ss is both the location of the upper surface and the ice thickness.An evolution equation of ss can be derived by integrating (3.22) from zs = 0 to zs = ss, andfrom ys = −yM to ys = yM , and using (3.24)3:∫ ss0∫ yM−yM∂us∂xsdys dzs +∫ yM−yMw|zs=ss dys +∫ zs0v|yM − v|−yM dzs = 0. (3.28)Using the fact that us is independent of zs and also vs is independent of zs, see (3.21), and usingthe kinematic boundary condition (3.23b) we can write:∫ yM−yMss∂us∂xsdys +∫ yM−yM(∂ss∂t+ us∂ss∂xs−sra)dys + 2v|yMss = 0, (3.29)46where we have also used v|yM = −v|−yM as the geometry is symmetric. Using the fact thatss is independent of ys, using the symmetry of the system and manipulating the derivativesaccording to Leibniz’ rule, we can rewrite (3.29) as∂∂ts(2yMss)− 2ss∂yM∂ts+∂∂xs∫ yM−yMssus dy − 2ssus|yM∂yM∂xs=sr∫ yM−yMa dy − 2v|yMss (3.30)or in more succinct form∂ (2yMss)∂ts+∂Q∂xs=sr∫ yM−yMa dy + qin (3.31a)whereQ(xs, ts) =∫ yM (xs,ts)−yM (xs,ts)ss(xs, ts)us(xs, ys, ts) dys (3.31b)andqin = 2ss(∂ym∂ts− vs(xs, yM , ts))+ 2ssus|yM∂yM∂xs. (3.31c)(3.27) together with (3.31) is the leading order model for the ice stream. In these equations,the necessary boundary conditions on us in (3.27) as well as qin and any changes in the marginposition yM in (3.31) represent the coupling between the ice stream and the surrounding regionsof the ice sheet. Note that the ratio s/r appears in (3.31a); it equals Ws/W and is alwaysless than 1. s/r is small if the ice stream width is much less than the domain width, in whichcase the first term on the right hand side of (3.31a), representing surface mass balance of theice stream, becomes negligible compared with the second term, which represents mass exchangebetween ice stream and the surrounding parts of the ice sheet.Determining qin and changes in yM forces us to consider the flow of the ice ridges and of theboundary layer separating the ridges from the stream.3.4 Ice ridgeIce ridges are regions of slowly flowing ice between ice streams. In our model we assume thatthe ice is frozen to the bed and that there is no slip. Ice therefore moves by vertical shear alonewith surface velocities that are much slower than in the adjacent ice stream. We also assume47that ice is thicker in the ice ridge than in the ice stream and the resulting gradient in surfaceelevation leads to a gravity driven inflow of ice into the ice stream.The dimensions of the ice ridge are as shown in figure 3.1. We restrict ourselves to the iceridge to the left at y < −ym (looking upstream). Consequently, we rescale the lateral coordinatey = [y]ryr−W with [y]r = W . The length scale [x]r = [x]s = L is the same as for the ice stream.Typical values are W = 50 km and L = 1000 km, with the aspect ratio r = W/L taking atypical value around 0.05. We also assume that the ice ridge is shallow, with its thickness muchsmaller than its horizontal extent.If the dominant motion in the ice ridge is through vertical shear, then vertical shear stressbalances the driving pressure gradient, so [τxz]r/[z]r = [p]r/L and [τyz]r/[z]r = [p]r/W , where[pr] = ρg[z]r. Flow is dominantly towards the ice stream, so [v]r/[z]r = A[τyz]nr while [u]r/[z]r =A[τyz]n−1r [τxz]r. Mass conservation (3.3) suggests [v]r/W = [w]r/[z]r with [w]r = [a] where[a] is the scale of the surface accumulation. Scales for the remaining stress components are[u]r/[x]r = A[τyz]n−1r [τxx]r, [v]r/[y]r = A[τyz]n−1r [τyy]r, [w]r/[z]r = A[τyz]n−1r [τzz]r, and [u]r/[y]r =A[τyz]n−1r [τxy]r. Lastly we assume that the ice stream sets the dominant evolution time scaleand put [t]r = [t]s as determined in the previous section. In terms of the known dimensions W ,L and the imposed accumulation rate scale [a], the remaining scales can be expressed as:[x]r = L, [y]r = W, [z]r =([a]W (n+1)A(ρg)n) 12n+2(3.32a)[u]r = r(A(ρg)nW n+1[a]2n+1) 12(n+1) (3.32b)[v]r =(A(ρg)nW n+1[a]2n+1) 12(n+1) = −1r [u]r, [w]r = [a] = δr[v]r = δr−1r [u]r, (3.32c)p =([a](ρg)n+2W n+1A) 12(n+1)= ρg[z]r, [τ ]r = [τyz]r =([a]ρgA) 1n+1= δr[p]r (3.32d)48[τxx]r = 2r([a]3A3(ρg)n−2W n+1) 12(n+1)= 2rδ2r [p]r, (3.32e)[τyy]r = [τzz]r =([a]3A3(ρg)n−2W n+1) 12(n+1)= δ2r [p]r, (3.32f)[τxy]r = r([a]3A3(ρg)n−2W n+1) 12(n+1)= rδ2r [p]r, [τxz]r = r ([a]ρg/A)1n+1 = rδr[p]r, (3.32g)[t]r = [t]s =1r(1A(ρg)n[a]nrs) 1n+1, [T ]r = [T ]s = Tm − T0. (3.32h)The ice ridge thickness [z]r allows us to define an aspect ratio for the ice ridge, δr = [z]r/W .Our leading order model below is based on δr  1, as well as r  1. Note however that δr isnot independent of the previously defined parameters s, r, δs, but can be expressed asδr = 2n+32n+2s − n+22n+2r δ12s . (3.33)We are therefore not only assuming that s, r, δs are small, but also that2n+3n+2s δn+1n+2s  r  1. (3.34)Such a parameter regime however clearly exists, and with the scales in table 3.1 we get δr = Leading order ice ridge modelWe put u = [u]ur, v = [v]vr, etc., and T = [T ]Tr + Tm, y = [y]r(yr − 1). With these, the stressbalance at leading order is∂τxz,r∂zr−∂pr∂xr=0, (3.35a)∂τyz,r∂zr−∂pr∂yr=0, (3.35b)−∂pr∂zr=1. (3.35c)49Here, terms of O(δ2r) were omitted. At leading order, the effective stress in our analysis is solelya function of the vertical shear stress τyz,r, i.e.,τr = |τyz,r| , (3.36)so the stress components are related to the velocity gradients at leading order through12(∂ur∂yr+∂vr∂xr)= |τyz,r|n−1 τxy,r,12∂ur∂zr= |τyz,r|n−1 τxz,r,12∂vr∂zr= |τyz,r|n−1 τyz,r,∂ur∂xr= |τyz,r|n−1 τxx,r,∂vr∂yr= |τyz,r|n−1 τyy,r,∂wr∂zr= |τyz,r|n−1 τzz,r.(3.37)Omitting terms of O(2r), mass conservation becomes∂vr∂yr+∂wr∂zr= 0. (3.38)The boundary conditions at the upper surface zr = sr areτxz,r = 0, τyz,r = 0, pr = 0, wr + a = vr∂sr∂yr+2rδr2sδs∂sr∂ts. (3.39a)The boundary conditions at the base of the ice areur = vr = wr = 0, (3.40a)where we have restricted ourselves to the region of no slip, yr < 1− (s/r)yM .The rescaled heat equation at leading order (omitting terms of O(δ2r), O(2r)) isPer(2rδr2sδs∂Tr∂ts+ vr∂Tr∂yr+ wr∂Tr∂zr)−∂2Tr∂z2r= αr|τyz,r|n+1 for 0 < zr ≤ sr, (3.41a)PerΓ2rδr2sδs∂Tr∂ts− κ∂2Tr∂z2r= 0 for zr ≤ 0. (3.41b)The ice ridge Peclet number is Per = ρc[a][z]r/k and αr = 2A[τyz]n+1r [z]2r/(k[T ]r) is a measure50for the magnitude of heat production relative to conductive cooling. Note that as with δr, Perand αr are not independent but can be expressed in terms of αs and Pes as well as s, r andδs as Per = (2n+3)/(2n+2)s −(n+2)/(2n+2)r δ−1/2s Pes and αr = (n+2)/(n+1)s −1/(n+1)r δ−1s αs. With thevalues in table 3.1, we get Per ≈ 7 and αr ≈ 4. Boundary conditions areTr = −1 at zr = sr, (3.42a)Tr < 0 and∂Tr∂zr∣∣∣∣+− κ∂Tr∂zr∣∣∣∣−= 0 at zr = 0, (3.42b)∂Tr∂zr→ −νr as zr → −∞, (3.42c)with νr = [z]rqgeo/kbed[T ] = 1/(2n+2)s n/(2n+2)r δ−1/2s νs. With the values in table 3.1 we getνr ≈ Depth integrationThe ice flow model (3.35a)–(3.40) takes the form of a standard ‘shallow ice’ model for flow inthe yr–direction (Fowler and Larson, 1978; Morland and Johnson, 1980; Hutter, 1983), and canbe integrated straightforwardly to show thatur = −2sn+1rn+ 1[1−(1−zrsr)n+1] ∣∣∣∣∂sr∂yr∣∣∣∣n−1 ∂sr∂x, (3.43a)vr = −2sn+1rn+ 1[1−(1−zrsr)n+1] ∣∣∣∣∂sr∂yr∣∣∣∣n−1 ∂sr∂y, (3.43b)from (3.35). The evolution of the ice ridge surface is governed by2rδr2sδs∂sr∂ts+∂qr∂yr= a, (3.44)whereqr =∫ sr0vr dzr = −2n+ 2sn+2r∣∣∣∣∂sr∂yr∣∣∣∣n−1 ∂sr∂yr. (3.45)(3.44) together with (3.45) is a non-linear diffusion problem which requires a boundary conditionat the ice stream margin yr = 1− (s/r)yM . This requires us to couple the ice stream problem51to the ice ridge problem, which leads us to consider a boundary layer separating the two.3.5 Boundary layer scalingIn the last two sections we have derived equations that describe the evolution of an ice stream,(3.27)–(3.31), and of an ice ridge adjacent to it, (3.44)–(3.45). The ice stream moves as a plugflow through a ‘channel’ formed by the ice ridges, where ice moves through vertical shear towardsthe ice stream. At present we can only speculate about the appropriate coupling between thetwo models. For instance, we might expect the inflow qin in (3.31c) to be related to the iceridge flux qr in (3.45) at the margin. However, due to the distinct nature of the velocity fieldsin stream and ridge – one is a plug flow, the other a shearing flow – coupling is not trivial anda boundary layer is required. Similar boundary layers are encountered in other problems inglaciology, such as marine ice sheets (Chugunov and Wilchinsky, 1996).The boundary layer is located between the ice stream and the ice ridge in the region aroundy = ±ym(x, t), which we refer to as the ice stream margin. As with the ice ridge, we focuson the margin to the left of the ice stream near y = −ym. Within the boundary layer thespatially rapid transition from the Poiseuille-type flow in the ice ridge to the plug flow in the icestream has to take place. The boundary layer cannot be shallow – otherwise we would expectone kind of thin film approximation or the other to hold – and to account for continuity of thesurface with the adjacent stream, we choose the scales [y]BL = [z]BL = [z]s. We rescale thelateral coordinate y = [y]BLY − [y]syM(xs, ts), so that the boundary layer coordinate is centeredaround the transition between a frozen bed and an unfrozen temperate bed that occurs in theright hand margin at Y = 0.To account for mass flux from the ridge through the margin into the ice stream, we put[v]BL[z]BL = [v]s[z]s = [a]W , or [v]BL = [v]s, while mass balance requires [v]BL = [w]BL. Theice stream imposes a lateral shear stress on the boundary layer, which must be balanced byvertical shear stresses in the boundary layer and leads us to put [τxz]BL = [τxy]BL = [τxy]s.Shearing due to transverse flow in the boundary layer leads to [v]BL/[z]BL = A[τ ]n−1BL [τyz]BL.There are now two possible cases: the second stress invariant τ in the boundary layer canbe dominated either by shear stresses due to transverse flow, so that [τ ]BL = [τyz]BL, or by52lateral shear stresses, so that [τ ]BL = [τxy]BL. In the first case, the stress invariant becomesT =√T 2xy + T 2xz + (s/δs)2 T 2yz + ..., and hence this scaling is only appropriate for s/δs . 1.Similarly, the second version is only appropriate for s/δs & 1, but to be definite when s ∼ δs,we define [τ ]BL = [τyz]BL only when s/δs  1, and put [τ ]BL = [τxy]BL otherwise.The scale of the downstream velocity is defined through [u]BL/[z]BL = A[τ ]n−1BL [τxy]BL, andthe normal components of the stress are determined through [u]BL/[x]BL = A[τ ]n−1BL [τxx]BL,[v]BL/[y]BL = A[τ ]n−1BL [τyy]BL and [w]BL/[z]BL = A[τ ]n−1BL [τzz]BL. We assume that the basal dragin the unfrozen part of the margin is comparable with basal drag in the ice stream, which forinstance would be appropriate if the bed of the margin is hydraulically well connected to theice stream bed, and we therefore put [τb]BL = [τb]s.In terms of the imposed geometrical parameters Ws, W , L, and the accumulation scale [a],this gives[x]BL = L, [y]BL = [z]BL = [z]s =LWs([a]A(ρg)nWWs) 1n+1=1s([a]A(ρg)nrs) 1n+1, (3.46a)[u]BL =WLW 2s[a] = r2s [a] forsδs. 1,r2/ns([a]A(ρg)nWWn+2s) 1−n1+n1n= rδ−1s λ−1 for sδs  1,(3.46b)[v]BL = [w]BL =WL(A(ρg)n[a]nW n+2sW)1/(n+1)= λ−1s [u]BL, (3.46c)[p]BL =LWs(ρg[a]AWWs) 1n+1= ρg[z]BL, (3.46d)[τ ]BL =(ρg[a]AWWs) 1n+1= s[p]BL for sδs . 1,(WsL)2/n(A(ρg)n[a]W21−nsW) 1−n1+n1n= λ[p]BL for sδs  1,(3.46e)[τxx]BL =([a]A(ρg)(n−1)/2WW (n+3)/2s)2/(n+1)= 2s[p]BL, (3.46f)[τyy]BL = [τzz]BL = [τyz]BL = λ[p]BL =ρgW 3s L−2 for sδs . 1,(WsL)2/n(A(ρg)n[a]W21−nsW) 1−n1+n1nfor sδs  1,(3.46g)53[τxy]BL = [τxz]BL =(ρg[a]AWWs) 1n+1= s[p]BL, (3.46h)[τb]BL = [τb]s = δs[τxy]BL =LWs([a]2A2(ρg)n−1W n+1sW 2W 2s) 1n+1(3.46i)[t]BL = [t]s = [t]r =LW(1A(ρg)n[a]nWWs) 1n+1, [T ]BL = [T ]s = [T ]r = Tm − T0. (3.46j)The leading order models for the two cases in s/δs are almost identical and we can treatthem both succinctly by defining a parameter λ throughλ =2s/δs forsδs. 1,(n+1s /δs)1/n for sδs  1.(3.47)It will be essential that λ  1 and 2s/λ  1, which is assured when s/δs . 1, while fors/δs  1, it leads to an additional constraint n+1s  δs  s. This is in addition to theconstraint (2n+3)/(n+2)s δ(n+1)/(n+2)s  r  1 in (3.34).3.5.1 Boundary layer modelWe put u = [u]BLU , v = [v]BLV , w = [w]BLW , τxy = [τxy]BLTxy, etc., and T = [T ]BLTBL + Tm,y = [y]BLY −[y]syM(xs, ts). In the dimensionless variables, we give the boundary layer equationsup to the order indicated. Momentum conservation becomes∂Txy∂Y+∂Txz∂Z−δs∂P∂X−∂yM∂xs∂P∂Y=0, (3.48a)λ∂Tyy∂Y+λ∂Tyz∂Z−∂P∂Y=0, (3.48b)λ∂Tyz∂Y+λ∂Tzz∂Z−∂P∂Z=1, (3.48c)omitting terms of O(2s). The effective stress to O(2s) is:T =√T 2xy + T 2xz +2sδ2sT 2yz +122sδ2sT 2yy +122sδ2sT 2zz forsδs. 1,√T 2yz +12T2yy +12T2zz forsδs 1.(3.49)54The stress-strain relationships in rescaled form are now:12(∂U∂Y+ λ∂yM∂xs∂V∂Y)= T n−1Txy,12(∂U∂Z+ λ∂yM∂xs∂W∂Y)= T n−1Txz,12(∂V∂Z+∂W∂Y)= T n−1Tyz,∂yM∂xs∂U∂Y+ δs∂U∂X= T n−1Txx,∂V∂Y= T n−1Tyy,∂W∂Z= T n−1Tzz,(3.50)omitting terms of O(2s) and O(δsλ). To an error of O(2sδs/λ), mass conservation can beexpressed as2sλ∂yM∂xs∂U∂Y+∂V∂Y+∂W∂Z= 0. (3.51)The boundary conditions at the upper surface at Z = S are, omitting terms of O(2sδs)− Txy∂S∂Y+ Txz+δsP∂S∂X+ P∂yM∂xs∂S∂Y=0, (3.52a)−λTyy∂S∂Y+λTyz+ P∂S∂Y=0, (3.52b)−λTyz∂S∂Y+λTzz− P =0, (3.52c)and, omitting terms of O(δ2s)W +srδsa = V∂S∂Y+∂yM∂ts∂S∂Y+2sλU∂yM∂xs∂S∂Y+ δs∂S∂ts. (3.53)At the base of the ice at Z = 0 the boundary conditions are, omitting terms of O(δs)Txz = 0 Tyz = 0 W = 0 for Y > 0, (3.54a)U = V = W = 0 for Y < 0. (3.54b)55Finally, the heat equation in rescaled form isPeBL(δs∂TBL∂ts+∂yM∂ts∂TBL∂Y+2sλU∂yM∂xs∂TBL∂Y+ V∂TBL∂Y+W∂TBL∂Z)−(∂2TBL∂Y 2+∂2TBL∂Z2)= αBLTn+1 for 0 < Z < S,(3.55a)PeBLΓ(δs∂TBL∂ts+∂yM∂ts∂TBL∂Y)− κ(∂2TBL∂Y 2+∂2TBL∂Z2)= 0 for Z < 0,(3.55b)omitting terms of O(δ2s), with PeBL = ρc[a]W/k and αBL = 2A[τBL]n+1[z]2BL/(k[T ]BL). As withPer and αr, PeBL and αBL are not independent of Pes, αs, s, r and δs, but can be written inthe formPeBL = δ−1s Pes, αBL =αs for sδs . 11+1/ns δ−1−1/ns αs for sδs  1.(3.56)PeBL ≈ 380 and αBL ≈ 6 with the typical values in table 3.1. The boundary conditions are thesame for both casesTBL = −1 at Z = S, (3.57a)−∂TBL∂Z= νs as Z → −∞, (3.57b)TBL < 0 and −∂TBL∂Z∣∣∣∣++ κ∂TBL∂Z∣∣∣∣−= 0 for Y < 0, Z = 0, (3.57c)TBL = 0 for Y > 0, Z = 0. (3.57d)3.5.2 Series expansionOur primary objective is to find leading order solutions for velocities U , V and W , surfaceelevation S and temperature TBL in the boundary layer. To do so we need to develop an56asymptotic expansion in λ:P (X, Y, Z, ts) = P(0)(X, Y, Z, ts) + λP(1)(X, Y, Z, ts) + o(λ),S(X, Y, ts) = S(0)(X, Y, ts) + λS(1)(X, Y, ts) + o(λ),U(X, Y, Z, ts) = U(0)(X, Y, Z, ts) + o(1), V (X, Y, Z, ts) = V(0)(X, Y, Z, ts) + o(1),W (X, Y, Z, ts) = W(0)(X, Y, Z, ts) + o(1).To evaluate the pressure at the surface Z = S = S(0) + λS(1) + ..., we can write P as a TaylorexpansionP (X, Y, S, ts) = P(0)(X, Y, S(0), ts) + λ∂P (0)∂ZS(1)(X, Y, ts) + λP(1)(X, Y, S(0), ts) + o(λ). (3.59)At zeroth order we get from (3.48b)–(3.48c)−∂P (0)∂Y= 0, −∂P (0)∂Z= 1, (3.60)with boundary condition P (0) = 0 at Z = S(0) from (3.52c). Integration of (3.60)2 givesP (0) = S(0) − Z and consequently (3.60)1 requires ∂S(0)/∂Y = 0. The upper surface is flat atlowest order, and the ice flow domain is a parallel-sided strip.Combining (3.48) and (3.50), the leading order problem for the downstream flow is∂∂Y(η∂U (0)∂Y)+∂∂Z(η∂U (0)∂Z)= 0 (3.61)where we have used the definition of the viscosityη =121/n[∣∣∣∂U(0)∂Y∣∣∣2+∣∣∣∂U(0)∂Z∣∣∣2+ 2sδ2s(∣∣∣∂V (0)∂Z +∂W (0)∂Y∣∣∣2+ 2∣∣∣∂V (0)∂Y∣∣∣2+ 2∣∣∣∂W (0)∂Z∣∣∣2)] 1−n2nfor sδs . 1,[∣∣∣∂V (0)∂Z +∂W (0)∂Y∣∣∣2+ 2∣∣∣∂V (0)∂Y∣∣∣2+ 2∣∣∣∂W (0)∂Z∣∣∣2] 1−n2nfor sδs  1.(3.62)57Boundary conditions at the ice surface Z = S(0)(X, ts) and base Z = 0 areη∂U (0)∂Z= 0 at Z = S(0), (3.63a)η∂U (0)∂Z= 0 at Z = 0, Y > 0, (3.63b)U (0) = 0 at Z = 0, Y < 0. (3.63c)The across-flow velocities in the (y, z)–direction perpendicular to the direction of bulk ice streamflow are solutions of∂∂Y(2η∂V (0)∂Y)+∂∂Z[η(∂V (0)∂Z+∂W (0)∂Y)]−∂P (1)∂Y= 0, (3.64a)∂∂Y[η(∂V (0)∂Z+∂W (0)∂Y)]+∂∂Z(2η∂V (0)∂Y)−∂P (1)∂Z= 0, (3.64b)∂V (0)∂Y+∂W (0)∂Z= 0, (3.64c)subject to the boundary conditionsW (0) = η(∂V (0)∂Z+∂W (0)∂Y)= 0, 2η∂W (0)∂Z− P (1) + S(1) = 0 at Z = S(0) (3.65)at the ice surface. Boundary conditions at the base areη(∂V (0)∂Z+∂W (0)∂Y)= W (0) = 0 for Y > 0, Z = 0, (3.66a)V (0) = W (0) = 0 for Y < 0, Z = 0. (3.66b)For n = 1, the transverse flow model (3.64)–(3.66) is similar to the model of an ice stream flow-ing across a no slip–free slip boundary in Barcilon and MacAyeal (1993), while the downstreamflow model (3.61)–(3.63) is the same as the model in Schoof (2004). We leave the temperatureproblem unchanged from (3.55) for now, deferring possible considerations of large Peclet num-bers PeBL until later. The domain however now simplifies as the upper surface is known to beflat.583.6 MatchingThe boundary layer solutions have to match the solutions in the outer regions, the ice streamand ridge. Even though we have motivated all our scales physically, the spatial coordinates anddependent variables used to describe stream, ridge and boundary layer can be related purelythrough rescalings using the independent aspect ratios s, r, δs, and the derived dimensionlessgroups δr and λ defined in terms of s, r and δs in (3.33) and (3.47), respectively. Theserescalings take the formxs = xr = X, ys =rs(yr − 1) = δsY − yM(xs, ts), zs =sδsrδrzr = Z, (3.67a)us =2sδsδrur =2sλU, vs =sδsrδrvr = V, ws =srwr =1δsW, (3.67b)ps =rδrsδspr = P, τs =rδ2r2sδsτr =T for sδs . 1,λsT for sδs  1,τxx,s =3rδ3rsδsτxx,r = Txx, (3.67c)τyy,s =rδ3r3sδsτyy,r =λ2sTyy, τzz,s =rδ3r3sδsτzz,r =λ2sTzz, τxy,s =2rδ3r2sδsτxy,r = Txy, (3.67d)τxz,s =2rδ2r2sδ2sτxz,r =rδrsδsTxz, τyz,s =rδ2r3sδ2sτyz,r =λ2sδsTyz, Ts = Tr = TBL. (3.67e)3.6.1 Matching the stream with the boundary layerThe matching region between stream and boundary layer corresponds to ys → −yM(xs, ts)and Y → ∞. We have shown that ss and S are independent at leading order of ys and Y ,respectively. Matching surface elevation therefore requires thatS(0)(xs, ts) = ss(xs, ts), (3.68)which we may understand to mean that the ice thickness in the boundary layer is dictated bythe thickness in the ice stream. Matching ice velocities in the matching region, we get us ∼ δsU59for s/δs . 1, us ∼ s(δs/s)1/nU for s/δs  1, vs ∼ V , and ws ∼ δ−1s W , or at leading orderus|y=−yM = 0, (3.69a)vs|y=−yM = limY→∞V (0). (3.69b)Because vs = vs(xs, ys, ts) is independent of zs, the second equation also implies that the velocityV (0) has no vertical shear as Y →∞, soW (0) → 0∂V (0)∂Z → 0for Y →∞. (3.70)Matching lateral shear stresses in the boundary layer requires τxy,s ∼ Txy, orlimY→∞η∂U (0)∂Y= τxy,s|ys=−yM =121/n∣∣∣∣∂us∂ys∣∣∣∣1/n−1 ∂us∂ys∣∣∣∣ys=−yM. (3.71)Matching temperatures and heat fluxes leads to TBL ∼ Ts, δs∂Ts/∂ys ∼ ∂TBL/∂Y , or at leadingorderTs|ys=−yM = limY→∞TBL, (3.72a)∂TBL∂Y→ 0 as Y →∞. (3.72b)These matching conditions provide boundary conditions on the stream and boundary layerproblems and we can understand them in terms of the equations to which they most readilyapply. For instance (3.69a) closes the ice stream flow problem (3.27), while (3.69b) couplesmass conservation of the stream (3.31c) to the boundary layer. Conversely, (3.71) providesboundary conditions for the longitudinal Stokes flow problem in the boundary layer (3.61),while (3.70) provides the necessary boundary conditions for the transverse Stokes flow problem(3.64). Assuming there is net inflow of ice into the ice stream, −vs(xs,−yM) − ∂yM/∂ts < 0,(3.72b) helps to close the heat equation for the boundary layer (3.55), while (3.72a) acts as aboundary condition for the temperature in the stream.603.6.2 Matching the ridge with the boundary layerThe matching region between ridge and boundary layer is at Y → −∞, yr → 1 − (s/r)yM .Recall that s/r by definition is the ice stream half width divided by the half width of thedomain, and must therefore be less than 2. There is the possibility that s/r is small, in whichcase the matching region at leading order is at a fixed position yr → 1, which however does notaffect the matching procedure below.When matching boundary layer and ice ridge, the ratioγ =[z]s[z]r=sδsrδr(3.73)turns out to be of central importance. As a last constraint on our independent parameters s,r and δs we will assume in the procedure below that γ . 1, meaning that the natural ice ridgethickness is not less than the natural ice stream thickness. This constraint can be expressed asδs . 1n+1s nn+1r . (3.74)This constraint is in addition to (2n+3)/(n+2)s δ(n+1)/(n+2)s  r  1 from (3.34) and eithers . δs  1 or n+1s  δs  s  1 from (3.47) and it is easy to see that these constraints arenot mutually incompatible. With the typical parameter values from table 3.1, γ ≈ 0.9, so ourassumption of γ . 1 is reasonable. That being said, it is possible in principle to deal with thecase γ  1 through a rescaling of the ice ridge equations which we consider in appendix B.The surface elevation in the boundary layer is flat at leading order, so that matching surfaceelevations leads tosr(xs, 1− (s/r)yM , ts) = γS(0)(xs, ts) = γss(xs, ts) (3.75)using (3.68). It is also possible to match surface slopes, which however requires the next ordercorrection S(1).Matching the velocities we either have U ∼ 2sδ−1r ur for s/δs . 1 or U ∼ δsδ−1r s(s/δs)1/nurfor s/δs  1, V ∼ γvr, and W ∼ γδrwr. To simplify the matching conditions for the transverse61components V and vr, we can combine (3.43b) and (3.45) intovr =(n+ 2)(n+ 1)qrsr[1−(1−zrsr)n+1]. (3.76)With sr ∼ γS, zr ∼ γZ, this means that γvr ∼ V in the matching region can be expressed aslimY→−∞V (0) =(n+ 2)(n+ 1)qr|yr=1−(s/r)yMS(0)(xs, ts)[1−(1−ZS(0)(xs, ts))n+1]. (3.77)The other two velocity conditions in the matching region becomeU (0) → 0, W (0) → 0 as Y → −∞. (3.78)(3.75) shows that the boundary value for ice thickness in the ice ridge is set by the icestream and provides one instance of direct coupling between stream and ridge. We can furthershow that inflow into the stream is directly given by the value of ice flux into the stream atthe boundary. We can integrate the continuity equation for the boundary layer (3.64c) to showthat0 =∫ Z=S(0)Z=0∂V (0)∂YdZ + W |Z=S(0)Z=0 =∫ Z=S(0)Z=0∂V (0)∂YdZ (3.79)using the boundary conditions W (0) = 0 at Z = 0 and Z = S(0), (3.65)1 and (3.66). Furtherintegration over the width of the boundary layer, assuming that all the relevant limits commute,yields0 =∫ ∞−∞(∫ Z=S(0)Z=0∂V (0)∂YdZ)dY =∫ Z=S(0)Z=0(∫ ∞−∞∂V (0)∂YdY)dZ=∫ Z=S(0)Z=0limY→∞V (0) dZ −∫ Z=S(0)Z=0limY→−∞V (0) dZ=∫ Z=S(0)Z=0vs|y=−yM dZ −∫ Z=S(0)Z=0(n+ 2)(n+ 1)qr|yr=1−(s/r)yMS(0)(xs, ts)[1−(1−ZS(0)(xs, ts))n+1]dZ= vs|y=−yMss − qr|yr=1−(s/r)yM (3.80a)where we used (3.60), (3.69b), (3.77), and (3.21). Together with (3.69a), this shows that the62inflow of ice into the stream qin in (3.31c) can be written in terms of the discharge of ice fromthe ice ridge qr|yr=1−s/ryM and the margin migration rate ∂yM/∂ts asqin = 2ss∂yM∂ts+ 2qr|yr=1−(s/r)yM . (3.81)Finally, we need to match temperatures,limY→−∞TBL(X, Y, Z, ts) = Tr(xr, 1− (s/r)yM , zr, ts). (3.82)This acts as an inflow boundary condition for the heat equation (3.55).3.7 DiscussionIn section 3.3 we have derived a mixed hyperbolic-elliptic stream model, (3.31a)–(3.31b) with(3.27), which is closed by (3.81) and (3.69a). The model for the ice ridge, derived in section3.4 and described by (3.44) and (3.45), is a diffusive shallow ice model and closed by (3.75).The evolution of the ice ridge, (3.44), also depends on the ratio (r/s)γ. With γ ∼ O(1) andr/s > 2, the time scale for changes in the ice ridge is therefore of the same order or slowerthan the time scale for changes in the ice stream. For r/s = W/Ws  1, the ice ridge is in aquasi-steady state.If yM was fixed, the outer problems could be solved without further recourse to the boundarylayer model: the inflow of ice into the stream is set through the ice ridge, while the ice thicknessof the stream provides the necessary boundary condition for the ice ridge. In order to find∂yM/∂ts we need to solve the boundary layer problem, in particular we need to solve for thetemperature TBL.Our scaling suggests PeBL  1, but to facilitate our analysis and for comparison with priorwork (Schoof, 2012), we start by treating Pe as an O(1) parameter and simultaneously assumeαBL ∼ O(1). The case PeBL  1 will be considered separately in section 4.5.The ratio s/δs = [τyz]BL/[τxy]BL = [Q]r/[Q]s describes the relative magnitude of verticaland lateral shear stress in the boundary layer, which equals the ratio between the scale for influx63of ice from the ridge, [Q]r = [v]r[z]r, and the scale of downstream discharge [Q]s = A[τxy]ns [z]2s.With the values in table 3.1 we obtain s/δs ≈ 0.7, but the Siple Coast ice streams drain abigger region than indicated by our choice of L = 1000 km, leading to even smaller values ofs/δs. We therefore limit ourselves to simulations for s/δs  1.The limit of small s/δs should in principle let us simplify viscosity and heat productionrates asη =121/n(∣∣∣∣∂u∂y∣∣∣∣2+∣∣∣∣∂u∂z∣∣∣∣2)1/n−1and a =121/n(∣∣∣∣∂u∂y∣∣∣∣2+∣∣∣∣∂u∂z∣∣∣∣2)1/n+1. (3.83)However, in the case of the viscosity, omitting the O(2s/δ2s) terms due to lateral inflow isproblematic as the derivatives of u vanish in the far field as Y → −∞. With n > 1, this leadsto an infinite viscosity. In principle this can be fixed through the introduction of yet another(rather trivial and passive) boundary layer in which all the velocity gradients are of similar size,but it is in practice (especially numerically) simpler to retain the additional velocity gradientsdue to lateral inflow in the viscosity.With these assumptions, the terms δsPeBL∂TBL/∂ts and δsU(∂yM/∂xs)∂TBL/∂Y in (3.55)are small and the leading order temperature model isPeBL(∂yM∂ts∂TBL∂Y+ V (0)∂TBL∂Y+W (0)∂TBL∂Z)−(∂2TBL∂Y 2+∂2TBL∂Z2)= αBLTn+1 for 0 < Z < S(0) (3.84a)PeBLΓ(∂yM∂ts∂TBL∂Y)− κ(∂2TBL∂Y 2+∂2TBL∂Z2)= 0 for Z < 0, (3.84b)withT =(∣∣∣∣∂U (0)∂Y∣∣∣∣2+∣∣∣∣∂U (0)∂Z∣∣∣∣2) 12n. (3.85)This is a travelling wave problem and we have arrived at a generalized version of the model inSchoof (2012) with the main novelty that we can account for lateral inflow of ice from the ridge.The assumptions above also imply that the ridge problem has a small Peclet number as Per =δrPeBL. The heating parameter is αr = (s/δs)(s/r)1/(n+1)αBL. The ratio (s/r)1/(n+1) =64[τyz]r/[τxy]s describes the relative magnitude of vertical shear stress in the ice ridge and lateralshear stress in the ice stream. Recall that s/r equals Ws/W and that it is always less than 1,i.e., s/r < 1. Then it also follows that αr  αs ∼ O(1) and αr  1, and for κ = kbed/k = 1we obtain a linear temperature profile Tr = νr(sr − zr)− 1 as solution of (3.41). Together with(3.82) this yieldslimY→−∞TBL(X, Y, Z, ts) = νs(S(0) − Z)− 1 for αBL ∼ O(1), PeBL ∼ O(1). (3.86)This simplification will not hold for a real ice sheet, where αBL  1 and PeBL  1. In section 4.5we turn to the solution of the boundary layer model for this case, which per the considerationsabove implies that Per  1 and αr  1 does not necessarily hold and we have to revisit ourlateral boundary condition for Y → −∞.For a widening ice stream, ∂yM/∂ts > 0, the work in Schoof (2012) suggests that we mightbe able to find a unique solution of the migration rate, but only if we supplement the boundarycondition (3.57d) with the condition that the stream-side part of the bed is not actively freezing−∂TBL∂Z∣∣∣∣++ κ∂TBL∂Z∣∣∣∣−≤ 0 for Y > 0, Z = 0. (3.87)If this condition were violated, then as the ice stream widens the newly unfrozen bed wouldbe subjected to net heat loss. At this point the bed should not yet contain any liquid waterto freeze and would therefore not remain unfrozen. This is in contradiction to our assumptionthat the ice stream is widening.The situation is different for a shrinking ice stream, ∂yM/∂ts < 0, where we expect freezingand therefore a net negative heat flux on the temperate side for Y > 0. This heat flux togetherwith the available amount of water stored in the bed then dictates the migration rate, as thiswater needs to be frozen before the ice stream margin migrates inward (see section 2.3 andthe appendix of Schoof, 2012). The inwards migration rate therefore depends not only on themechanical and thermal processes in the ice stream and at the bed, but potentially also on thehydraulic properties of the bed.In the remainder of the thesis we focus on ice stream widening and pursue the numerical65solution of (3.61)–(3.66) with (3.70)–(3.71), (3.77)–(3.78), (3.84), (3.57) and (3.87) next.664 | Solution of the boundary layer modelIn chapter 3 we developed a boundary layer model that describes the physics of the marginbetween an ice stream and an ice ridge. For dynamically evolving geometries, the solution ofthis model provides the necessary coupling of stream and ridge. In particular, it is necessary todetermine the term ∂yM/∂ts, which describes the evolution of the ice stream margin and whichis required in the boundary conditions of the ice stream model (3.27)–(3.31).The goal of this chapter is to understand which physical processes determine margin mi-gration and to derive a parameterization of the migration velocity in terms of a minimal setof parameters. We will address both objectives by using numerical methods. Before we turnto the solution of the boundary layer model, we express the model in the most succinct formpossible. To do so we introduce another rescaling which allows us to absorb quantities such asthe dimensionless ice thickness, dimensionless inflow rate, and dimensionless lateral shear stressin the ice stream margin into a minimal number of dimensionless groups. We defineZ∗ =Zss, P ∗ =P (1)ss, U∗ =U (0)τnxy,s|ys=−yMss, (4.1a)(V ∗,W ∗) =(n+ 1)(n+ 2)ssqr|yr=1−(s/r)yM(V (0),W (0)). (4.1b)The velocity in the downstream direction satisfies∂∂Y(η∂U∂Y)+∂∂Z(η∂U∂Z)= 0, (4.2a)η∂U∂Z= 0 at Z = 0, Y > 0, U = 0 at Z = 0, Y < 0, (4.2b)∂U∂Z= 0 at Z = 1, U → 0 for Y → −∞, η∂U∂Y→ 1 for Y →∞, (4.2c)where we have dispensed with asterisks on the rescaled variables. The viscosityη =121/n[∣∣∣∣∂U∂Y∣∣∣∣2+∣∣∣∣∂U∂Z∣∣∣∣2+ χ2(∣∣∣∣∂V∂Z+∂W∂Y∣∣∣∣2+ 2∣∣∣∣∂V∂Y∣∣∣∣2+ 2∣∣∣∣∂W∂Z∣∣∣∣2)] 1−n2n(4.3)67couples downstream and across-stream flow. Hereχ =n+ 2n+ 1QrAτns h2s, (4.4)where hs is the dimensional ice stream thickness, τs = τxy,s|ys=−yM [τxy]s is the lateral shear stressat the ice stream margin, and Qr = [v]r[zr]qr is the ice flux from the ridge into the margin,evaluated at the margin with the ice stream. In this formulation, the ice flow problem dependsonly on n and χ.As discussed in the previous chapter, we retain χ ∼ s/δs in equation (4.3) even thoughs/δs  1. This ensures that the viscosity remains finite in the far field when Y → −∞. In thenumerical solutions presented in this chapter, we typically use χ = 0.01. We have confirmedthat smaller values of χ do not change our results.The across-stream flow is described by∂∂Y(2η∂V∂Y)+∂∂Z[η(∂V∂Z+∂W∂Y)]−∂P∂Y= 0, (4.5a)∂∂Y[η(∂V∂Z+∂W∂Y)]+∂∂Z(2η∂V∂Y)−∂P∂Z= 0, (4.5b)∂V∂Y+∂W∂Z= 0, (4.5c)with boundary conditionsW = η(∂V∂Z+∂W∂Y)= 0, 2η∂W∂Z− P + S ′ = 0 at Z = 1, (4.5d)V = W = 0 for Y < 0, Z = 0, η∂V∂Z= W = 0 for Y > 0, Z = 0, (4.5e)W → 0,∂V∂Z→ 0 for Y →∞, (4.5f)V → 1− (1− Z)n+1 , W → 0 for Y → −∞, (4.5g)where S ′ is related to S(0) and S(1) through S ′ = S(1)/S(0).To simplify the temperature model, we define a reduced temperatureΘ = −TBLTR, (4.6)68relative to the dimensionless far field bed temperature in the ice ridge TR,TR = limY→−∞TBL(xs, Y, 0, ts). (4.7)TR is measured relative to the melting point and we require TR < 0. With the assumed lineartemperature profile in the far field introduced at the end of chapter 3,limY→−∞TBL = −νs(Z − 1)− 1, (4.8)we have TR = −(1− νs).Omitting the asterisks on V ∗, W ∗, Z∗ once more, the rescaled heat equation takes the formvm∂Θ∂Y+ Pe(V∂Θ∂Y+W∂Θ∂Z)−(∂2Θ∂Y 2+∂2Θ∂Z2)= αa for 0 < Z < 1, (4.9a)vmΓ∂Θ∂Y− κ(∂2Θ∂Y 2+∂2Θ∂Z2)= 0 for Z < 0, (4.9b)with the heat production terma =121+1/n(∣∣∣∣∂U∂Y∣∣∣∣2+∣∣∣∣∂U∂Z∣∣∣∣2) 1+n2n, (4.10)where we have omitted an O(χ2) term, analogous to that in (4.3). We obtain the new dimen-sionless groupsPe = PeBL(n+ 2)(n+ 1)qr =(n+ 2)(n+ 1)ρcQrk, ν =qgeohskbed(Tm − T0)=(Tb − T0)(Tm − T0), (4.11a)α =αBL−TRτxy,s|n+1y=−yms2s = 2Aτn+1sk(Tm − Tb)h2s, (4.11b)where Tb = [T ]TR + Tm is the dimensional bed temperature in the far field of the ridge. Wehave also definedvm = PeBLss∂yM∂ts(4.12)as the dimensionless speed with which the ice stream margin migrates outwards. This is related69to the dimensional migration speed asdymdt=kρchsvm. (4.13)The boundary conditions are:Θ =1TRat Z = 1,∂Θ∂Z→ 0 as Z → −∞, (4.14a)Θ→ ΘR(Z) for Y → −∞,∂Θ∂Y→ 0 for Y →∞, (4.14b)Θ < 0 and −∂Θ∂Z∣∣∣∣++ κ∂Θ∂Z∣∣∣∣−= 0, for Y < 0, Z = 0, (4.14c)Θ = 0 and −∂Θ∂Z∣∣∣∣++ κ∂Θ∂Z∣∣∣∣−≤ 0, for Y > 0, Z = 0, (4.14d)where ΘR is the scaled ridge temperature with ΘR(0) = −1.We have now arrived at a model for the margin migration velocity vm which depends onthree dimensionless groups (Pe, α and ν) and on three material constants (κ, Γ and n). Thedimensionless groups depend on parameters that are determined by the outer problems in the iceridge and the ice stream. They therefore describe the coupling of the boundary layer dynamics,and consequently the migration velocity, to these outer regions.For comparison with previous work (Schoof, 2012), we assume κ = 1 and Γ = 1, so that vmcan only depend on Pe, α, ν, and n. As outlined in section 3.7, we start by treating α and Peas O(1) parameters in section 4.3 and then turn to the asymptotic limit of α 1 and Pe 1in section 4.5. In the next section we describe the numerical solution of the model.4.1 The finite-element solver Elmer/IceTo solve the boundary layer model (4.2)–(4.14), we use the finite element package Elmer/Ice(Gagliardini et al., 2013; Favier et al., 2012; Zwinger et al., 2007). Specifically, we use the Stokesand heat solvers implemented in Elmer/Ice. For a non-linear rheology (n > 1 in equation (4.3)),the ice flow model (4.2)–(4.5) is solved using a combination of Picard and Newton iterations.A two-dimensional unstructured triangular mesh is created with the software Gmsh (Geuzaine70and Remacle, 2009). Even though our geometry is two-dimensional, the mesh is subsequentlyextruded to a three-dimensional tetrahedral mesh in order to solve for all three velocity compo-nents. The boundary layer domain is a strip of ice on of a half space of bedrock. The numericalmodel must be solved on a finite domain. We use the rectangle −5 < Y < 5, −9 < Z < 1 (seefigure 4.1) for most of our computations. The far field boundary conditions (4.2c)2, (4.5f) and(4.14b)1 are applied at Y = −5, (4.2c)3, (4.5g) and (4.14b)1 are applied at Y = 5, and (4.14a)2is applied at Z = −9. High resolution is required at the transition between no slip and freeslip at (Y, Z) = (0, 0). In this region we locally refine the mesh with element edge lengths ofapproximately 2.5 · 10−6.4.2 Ice flow and heat productionWe begin with the solutions to the ice flow problem (4.2)–(4.5), which are necessary for solvingfor the temperature field. The velocity U in the downstream direction determines the heatproduction, while the velocities V and W in the across-stream and vertical directions determinethe advection. The governing equations for the velocities do not depend on any parameters apart−5 0 5−401YZFigure 4.1: Example of a mesh of the computational domain with refinement near the originand along the bed. The resolution in the center at (Y, Z) = (0, 0) is 2.5 · 10−6 hs, where hs isthe ice thickness. For hs = 1000 m, this corresponds to 2.5 mm.71↓e1+λS′dridge streamsedimentYZ−1 0 1210 0 1−0.03−0.02−0.0100.01 S′YS′e−1 0 100.51YPfYb−1 0 1|(V,W)|01Yc−1 0 1a03YZ a−1 0 101U02Figure 4.2: Panel a: Contours of downstream velocity component U at contour intervals of 0.1.Panel b: transverse velocity field (V,W ) shown as arrows, shading indicates magnitude of thetransverse velocity. Panel c: contours of heat production a at contour intervals of 0.1, a = 1shown as a bold line. Note that contours of a are less smooth than those of U because theerror in ∇U is larger than the error in U itself. Panel d: First order corrected surface elevation1 +λS ′ with λ = 0.1 for illustrative purposes. Note that the calculation used χ = 0.01. Panel e:Close up of first order surface correction near the origin S ′. Panel f: Pressure with hydrostaticcontribution removed at the upper surface Z = 1. Calculations were done with n = 3.from n. They are also uncoupled from the temperature problem and are therefore identical forall choices of Pe, α and ν.The downstream velocity U exhibits no vertical shear in the ice stream far field, but insteadhas a prescribed lateral gradient of η∂U/∂Y → 1, as shown in figure 4.2a, where contours areplotted in the (Y, Z) plane. There is intense shear around the transition from no slip to freeslip at the origin, and the downstream velocity goes to zero in the ice ridge far field. Theintense shear around the origin translates to very high dissipation rates, a = 2−1−1/n|∇U |1+1/n,as shown in figure 4.2c (∇ refers to the gradient operator in the (X, Y ) plane). In fact, as wewill show in section 4.5, there is a singularity in strain rate and therefore in heat production atthe origin.The lateral velocity field V transitions from a shearing flow in the ice ridge to a plug flowin the ice stream, as shown in figure 4.2b. The vertical velocity component W is zero at allboundaries, but is negative near the origin, corresponding to downward motion of ice towardsthe bed. The near-field behaviour close to the transition point is analyzed in more detail in72section 4.5, where we compare our results with an asymptotic solution in the region around(Y, Z) = (0, 0).Although it is not required in the solution for the migration velocity, we can also determinethe first order correction to the ice surface in the boundary layer, which is flat at leading order.Figure 4.2d shows the corrected upper surface 1+λS ′ as a blue line relative to the leading orderboundary layer geometry. For legibility, we display the corrected surface as if λ = 2s/δs = 0.1,even though the calculations were done with χ = 0.01, where χ scales as s/δs  λ.Panel e shows a close up of S ′. The corrected surface approaches a uniform upward slopeinto the ridge. Towards the ice stream (Y → ∞) the surface slope asymptotes to zero, as isto be expected from the behaviour of the outer problems. Around the transition from free slipto no slip, there is a surface dip with a local minimum, which is characteristic for such sharptransitions in sliding (e.g., Barcilon and MacAyeal, 1993; Sergienko et al., 2007). Figure 4.2fshows the first order correction in the pressure at Z = 1, which also increases towards the iceridge due to the increase in ice thickness.4.3 Temperature field and migration velocityWe use the ice flow solutions above to determine the temperature field, which is the solutionof (4.9)-(4.14). Solving for the temperature field requires that we determine the outwardsmigration speed of the ice stream margin, vm, which appears as a parameter in (4.9). If wedisregard the inequality constraints (4.14c)1 (requiring that the temperature is less than zerofor Z < 0) and (4.14d)2 (requiring no freezing for Z > 0), the temperature equation can besolved with any choice of vm. However, with an arbitrary choice of vm these inequalities are notguaranteed to be satisfied.The work in Schoof (2012) shows that for Pe = 0 and n = 1, there is a unique value of vmfor a given value of α and ν such that both inequality constraints are satisfied simultaneously.vm effectively acts as an extra advection term that ‘transports’ cold ice from the ice ridge intothe boundary layer and thereby counteracts the heating term αa.In Schoof (2012), migration velocities that are too small correspond to insufficient transportof cold ice into the ice stream and lead to a build-up of excess heat. Temperatures at the bed73on the ice ridge side then rise above the melting point, which is in violation of the inequalityconstraint (4.14c)1. Conversely, migration velocities that are too large lead to excess transportof cold ice and therefore to freezing inside the stream bed, with singular freezing rates near theslip to no-slip transition. In this case, the constraint (4.14d)2 is violated.It is important to realize that the second constraint (4.14d)2 does not need to hold uncondi-tionally and only applies here because we assume a widening ice stream, so that vm is positive.More generally, in order for sliding to occur the bed on the unfrozen side of the no-slip to sliptransition needs to have a non-zero water content. As the transition moves into the ice ridge,the newly unfrozen bed starts with zero water content. We therefore require a non-negative heatflux into the bed, at least close to the transition point, to ensure that the water content increasesrather than decreases. Hence the constraint (4.14d)2 must hold at least in the neighborhood ofthe transition point.In the reverse case of a narrowing ice stream, we expect a negative heat flux into the bedto cause freezing and remove water from the bed. A singular freezing rate can be shown tocorrespond to a finite migration rate into the ice stream, see appendix B in Schoof (2012).We test whether the behaviour described above (values of vm that are too small violate(4.14c)1, values of vm that are too large violate (4.14d)2 and only a single value of vm fulfillsboth constraints) also holds for n > 1 and Pe > 0. In figure 4.3 we show solutions to the heatflow problem without the inequality constraints imposed using fixed values of α = 20, ν = 0.5,n = 3, and Pe = 10. Here and in all other plots of the temperature field in this chapter, weshow the non-dimensional boundary layer temperature for ease of interpretation, i.e., T = TBL.To the left (column a) we show a case where vm is too small, and temperatures for a smallregion of Y < 0 exceed the melting point (panel a2). In the middle (column b) we show vm toolarge, which leads to freezing for Y > 0 with an apparently singular heat flux at the origin (panelb3). In column c, we show results with a value of vm for which both inequality constraints arefulfilled up to the accuracy of the numerical method. Numerically we can only determine thatthe inequalities are satisfied at nodes and edges of the grid. We can therefore only test whethera certain vm fulfills the inequalities up to the discretization error. It is possible for more thanone vm to do this and for the shown example with a resolution of 2.5 · 10−6, viable values ofvm were found in the range 6.0 < vm < 6.3. This range increases with coarser resolution. We74Z0a1−1 0 1010b1−1 0 10c1  −1 0 1T−10−0.0500.05Ta2−0.25 0.25−50050Y∂T ∂z|+-∂T ∂z|− a3b2−0.25 0.25Yb3c2−0.25 0.25Yc3Figure 4.3: Numerical scheme to determine vm. Each column (a-c) shows the temperature fieldat contour intervals of 0.1 (row 1), temperature at the bed (row 2) and net heat flux into thebed (row 3) for α = 20, ν = 0.5, Pe = 10, and n = 3. Column a shows results for vm = 3.2 withconstraint (4.14c)1 violated in panel a2. Column b shows results for vm = 9.6 with constraint(4.14d)2 violated in panel b3 and an apparently singular heat flux at the origin. Column cshows results for vm = 6.2, satisfying both constraints. Note that the results in rows 2 and3 are plotted for a narrow range of Y around the origin. The region in which the inequalityconstraints are violated can be quite small even for substantially incorrect values of vm. Thisunderlines the need for a high grid resolution around the origin in our computations. Note thatT = TBL is the boundary layer temperature.assume that there is a unique value of vm in the continuum limit.We determine the migration speed vm iteratively by adapting the bisection method. Theupper limit of the search interval is a migration velocity that is too big and therefore violates(4.14d)2, the lower limit of the search interval violates (4.14c)1. As with a standard bisectionmethod we halve the search interval at every iteration. We determine in which interval wecontinue the search based on which inequality constraint is violated at the midpoint: if (4.14c)1is violated we continue in the upper half, otherwise in the lower half.Note that while only the case of vm too small violates the temperature constraint at the bed,75all three temperature fields in figure 4.3 have solutions with T > 0 in the ice. This is a featureof all the temperature fields corresponding to margin migration velocities calculated here. Itimplies that temperate ice (i.e., ice at the melting point with a non-zero melt water fractionin the ice matrix) forms in the margins of a widening ice stream. Our model at this pointassumes that the same heat equation can be used to model temperate ice as cold ice, whichis a very special and not necessarily realistic case of an enthalpy gradient model (Aschwandenet al., 2012). We continue with this model for now to facilitate comparison with prior work. Insection 4.7, we extend the theory to a more sophisticated model for temperate ice.4.4 Migration velocity as a function of forcing parame-tersWe now investigate how the migration velocity depends on the parameters Pe, α and ν, whichcouple the boundary layer to the outer problems of the ridge and stream flow.Figure 4.4a shows that the migration speed vm increases monotonically with the amount ofshear heating, provided α exceeds a minimum value. The migration term vm∂T/∂Y effectivelyacts as a cooling term that balances heating, so that larger α requires larger vm. In fact, itappears that vm is an almost linear function of α, an observation to which we return later.In the absence of lateral advection (Pe = 0) and with a constant viscosity (n = 1), Schoof(2012) solves the boundary layer model, (4.2)–(4.3) and (4.9)–(4.14) semi-analytically with aWiener-Hopf method. We use this solution to test our numerical method. When Pe = 0 it isstraightforward to see from (4.9)–(4.14) that vm depends only on α by defining the temperaturefield in the form Θ˜ = Θ + 1 − (1 + 1/TR)Z (Schoof, 2012). In figure 4.4a we show numericalsolutions of vm as a function of α for n = 1 and Pe = 0, computed at different mesh sizes nearthe origin. We compare the numerical results with the solutions in Schoof (2012). Only for thehighest resolution of 2.5 · 10−6 are we able to match the benchmark solution to an error of lessthan 1%. The need for this high resolution can probably be ascribed to the singular behaviourin the heating term a at the origin, which is illustrated in figure 4.2c (see also section 4.5).Adding advection to the heat balance shifts the migration speed towards smaller values.760 20 40510152025αmargin migration speed v m  aanalytical10−55 ⋅ 10−62.5 ⋅ 10−60 20 40510152025αn = 1  bPe = 0Pe = 10Pe = 30Pe = 500 20 40510152025α  n = 1, Pe = 10cν = 0.1ν = 0.5ν = 0.90 20 40510152025α  dn = 1, Pe = 0n = 3, Pe = 0n = 3, Pe = 13n = 3, Pe = 26Figure 4.4: Migration velocity as function of forcing parameters. Panel a: Comparison ofnumerical results and semi-analytical results for a constant viscosity and no advection. Panelb: Effect of adding advection. Panel c: Solutions for n = 1,Pe = 10 and different values of ν.Panel d: Runs with Pe = 0, Pe = 13 and Pe = 26 for n = 3.Figure 4.4b shows that increasing values of Pe reduces the corresponding vm. This is becausethe cooling required to balance heat production no longer comes purely from migration into coldridge ice (represented by vm), but can also be caused by the flow of cold ice into the boundarylayer (represented by Pe).For non-zero advection (Pe > 0) it is not obvious from (4.9)–(4.14) that the migration speedshould not depend on ν. Figure 4.4c shows solutions of vm for values of ν between 0.1 and 0.9for Pe = 10. For increasing ν, vm is shifted towards smaller values. This shift appears to beconstant with increasing α, so that relative discrepancies in the migration velocities vanish forlarge α.Any additional dependence on ν is small enough that it is not visible. For n = 3, we see aqualitatively similar picture of a migration speed that increases with α and decreases with Pe(figure 4.4d). The calculated dimensionless migration velocity is smaller than for n = 1, butthe two cases cannot be compared directly in a meaningful manner.The temperature fields in the ice stream margin will depend on the choice of forcing parame-ters α and Pe and the migration velocity vm. Typical temperature fields for migration velocitiesthat fulfill the inequality constraints (4.14c)–(4.14d) are shown in figure 4.5, using the sameplotting scheme as in figure 4.2.77Z  a01  b  cT−10YZ  d−1 0 101Y  e−1 0 1 Y  f−1 0 1T−10Figure 4.5: The effect of shear heating and advection on the temperature field, same plottingscheme as in row 1 of figure 4.2. Top row has Pe = 0 and α = 5 (panel a), α = 7 (panel b),α = 9 (panel c). Bottom row has α = 8 and Pe = 3.2 (panel d), Pe = 6.4 (panel e), Pe = 12.8(panel f). All calculations were done with n = 3 and ν = 0.5.The top row shows how temperature solutions change as the heating parameter α increaseswhile the advection parameter Pe is held constant. The immediate effect of increasing α isto heat more of the ice around the no-slip to free slip transition and inside the ice stream.This however leads to a simultaneous increase in the migration velocity, which leads to lowertemperatures in the bed and to the left of the transition point, because these regions have lesstime to heat up.Conversely, increasing the lateral advection rate Pe while keeping α constant has the oppositeeffect. There is a reduction in the migration velocity vm, giving the bed and parts of theice to the left of the no-slip to free slip transition more time to warm up (figure 4.5d–f).Advection velocities are greater in the upper layers of the ice (see figure 4.2b), which decreasesthe temperatures there for larger Pe.4.5 Asymptotic solutionsSo far we have studied the dependence of vm on a few combinations of α, ν and Pe, but we stillneed to quantify this dependence in a useful manner. In particular, it would be useful to havea formula that allows the implemention of margin migration in large scale models of the outer78regions, the stream and the ridge. Deriving such a formula is difficult in general, but the aboveresults strongly suggest that vm increases linearly with α when α is large. Moreover, we alreadyknow that Pe and α are large parameters.In this section we will use these properties to simplify the relationship between vm and Pe,α and ν. First, we will show that vm indeed increases linearly with α. This is based on ananalysis of the case α → ∞ for fixed Pe. Subsequently we will investigate what happens if αand Pe are both allowed to be large (the ‘simultaneous’ limit α→∞, Pe→∞).In both cases (α large, and Pe and α simultaneously large), the migration velocity is deter-mined by processes that occur close to the transition from no slip to slip. This motivates us tostudy the local behaviour of the velocity and temperature field in this region in more detail.4.5.1 The velocity field close to the transition pointNear the origin of our geometry, i.e. for R = (Y 2 + Z2)1/2 → 0, the equation for the down-stream velocity (4.2a) with boundary conditions (4.2b) is identical to the model for a crack-tipconsidered in Rice (1967, 1968). In polar coordinates, the velocity solution close to the transitionfrom free slip to no slip is of the formU ∼ CuR1n+1√2nn+ 1A2n+1θ + cos θA1−n1+nθ for R→ 0, (4.15)where R =√Y 2 + Z2, cos θ = Y/R, Cu is a constant that depends on the far field conditionsandAθ =n2 − 14ncos θ +√(n2 − 14n)2cos2 θ +(n+ 1)24n. (4.16)Figure 4.6a confirms that our numerical solution reproduces this behaviour. From (4.15) theasymptotic behaviour of the heat production rates (4.10) isa ∼(Cu2)1+1/nR−1A−1θ . (4.17)The important feature of this result is that the heat production is singular, behaving as R−1near the transition point. This should not come as a surprise: a similar behaviour for n = 17910−3 10−1100RU  aElmer/Ice0.6 R1/(n+1)10−3 10−110−210−1RV  bElmer/Ice3.5⋅10−2 Rβ10−3 10−1−10−2−10−3RW  cElmer/Ice−3.8⋅10−3 RβFigure 4.6: Comparison of numerical velocity solutions with asymptotic solutions from Rice(1967) and the solutions of the boundary value problem (4.20) for θ = pi/8 (see appendix C.1).Panel a shows solutions of the downstream solution U , panel b shows solution of the across-stream solution V and panel c shows solutions of the vertical velocity W . n = 3 in all threecases.appears in Schoof (2004, 2012) and for n = 3 in Suckale et al. (2014). For the frequently usedspecial cases of n = 1 and n = 3, a can alternatively be written asa ∼ CaR−1 ×const. for n = 1,(√3 + cos2 θ + cos θ)−1for n = 3.(4.18)The local behaviour of the across-stream velocities (V,W ) is more difficult to determine. Fora constant viscosity (n = 1), Barcilon and MacAyeal (1993) show thatV ∼ R1/2(11 cosθ2− cos3θ2), W ∼ −R1/2(sinθ2+ sin3θ2). (4.19)For n 6= 1, the problem of finding the local behaviour of V and W is complicated by the factthat the viscosity is determined by |∇U |, where the local behaviour of U is given by (4.15). Tofind a generalisation of (4.19) for n 6= 1, we assumeV ∼ RβVθ(θ), W ∼ RβWθ(θ) as R→ 0, (4.20)and solve a local version of (4.5a)–(4.5g) as an eigenvalue problem for β, see appendix C.1for details. With n = 3 we obtain β ≈ 0.271, so that the gradients in transverse velocityare slightly less singular than the gradients in downstream velocity with U ∼ R1/(n+1), where1/(n + 1) = 0.25. Once again we find that our numerical solutions reproduce this behaviour,see figure 4.6b-c.80Equipped with (4.17) and (4.20), we now turn to the behaviour of the temperature fieldclose to R = 0. We begin by considering the limiting behaviour for α Temperature solutions and migration velocity for α 1For large α with constant Pe, there is strong heating in the margin. To achieve sufficient coolingto balance this heating rate, the margin migrates rapidly into the ice ridge. The conductivetransport of heat into this cold ice can only occur on a short length scale, meaning that con-duction is confined to a small, nested boundary layer around the transition point ( by ‘nested’we mean that it is embedded in another boundary layer).We find that the relevant conductive length scale is α−1. We set (Y, Z) = α−1(Y˜ , Z˜), a = αa˜and (V,W ) = α−β(V˜ , W˜ ) using (4.17) and (4.20). The heat equation (4.9) can then be writtenasvmα−1∂Θ∂Y˜+ α−(1+β)Pe(V˜ , W˜ ) ·(∂Θ∂Y˜,∂Θ∂Z˜)−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= a˜ for 0 < Z˜, (4.21a)vmα−1∂Θ∂Y˜−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= 0 for Z˜ < 0, (4.21b)with boundary conditions near the transition pointΘ < 0 and −∂Θ∂Z˜∣∣∣∣++ κ∂Θ∂Z˜∣∣∣∣−= 0, for Y˜ < 0, Z˜ = 0, (4.22a)Θ = 0 and −∂Θ∂Z˜∣∣∣∣++ κ∂Θ∂Z˜∣∣∣∣−≤ 0, for Y˜ > 0, Z˜ = 0. (4.22b)(4.21a)–(4.21b) show that the migration speed increases linearly with α for α → ∞, i.e.,vm ∼ α. We therefore set vm = αVm. The leading order temperature model (omitting terms ofO(α−1)) is thenVm∂Θ∂Y˜−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= a˜ for 0 < Z˜, (4.23a)Vm∂Θ∂Y˜−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= 0 for Z˜ < 0, (4.23b)while the leading order outer problem for temperature (away from the conductive boundary81Figure 4.7: Boundary layer structure for α  1 and Pe  1. The conductive boundary layeris of extent Lα = α−1, the vertical length scale of the advective boundary layer is Lp = Lβα =Pe−β/(1+β).layer) is simplyVm∂Θ∂Y= a, (4.24)where a is the heat production rate in the original boundary layer coordinates Y , Z (recallthat a depends only on Y , Z and n but is independent of α). If we continue to assume thatthere is a unique solution for the migration velocity vm, then the fact that (4.23) and (4.24) areindependent of α and ν implies that Vm is simply a constant Cα, so thatvm ∼ Cαα for α→∞ (4.25)at constant Pe. This confirms that the apparent linear relationship between vm and α infigure 4.4a–d reflects a robust asymptotic behaviour. Note that Cα is independent of Pe, whichis consistent with figure 4.4a–d: the slope of the curves of vm against α is always the same forlarge α for a given value of n.Panel a of figure 4.8 shows this limiting behaviour, plotting Vm = vm/α against α at constantPe. The calculations were done for n = 3. Regardless of the value of Pe chosen, Vm asymptotesto a value of ≈ 0.99.We have arrived at (4.25) without computing and matching the solutions of the conductiveboundary layer and the outer problems. To do this, we would need the temperatures in theouter region close to the bed to match with the inner solution, because the conductive boundarylayer remains close to the bed. As (4.24) is an advective problem with transport purely in the82positive Y –direction, temperatures near the bed are determined by the far field temperaturesat the bed for Y → −∞ and the integrated heat production rate along the bed. This is therationale for defining Θ in terms of the far field bed temperature TR; as we will see later thisholds even when Pe is large as well.The limiting behaviour (4.25) above holds for any constant Pe. How quickly we approach thislimit depends on the value of Pe: α must be sufficiently large so that advection terms becomenegligible in the conductive boundary layer problem (4.23) and the outer problem (4.24). Inthe next section we consider the distinguished limit in α → ∞ and Pe → ∞. In this case wecannot neglect the effect of advection.4.5.3 Temperature solutions and migration velocity for α  1 andPe 1Restoring advection terms in the boundary layer problem (4.21) requires the product α−(1+β)Pein (4.21a) to be of O(1), i.e., the Peclet number must scale as Pe ∼ α1+β. We therefore putΛ = Pe−11+βα (4.26)and consider the case of Λ strictly O(1) in (4.23) which becomes (neglecting the O(α−1) =O(Pe−1/(1+β)) term):Vm∂Θ∂Y˜+ Λ−(1+β)(V˜ , W˜ ) ·(∂Θ∂Y˜,∂Θ∂Z˜)−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= a˜ for 0 < Z˜, (4.27a)Vm∂Θ∂Y˜−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= 0 for Z˜ < 0. (4.27b)In addition, we need to consider the equivalent of (4.24), which forms the outer problem tothis conductive boundary layer. The complication in adding advection to the outer problem isthat the velocity (V,W ) is partly a shearing velocity, with V = W = 0 at the bed for Y < 0.However, as discussed above, the far field of the conductive boundary layer must match withtemperatures in the outer region near the bed, which is precisely why we need to consideradvection near the bed.83To identify leading order terms in the advective boundary layer, we first need to understandthe transverse velocity field near the bed. Note that V = 0 implies that ∂V/∂Y = 0 at the bed,so ∂W/∂Z = 0 by mass conservation. By a Taylor expansion, we obtain V ∼ Z, W ∼ Z2.Near-bed advection in the outer problem is captured by considering a thin region of verticalextent Pe−β/(1+β)  1 relative to ice thickness, labelled the ‘advective boundary layer’ in fig-ure 4.7. Within this region, V ∼ Pe−β/(1+β), W ∼ Pe−2β/(1+β), and we rescale V = Pe−β/(1+β)Vˆ ,W = Pe−2β/(1+β)Wˆ , Z = Pe−β/(1+β)Zˆ.The vertical coordinate in the advective region is related to the vertical coordinate in theconductive boundary layer through Z˜ = ΛPe(1−β)/(1+β)Zˆ. For β < 1, Zˆ = O(1) implies thatZ˜  1. For n = 1, the exponent β equals 1/2, and for n = 3 we have β ≈ 0.27. Thereforethe near-bed advective layer is a viable outer region to the conductive boundary layer. Theadvective layer has a much larger vertical and horizontal extent than the conductive boundarylayer.The rescaled advective model isVm∂Θ∂Y+ Λ−1(Vˆ , Wˆ ) ·(∂Θ∂Y,∂Θ∂Zˆ)− Pe2β−11+β Λ−1∂2Θ∂Zˆ2− Pe−1/(1+β)Λ−1∂2Θ∂Y 2= a for 0 < Zˆ,(4.28a)Vm∂Θ∂Y− Pe2β−11+β Λ−1∂2Θ∂Zˆ2− Pe−1/(1+β)Λ−1∂2Θ∂Y 2= 0 for Zˆ < 0.(4.28b)The leading order outer problem for n = 1 is, omitting terms of O(Pe−1/(1+β))Vm∂Θ∂Y+ Λ−1(Vˆ , Wˆ ) ·(∂Θ∂Y,∂Θ∂Zˆ)− Λ−1∂2Θ∂Zˆ2= a for 0 < Zˆ, (4.29a)Vm∂Θ∂Y−∂2Θ∂Zˆ2= 0 for Zˆ < 0. (4.29b)For n = 3 (or generally for n > 1 and β < 1/2), the outer problem is,Vm∂Θ∂Y+ Λ−1(Vˆ , Wˆ ) ·(∂Θ∂Y,∂Θ∂Zˆ)= a for 0 < Zˆ, (4.30a)Vm∂Θ∂Y= 0 for Zˆ < 0, (4.30b)84101 103 105 10700.51PeV m  bΛ = 100Λ = 10Λ = 1Λ = 0.5102 104 10600.51αV m  aPe = 103Pe = 104Pe = 105Pe = 106Figure 4.8: Asymptotic behaviour of Vm = vm/α. Panel a: Vm against α for different valuesof Pe = 103 (blue circles), Pe = 104 (magenta triangles), Pe = 105 (green stars) and Pe = 106(black dots). Panel b: Vm against Pe for Λ = 100 (blue circles), Λ = 10 (black stars), Λ = 1(red triangles) and Λ = 0.5 (green squares). Note that holding Λ constant as Pe changes impliesthat α changes proportional to Pe1/(1+β).to an error of O(Pe(2β−1)/(1+β)). Once more, we are considering an outer problem that describesa slender region near the bed. Our choice of reduced temperature Θ in terms of the far field bedtemperature TR means that the relevant boundary condition is Θ(Z = 0)→ −1 as Y →∞.From the conductive boundary layer problem (4.27) and the near-bed advective layer (4.29)or (4.30), we can conclude that Vm is a function of Λ alone at leading order. Vm depends on αand Pe only through the ratio Λ = αPe−1/(1+β) and is independent of the far field temperaturegradient ν:Vm ∼ g(Λ) = Λ−1f (Λ) , (4.31)if f(Λ) = Λg(Λ). For n = 1 this appears to be true to the order of O(Pe−1/(1+β)), for n = 3 theerror is substantially larger at O(Pe(2β−1)/(1+β)); with β = 0.27, this is Pe−0.36.Our goal now is to confirm this relationship numerically and to find the approximate formof f (or equivalently, g). The function f now depends on a single variable Λ and is thereforeeasier to find than a function that depends on multiple variables (in this case α, Pe and ν).From the previous section, we already know that when Λ is large, which corresponds to α beinglarge in comparison to Pe1/(1+β), then g(Λ)→ const. or f(Λ) ∼ Λ. This behaviour is confirmedby our numerical results displayed in figure 4.8a, which shows Vm against α for different valuesof Pe.Figure 4.8b shows the convergence of Vm to its limiting form g(Λ) as we make Pe large while850 10 20 30 40 5002040ΛV mΛ  aP e = 107P e = 105P e = 1030.99Λ− 0.320 2 4 6 8 100510ΛV mΛbFigure 4.9: Relation between α, Pe and Vm. Numerical results of Vm = vmPe−1/(1+β) are plottedagainst Λ = Pe−1/(1+β)α. Dashed line: fitted curve Vm = Cp +CαΛ. B: For smaller values of Pethe same qualitative behaviour can be observed, but with a smaller slope.holding Λ fixed. Note that holding Λ fixed means that α must grow in lock step with Pe1/(1+β),meaning that the limit Pe→∞ corresponds to α→∞. Results are shown for n = 3, and theapproach to the limit is relatively slow, with convergence requiring values of Pe ' 106.We already know that the limiting behaviour of g(Λ) and f(Λ) is g(Λ)→ Cα and f(Λ) ∼ CαΛas Λ→∞, where Cα = 0.99... for n = 3. Figure 4.9 shows VmΛ = vm/Pe1/(1+β) plotted againstΛ. Panel a shows results for 0 < Λ < 50 at different values of Pe. In the limit of Pe→∞ thisplot represents f(Λ). The expected linear behaviour of f(Λ) for large Λ is apparent, althoughfor VmΛ to approach f(Λ) requires fairly large values of Pe around 105.It is notable that f(Λ) appears to be approximately affine (an approximately linear functionwith a non-zero intercept) down to small values of Λ. This is shown in panel b, where we zoomin to the range 0 < Λ < 10. As discussed earlier and also in Schoof (2012), a positive migrationvelocity requires a minimum value of α to be exceeded. Here this translates to f(Λ) beingpositive only if Λ exceeds a minimum value Λmin ≈ 0.2 for n = 3. A best fit to f(Λ) in therange 0.2 < Λ < 50 that preserves the limiting slope isf(Λ) = Cp + CαΛ, (4.32)with Cp ≈ −0.3. This gives the limiting behaviour of the migration velocity asvm = CpPe1/(1+β) + Cαα, (4.33)for α 1 and Pe 1.86−1 −0.9 −0.8 −0.7 −0.6 −0.500.51TZ  aT1 = (ν − 1) − νZT2 = (1− ν)Z2 − Z + (ν − 1)T3 = ν(Z − 1)4 − 10 10 20 30 40 5001020304050ΛV mΛ  bT1 = (ν − 1) − νZT2 = (1− ν)Z2 − Z + (ν − 1)T3 = ν(Z − 1)4 − 1Figure 4.10: Influence of temperature field in the ridge on the migration rate. Panel a: tem-perature fields prescribed towards the ice ridge, Y → −∞. Panel b: Vm against Λ for the threedifferent temperature solutions shown in panel a.In the derivation of the boundary layer model we have assumed that Pe is not large and thatwe have a linear temperature profile towards the ridge. We have now investigated the limit oflarge Pe within this reduced model. This raises the question whether there is a limit to howlarge we can make Pe and whether a different shape of the temperature profile in the far fieldtowards the ridge (which is permissible in the large Pe limit) would change our results.If we had assumed a priori that Pe is large in our derivation of the boundary layer model andrelaxed the assumption of a linear temperature profile in the ridge, we would have found that theresults of this chapter still hold as long as δs  Pe−β/(1+β). δs is a the aspect ratio of ice streamthickness to ice stream width with typical values of δs ≈ 0.02. With −β/(1 + β) ≈ −0.2 this isnot a very restrictive requirement. We show this in detail in appendix C.2. In appendix C.2 wealso confirm that the only relevant part of the far field boundary condition is the temperature atthe bed in the far field, which has already been incorporated in our definition of α. Here we onlyshow numerically that different far field temperature profiles with the same bed temperaturesyield the same relationship between Vm and Λ in the limit of large Pe, see figure Sub-temperate slidingThe assumed abrupt transition from no slip to free slip at (Y, Z) = (0, 0) is an idealisation thatis unlikely to hold in a real ice stream margin (Fowler, 2013). This transition gives rise to thesingular stresses and heat production rates described in the previous section, which imply thatno sliding is possible in the ridge for any amount of stress. In reality we would expect slip tooccur, either due to mechanical failure of the ice–bed contact, or due to a residual premelted87YZa2−1 0 101Yb2−1 0 1Yc2−1 0 1Za101 b1 c1a03T−10Figure 4.11: Heat production rates and temperature solutions, same plotting scheme as infigure 4.2. Parameter values are n = 3, α = 12, Pe = 10 and Tc = 3.0 (column a), Tc = 0.9(column b), Tc = 0.5 (column c). Row 1 shows heat production rates a at contour intervals of0.1 and row 2 shows temperature contours at intervals of 0.1.water film at the ice–bed contact, leading to sub-temperate sliding (Fowler, 1986; Echelmeyerand Zhongxiang, 1987; Cuffey et al., 1999). In addition, there may be mechanical failure ordamage production in the ice, helping to alleviate the high stress concentrations (Pralong andFunk, 2005).To investigate how a continuous, rather than an abrupt, transition between free slip and noslip affects the solutions presented here, we consider the simplest case of a plastically deformingbed with a dimensionless yield stress Tc (Schoof, 2004, 2010). Instead of the boundary conditions(4.5e)1, we imposeTxz = Tc U|U | , Tyz = TcV|U | , |U | > 0, W = 0|Txz| ≤ Tc, |U | = 0, W = 0for Y < 0, Z = 0 (4.34)(see appendix C.3 for details). Sliding along the bed dissipates heat, and we replace (4.14c),describing the energy balance of the frozen portion of the bed, withT < 0 and∂T∂Z∣∣∣∣+− κ∂T∂Z∣∣∣∣−= αTc|U | for Y < 0, Z = 0. (4.35)Figure 4.11 shows examples of heat production (row 1) and temperature field (row 2) fora variety of values of Tc. While values of Tc above about 10 have no discernable effect on the88100 101050100TcV m  bPe=5x104Pe=105Pe=5x105100 1010510Tcv m  aabcPe=0Pe=1Pe=10Figure 4.12: Migration speed against Tc for different values of Pe. Panel a: vm against Tc atthe values of Pe indicated on the right, with n = 3, α = 12 and ν = 0.5. Points a, b, c and dcorrespond to the equivalently labelled columns in figure 4.11. Panel b: Vm = vm/α against Tcat the values of Pe indicated, with Λ = αPe−1/(1+β) = 1 fixed, n = 3 and ν = 0.5.heat production and temperature field when viewed at the ice thickness scale as in figure 4.11,smaller values of Tc allow significant slip for Y < 0. As panels b1 and c1 in figure 4.11 illustrate,heat dissipation is then shifted away from the frozen–unfrozen bed transition at Y = 0. Instead,we observe dissipation due to sliding along the bed, as well as shearing in the ice towards thetransition between no slip and slip, which is now at Y < 0 (for instance at Y ≈ −0.3 in panelb1). The location of this transition can also be identified as the point where the heat productionrate at the bed, Tc|U |, goes to zero. It is outside the displayed part of the domain in panel c2of figure 4.11.Figure 4.12a shows how the migration velocity vm is affected by the choice of Tc. Finite Tclarger than about 10 has a negligible effect on the margin migration velocity. For smaller valuesof Tc, the migration speed increases as heat production becomes less concentrated around thefrozen-unfrozen transition point and more heat is produced along the ice–bed interface.Thus far, we have seen that for sufficiently large Tc, regularization of the no-slip to free sliptransition in (4.34) has an insignificant effect on the migration velocity vm. This is presumablybecause heat production and advection change significantly only in a small neighbourhood of theunfrozen to frozen bed transition. Previously, we have found that for large advection (Pe 1)and large heating (α  1), vm is determined by the behaviour of the temperature field in aboundary layer around the unfrozen to frozen bed transition.We now consider whether sub-temperate sliding changes the heat production and advection89in a region that is of comparable size to the boundary layer. In the limit of infinite Tc and forsufficiently large Pe, we expect that vm/α is only a function of Λ. The size of the advectiveboundary layer shrinks as Pe increases, while the region that is affected by the regularization(4.34) shrinks as Tc increases. We might therefore naively expect that for larger values of Pe,larger values of Tc are also required in order to preserve the asymptotic behaviour of vm. In otherwords, for fixed Tc we would expect agreement with our asymptotic results from the previoussection to worsen as Pe increases beyond a certain point at which the thermal boundary layerand the region affected by the regularization are of similar extent.We investigate this by once more calculating migration rates for a fixed choice of Λ = 1at various values of Pe. Our results in 4.12b show that Vm = vm/α is only affected by Pe forTc . 1. If Tc ≥ 10, then Vm appears to be essentially identical to the value calculated in theasymptotic limit. If Tc is decreased below approximately 1, there is an increase in migrationvelocities Vm for all cases. In agreement with our argument above, this increase is strongest forthe largest value of Pe considered, Pe = 5× 105.4.7 Temperate iceAn obvious shortcoming of our model is that it permits temperatures to exceed the meltingpoint in the ice. These elevated temperatures occur in all solutions we have computed so far.In reality, once the melting point is reached, further supply of heat does not lead to an increasein temperature but to the production of melt water. The resulting mixture of water and ice isgenerally referred to as temperate ice. In contrast we will refer to ice below the melting pointas ‘cold ice’ below.In temperate ice, the local volume fraction of melt water is the relevant measure of internalenergy density, just as ρcT is the internal energy density in cold ice. By using the same heatequation everywhere in our domain, we have effectively assumed that melt water in temperateice flows down gradients in volume fraction of melt water with the same diffusivity as that whichapplies to heat flow in cold ice.The assumption that melt water flow in temperate ice can be described by a diffusive modelis the basis of so-called enthalpy gradient methods (Hutter et al., 1988; Aschwanden and Blatter,902009; Aschwanden et al., 2012). These methods usually assume a much smaller diffusivity formelt water content in temperate ice than for heat transport in cold ice. The physical basisof these methods is, however, questionable as they neglect the compaction and gravity-drivendrainage processes that are likely to occur in temperate ice (Fowler, 2001).In this section, we pursue an alternative approach used previously in Zwinger et al. (2007)and Suckale et al. (2014). We revert to treating T as temperature only (whereas T effectivelyserves as a proxy for volume fraction of melt water if we interpret our previous results in thecontext of enthalpy gradient methods when T > 0) and divide the part of the domain occupiedby ice into a cold part Ωc and a temperate part Ωt, with a free boundary Γct separating thetwo, see figure 4.13. In Ω− the heat equation (3.84) continues to hold, supplemented with theconstraint T < 0.Given that T satisfies an advection-diffusion equation, we expect that two boundary condi-tions are required to locate the free boundary Γct. We take these to beT = 0 and −∇T · n = 0 on Γct. (4.36)The first condition, T = 0 on the free boundary, follows from requiring continuity of temperatureacross the boundary. The second condition equation (4.36)2 states that there is no conductiveheat flux from the free boundary into the cold region. More generally, we would expect a Stefancondition to hold: the difference between the normal components of heat flux on either side ofthe boundary must account for the rate at which latent heat is absorbed or released by freezingand melting at the free boundary.A zero heat flux into the cold part of the domain as assumed in (4.36) is then an appropriatecondition where ice moves across the cold-temperate boundary from the cold region into thetemperate region. The temperature is below the melting point (T < 0) in the cold region,which prevents the temperature gradient −∇T · n from being negative, as the temperaturemust increase to T = 0 towards the free boundary. In other words, this scenario can only leadto heat loss from the temperate region and to freezing at the boundary.However, where cold ice flows into the temperate region, we expect it to have initially zerowater content ϕ (solid line in figure 4.13). Consequently, there is no water to freeze at these parts91xy zsheary i shenrgizshear hflozwar wvrp hehuh∇ n arhh(,hwarp hehuh∇ n trhh(,hwtrhnFigure 4.13: Geometry for the boundary conditions of temperate ice. Ωc: ice below the meltingpoint, Ωt: temperate ice, Γct: cold/temperate boundary, ϕ: melt water content.of the cold-temperate boundary, leading to the conclusion that we cannot have −∇T · n > 0.This only leaves the possibility −∇T · n = 0.Conversely, where ice moves across Γct from temperate to cold (dashed line in figure 4.13),it is conceivable that the temperate ice has non-zero water content at the free boundary anda finite temperature gradient (−∇T · n > 0 if n points into the cold region) is required tofreeze this water. Our choice of boundary condition (4.36)2 is therefore legimitate if there areno parts of Γct where ice flows from temperate to cold. This turns out to be the case for all thecalculations we report in this section.We can go further and view (4.36)2 as appropriate if the amount of water that needs to befrozen at the boundary Γct can be minimised to a negligible residual by very efficient gravitydriven drainage. This works if parts of Γct where temperate ice flows into the cold subdomainnever have cold ice underneath them. This is the justification for using the boundary condition(4.36)2 in Suckale et al. (2014).We solve the heat equation (4.9) with boundary conditions (4.14a)–(4.14d) and (4.36) usingthe active set method described in Zwinger et al. (2007). Note that the free boundary conditions(4.36) are not spelled out in Zwinger et al. (2007) but are implicit in their formulation as adiscretized obstacle problem (Hintermu¨ller et al., 2002).One of the advantages of this formulation is that it allows us to solve for the temperature andthe location of the free boundary without having to know the details of drainage and compactionin the temperate region. Dissipation rates are still positive in the temperate region and there92is no transport of sensible heat as the region is isothermal. As discussed earlier, we expect thatthis heat leads to the production of melt which is then likely to drain under gravity and pressuregradients (Fowler, 2001), but we do not model these processes here.Before we present numerical solutions of this adapted model, we derive a local solution forthe temperature field close to the origin. This is analogous to the proceduce in Schoof (2012)and allows us to check our numerical solutions in sections 4.7.3 and 4.7.4. As we will see, thelocal solution is particularly useful as the temperate ice model aggravates resolution challengesin the heat equation around the no slip to free slip transition.4.7.1 Local solution and numerical solutionNear the origin, the heat equation (4.9) is dominated by the diffusive terms and is thereforeapproximately−(∂2Θ∂Y 2+∂2Θ∂Z2)= αa for 0 < Z, (4.37a)−κ(∂2Θ∂Y 2+∂2Θ∂Z2)= 0 for Z < 0, (4.37b)with boundary conditions (4.14c)–(4.14d):Θ < 0 and −∂Θ∂Z∣∣∣∣++ κ∂Θ∂Z∣∣∣∣−= 0, for Y < 0, Z = 0, (4.38a)Θ = 0 and −∂Θ∂Z∣∣∣∣++ κ∂Θ∂Z∣∣∣∣−≤ 0, for Y > 0, Z = 0. (4.38b)For R = (Y 2 + Z2)1/2 and θ the angle between the Y –axis and R (see figure 4.14), the heatproduction rate is given by equation (4.18),a ∼ CaR−1 ×const. for n = 1,(√3 + cos2 θ + cos θ)−1for n = 3.(4.39)From (4.38b) and the results in the previous sections, we expect a temperate ice regionto form above Z = 0 at Y > 0. Assuming the boundary of this region can be approximated93xyzs y s h eabedtemperate iceris zFigure 4.14: Geometry and definitions for the local solution. Ωt: temperate ice, Ωc: ice belowthe melting point, Ωb: bed, Γct: cold/temperate boundary.locally by a linear function of distance from the origin (Z ∼ R cos(φ)) allows only solutions withcontact angle φ = pi (Schoof, personal communication, included in appendix C.4). To get finiteseparation between the cold–temperate boundary and the bed as required by (4.38a) requiresus to go to higher order in the shape of the boundary. We therefore assume a power-law shapedboundary Z = A|Y |m with m > 1, as shown by the curve Γct in figure 4.14. A power-law shapeof the boundary implies that there is only a very slender part of the cold region separating Γctfrom the bed. In this slender region, vertical heat conduction is dominant, and we thereforehave to solve−∂2Θˆ∂Zˆ2= a0|Yˆ |−1. (4.40)The right hand side represents the heat production rate (4.39) near the bed. With the boundaryconditions Θ = 0 and ∂Θ/∂Z = 0 at Z = A|Y |m, this has solutionΘ = −a02|Y |−1 (Z − A|Y |m)2 , (4.41)which gives temperature and heat flux at Z = 0:Θ = −a02A|Y |2m−1 and −∂Θ∂Z= a0A|Y |m−1 = a0ARm−1 at Z = 0, Y < 0. (4.42)A cold–temperate boundary above the bed (A > 0) automatically ensures that the temperatureat the bed is below freezing through (4.42)1. In the bed, we solve (4.37b) in polar coordinates:1R∂∂R(R∂Θ∂R)+1R2∂2Θ∂θ2= 0 for pi < θ < 2pi (4.43)94subject to (4.42) at θ = pi and Θ˜ = 0 at θ = 2pi. We make the ansatz Θ ∼ Rνf(θ), leading tof = c1 sin(ν(θ − 2pi)) + c2 cos(ν(θ − 2pi)). (4.44)Applying the boundary condition (4.38b), Θ˜ = 0 at θ = 0, gives c2 = 0. The boundary condition(4.42)1 at θ = pi givesa02AR2m−ν−1 ∼ c1 sin(νpi). (4.45)The heat flux at Z = 0, Y < 0 is−∂Θ∂Z∣∣∣∣Z=0= −1R∂Θ∂θ∣∣∣∣θ=pi∼ −c1νRν−1 cos(νpi), (4.46)so thata0ARm−ν ∼ −c1ν cos(νpi). (4.47)With m > 1, R2m−ν−1  Rm−ν near the origin, the appropriate leading order forms of (4.45)and (4.47) arec1 sin(νpi) = 0 and − c1ν cos(νpi) = aoA (4.48)with m = ν. The first of these requires ν ∈ N and with m > 1, the smallest plausible choice ism = ν = 2, so that c1 = −a0A/2. Therefore, below the bedΘ = −a0A2R sin(2θ). (4.49)For diagnostic purposes the vertical heat flux is useful,−dΘdZ= a0AR cos θ. (4.50)Both a0 and A depend on the far field conditions. The factor a0 can be determined from theheat production alone, A must come from a local fit of the temperature field or the heat fluxto the local solution. Note that the constraint (4.14c)1 requires the cold–temperate boundaryto be above the bed, so A > 0, in which case (4.50) automatically fulfills the flux constraint(4.14d)2,95xyz xyxz sheasrihsng hflsowvpuewmargriFigure 4.15: Sample solution with additional boundary condition (4.36) for temperate ice andα = 20, ν = 0.5, Pe = 10. Panel a: temperature field, with contours at intervals of 0.1. Thebold red line marks T = 0. Panel b: close up of panel a, with contours at intervals of 0.01.The dashed green line is a parabolic fit of the T = 0 isotherm. Panel c: Comparison of thenumerical solution for temperature along the half circle in panel b (red circles) and the localsolution (4.49) (black line).Figure 4.15a shows a sample solution of the heat equation with a free cold–temperate bound-ary (4.36) for n = 3, α = 20, ν = 0.5 and Pe = 10. These are the same values as used in figure4.3c. A close up of the region around the origin is shown in panel b, where the green dashed linedescribes a parabolic fit Z = A|Y |2 to the cold–temperate boundary with A = 1.5. Solutionsalong the black half-circle are shown in panel c (red circles: numerical solutions, black line:local solution).With a parabolic cold–temperate boundary, resolving the slender cold ice region betweenthe cold–temperate boundary and the bed requires either progressively more flattened elementsin the mesh or element sizes shrinking with Y 2 towards the origin. With a finite numberof elements in the mesh and the requirement that the interior angles of the elements cannotshrink to zero, it is not possible to resolve the cold region up to the origin: the numericallyinferred cold–temperate boundary will reach the bed at finite distance to the origin, causing theconstraint (4.14c)1 to be violated.This is confirmed by our numerical results, see figure 4.16. Close to the origin for Y < 0, thevertical distance between the temperate region and the bed goes to zero as Z = A|Y |2. Thisline is plotted on top of a close up of the numerical mesh close to the origin in figure 4.16a. Thelowermost nodes with Z > 0 are marked by a purple bold line. Where this line lies above theZ = A|Y |2 curve of the temperate ice profile, we expect not to be able to resolve the slendercold region between Z = 0 and Z = A|Y |2. This is confirmed in figure 4.16b, which shows thattemperatures at the bed are at the melting point up to Y ≈ −0.009, thereby violating (4.14c)1.96b Z~|Y|2ax 10-3Figure 4.16: Numerical mesh and temperature along the bed for the solution shown in figure4.15. Panel a: Numerical mesh close to the origin. The dashed green line is the same as infigure 4.15b, the magenta solid line connects the lowermost nodes above the bed. Panel b:temperatures along the ice–bed interface at Z = 0 for the same range of Y as in panel a. Notethat the temperature constraint T < 0 for Y < 0 is typically violated where the dashed greenline of the local solution lies below the purple line and the cold region is not at least one elementwide.In the numerical determination of the migration velocity vm we therefore relax our require-ment (4.14c)1 that T < 0 at the bed for Y < 0. We still require numerically that the heat fluxcondition (4.14d)2 is not violated at the bed for Y > 0, and subject to this we choose vm so asto minimize the number of nodes with T = 0 at the bed for Y < 0. There is no guarantee thatthis procedure gives a solution, or if it does that this is solution is unique. However, minimizingthe number of nodes that violate the T < 0 constraint reproduces the local solution close to theorigin, while other choices of vm do not. It is important to note that relaxing the temperatureconstraint is necessitated by the limitations of the numerical model and that we assume thatboth inequality constraints still hold in the continuum limit.Figure 4.17 shows migration rates vm for n = 1 and different values of Pe. In our solutionswithout the constraint on temperature we observed a quasi-linear increase in the migrationvelocity with α, see figure 4.4b. Now it appears that the migration velocity increases sub-linearly with α.Sample temperature profiles are shown in figure 4.18, where the top row shows variationsin α and the bottom row shows variations in Pe. The most obvious feature of these solutions970 20 40 6000.511.522.5αv m  Pe=0 without (4.36)Pe=0 with (4.36)Pe=10 with (4.36)Pe=20 with (4.36)Figure 4.17: Margin migration rates vm with additional boundary condition (4.36) for n = 1,ν = 0.5 and the Pe values indicated. See figure 4.4b for comparison with solutions without theboundary condition (4.36)is that the temperate region expands as α increases. The temperature and heat flux conditions(4.36) we impose at the cold-temperature boundary do not allow any heat flow out of thetemperate region. Therefore the expansion of the temperate region limits the additional heatdissipation available to heat the bed when α is increased, compared with the previous casewith no temperature constraint applied inside the ice. This presumably explains why migrationvelocities are smaller now.As before, we are now interested in the asymptotic behaviour for large Peclet numbers Peand large heating parameters α. As we have changed the model through the free boundarycondition (4.36), we cannot expect the behaviour found in section 4.5 to hold any longer. Wewill therefore once again go through the exercise of identifying the asymptotic behaviour of vmfor α 1 and Pe 1. First we derive an auxiliary result.4.7.2 Heat production at large distances from the originFrom the profiles in figure 4.18 it is apparent that an increase in α leads to the bulging of alarge temperate area into the ice domain at Y < 0. Heat dissipation in the cold ice subdomaintherefore primarily happens at distances of O(1) from the slip–to–no–slip transition. Antici-pating that these distances become large in the limit of large α, we consider the behaviour ofa(Y, Z) for large negative Y in this section, limiting ourselves to the case of n = 1.98Z  a01  b  cT−10YZ  d−1 0 101Y  e−1 0 1 Y  f−1 0 1Figure 4.18: Sample profiles of temperate ice solutions for α and Pe of O(1). Panel a: α = 10,Pe = 0, panel b: α = 20, Pe = 0, panel c: α = 40, Pe = 0, panel d: α = 40, Pe = 10, panel e:α = 40, Pe = 20, panel f: α = 40, Pe = 40. ν = 0.5 and n = 3 in all simulations.The heat production a is determined from the solution of the mechanical boundary layermodel, equations (4.2) and (4.10). For n = 1, Schoof (2012) gives an analytical solution for a:a =exp(piY2)2√sinh2(piY2)+ sin2(piZ2) for n = 1. (4.51)When −Y  1, this reduces toa ∼ exp (−pi|Y |) (4.52)as shown figure 4.19a.For n > 1, no analytical solution of (4.2) and (4.10) exists to our knowledge. Solving (4.2)with χ = 0 in (4.3), we can find an approximate exponential fit to the resulting dissipation rateasa ∼ exp (−bpi|Y |) (4.53)with b ≈ 1.375, plotted as dashed black line in figure 4.19b.When χ is not zero and the omitted O(χ2) terms in (4.10) are taken into account, then thisasymptotic result holds only up to limited values of Y , beyond which heat dissipation due to thelateral inflow becomes dominant, see figure 4.19. Ordinarily, we neglect these very small residualterms, but they may become important in the physically unrealistic case that the temperateregion extends sufficiently far beyond the slip to no slip transition.99−10 −8 −6 −4 −2 0 2 410−1510−1010−5100Yheat production a  aχ=10−2, n=1χ=10−3, n=1χ=0, n=1analyticalexp(piY)−10 −8 −6 −4 −2 0 2 410−1510−1010−5100Yheat production a  bχ=10−2, n=3χ=10−3, n=3χ=0, n=3exp(1.375*piY)Figure 4.19: Heat production rates for n = 1 (panel a) and n = 3 (panel b). Numericallyobtained heat production rates are shown for different scales of the across-stream velocity V .For n = 1, the analytical solution (4.51) is plotted for comparison (red line in panel a).To avoid any sensitivity to the choice of χ, we prescribe the heat production in the form(4.52) for the case of large advection in section 4.7.4. The obvious caveat of this choice is thatour numerical results in this section are limited to n = 1. First we again consider the asymptoticbehaviour for large shear heating (α  1) and none or moderate advection (Pe  α) in thenext section.4.7.3 Asymptotic behaviour for α 1When the shear heating parameter α is much larger than the Peclet number Pe, we observe thetemperate region to extend significantly into the no slip region at Y < 0. This is illustrated infigure 4.20b, which shows a solution for α = 2.5× 106, Pe = 0, ν = 0.8, and n = 3. The red lineshows the T = 0 isotherm. We can use this result to guide our understanding of the temperaturefield and the migration velocity solution for large α. Consider the simplified geometry in figure4.20b. The horizontal position of the cold–temperate boundary is denoted by Yct(Z). Thenumbers in black circles indicate three distinct regions in which different leading order regimesapply. Region 1 corresponds to the region around the origin in which the local solution derivedpreviously applies. In region 2 the pseudo-advection term vm∂Θ/∂Y , that is introduced by themigration of the margin into the cold ice, dominates, and in this region temperature increaseswith increasing Y due to dissipation. In region 3 vertical conduction into the bed balances thesame pseudo-advection term as above. There are additional boundary layers that we do notdiscuss in detail below.100xyzxszxhzearY ct= vmπ (1−νν )21−ννxszing flnx onnsnzwvpu−1−ννac(b)eea TTFigure 4.20: Boundary layer geometry for temperate ice with small advection. Panel a: Bound-ary layer geometry for small advection, region 1 marks the position of the local solution, region2 marks the position of the outer solution, region 3 marks region where advection and verticaldiffusion are of leading order. Panel b: Solution with Elmer/Ice for α = 2.5 × 106, ν = 0.8,Pe = 0 and n = 3. Panel c: Analytic solution (4.58) in region 3.Our numerical results have indicated that vm increases with α, albeit slowly. We thereforeassume that at large α, vm is large and consequently omit the diffusive term in the heat equationfor region 2. The leading order version of (4.9) is thenvm∂Θ∂Yˆ= α exp (bpiYct) exp(bpiYˆ)for α 1, Yˆ < 0 (4.54)where Yˆ = Y + Yct and we have also used (4.52). Towards the ridge Θ still needs to satisfy(4.14b)1. Having omitted diffusive terms, we expect to be able to satisfy only the temperaturecondition Θ = 0 at the free boundary Yct. In order to satisfy the flux condition (4.36) a diffusiveboundary layer of thickness 1/vm is required.It is straightforward to solve equation (4.54) and we get the temperature fieldΘ =αbpivmexp(bpiYˆ)exp (bpiYct) + ΘR. (4.55)Putting Θ = 0 at Yˆ = 0 givesYct(Z) = −1bpilogα +1bpilog vm +1bpilog(bpi)−1bpilog(−ΘR(Z)). (4.56)101The right hand side of this expression still contains the migration velocity vm as an unknown.To determine it, we will now turn to the temperature field below the temperate ice region, whichwill give us another expression of Yct in terms of vm, α, ΘR, and b. Together with equation(4.56) this allows us to determine vm.Note that the argument above does not hold at small distances to the bed. The proximitybetween the ice, where heat is dissipated and the bed, where no heat is dissipated means thata diffusive boundary layer must form near the ice–bed interface, and Θ will not be given by(4.55) in this boundary layer. Although (4.55) gives an expression for Yct(0), we therefore do notexpect temperate ice to extend all the way to the bed at this location. Instead the temperateice is separated from the bed by a thin region of cold ice.By analogy with the above, we expect the temperature field below the bed and for Y <Yct(0) to be dominated by advection. With no heat production this implies that the far fieldtemperature field ΘR is simply advected without change. This can no longer hold in region 3,where there is temperate ice immediately above the bed, meaning for Y > Yct(0). From thelocal solution for region 1, we know that a Θ = 0–isotherm in the bed must connect to theorigin, which is not possible if we simply advect ΘR. Conduction of heat must therefore balanceadvection in region 3.Balancing advection and vertical diffusion, and assuming region 3 to be much longer thanit is deep, yieldsvm∂Θ∂Y−∂2Θ∂Z2= 0 for Z < 0. (4.57)As boundary conditions we assume that the proximity of temperate ice to the bed imposes afixed Θ = 0 at Z = 0 on region 3, while at the left hand end of region 3 (at Y = Yct(0)) weassume Θ to be given by Θ = ΘR(Z) = −1 − ν/(1 − ν)Z. This problems admits a similaritysolutionΘ = −erf(12√vmZ√Y − Yct(0))−ν1− νZ. (4.58)Having a Θ = 0–isotherm connect to the origin as required by the local solution in region 1implies ∂Θ/∂Z = 0 at Y = 0, which allows us to determine the extent to which temperate ice102xyzshye xyashyem aFigure 4.21: Asymptotic results for temperate ice with small advection. Panel a: n = 1, panelb: n = 3.protrudes towards the ridge:Yct(0) = −vmpi(1− νν)2. (4.59)Together with (4.56) we get an implicit equation for vm as function of α:logα = bvm(1− νν)2+ log(vm) + log bpi ≡ f(vm). (4.60)We thus expect a logarithmic increase of the margin migration rate vm and temperate ice extentYct(0) with α. This incidentally justifies our assumption of large vm and a shallow region 3,although in both cases the respective limit is approached logarithmically.Figure 4.21 shows corresponding results with Elmer/Ice for n = 1 (panel a) and n = 3 (panelb) plotted as lnα − f(vm) for different values of ν and Pe = 0. Convergence to zero indicatesagreement between the numerical solutions and the asymptotic solution. Note that in the limitof large α and vm we can approximately writevm ∼1b(ν1− ν)2logα. (4.61)4.7.4 Asymptotic behaviour for α 1 and Pe 1Next we extend the theory of the previous section to encompass large advection with Pe  1.If we compare 4.20b with 4.22b, it is clear that the introduction of advection reduces the extentof temperate ice, especially at greater heights above the bed. We attribute this to PeV ∂Θ/∂Y103being the dominant advection term, rather than vm∂Θ/∂Y , as V increases with height abovethe bed.In this case we have to consider four, rather than three distinct regions of heat transport,as depicted in figure 4.22a: it turns out that we can no longer legitimately ignore the diffusiveboundary layer. Again, region 1 is where we expect the local solution, given by equation(4.49), to hold. Region 2’ is again advection dominated, but in this case we assume that thedominant advection term is PeV ∂Θ/∂Y . The temperature field in the bed immediately belowthe temperate region should again by determined by vertical diffusion and advection due tomargin migration, vm∂Θ/∂Y . This is region 3. The region marked 4 is the diffusive boundarylayer near the bed.The leading order balance in the heat equation (4.9) in region 2’ isPeV∂Θ∂Y∼ α exp (−bpi|Y |) . (4.62)At the cold–temperate boundary the ice is at the melting point, Θ = 0, so thatYct = −1bpilogα +1bpilog(PeV (Z)) +1bpilog(bpi), (4.63)where we have assumed that Yct is large, so that V is given by the far field velocity profile andtherefore depends on Z only; this also justifies our neglect of vertical advection. Figure 4.22bshows a sample solution of the model and Yct from (4.63) as bold dashed line. Close to the bed,V is approximately V ∼ Z. This suggest that Yct goes to −∞ logarithmically as Z → 0.Yct does not actually go to infinity, because close to the bed vertical diffusion plays a leadingorder role. If we scale (4.9) with V = Pe−1/3Vˆ , Z = Pe−1/3Zˆ, we get the leading order balanceVˆ∂Θ∂Y−∂2Θ∂Zˆ2= αPe−2/3 exp (−bpi|Y |) , (4.64)assuming that vm  Pe1/3 (this effectively tells us how large we need to make the Pecletnumber before advection affects the solution at leading order). We do not attempt to solve(4.64) explicitly, but observe that we now have a scale for the distance above the bed at which(4.63) breaks down. Insisting that the right hand side of (4.64) is O(1), when Y ∼ Yct,max, the104xyzsyhsyesarY ct= vmπ (1−νν )21−ννyhsing flny onnhnswvpu( ,rtw )) mdWacc(bFigure 4.22: Boundary layer geometry for temperate ice with large advection. Panel a: Sketch ofthe geometry. Panel b: Sample solution with α = 105, ν = 0.9 and Pe = 104. The dashed blackline is the shape of the cold–temperate ice boundary from equation (4.63). Panel c: Comparisonof numerical results and asymptotic results for α = 105, ν = 0.9, and n = 1.maximum extent of temperate ice, we findYct,max ∼1bpilogα−1bpilog(Pe2/3). (4.65)The diffusive boundary layer has a counterpart in the bed, which may be considerablythicker, scaling in thickness as 1/√vm. Below this boundary layer, heat advection by vm∂Θ/∂Yis dominant and the far field profile ΘR is simply advected to the right until we reach Y = Yct,max.In the bed immediately below the temperate region to the right of Y = Yct,max, we expectthe leading temperature field to be again determined by (4.57), so that we can use (4.59), whichtogether with (4.65) gives the migration velocity as function of Pe, α and ν only:vm ∼1b(ν1− ν)2(logα−23log Pe), (4.66)in agreement with (4.61). Figure 4.22c shows the numerically calculated behaviour of vm for n =1, ν = 0.9 and a given by (4.52) for different values of α, together with the asymptotic solution(4.66). Note that we expect the best agreement between numerically calculated migration ratesand the asymptotic solution in the limit where α, Pe, and vm are all large.1054.8 DiscussionChapter 3 provided us with a model for the velocity and temperature fields in an ice streammargin, with the mechanical model largely decoupled from the thermal model. The mainpurpose of the thermal model is to provide a migration velocity for the position of the icestream margin. The rate of migration is dictated by the rate at which the temperature of thebed outside of the ice stream can be raised to the melting point by mechanical dissipation aswell as by conduction and advection of heat. Once the melting point is reached in the bed,sliding commences.Our boundary layer model formulates this as a travelling wave problem in which the ‘wavespeed’ needs to be chosen so that the bed does not reach the melting point prematurely outsideof the ice stream (in which case the migration is too slow) and does not freeze inside the icestream (in which case the migration is too fast). Mathematically, this is implemented throughan inequality constraint on the temperature on the ridge side of the margin and on heat flux onthe stream side of the margin (Schoof, 2012).Different aspects of this model can be found in the literature. The basic idea of solving athermomechanical model for an ice stream margin goes back to Raymond (1996) and Jacobsonand Raymond (1998), though with different boundary conditions and a different transversevelocity field from the one used here. Our model for the transverse flow field is a generalizedversion of the model in Barcilon and MacAyeal (1993), slightly modified due to different farfield boundary conditions. It describes the transition of the across-stream velocities from ashearing flow to a plug flow, and these velocities account for advection in the heat equation inthe shear margin. The velocity field in the downstream direction, whose gradient provides theheat production rate a, is a generalized version of the model in Schoof (2004) and Schoof (2012)and is very similar to that in Suckale et al. (2014).The novelty here lies in bringing these different ingredients together. In this chapter we havesolved the mechanical and thermal models numerically using the methods described in sections4.1–4.3. Using dimensional analysis (at the beginning of this chapter) has also allowed us to106show that the migration velocity depends on the far field conditions in stream and ridge throughdymdt=kρchsf(2Aτn+1s h2sk(Tm − Tb),(n+ 2)(n+ 1)ρcQrk,(Tb − T0)(Tm − T0)); (4.67)this is the dimensional version of the statement vm = f(α,Pe, ν). Our numerical work hasmostly been aimed at finding a representation of the function f .In 4.67, the first argument of f is α, a dimensionless strength of heat dissipation, drivenby the lateral shear stress τs imposed by the ice stream. hs is the thickness of the ice streamand Tm − Tb is the difference between the melting point and the far field bed temperaturein the ice ridge. The second argument is the Peclet number Pe, which can be understood asa dimensionless version of the ice flux Qr across the shear margin. The third argument, ν,represents the cooling of the bed through the surface temperature T0; alternatively one can alsounderstand this as a dimensionless geothermal heat flux.To determine how vm depends on α and Pe we made use of the fact that both are typicallylarge. This has led to (4.33) as a leading order version of (4.67), which gives the dimensionalrate at which an ice stream margin migrates outwards asdymdt= max(0,kρchs(Cαα + CpPe1/(1+β)))(4.68)with Cα ≈ 1, Cp ≈ −0.3 and β ≈ 0.27 for n = 3. This formula can only be applied to determinethe rate of ice stream widening, which should occur in the range where (4.68) is positive. Whenthe second argument in (4.68) is negative, we expect the ice stream to narrow. In this casedifferent physics apply, as discussed in section 3.7: inwards migration requires freezing of apotentially finite amount of water, which therefore dictates the time scale of this process.In deriving (4.68), we have treated regions of temperate ice with an extremely simplifiedenthalpy gradient method that assumes the same diffusivity throughout the entire domain (As-chwanden et al., 2012). There is no physical justification for this approach. We thereforeadapted our model and treated temperate ice in a more sophisticated manner in section 4.7.There we argued that the temperate and cold regions are separated by a free boundary, wherenot only the temperature is constrained to be at the melting point, but also the conductive heat1070 20 40 600204060Ws in kmdy m/dt in m/yeara0 20 40 6000.20 500 1000 15000204060Hs in mdy m/dt in m/year  bequation (4.68)equation (4.69)Figure 4.23: Dimensional migration velocities against ice stream width (panel a) and ice streamthickness (panel b). Inset in panel a shows a blow up for smaller values of ∂ym/∂t. We usedthe parameters in table 3.1 and τs = ρgHsWS/L, Qr = [a]W and Tb = 265 K.flux across the boundary must vanish. The resulting formulation is equivalent to the obstacleproblem in Zwinger et al. (2007).Going through the same exercise of finding the asymptotic behaviour of vm, we found analternative migration rate in the formdymdt=1bkρchsmax(0, logα−23log Pe)(4.69)with b ≈ 1.375 for n = 3 and b = 1 for n = 1 applicable to our temperate ice model.Figure 4.23 shows how (4.68) and (4.69) compare when varying the ice stream width (panela) or the ice stream thickness (panel b). The migration rates with the additional free boundarygiven by (4.36) are typically two orders of magnitude smaller than the migration rates without.This is due to the fact that the heat flux condition at the cold–temperate free boundary (−∇T ·n = 0) prohibits any heat transport from the temperate region to the cold region. Thereforeany heat that is dissipated in the temperate region cannot be used warm the bed on the ridgeside.Additionally, there is a subtle difference in the local behaviour around the origin betweenthe case with a free cold–temperate boundary and the case without. The local solution with afree boundary in section 4.7.1 requires substantial parts of the bed to reach the melting pointat the slip to no-slip transition. Specifically, the local solution (4.49) predicts that a meltingpoint isotherm emerges from the origin at an angle to the bed of θ = 3pi/2, see figure 4.15b.The melting point isotherm connects the ridge far field with the origin below the ice, which is108evident from figure 4.20b. This requires a substantial amount of heat to be provided for heatinglarge regions of the bed domain at Z < 0. There is no analogue to this requirement in the casewithout a free cold–temperate boundary in the ice domain (see section 4.3).Naively, one would assume that (4.69) gives better results than (4.68) in describing realworld migration rates, but observational data suggest that the opposite is the case: Echelmeyerand Harrison (1999) found Whillans ice stream to widen at a rate of 9.7 ± 1.1 m/year. Awidening of tributary B2 of the same ice stream has been reported at rates of 7 m/year byHarrison et al. (1998), while Hamilton et al. (1998) report a widening of 17 m/year for the sametributary. While the values we obtain with (4.68) lie in the range of these observations, thevalues obtained with (4.69) are clearly too small.We argue that this underestimation of migration velocities by (4.69) is probably due to theneglect of sub-temperate basal sliding (Cuffey et al., 2000; Fowler, 2013). The temperatureprofiles in 4.7 show that we expect the temperate ice region to intrude far into the no-slipregion at Y < 0 with an extremely thin region of cold ice separating the temperate ice and thesub-temperate bed, see for instance figures 4.20b and 4.22b. The proximity of the temperate iceto the bed implies that the temperatures at the bed are very close to the melting point, but westill assume that there is no slip at the bed. This is hardly realistic, and a regularized transitionfrom free slip to no slip as explored in section 4.6 is called for (Fowler, 2013). Combining ourtreatment of sub-temperate slip with our temperate ice model is left to future work. We donote, however, that while a temperate ice region can ‘hide’ much of the heat dissipation in theice from the bed, this is not true for heat dissipated through sliding at the bed. We anticipatethat the latter gives rise to migration velocities much closer to those observed.Another simplification we have made in this chapter is to assume that the viscosity is effec-tively temperature-independent. This allowed us to solve the mechanical model independently ofthe temperature model, but neglected the potential coupling between a temperature-dependentviscosity and the mechanical model, in particular the heat production. This effect is believedto be significant in ice stream margins, where steep temperature gradients exist (Truffer andEchelmeyer, 2003).The viscosity decreases exponentially with temperature (Paterson, 1994), which should leadto a relative hardening of the ice in the cold upper parts of the ice stream margin and further109YZa−1 0 101Y  b−5 0 5 10 15 20 25T−10Figure 4.24: Temperature field in the boundary layer and the adjacent ice stream. Panel a:Close up of the margin, panel b: ice stream domain. We used α = 20, Pe = 0, n = 3, ν = 0.5,and f = 1/25.softening in the temperate ice regions. However, incorporating these effects into our modelwould likely only affect the results for α and Pe of O(1). For α  1, the migration velocityis controlled by processes in a small diffusive boundary layer, so that using the local value ofA(T ) in equation (4.68) should be a good first approximation of the migration velocity even ifthermo-mechanical coupling is taken into account.A further effect neglected here is the weakening of ice through rescrystallization under largestrain (Kamb, 1973; Budd and Jacka, 1989). To parameterize such effects, which occur on thegrain-size level, so-called ‘enhancement factors’ are commonly introduced in the constitutiverelationship (3.2a) (Ma et al., 2010). While it would be straightforward to extend our model tosuch a parameterization, an in-depth treatment of these effects would be more desirable, but isbeyond the scope of this work.The temperature profiles we have shown above, for instance in figures 4.5 and 4.18, all lookvery warm, with large regions of temperate ice extending towards Y → ∞. The reason forthis is that we have prescribed a constant lateral shear stress towards the ice stream, and theheat production becomes constant in this limit (see figure 4.2c). In a real ice stream the lateralshear stress would decrease quickly away from the margin, and consequently so would the heatproduction. We therefore expect the temperatures to decrease towards the ice stream center.To illustrate this behaviour, we change the boundary layer model for the downstream flowand replace (4.2a) with∂∂Y(η∂U∂Y)+∂∂Z(η∂U∂Z)= −f (4.70)and discard (4.2c)3. By choosing f = 1/Ws with Ws the half width of the ice stream, we get a110lateral shear stress τxy = 1 in the ice stream margin, which corresponds to the value of τs = 1used in the other examples shown in this chapter. Figure 4.24 shows a temperature profilecalculated with this altered model. Around the margin at Y = 0 the temperature field looksvery similar to those calculated above, but towards the ice stream center (in this case at Y = 25)temperatures drop and a uniform temperature gradient is established.1115 | ConclusionsThe aim of this thesis has been to achieve a better understanding of the flow of ice streams andthe migration of their margins. In the introduction, we identified two specific research goals:(i) to better understand the physical processes that lead to the migration of ice stream mar-gins, and(ii) to derive a parametrization of these processes that can be used to better quantify icestream dynamics in continental-scale ice sheet models.This thesis focuses almost exclusively on the outwards migration of ice stream margins. Inthe next section, we summarize our results on ice stream widening. We comment on inwardsmigration in section 5.2, and we conclude with suggestions for future work in section Outwards migration of ice stream marginsWe have considered two different mechanisms for ice stream widening: drainage-driven wideningand widening through heat dissipation in the ice stream margin. In chapter 2, we considereddrainage-driven widening and how it is affected by basal freeze-on. We found that subglacialdrainage can drive ice stream widening, but only if the bed in the surrounding, slowly movingregions is at least partially unfrozen.Due to a positive feedback between fast ice flow, heat dissipation, and melt water production,the ice stream bed is well lubricated with high subglacial water pressures, while water pressuresare low at the bed of the slowly moving ice ridges bordering the ice stream. High water pressuregradients are therefore expected across the ice stream margin, and it is conceivable that thesegradients can drive water flow along the bed from the ice stream towards the ice ridge. Usinga simple model for ice stream flow over a bed that allows for the lateral movement of water, aswell as freezing and melting of the ice stream bed, we found that water seeping into the regionsof low water pressure can indeed lubricate the bed in the regions of slow flow. This can lead toan increase in ice velocities, and thereby cause the ice stream to widen.112However, because heat dissipation is low in the ice ridge, we also expect freezing of the bedin these regions. This freezing has profound consequences: a completely frozen sediment consti-tutes a thermal barrier through which we assume subglacial water cannot penetrate. This mech-anism operates independently of other possible processes that might stop ice stream widening.For example, dynamic thinning of the ice stream can lead to the formation of ice ridges. Thiscreates inwards gradients in normal stresses at the bed that can suppress outwards drainage,thereby preventing further widening (Kyrke-Smith et al., 2014).Observations suggest that the beds of ice ridges are indeed often frozen (Bentley et al., 1998;Catania et al., 2003), but that widening of ice streams into these regions is nevertheless possible(Stephenson and Bindschadler, 1988; Fahnestock et al., 2000; Conway et al., 2002; Cataniaet al., 2012). Drainage-driven ice stream widening therefore likely plays only a minor role inthe dynamics of the Siple Coast ice streams.An alternative mechanism for ice stream widening is through heat dissipation in the icestream margins. There the fast ice stream flow slows down over relatively short distances to thespeed of the surrounding ice ridges. This results in high strain rates and high heat dissipationrates in the margins, and through lateral heat transport this can melt the frozen bed adjacentto the ice stream and thereby cause widening (Jacobson and Raymond, 1998; Schoof, 2012).A complication to this picture is introduced by lateral advection of ice from the ridges to thestream, which leads to a cooling of the ice stream margin. Studies of stationary ice streammargins have identified this as a potentially significant process (Jacobson and Raymond, 1998;Suckale et al., 2014). However, these models used ad hoc parameterizations of the advectiveflow field.We show in chapter 2 that widening due to heat dissipation in the ice stream margin cannotbe modelled with the thin film models typically employed in continental-scale ice sheet simula-tions. To capture the relevant physics, processes at the length scale of the ice stream thicknessmust be accounted for. To address this issue, we derived a boundary layer model for the icestream margin in chapter 3.The boundary layer model has two functions: it provides the coupling between the ice streamand the ice ridge, and its solution describes how the boundary between these two regions evolvesin time (in other words, how fast the ice stream margin migrates). The mechanical component113of this model describes the transition of the shearing flow in the ice ridge to the plug flow inthe ice stream, and it provides the heat production rates and advection velocities necessary tosolve the thermal component of model. The latter takes the form of a travelling wave problemin which the migration of the margin appears as an advection rate. This allows us to treatthe rate vm at which the margin migrates as a scalar parameter that needs to be determined,somewhat like an eigenvalue.The temperature field in the ice stream margin is primarily influenced by heat dissipationthrough lateral shear and by the inflow of cold ice from the margins. The effect of both processescan be described through two dimensionless groups (α and Pe), which are linked to the far fieldproperties in the ice ridge and the ice stream, so that they can be calculated without knowledgeof the boundary layer dynamics.The outwards migration of an ice stream margin effectively introduces an advective coolingterm into the heat equation, as the margin moves into the colder ice of the ridge. The velocityat which the margin migrates outwards must therefore be slow enough to prevent freezing onthe ice stream side of the bed, and fast enough to prevent temperatures above the meltingpoint on the ice ridge side. These constraints were first introduced by Schoof (2012) in asimplified model that treats ice as a Newtonian fluid and neglects advection. Our numericalresults suggest that these constraints still uniquely determine the migration speed when ice istreated as a non-Newtonian fluid and advection is taken into account.We reproduce the original results by Schoof (2012), but only if the numerical mesh aroundthe transition between no slip and free slip is highly refined. It is therefore unlikely that theprocesses that lead to margin migration can be captured with continental-scale models, even ifthey are ‘non-shallow’. Instead, a parameterization is necessary to model ice stream dynamicsin continental-scale simulations: we need vm as a function of the shear heating parameter αand the advection parameter Pe. This function is easiest to compute in the limit of large α andPe, which also happens to be physically realistic. By considering the local behaviour of velocityand temperature near the no-slip to free slip transition, we were able to derive equation (4.68)for vm as a simple increasing function of α and decreasing function of Pe.Equation (4.68) can be used as a parameterization of the migration velocity to model themigration of ice stream margins in large scale models. It is similar to a Stefan condition114describing the motion of the free boundary between ice stream and ice ridge.5.2 Inwards migration of ice stream marginsMost of the work in this thesis has focused on ice stream widening and equation (4.68) onlydescribes the rate at which an ice stream margin migrates outwards, but not inwards. Inwardsmigration of ice stream margins involves different physics, as our results in chapter 2 illustrate:freezing of the entire vertical water column in the ice stream bed is necessary, while we assumethat widening only requires an infinitesimal amount of melting (in other words heating the ice–bed interface to the melting point). The amount of water that needs to be frozen depends notonly on the local energy balance, but also on the supply of melt water from adjacent regions.The boundary layer model may not be necessary to capture inwards migration, since ourresults in chapter 2 show that it is possible to model the inwards migration and shutdown of anice stream with a depth-integrated model. We found that the inwards migration of ice streammargins, which in our model inevitably leads to a shutdown of the ice stream, can result frominsufficient heat dissipation over the entire ice stream. Owing to the positive feedback describedabove and in the introduction, heat dissipation in the ice stream margin and at the ice streambed both increase with ice stream velocity. If an ice stream is not wide enough or thick enoughto dissipate sufficient heat to balance conductive cooling of the bed, the ice stream bed willeventually experience freezing everywhere. Once the entire bed is frozen, no more sliding ispossible and the ice stream shuts down.5.3 Future workThe thin film model in chapter 2 has made a number of simplifying assumptions, and futurework should explore the implications of relaxing these assumptions. In particular, we foundthat no steady state velocity solutions exist on a bed that is impermeable to subglacial waterflow. However, we cannot guarantee that this result will hold in a more realistic representationof heat and ice fluxes. For example, we have considered an infinitely long domain, so that thereis no net transport of ice or subglacial water in the axial (downstream) direction. This excludes115the possibility of forming ice ridges through dynamic thinning of the ice stream, so that thereis no transverse transport of ice.Incorporating ice stream widening with equation (4.68) in large scale models is not necessar-ily straightforward and two main challenges have to be overcome to use this parameterization.First, (4.68) gives the margin migration rate perpendicular to the direction of the ice streambulk flow, with the parameters α and Pe evaluated at the boundary between ice stream and iceridge. This boundary does not necessarily align with the (often fixed) grid or mesh of continen-tal scale models. Secondly, migration rates of a few meters per year will likely result in changesin the ice stream width over a single time step that are smaller than typical mesh sizes. Use ofequation (4.68) will therefore require methods for the treatment of moving boundaries (see e.g.,Shyy et al., 2012).A more immediate direction for future research is the interaction between sub-temperatesliding and the formation of temperate ice. In deriving equation (4.68) we used a simplifiedenthalpy gradient model for ice whose temperature is above the melting point. This modeldoes not have a rigorous physical justification. A more sophisticated treatment of temperate iceincorporating a free boundary between temperate and cold ice has been used in section 4.7. Atthis boundary, the temperature is constrained to the melting point and there are no conductiveheat fluxes across the boundary. Determining the margin migration speed in the limit whereboth α and Pe are large leads to equation (4.69). The formula therein predicts migration ratesthat are much smaller than those observed in reality. It is likely that this is because we neglectsub-temperate sliding on the ice ridge side of the margin, and future work should take this intoaccount.To further our understanding of ice stream dynamics, some open questions still need tobe addressed in future research. For instance, we have only considered the margin migrationin a fully evolved ice stream, but the physics of ice stream onset at their tributaries are stilllargely unknown, which also poses questions about the initial formation of ice streams througha putative instability mechanism. Observations suggest that geological factors may play a rolein the locations of the tributaries. Peters et al. (2006) found the onsets of two tributaries ofthe Siple Coast ice streams to coincide with patches of soft sediment at the bed, while thesurrounding, slower moving ice is underlain by harder bedrock. However, this does not mean116that we can completely disregard thermal transitions at the bed, although it is possible thatthe geologic and thermal boundaries do not coincide. If a thermal transition is present in anice stream onset region, it is conceivable that a similar boundary layer problem to the one westudied here might hold the key to study the evolution of ice stream onset regions. An analysisof this situation is however beyond the scope of this thesis and should be addressed in futurework.117BibliographyAlley, R. B., Blankenship, D., Bentley, C., and Rooney, S. (1987). Till beneath ice stream B, 3,Till deformation: Evidence and implications. Journal of Geophysical Research, 92:8921–8930,doi:10.1029/JB092iB09p08921.Anandakrishnan, S. and Alley, R. B. (1997). Stagnation of ice stream C, West Antarctica bywater piracy. Geophysical Research Letters, 24(3):265–268, doi:10.1029/96GL04016.Anandakrishnan, S., Blankenship, D. D., Alley, R. B., and Stoffa, P. L. (1998). Influence ofsubglacial geology on the position of a West Antarctic ice stream from seismic observations.Nature, 394:62–65, doi:10.1038/27889.Aschwanden, A. and Blatter, H. (2009). Mathematical modeling and numerical sim-ulation of polythermal glaciers. Journal of Geophysical Research, 114(F1):F01027,doi:10.1029/2008JF001028.Aschwanden, A., Bueler, E., Khroulev, C., and Blatter, H. (2012). An enthalpy formulation forglaciers and ice sheets. Journal of Glaciology, 58:441–457.Bamber, J. L., Vaughan, D. G., and Joughin, I. (2000). Widespread Complex Flow in the Interiorof the Antarctic Ice Sheet. Science, 287:1248–1250, doi:10.1126/science.287.5456.1248.Barcilon, V. and MacAyeal, D. R. (1993). Steady flow of a viscous ice stream across a no-slip/free-slip transition at the bed. Journal of Glaciology, 39:167–185.Bennett, M. R. (2003). Ice streams as the arteries of an ice sheet: their mechanics, stabilityand significance. Earth Science Reviews, 61:309–339, doi:10.1016/S0012–8252(02)00130–7.Bentley, C. R., Lord, N., and Liu, C. (1998). Radar reflections reveal a wet bed beneath stagnantIce Stream C and a frozen bed beneath ridge BC, West Antarctica. Journal of Glaciology,44:149–156.118Bindschadler, R., Chen, X., and Vornberger, P. (2000). The onset area of Ice Stream D, WestAntarctica. Journal of Glaciology, 46:95–101.Bindschadler, R. and Vornberger, P. (1998). Changes in the West Antarctic IceSheet Since 1963 from Declassified Satellite Photography. Science, 279(5351):689–692,doi:10.1126/science.279.5351.689.Bindschadler, R., Vornberger, P., Blankenship, D., Scambos, T., and Jacobel, R. (1996). Surfacevelocity and mass balance of Ice Streams D and E, West Antarctica. Journal of Glaciology,42:461–475.Bindschadler, R. A. and Scambos, T. A. (1991). Satellite-image-derived velocity field of anAntarctic ice stream. Science, 252:242–246, doi:10.1126/science.252.5003.242.Bindschadler, R. A., Stephenson, S. N., MacAyeal, D. R., and Shabtaie, S. (1987). Ice dynamicsat the mouth of ice stream B, Antarctica. Journal of Geophysical Research, 92:8885–8894,doi:10.1029/JB092iB09p08885.Blankenship, D. D., Bentley, C. R., Rooney, S. T., and Alley, R. B. (1986). Seismic mea-surements reveal a saturated porous layer beneath an active Antarctic ice stream. Nature,322:54–57, doi:10.1038/322054a0.Blankenship, D. D., Bentley, C. R., Rooney, S. T., and Alley, R. B. (1987). Till beneath icestream B. I - Properties derived from seismic travel times. II - Structure and continuity. III -Till deformation - Evidence and implications. IV - A coupled ice-till flow model. Journal ofGeophysical Research, 92:8903–8940, doi:10.1029/JB092iB09p08903.Bond, G., Heinrich, H., Broecker, W., Labeyrie, L., McManus, J., Andrews, J., Huon, S.,Jantschik, R., Clasen, S., Simet, C., Tedesco, K., Klas, M., Bonani, G., and Ivy, S. (1992).Evidence for massive discharge of icebergs into the North Atlantic Ocean during the lastglacial period. Nature, 360:245–249, doi:10.1038/360245a0.Bougamont, M., Tulaczyk, S., and Joughin, I. (2003a). Numerical investigations of the slow-down of Whillans Ice Stream, West Antarctica: is it shutting down like Ice Stream C? Annalsof Glaciology, 37:239–246.119Bougamont, M., Tulaczyk, S., and Joughin, I. (2003b). Response of subglacial sediments tobasal freeze-on 2. Application in numerical modeling of the recent stoppage of Ice Stream C,West Antarctica. Journal of Geophysical Research, 108(B4):2223, doi:10.1029/2002JB001936.Boulton, G. S. and Hindmarsh, R. C. A. (1987). Sediment deformation beneath glaciers:rheology and geological consequeces. Journal of Geophysical Research, 92:9059–9082,doi:10.1029/JB092iB09p09059.Broecker, W., Bond, G., Klas, M., Clark, E., and McManus, J. (1992). Origin of the northernAtlantic’s Heinrich events. Climate Dynamics, 6:265–273.Budd, W. and Jacka, T. (1989). A review of ice rheology for ice sheet modelling. Cold RegionsScience and Technology, 16(2):107–144.Bueler, E. and Brown, J. (2009). The shallow shelf approximation as a sliding law in athermomechanically coupled ice sheet model. Journal of Geophysical Research, 114:F03008,doi:10.1029/2008JF001179.Calov, R., Greve, R., Abe-Ouchi, A., Bueler, E., Huybrechts, P., Johnson, J. V., Pattyn, F.,Pollard, D., Ritz, C., Saito, F., et al. (2010). Results from the Ice-Sheet Model Intercom-parison Project–Heinrich Event INtercOmparison (ISMIP HEINO). Journal of Glaciology,56(197):371–383.Catania, G., Hulbe, C., Conway, H., Scambos, T. A., and Raymond, C. F. (2012). Variabilityin the mass flux of the Ross ice streams, West Antarctica, over the last millennium. Journalof Glaciology, 58:741–752.Catania, G. A., Conway, H. B., Gades, A. M., Raymond, C. F., and Engelhardt, H. (2003). Bedreflectivity beneath inactive ice streams in West Antarctica. Annals of Glaciology, 36(1):287–291.Catania, G. A., Scambos, T. A., Conway, H., and Raymond, C. F. (2006). Sequen-tial stagnation of Kamb Ice Stream, West Antarctica. Geophysical Research Letters,33:L14502,10.1029/2006GL026430.120Christoffersen, P., Bougamont, M., Carter, S. P., Fricker, H. A., and Tulaczyk, S. (2014).Significant groundwater contribution to Antarctic ice streams hydrologic budget. GeophysicalResearch Letters, 41(6):2003–2010, doi:10.1002/2014GL059250.Christoffersen, P. and Tulaczyk, S. (2003). Response of subglacial sediments to basal freeze-on 1.Theory and comparison to observations from beneath the West Antarctic Ice Sheet. Journalof Geophysical Research, 108:2222, doi:10.1029/2002JB001935.Chugunov, V. A. and Wilchinsky, A. V. (1996). Modelling of a marine glacier and ice-sheet-ice-shelf transition zone based on asymptotic analysis. Annals of Glaciology, 23:59–67.Church, J., Clark, P., Cazenave, A., Gregory, J., Jevrejeva, S., Levermann, A., Merrifield, M.,Milne, G., Nerem, R., Nunn, P., Payne, A., Pfeffer, W., Stammer, D., and Unnikrishnan, A.(2013). Climate Change 2013: The Physical Science Basis. Contribution of Working Group Ito the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, chapterSea level change. Cambridge University Press, Cambridge, United Kingdom and New York,NY, USA.Clarke, G. K. C., Nitsan, U., and Paterson, W. S. B. (1977). Strain heating andcreep instability in glaciers and ice sheets. Reviews of Geophysics, 15(2):235–247,doi:10.1029/RG015i002p00235.Conway, H., Catania, G., Raymond, C. F., Gades, A. M., Scambos, T. A., and Engelhardt,H. (2002). Switch of flow direction in an Antarctic ice stream. Nature, 419:465–467,doi:10.1038/nature01081.Creyts, T. T. and Schoof, C. G. (2009). Drainage through subglacial water sheets. Journal ofGeophysical Research, 114:F04008, doi:10.1029/2008JF001215.Cuffey, K., Conway, H., Hallet, B., Gades, A., and Raymond, C. (1999). Interfacial water inpolar glaciers and glacier sliding at −17 oC. Geophysical Research Letters, 26(6):751–754.Cuffey, K. M., Conway, H., Gades, A. M., Hallet, B., Lorrain, R., Severinghaus, J. P., Steig,E., Vaughn, B., and White, J. W. C. (2000). Entrainment at cold glacier beds. Geology,28(4):351–354.121Echelmeyer, K. and Zhongxiang, W. (1987). Direct observation of basal sliding and deformationof basal drift at sub-freezing temperatures. Journal of Glaciology, 33(113):83–98.Echelmeyer, K. A. and Harrison, W. D. (1999). Ongoing margin migration of Ice Stream B,Antarctica. Journal of Glaciology, 45:361–369.Echelmeyer, K. A., Harrison, W. D., Larsen, C., and Mitchell, J. E. (1994). The role of themargins in the dynamics of an active ice stream. Journal of Glaciology, 40:527–538.Engelhardt, H. and Kamb, B. (1997). Basal hydraulic system of a West Antarctic ice stream:constraints from borehole observations. Journal of Glaciology, 43:207–230.Engelhardt, H. and Kamb, B. (1998). Basal sliding of Ice Stream B, West Antarctica. Journalof Glaciology, 44:223–230.England, A. (1971). Complex Variable Methods in Elasticity. J. Wiley & Sons, Ltd., London.Fahnestock, M. A., Scambos, T. A., Bindschadler, R. A., and Kvaran, G. (2000). A millenniumof variable ice flow recorded by the Ross Ice Shelf, Antarctica. Journal of Glaciology, 46:652–664.Favier, L., Gagliardini, O., Durand, G., Zwinger, T., et al. (2012). A three-dimensional fullstokes model of the grounding line dynamics: effect of a pinning point beneath the ice shelf.The Cryosphere, 6:101–112.Fowler, A. C. (1986). Sub-temperate basal sliding. Journal of Glaciology, 32:3–5.Fowler, A. C. (1987). Sliding with cavity formation. Journal of Glaciology, 33:255–267.Fowler, A. C. (2001). Modelling the flow of glaciers and ice sheets. Continuum Mechanics andApplications in Geophysics and the Environment, Springer, 2001, pages 201–221.Fowler, A. C. (2013). The motion of ice stream margins. Journal of Fluid Mechanics, 714:1–4.Fowler, A. C. and Johnson, C. (1996). Ice-sheet surging and ice-stream formation. Annals ofGlaciology, 23:68–73.122Fowler, A. C. and Larson, D. A. (1978). On the Flow of Polythermal Glaciers. I. Model andPreliminary Analysis. Proceedings of the Royal Society of London. Series A, Mathematicaland Physical Sciences, 363(1713):217–242.Fowler, A. C. and Schiavi, E. (1998). A theory of ice-sheet surges. Journal of Glaciology,44(146):104–118.Fricker, H. A., Scambos, T., Bindschadler, R., and Padman, L. (2007). An Active Sub-glacial Water System in West Antarctica Mapped from Space. Science, 315:1544–1548,doi:10.1126/science.1136897.Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., Fleurian, B. d., Greve,R., Malinen, M., Mart´ın, C., R˚aback, P., et al. (2013). Capabilities and performance ofElmer/Ice, a new generation ice-sheet model. Geoscientific Model Development, 6:1299–1318,doi:10.5194/gmd–6–1299–2013.Geuzaine, C. and Remacle, J.-F. (2009). Gmsh: A 3-d finite element mesh generator withbuilt-in pre-and post-processing facilities. International Journal for Numerical Methods inEngineering, 79(11):1309–1331.Gray, L., Joughin, I., Tulaczyk, S., Spikes, V. B., Bindschadler, R., and Jezek, K. (2005).Evidence for subglacial water transport in the West Antarctic Ice Sheet through three-dimensional satellite radar interferometry. Geophysical Research Letters, 32(3):L03501,doi:10.1029/2004GL021387.Hamilton, G. S., Whillans, I. M., and Morgan, P. J. (1998). First point measurements ofice-sheet thickness change in Antarctica. Annals of Glaciology, 27:125–129.Harrison, W. D., Echelmeyer, K. A., and Larsen, C. F. (1998). Measurement of temperaturein a margin of Ice Stream B, Antarctica: implications for margin migration and lateral drag.Journal of Glaciology, 44:615–624.Heinrich, H. (1988). Origin and consequences of cyclic ice rafting in the Northeast AtlanticOcean during the past 130,000 years. Quaternary Research, 29:142–152.123Hewitt, I. J. and Fowler, A. C. (2008). Partial melting in an upwelling mantle column. Proceed-ings of the Royal Society A: Mathematical, Physical and Engineering Science, 464(2097):2467–2491.Hindmarsh, R. (1997). Deforming beds: Viscous and plastic scales of deformation. QuaternaryScience Reviews, 16:1039–1056.Hindmarsh, R. C. A. (2006). Stress gradient damping of thermoviscous ice flow instabilities.Journal of Geophysical Research, 111:B10, doi:10.1029/2005JB004019.Hindmarsh, R. C. A. (2009). Consistent generation of ice-streams via thermo-viscous in-stabilities modulated by membrane stresses. Geophysical Research Letters, 36:L06502,doi:10.1029/2008GL036877.Hintermu¨ller, M., Ito, K., and Kunisch, K. (2002). The Primal-Dual Active Set Strategy as aSemismooth Newton Method. SIAM Journal on Optimization, 13(3):865–888.Hooke, R. L., Hanson, B., Iverson, N. R., Jansson, P., and Fischer, U. H. (1997). Rheology oftill beneath Storglacia¨ren, Sweden. Journal of Glaciology, 43:172–179.Hulbe, C. and Fahnestock, M. (2007). Century-scale discharge stagnation and reactivationof the Ross ice streams, West Antarctica. Journal of Geophysical Research, 112:F03S27,doi:10.1029/2006JF000603.Hulbe, C. L. and Whillans, I. M. (1997). Weak bands within Ice Stream B, West Antarctica.Journal of Glaciology, 43:377–386.Hutter, K. (1983). Theoretical Glaciology. D. Reidel Publishing Company/ Terra ScientificPublishing Company.Hutter, K., Blatter, H., and Funk, M. (1988). A model computation of moisture con-tent in polythermal glaciers. Journal of Geophysical Research, 93(B10):12205–12214,doi:10.1029/JB093iB10p12205.Hutter, K. and Olunloyo, V. O. S. (1980). On the distribution of stress and velocity in an icestrip, which is partly sliding over and partly adhering to its bed, by using a Newtonian viscous124approximation. Proceedings of the Royal Society of London. A. Mathematical and PhysicalSciences, 373(1754):385–403, doi:10.1098/rspa.1980.0155.Iken, A. and Bindschadler, R. A. (1986). Combined measurements of subglacial water pressureand surface velocity of Findelengletscher, Switzerland: conclusions about drainage systemand sliding mechanism. Journal of Glaciology, 32:101–119.Iverson, N. R., Baker, R. W., LeB Hooke, R., Hanson, B., and Jansson, P. (1999). Couplingbetween a glacier and a soft bed: I. A relation between effective pressure and local shearstress determined from till elasticity. Journal of Glaciology, 45:31–40.Iverson, N. R., Hooyer, T. S., and Baker, R. W. (1998). Ring-shear studies of till deforma-tion: Coulomb-plastic behavior and distributed strain in glacier beds. Journal of Glaciology,44:634–642.Jackson, M. and Kamb, B. (1997). The marginal shear stress of Ice Stream B, West Antarctica.Journal of Glaciology, 43:415–426.Jacobson, H. P. and Raymond, C. F. (1998). Thermal effects on the location of ice streammargins. Journal of Geophysical Research, 103:12111–12122, doi:10.1029/98JB00574.Joughin, I., Tulaczyk, S., Bindschadler, R., and Price, S. F. (2002). Changes in west Antarcticice stream velocities: Observation and analysis. Journal of Geophysical Research, 107:2289,doi:10.1029/2001JB001029.Kamb, B. (1973). Experimental recrystallization of ice under stress.Kamb, B. (1991). Rheological nonlinearity and flow instability in the deforming bedmechanism of ice stream motion. Journal of Geophysical Research, 96:16585–16595,doi:10.1029/91JB00946.Kamb, B. (2001). Basal zone of the West Antarctic ice streams and its role in lubrication oftheir rapid motion. In Alley, R. B. and Bindschadler, R. A., editors, The West AntarcticIce Sheet: Behaviour and Environment, Antarctic Research Series vol. 77, pages 157–199.American Geophysical Union.125Kyrke-Smith, T. M., Katz, R. F., and Fowler, A. C. (2014). Subglacial hydrology and theformation of ice streams. Proceedings of the Royal Society A: Mathematical, Physical andEngineering Science, 470(2161):doi:10.1098/rspa.2013.0494.Le Brocq, A. M., Payne, A. J., Siegert, M. J., and Alley, R. B. (2009). A subglacial water-flowmodel for West Antarctica. Journal of Glaciology, 55:879–888.Ma, Y., Gagliardini, O., Ritz, C., Gillet-Chaulet, F., Durand, G., and Montagnat, M. (2010).Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flowmodel. Journal of Glaciology, 56(199):805–812.MacAyeal, D. R. (1989). Large-scale ice flow over a viscous basal sediment - Theory andapplication to ice stream B, Antarctica. Journal of Geophysical Research, 94:4071–4087,doi:10.1029/JB094iB04p04071.MacAyeal, D. R. (1993). Binge/purge oscillations of the Laurentide Ice Sheet as a cause of theNorth Atlantic’s Heinrich events. Paleoceanography,, 8(6):775–784, doi:10.1029/93PA02200.Moore, P. L., Iverson, N. R., and Cohen, D. (2010). Conditions for thrust faulting in a glacier.Journal of Geophysical Research, 115:F02005, doi:10.1029/2009JF001307.Morland, L. W. and Johnson, I. R. (1980). Steady motion of ice sheets. Journal of Glaciology,25:229–246.Muszynski, I. and Birchfield, G. E. (1987). A coupled marine ice-stream-ice-shelf model. Journalof Glaciology, 33:3–15.Nereson, N. A. and Raymond, C. F. (2007). Radar investigations of Antarctic ice streammargins, Siple Dome, 1998. Technical report, Boulder, Colorado USA: National Snow andIce Data Center. http://dx.doi.org/10.7265/N52B8VZP.Paterson, W. S. B. (1994). The Physics of Glaciers. Elsevier, Oxford.Payne, A. J. and Dongelmans, P. W. (1997). Self-organization in the thermomechanical flow ofice sheets. Journal of Geophysical Research, 102:12219–12234, doi:10.1029/97JB00513.126Peters, L. E., Anandakrishnan, S., Alley, R. B., Winberry, J. P., Voigt, D. E., Smith, A. M.,and Morse, D. L. (2006). Subglacial sediments as a control on the onset and location oftwo Siple Coast ice streams, West Antarctica. Journal of Geophysical Research, 111:B01302,doi:10.1029/2005JB003766.Pralong, A. and Funk, M. (2005). Dynamic damage model of crevasse opening and application toglacier calving. Journal of Geophysical Research, 110(B1):B01309, doi:10.1029/2004JB003104.Raymond, C. (1996). Shear margins in glaciers and ice sheets. Journal of Glaciology, 42:90–102.Raymond, C. F. (2000). Energy balance of ice streams. Journal of Glaciology, 46:665–674.Rempel, A. W. (2009). Effective stress profiles and seepage flows beneath glaciers and ice sheets.Journal of Glaciology, 55(191):431–443.Retzlaff, R. and Bentley, C. R. (1993). Timing of stagnation of Ice Stream C, West Antarctica,from short-pulse radar studies of buried surface crevasses. Journal of Glaciology, 39:553–561.Retzlaff, R., Lord, N., and Bentley, C. R. (1993). Airborne-radar studies: Ice Streams A, B andC, West Antarctica. Journal of Glaciology, 39:495–506.Rice, J. (1967). Stresses due to a sharp notch in a work-hardening elastic-plastic material loadedby longitudinal shear. Journal of Applied Mechanics, 34(2):287–298.Rice, J. R. (1968). A path independent integral and the approximate analysis of strain concen-tration by notches and cracks. Journal of Applied Mechanics, 35(2):379–386.Rignot, E., Mouginot, J., and Scheuchl, B. (2011). Ice Flow of the Antarctic Ice Sheet. Science,333(6048):1427–1430, doi:10.1126/science.1208336.Robel, A., DeGiuli, E., Schoof, C., and Tziperman, E. (2013). Dynamics of ice stream temporalvariability: Modes, scales, and hysteresis. Journal of Geophysical Research, 118(2):925–936,doi:10.1002/jgrf.20072.Robel, A., Schoof, C., and Tziperman, E. (2014). Rapid grounding line migration inducedby internal ice stream variability. Journal of Geophysical Research, 119(11):2430–2447,doi:10.1002/2014JF003251.127Sayag, R. and Tziperman, E. (2008). Spontaneous generation of pure ice streams via flowinstability: Role of longitudinal shear stresses and subglacial till. Journal of GeophysicalResearch, 113:B05411, doi:10.1029/2007JB005228.Sayag, R. and Tziperman, E. (2011). Interaction and variability of ice streams under a triple-valued sliding law and non-Newtonian rheology. Journal of Geophysical Research, 116:F01009,doi:10.1029/2010JF001839.Schoof, C. (2004). On the mechanics of ice-stream shear margins. Journal of Glaciology, 50:208–218.Schoof, C. (2010). Coulomb friction and other sliding laws in a higher-order glacier flow model.Mathematical Models and Methods in Applied Sciences, 20(01):157–189.Schoof, C. (2012). Thermally driven migration of ice-stream shear margins. Journal of FluidMechanics, 712:552–578.Schoof, C. and Hewitt, I. (2013). Ice-sheet dynamics. Annual Review of Fluid Mechanics,45:217–239.Schoof, C., Hewitt, I. J., and Werder, M. A. (2012). Flotation and free surface flow in a model forsubglacial drainage. Part 1. Distributed drainage. Journal of Fluid Mechanics, 702:126–156.Sergienko, O. V., MacAyeal, D. R., and Bindschadler, R. A. (2007). Causes of sud-den, short-term changes in ice-stream surface elevation. Geophysical Research Letters,34(22):doi:10.1029/2007GL031775.Shabtaie, S. and Bentley, C. R. (1987). West Antarctic ice streams draining into the RossIce Shelf: Configuration and mass balance. Journal of Geophysical Research, 92:1311–1336,doi:10.1029/JB092iB02p01311.Shabtaie, S. and Bentley, C. R. (1988). Ice-thickness map of the West Antarctic ice streams byradar sounding. Annals of Glaciology, 11:126–136.128Shyy, W., Udaykumar, H., Rao, M. M., and Smith, R. W. (2012). Computational fluid dynamicswith moving boundaries. Dover Publications, East 2nd Street, Mineola, New York, NY 11501,USA.Stearns, L. A., Jezek, K. C., and Van der Veen, C. (2005). Decadal-scale variations in iceflow along Whillans Ice Stream and its tributaries, West Antarctica. Journal of Glaciology,51(172):147–157.Stephenson, S. N. and Bindschadler, R. A. (1988). Observed velocity fluctuations on a majorAntarctic ice stream. Nature, 334:695–697, doi:10.1038/334695a0.Stokes, C. R. and Clark, C. D. (2001). Palaeo-ice streams. Quaternary Science Reviews,20(13):1437 – 1457.Suckale, J., Platt, J. D., Perol, T., and Rice, J. R. (2014). Deformation-induced melting in themargins of the West Antarctic ice streams. Journal of Geophysical Research, 119:1004–1025,doi:10.1002/2013JF003008.Truffer, M. and Echelmeyer, K. A. (2003). Of isbrae and ice streams. Annals of Glaciology,36(1):66–72.Tulaczyk, S., Kamb, W. B., and Engelhardt, H. F. (2000a). Basal mechanics of Ice StreamB, West Antarctica I. Till mechanics. Journal of Geophysical Research, 105:463–482,doi:10.1029/1999JB900329.Tulaczyk, S., Kamb, W. B., and Engelhardt, H. F. (2000b). Basal mechanics of Ice Stream B,West Antarctica II. Plastic-undrained-bed-model. Journal of Geophysical Research, 105:463–482.Valle´e, O. and Soares, M. (2004). Airy functions and applications to physics. World Scientific.Van der Veen, C. J. and Whillans, I. M. (1996). Model experiments on the evolution andstability of ice streams. Annals of Glaciology, 23:129–137.129van der Wel, N., Christoffersen, P., and Bougamont, M. (2013). The influence of subglacialhydrology on the flow of Kamb Ice Stream, West Antarctica. Journal of Geophysical Research,118(1):97–110, doi:10.1029/2012JF002570.Walder, J. S. and Fowler, A. (1994). Channelized subglacial drainage over a deformable bed.Journal of Glaciology, 40:3–15.Whillans, I. M., Bolzan, J., and Shabtaie, S. (1987). Velocity of ice streams B and C, Antarctica.Journal of Geophysical Research, 92:8895–8902, doi:10.1029/JB092iB09p08895.Whillans, I. M., Jackson, M., and Tseng, Y.-H. (1993). Velocity pattern in a transect across IceStream B, Antarctica. Journal of Glaciology, 39:562–572.Winberry, J. P., Anandakrishnan, S., and Alley, R. B. (2009). Seismic observations of transientsubglacial water-flow beneath MacAyeal Ice Stream, West Antarctica. Geophysical ResearchLetters, 36(11):L11502, doi:10.1029/2009GL037730.Zwinger, T., Greve, R., Gagliardini, O., Shiraiwa, T., and Lyly, M. (2007). A full Stokes-flowthermo-mechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka.Annals of Glaciology, 45:29–37.130A | Appendix to chapter 2A.1 Integrated energy balanceWe start by considering a simplified energy balance for an ice stream (corresponding to (3.25)with time-dependent terms and advection terms neglected)− k∂2T∂z2= 2xyτxy, (A.1)with xyτxy the heat dissipation in the ice, which we assume to be through lateral shear heatingonlyxyτxy =12η∣∣∣∣∂u∂y∣∣∣∣2, (A.2)and with η given through (2.2). The heat dissipation is constant throughout the ice becausethe velocity is uniform through the ice thickness. Solution of (A.1) is straightforward:− kT =12η∣∣∣∣∂u∂y∣∣∣∣2z2 + c1z + c2, (A.3)with c1 and c2 constants of integration. At the surface, the temperature is given through thesurface temperature Ts, at the bed we assume the ice to be at the melting point Tm, whichdetermines c1 and c2:− kT =12η∣∣∣∣∂u∂y∣∣∣∣2(z2 − zH)−k(Ts − Tm)Hz, (A.4)so that the heat flux is− k∂T∂z=12η∣∣∣∣∂u∂y∣∣∣∣2(2z −H)−k(Ts − Tm)H. (A.5)131At the bed, we assume that heat is dissipated through basal friction and that melting occurs,so that− k∂T∂z∣∣∣∣+= − k∂T∂z∣∣∣∣−+ τbu− ρwLhm = qgeo + τbu− ρwLhm. (A.6)with the superscript indicating the limit taken from above and below, and −k∂T/∂z|− = qgeothe geothermal heat flux. −k∂T/∂z|+ is the heat flux into the ice, given through (A.5), so thatwe getρwLhm = qgeo − k∆TH+ uτb +12ηH∣∣∣∣∂u∂y∣∣∣∣2. (A.7)with ∆T = −(Ts−Tb). Introducing a constant 2σ with 0 ≤ σ ≤ 1/2 in front of the lateral shearterm then gives (2.11).A.2 The upper bound on ice stream velocityIf the effective pressure reaches N = 0 everywhere, the ice stream and the bed become discon-nected by a water film and we consider the bed to be flooded. By integrating (2.1) twice, thevelocity can be calculated asu = 2(fB¯H)n (W/2)n+1 − |y|n+1n+ 1. (A.8)This limit is marked with a red line in figures 2.5 and 2.6.A.3 Ice stream flow on a fully permeable bedFor a fully permeable bed (i.e. K →∞), (2.14a)–(2.14b) together with (2.13a)–(2.13b) requirethe effective pressure to be constant at N = N0 everywhere. For D0 = u0 = σ = 0 and τb = µN0(i.e., a Coulomb-plastic friction law), the velocity can be straightforwardly calculated asu = 2(f − µN0B¯H)n (W/2)n+1 − |y|n+1n+ 1. (A.9)132The steady state condition∫W/20 m dy = 0 together with (A.9) then leads to an algebraicequation for N0:0 = 2µN0(f − µN0B¯H)n(W2)n+1+ (n+ 2)(qgeo −k∆TH). (A.10)For n = 3 this can solved by use of the general formula for roots for quartic equations.A.4 Solutions with small lateral shear stressesTo identify the limit in which lateral shear stresses are small compared to the basal shear stress,we non-dimensionalise the steady state model with (2.13a)–(2.14b) with a power law frictionlaw. We introduce [y] = W/2, [u] = q0/f , and [N ] = c−1f 1+rq−r0 with q0 = k∆T/H − qgeo thebackground heat flux. We set y = [y]y∗, N = [N ]N∗, and u = [u]u∗, drop the asterisks, andobtain∂∂y(∣∣∣∣∂u∂y∣∣∣∣1/n−1 ∂u∂y)+ 1−Nur = 0, (A.11a)δ∂∂y(N−a∂N∂y)+ 1−Nur+1 = 0. (A.11b)The boundary conditions are∂N∂y= u = 0 at y = ±1, and∂u∂y=∂N∂y= 0 at y = 0, (A.12)and we have two non-dimensional parameters: =B¯H21/n[u]1/n[y]1/n+11f, δ = KρwLhq0[N ]1−a[y]2. (A.13)Note that both  and δ depend on the ice stream width. The parameter  measures the relativestrength of lateral shear stress to basal shear stress. For  → 0, but δ = O(1), we are in alimit where lateral shear stresses are small but drainage remains significant. These equationscorrespond to the lateral version of the equations used by Fowler & Johnson (1996).1332050200100101102103102103  ←cW in km←abK in m3/s Pa u C in m/yearKminstableunstableFigure A.1: Location of asymptotic solutions. Black solid and dashed lines are numericalsolutions for a = 2, H = 900 m, and the values of K and W displayed. The upper, fast branchcorresponds to the upper, fast branch in figure 2.5, and similarly the lower, slow branch herecorresponds to the lower, slow branch in figure 2.5. The magenta lines show the position of theanalytic and semianalytic solutions relative to these. ‘a’ marks location of solutions (A.15), ‘b’marks solutions (A.18) and ‘c’ marks solutions (A.29b). See figure A.2 for sample plots.For  = 0 the second order differential equation for the velocity (A.11a) turns into analgebraic expression for the velocity as a function of effective pressure1−Nur = 0. (A.14)This means that we can fulfill only two of the four boundary conditions above. We chooseto keep the boundary conditions on the effective pressure, as we expect that the melt rateintegrated over the ice stream width should still be zero in a steady state. (A.11b) can then beintegrated twice to give an integral equation of ice stream width as function of effective pressureat the margin Nm = N(Y ± 1) and at the centre Nc = N(Y = 0)1 =√δ(a− 1)2∫ NcNmN−a dN√N1−a − r(1−a)(3−a)r+1N3−a+1/r − C. (A.15)The integration constant C(Nm, Nc, a, r) is determined through the boundary conditions (A.12)134−50 0 50681012y in kmN in kPa  anumericalsemi−analytical−200 0 200681012y in kmN in kPa  bnumericalanalytical−10 0 1000.51y in kmN in kPa  cnumericalanalyticalFigure A.2: Sample plots of semianalytical solutions in the limit of low and high lateral shear.Panel a shows solutions (A.15), panel b shows solutions (A.18), and panel c shows solutions(A.29b).above which require that ∂N/∂y = 0 at the margin and in the centre, i.e. from solving0 = N1−a −r(1− a)(3− a)r + 1N3−a+1/r − C. (A.16)For C ∈ [0, Cmax(r, a)] this equation has two solutions N1,2 = Nm,c, which define the effectivepressures in the ice stream centre and in the margin. The fact that solutions are possible only fora range of effective pressures translates to a range of permeabilities K per ice stream width Wfor which solutions are possible, see the thin strip marked as ‘a’ in figure A.1. Sample profiles ofthe effective pressure in this limit are shown in figure A.2a (dashed line: semi-analytic solutionfrom above, solid lines: numerical solutions).A second limit can be identified for → 0 and δ → 0. Then (A.11a)–(A.11b) simplify toN = u−r, N = u−r−1. (A.17)The solution of these equations isu = 1, N = 1, (A.18)which fulfills the boundary conditions (A.12)2 in the center. In the proximity of the margintowards y ± 1 (or equivalently y = W/2), a passive boundary layer forms which ensures thatthe boundary conditions on this side are also fulfilled. Figure A.2b shows the effective pressureprofiles of the analytical solution (dashed line) and the numerical solutions, which convergetowards the analytic solution for W →∞.135A.5 Solutions with large lateral shear stresses and smalldrainageTo obtain solutions in the limit where lateral shear is dominant, we perform an alternative scalingof u and N . We define [y] = W/2, [u] = 2[y]n+1fn/(BH)n, and [N ] = (k∆T/H − qgeo)µ−1[u]−1,and set y = [y]y˜, u = [u]u˜, and N = [N ]N˜ . This gives the scaled equations∂∂y˜(∣∣∣∣∂u˜∂y˜∣∣∣∣1/n−1 ∂u˜∂y˜)+ 1− ˜N˜ = 0, (A.19a)δ˜∂2N˜∂y˜2+ 1− N˜ u˜ = 0, (A.19b)with boundary conditionsu˜ =∂N˜∂y˜= 0 at y˜ = ±1 and∂u˜∂y˜=∂N˜∂y˜= 0 at y˜ = 0. (A.20)The limit in which lateral shear stresses are large corresponds to ˜→ 0, and the solution for thecase ˜ = 0 is already provided in dimensional form as equation (A.8) above and marked with ared line in figures 2.5 and 2.6.We now consider solutions for both δ and  small, but larger than zero. At leading order,the solution to (A.19) isu˜ =1− |y˜|n+1n+ 1(A.21a)N˜ =1u˜=n+ 11− |y˜|n+1. (A.21b)While (A.21) satisfies the boundary conditions at y˜ = 0, it does not at y˜ = ±1. Here rescaling isnecessary in order to bring back the diffusive term in (A.19b). With Y = (1−|y˜|)δ˜−x), u˜ = yU ,136N˜ = zN † the inner problem isδ˜−x1+nn +yn∂∂Y(∣∣∣∣∂U∂Y∣∣∣∣1n−1 ∂U∂Y)+ 1− ˜δ˜−zN † = 0 (A.22)δ˜1+z−2x∂2N †∂Y 2+ 1− δ˜z+yN †U = 0, (A.23)with boundary conditionsU = 0 at Y = 0,∂N †∂Y= 0 at Y = 0. (A.24)To bring back the diffusive term for the effective pressure we get the conditions1 + z − 2x = 0 and y = −z,which simplifies (A.22) to∂∂Y(∣∣∣∣∂U∂Y∣∣∣∣1n−1 ∂U∂Y)+ δ˜1n (3x−1)+x − ˜δ˜x(3n−1)+(1−1n )N † = 0, (A.25a)∂2N †∂Y 2+ 1−N †U = 0. (A.25b)We require the exponents of the δ˜–terms in front of the driving stress and the basal shear stressto be positive in order to be in the limit of large lateral shear. It is easy to see that this requiresx > 1/(3 + n). Then the leading order inner solution of the velocity isU ∼ Y,which turns the equation for the effective pressure into an inhomogeneous Airy equation,∂2N †∂Y 2= N †Y − 1. (A.26)Equation (A.26) can be solved in terms of homogeneous and inhomogeneous Airy functions137(Valle´e and Soares, 2004). Hence the inner solution at leading order isN † ∼ CAi (Y ) +[C√3−pi3]Bi (Y ) + piGi (Y ) . (A.27)Matching the inner solutions with the outer velocity solutions gives y ∼ δy−xy, i.e., x = 1/3.It is easy to verify that this choice of x does indeed render the exponents of the δ˜–terms in(A.25a) positive for n = 1 and n = 3.Asymptotically, for |y˜| → 1, the outer solution goes as N˜ ∼ 1/(1− |y˜|). Equation (A.27) isdivergent for Y →∞ due to the contribution of Bi. Therefore, in order to match the inner andouter solution, the integration constant C has to be chosen in such a way that the two terms ofthe Bi(x)-function cancel in (A.27), and we obtain for the inner solutionN † ∼[1√3Ai(1− |y˜|δ˜1/3)+ Gi(1− |y˜|δ˜1/3)]. (A.28)Now, for Y →∞, this goes as N † ∼ Y −1 ∼ 1/3/(1−|y∗|), so that the inner and outer solutionsof the effective pressure match as well.The composite solution is obtained by adding the inner and outer solutions and subtractingthe overlapping term:N˜ ={n+ 11− |y˜|n+1+piδ˜1/3[1√3Ai(1− |y˜|δ˜1/3)+ Gi(1− |y˜|δ˜1/3)]−11− |y˜|}, (A.29a)u˜ =1− |y˜|n+1n+ 1. (A.29b)An example of this solution is shown in figure A.2c.138A.6 The role of subglacial drainageA.6.1 Ice streams without subglacial drainage and σ > 0In the absence of lateral drainage the steady state model isB¯H21/n∂∂y(∣∣∣∣∂u∂y∣∣∣∣1/n−1 ∂u∂y)+ f − τb = 0, (A.30a)qgeo −k∆TH+ τbu+ σB¯H21/n∣∣∣∣∂u∂y∣∣∣∣1/n+1= 0, (A.30b)with boundary conditionsu = 0 at y = ±W/2 and∂u∂y= 0 at y = 0. (A.31)From (A.30b) we can identify τb asτb =q0 − σ B¯H21/n∣∣∣∂u∂y∣∣∣1/n+1u(A.32)where we introduced the abbreviation q0 = k∆T/H − qgeo. Use of (A.32) in (A.30a) givesB¯H21/n∂∂y(∣∣∣∣∂u∂y∣∣∣∣1/n−1 ∂u∂y)+ f −q0u− σB¯H21/n∣∣∣∣∂u∂y∣∣∣∣1/n+1 1u= 0. (A.33)By setting τ = (B¯H2−1/n)1/(n+1)|∂u/∂y|1/n−1∂u/∂y we getd|τ |n+1du+ (n+ 1)σ|τ |n+1u= (n+ 1)q0u− (n+ 1)f, (A.34)which we solve by multiplication with the integrating factor uσ(n+1). We obtain|τ |n+1 =q0σ−(n+ 1)σ(n+ 1) + 1fu+ Cu−σ(n+1). (A.35)139Bearing in mind that u→ 0 in the margins, we have the constraint C ≥ 0. Substituting (A.35)back into the basal shear stress in (A.32) givesτb =σ(n+ 1)f(n+ 1)σ + 1− Cσu−(n+1)σ−1 (A.36)which adds the condition C ≤ 0 to ensure that the basal shear stress remains non-negative.Therefore C = 0 andτb =σ(n+ 1)f(n+ 1)σ + 1. (A.37)We can use this in (A.30a)B¯H21/n∂∂y(∣∣∣∣∂u∂y∣∣∣∣1/n−1 ∂u∂y)+f(n+ 1)σ + 1= 0 (A.38)to obtain the velocity solutionu =1n+ 1(21/nfB¯H [(n+ 1)σ + 1])n[(W2)n+1− |y|n+1]. (A.39)We now have two equations for ∂u/∂y: one from (A.35) and one from (A.39). Those twoobviously need to agree, which gives the conditionW = 2[(B¯H)n2σ(k∆TH− qgeo)] 1n+1 σ(n+ 1) + 1f. (A.40)Figure A.40 shows these solutions as W against H for different choices of σ. In the absence ofdrainage but with dissipation through lateral shear, physically acceptable steady state solutionsexist only for a single combination of H and W and there is no reason why this combinationshould be achieved in practice. If ice streams are too narrow, no viable steady state solutionsexist (this is the case of C < 0, so |τ | cannot remain positive as u→ 0). For an ice stream thatis too wide, energy balance in the margins requires that the excessive heat dissipation throughlateral shear be offset by freezing at the bed, which would require a negative τb (this is thecase C > 0). It is possible to establish this interpretation by separating variables in (A.35) and1400 500 1000 1500 2000 2500020406080100H in mW in km  σ=1/8σ=1/4σ=1/2Figure A.3: Solutions (A.40) for K = 0 and σ > 0. We used the values listed in table 2.1 andf = ρgH2/L.integrating.A.6.2 Bifurcation diagrams for subglacial drainage parameter KThe bifurcation diagram in figure 2.7a shows the simplest possible behaviour of the steady statesolutions under variations of the drainage parameter K. At larger ice stream widths and icethicknesses, the upper, fast branch and the lower, stable branch of steady states do not connectin the range of physical solutions. Instead, physically plausible solutions are only possible upto a minimum Kmin below which the heat dissipation term would need to become negative.Figures A.4 and A.5 illustrate this effect for the high and low velocity branches, respectively.On the upper branch, for decreasing values of the permeability, the effective pressure gradientsbecome steeper until eventually the effective pressure becomes negative in the centre, see figureA.4c. Despite the drastic changes in the subglacial water system, no significant speed-up orincrease in discharge can be observed.On the lower branch, velocities towards the margin become negative, see figure A.5b. Notethat this is only possible because we chose the the regularization constant u0 to be zero in thiscase. If the regularization constant u0 is greater than zero, then the lower branch ends whereτb = 0 in the margin. This is because the velocities enter as u2 in the heat dissipation term,which is always positive.14110−6 10−4 10−28687K in m3 Pa−1 s−1u c in m/yeara−10 −5 0 5 10050100u in m/year b−10 −5 0 5 10−505y in kmN in kPacFigure A.4: Unphysical solutions at small permeabilities along the fast, upper branch. Panela shows the ice stream center line velocity against permeability parameter K. Solutions witha negative effective pressure are marked by the dashed-dotted line. Panel b: Sample profilesof velocity at the points marked with filled circles in panel a. Panel c: Corresponding effectivepressures. We used σ = 0, n = 3, a = 0, a Coulomb-plastic basal shear stress with u0 = 0,W = 27 km, and H = 2125 m.10−9 10−71020304050K in m3 Pa−1 s−1u c in m/yeara−10 −5 0 5 100204060u in m/year b−10 −5 0 5 1005001000y in kmN in kPacFigure A.5: Unphysical solutions at small permeabilities along the slow, lower branch. Panel ashows the ice stream center line velocity against permeability parameter K. Solutions with anegative velocity are marked by the dashed-dotted line. Panel b: Sample profiles of velocity atthe points marked with filled circles in panel a. Panel c: Corresponding effective pressures. Weused σ = 0, n = 3, a = 1, a Coulomb-plastic basal shear stress with u0 = 0, W = 27 km, andH = 2200 m.1420 1000 2000100102104H in mu c in m/yearaupper limit←−20 −10 0 10 201.822.2y in kmN in kPa  c050010001500u in m/year ba=0a=1/4a=1/3a=1/2Figure A.6: Influence of drainage parameter a on ice stream steady states. Panel a: centerline velocity against ice stream thickness H for different choices of drainage parameter a andK = 0.5 Pa−1 m3 s−1. Panels b and c: corresponding velocity (panel b) and effective pressureprofiles (panel c) at the point marked with a red cross in panel a. We used a Coulomb frictionlaw and σ = 0 for all calculations.Obviously, neither of these solutions is physical, therefore we do not track the behaviourof the steady state curves in these unphysical regimes. Unphysical limits also lead to thediscontinuous curves in figure A.1, which end where the heat production term would need tobecome negative (limit marked as Kmin in figure A.1).A.6.3 The effect of the drainage parameter aTo investigate the role of the drainage parameter a in equation (2.8), we varied a while keepingthe value of K fixed, see figure A.6. Panel a of figure A.6 shows that variations in this parameterhave little to no influence on the center line velocity of the ice stream, so that the velocity profilesin figure A.6b are indistinguishable. However, larger values of a lead to larger effective pressuregradients under the ice, as shown in figure A.6c.143B | Appendix to chapter 3Here we address the situation in which the ice stream thickness scale [z]s is much larger thanthe characteristic ice ridge thickness [z]r that we identified in (3.32a)3. This is the limit γ  1.Physically, we do not expect that the ridge will be thinner than the stream, as we expect iceridge flow to generally supply mass to the stream, not vice versa. As a result, we expect thatthe ice ridge thickness is determined by the ice stream thickness when γ is large, possibly atthe expense of smaller lateral surface slopes in the ridge. We investigate the relevant rescalinghere.It turns out that a qualitatively different behaviour in the ridge from that presented inchapter 3 occurs only when γ & (1−n)/(2n+2)r , which we assume to be the case. There are thentwo distinct regimes: If γ1/2r . 1, the ice stream sets the thickness scale for the ice ridge, butice ridge flow is no longer a simple shearing flow directed predominantly towards the ice stream.Instead there is non-negligible transport in the ridge parallel to the ice stream flow. This isthe case we investigate here. In addition there is the rather exotic possibility that γ1/2r  1,which can be shown to be equivalent to W/Ws = r/s  δ−(n+1)s . In that case the ridge iswide enough that the bulk of mass accumulation over the domain is evacuated by the flow ofthe ridge parallel to the ice stream. The ice stream then plays a minor role in overall massbalance and the thickness scale everywhere in the domain is set by the ridge. Our aim here isto investigate ice streams that play a significant role in ice sheet mass balance, and we thereforedo not consider the latter case.We rescale zr = γz∗r , pr = γp∗r, τxz,r = γ2τ ∗xz,r, τyz,r = γ−2n−n+1r τ∗yz,r, τxy,r = γ3τ ∗xy,r,τxx,r = γ3τ ∗xx,r, τyy,r = γ2n+1n−1r τ∗yy,r, τzz,r = γ2n+1n−1r τ∗zz,r, ur = γ2n+1n−1r u∗r, vr = γ−1v∗r ,wr = w∗r , and drop the asterisks immediately. The force balance in the ridge is then (omitting144terms of O(rδ2r))γ2δ2r∂τxy,r∂yr+∂τxz,r∂zr−∂pr∂xr= 0, (B.1a)γ−2(n+1)1−nr∂τyz,r∂zr−∂pr∂yr= 0, (B.1b)−∂pr∂zr= 1. (B.1c)The stress-strain relationships (omitting terms of O(γ−2n1−nr δ2r) and O(γ2δ2r)) are12(∂ur∂yr+ γ−2(n+1)1−nr∂vr∂xr)= |τr|n−1 τxy,r,∂ur∂xr= |τr|n−1 τxx,r, (B.2a)12∂ur∂zr= |τr|n−1 τxz,r,∂vr∂yr= |τr|n−1 τyy,r, (B.2b)12∂vr∂zr= |τr|n−1 τyz,r,∂wr∂zr= |τr|n−1 τzz,r, (B.2c)with τr = |τ 2xz,r + (γ−2(n+1)−nr )2τ 2yz,r|1/2, and the mass balance becomes(γ1/2r )2n+2∂ur∂xr+∂vr∂yr+∂wr∂zr= 0. (B.3)The boundary conditions at the base at zr = 0 remain unchanged from (3.40):ur = vr = wr = 0. (B.4)At the surface zr = sr the boundary conditions are now (neglecting terms of O(rδ2r))γ2δ2rτxy,r∂sr∂yr+ τxz,r + pr∂sr∂xr= 0, (B.5a)γ−2(n+1)1−nr τyz,r + pr∂sr∂yr= 0, (B.5b)pr = 0, (B.5c)(γ1/2r )2n+2ur∂sr∂xr+ vr∂sr∂yr+rs∂sr∂ts= wr + a. (B.5d)When γ ∼ (1−n)/(2n+2)r , we obtain a very similar model to the diffusion model for the ridge145thickness in chapter 3. We therefore consider the case γ  (1−n)/(2n+2)r , and introduce a seriesexpansion:pr = p(0)r + γ−2(n+1)1−nr p(1)r + o(γ−2(n+1)1−nr ), (B.6a)sr = s(0)r + γ−2(n+1)1−nr s(1)r + o(γ−2(n+1)1−nr ), (B.6b)ur = u(0)r + o(1), vr = v(0)r + o(1), wr = w(0)r + o(1). (B.6c)Withτ (0)xz,r =121/n∣∣∣∣∣∂u(0)r∂zr∣∣∣∣∣1/n−1∂u(0)r∂zr(B.7)(B.1) at leading order becomes∂τ (0)xz,r∂zr−∂p(0)r∂xr= 0, (B.8a)∂p(0)r∂yr= 0, (B.8b)−∂p(0)r∂zr= 1, (B.8c)with leading order surface boundary conditionsτ (0)xz,r + p(0)r∂s(0)r∂xr= 0, p(0)r∂s(0)r∂yr= 0, p(0)r = 0, at zr = s(0)r . (B.9)Integrating (B.8c) with (B.9)3 gives p(0)r = s(0)r − zr, and (B.8a) then implies that ∂s(0)r /∂yr = 0.At leading order the ridge surface has no transverse slope. Integration of (B.8a) yields thevelocity fieldu(0)r =2n+ 1[(s(0)r )n+1 − (s(0)r − zr)n+1]∣∣∣∣∣∂s(0)r∂xr∣∣∣∣∣n−1∂s(0)r∂xr, (B.10)which is independent of the transverse coordinate yr.146At first order, (B.1b)–(B.1c) are∂τ (0)yz,r∂zr−∂p(1)r∂yr= 0, (B.11a)−∂p(1)r∂zr= 0, (B.11b)withτ (0)yz,r =121/n∣∣∣∣∣∂u(0)r∂zr∣∣∣∣∣1/n−1∂v(0)r∂zr(B.12)and boundary conditionsτ (0)yz,r + p(0)r∂s(1)r∂yr= 0, p(1)r +∂p(0)r∂zrs(1)r = 0 at zr = s(0)r , (B.13a)w(0)r + a = (γ1/2r )2n+2u(0)r∂s(0)r∂xr+ v(0)r∂s(0)r∂yr+rs∂s(0)r∂tsat zr = s(0)r . (B.13b)We can integrate (B.11a) to obtainv(0)r =2n+ 1[(s(0)r )n+1 − (s(0)r − zr)n+1]∣∣∣∣∣∂s(0)r∂xr∣∣∣∣∣n−1∂s(1)r∂yr. (B.14)Depth integration of (B.3) yieldsrs∂s(0)r∂ts+ (1/2r γ)2n+2 ∂qx∂xr+∂qy∂yr= a (B.15)whereqx =∫ s(0)r0u(0)r dz = −2n+ 2[s(0)]n+2∣∣∣∣∣∂s(0)r∂xr∣∣∣∣∣n−1∂s(0)r∂xr, (B.16)qy =∫ s(0)r0v(0)r dz = −2n+ 2[s(0)]n+2∣∣∣∣∣∂s(0)r∂xr∣∣∣∣∣n−1∂s(1)r∂yr. (B.17)147We can use the fact that sr is independent of yr to integrate over the ridge as well(1−sryM)(rs∂s(0)r∂ts+ (1/2r γ)2n+2 ∂qx∂xr)+ qin =∫ 1− sr yM0a dyr (B.18a)withqin = qy|y=s/ryM . (B.18b)Matching the surface elevation yieldsss(xs, ts) = S(0)(xs, ts) = s(0)r (xs, ts). (B.19)Matching with the boundary layer gives the far field velocitieslimY→−∞U (0) = 0, limY→−∞V (0) =qinS(0)[1−(1−ZS(0))n+1], limY→−∞W (0) = 0, (B.20)which are identical to the far field conditions on the boundary layer in chapter 3.qin can therefore be identified with the inflow in equation (3.31a). Given ss = s(0)r , we caneliminate qin from (B.18a) and (3.31a) to find a single evolution equation for ice thickness ss thatinvolves fluxes in both ice stream and ice ridge (Q and qx). Conversely, to find the rate of inflowqin into the stream (which affects the migration rate ∂yM/∂ts), we can eliminate ∂ss/∂ts from(B.18a) and (3.31a). In this case qin can be expressed purely in terms of quantities determinedby current geometric and basal conditions of the ice stream (the latter through τb).148C | Appendix to chapter 4C.1 Local behaviour for the transverse velocityIn component-free notation, the model for the across-flow components of the velocity V =V ey +Wez, (4.5), can be written as−∇P +∇ · τ = 0, (C.1a)∇·V = 0, (C.1b)where τ = 2η(∇V + (∇V)T ) and ∇ is the gradient operator in the (Y, Z)-plane. In polarcoordinates (R, θ) this is−∂P∂R+1R∂∂R(RτRR) +1R∂τθR∂θ−τθθR= 0, (C.2a)−1R∂P∂θ+1R2∂∂R(R2τθR)+1R∂τθθ∂θ= 0, (C.2b)1R∂∂R(RvR) +1R∂vθ∂θ= 0. (C.2c)Here vR and vθ are the radial and angular velocity components, respectively, so V = vReR+vθeθ.The constitutive relations for τ in polar coordinates are:τRR = η∂vR∂R, τθθ = η1R(∂vθ∂θ+ vR), τθR =12η(1R∂vR∂θ+∂vθ∂R−vθR). (C.3)The boundary conditions (4.5e) at the base becomevθ = η1R∂vR∂θ= 0 for θ = 0, vθ = vR = 0 for θ = pi. (C.4)The downstream velocity U , given byU ∼ R1n+1√2nn+ 1F2n+1 + cos θF1−n1+n for R→ 0, (C.5)1490 0.5 100.51θ/piV  b0 0.5 1−1−0.50θ/piWc  0 0.5 100.51θ/piU  aRice (1967)Elmer/Ice BVPElmer/IceBVPElmer/IceFigure C.1: Comparison of numerical velocity solutions with asymptotic solutions from Rice(1967) and the solutions of the boundary value problem (4.20) for r = 0.01. Panel a showssolutions of the downstream solution U , panel b shows solution of the across-stream solution Vand panel c shows solutions of the vertical velocity W . n = 3 in all three cases.withF (θ) =n2 − 14ncos θ +√(n2 − 14n)2cos2 θ +(n+ 1)24n(C.6)(Rice, 1967) determines the viscosity η throughη ∼ R1−n1+nN with N(θ) = F (θ)n−1n+1 . (C.7)We put η = R1−n1+nN and make the ansatz (vR, vθ) = Rβ(v¯R(θ), v¯θ(θ)) and P = Rβ−2/(n+1)Pθ(θ),which gives in (C.1b)v¯R +1β + 1v¯′θ = 0. (C.8)Here a prime denotes an ordinary derivative with respect to θ, so v¯′θ = dv¯θ/ dθ. (C.1a) becomes−a0Pθ − a1Nv¯′θ + (a2Nv¯θ − a3Nv¯′′θ )′ = 0, (C.9a)−P ′θ + b1Nv¯θ − b2Nv¯′′θ + b3 (N′v¯′θ +Nv¯′′θ ) = 0, (C.9b)wherea0 =[β −2n+ 1], a1 =ββ + 1[β +2nn+ 1], a2 =(β − 1)2, a3 =121β + 1,(C.10a)b1 =12(β +2nn+ 1)(β − 1), b2 =12(β +2nn+ 1)1β + 1, b3 =ββ + 1. (C.10b)150The angular part of pressure satisfiesP ′θ = b1Nv¯θ − b2Nv¯′′θ + b3 (N′v¯′θ +Nv¯′′θ ) . (C.11)Elimination of the pressure in (C.9) leads to a fourth order homogenous differential equationfor v¯θ with non-constant coefficients0 =(b1 + c5N ′′N)v¯θ +(c4N ′N− c2)v′θ +(c3 −N ′′N)v¯′′θ − 2N ′Nv¯′′′θ − v¯′′′′θ (C.12)wherec1 =a0a3b1, c2 =a1a3, c3 =a0a3(b2 − b3) +a2a3, c4 =[2a2a3−a1a3−a0a3b3], c5 =a2a3, (C.13)and N is given by equation (C.7). The boundary conditions equation (C.4) are likewise ho-mogenous,v¯θ = v¯′′θ = 0 for θ = 0, v¯θ = v¯′θ = 0 for θ = pi, (C.14)and we have a generalized eigenvalue problem in which the eigenvalue β is somewhat uncon-ventionally hidden in the coefficients (C.13). We solve this problem using a shooting method,which gives β = 0.271... as the lowest positive eigenvalue for n = 3. Note that β is greater than1/(1 + n), so that the viscosity is indeed dominated by gradients in the downstream velocityU . The shooting method also gives us vθ, from which vR can be calculated through equation(C.8). The velocity components (V,W ) in cartesian coordinates can be calculated from (v¯R, v¯θ)throughV = Rβ(v¯R cos θ − v¯θ sin θ), W = Rβ(v¯R sin θ + v¯θ cos θ). (C.15)The angular dependence of U , V and W is shown in figure C.1.C.2 Boundary layer model for α 1 and Pe 1The model we solved in chapter 4 was derived for the case of α ∼ O(1) and Pe ∼ O(1). Wethen considered the asymptotic behaviour of vm for large α and Pe within this model. However,151in deriving the model equations for α ∼ O(1) and Pe ∼ O(1), we have neglected terms thatdepend on the migration velocity and Pe and which could become relevant in the limit of largeα and Pe.Here we show that neglecting these additional terms is consistent with the original boundarylayer model, provided δs  Pe−β/(1+β). To do so, we reverse the approach taken in chapter 3:we now first consider the leading order boundary layer model of chapter 3 for the case that thetwo dimensionless parameters α and Pe are large, and then confirm that the asymptotic resultsderived in section 4.5 still hold.Starting with the results in chapter 3, we follow the same steps as in chapter 4 to obtaina velocity and temperature model that depends on a minimum set of parameters. Specifically,we rescale the mechanical model equations according to (4.1). The model for the downstreamvelocity U is unaffected by the choice of α and Pe and retains its form of (4.2a)–(4.2c). Themodel equations for the across-stream velocity are also unchanged from (4.5a)–(4.5c), but weobtain an extra term in the surface boundary conditionW = λPe−1vm∂S ′∂Yat Z = 1, (C.16)omitting terms of O(δs) and O(λ). Equation (C.16) replaces (4.5d)1. The term Pe−1vm is smallfor vm increasing not more than linearly with Pe. Our results below in section C.2.2 showthat vm again increases with Pe1/(1+β) with 1/(1 + β) ≈ 0.8 for δs  Pe−β/(1+β). The otherconditions (4.5d)2–(4.5g) remain valid. Also the near-field behaviour of the heat productionand transverse velocities still satisfiesa ∼ R−1, V ∼ Rβ, W ∼ Rβ, (C.17)with β = 0.271... for n = 3.We introduce the same reduced temperature as before in (4.6), but relax the assumptionof a linear temperature profile towards the ice ridge, because αr = (s/δs)(s/r)1/(n+1)αBLand Per = δrPeBL are not necessarily small now that αBL  1 and PeBL  1. Then the152temperature equation at leading order isvm∂Θ∂Y+ Peδs∂Θ∂ts+ Pe · (V,W ) ·(∂Θ∂Y,∂Θ∂Z+1 + TRTR)−(∂2Θ∂Y 2+∂2Θ∂Z2)= αa for 0 < Z < 1, (C.18a)vm∂Θ∂Y+ Peδs∂Θ∂ts− κ(∂2Θ∂Y 2+∂2Θ∂Z2)= 0 for Z < 0. (C.18b)The main difference to the model considered before, equation (4.9), is the extra term Peδs∂Θ/∂ts,which describes the time dependence of the temperature solution. Moreover the scaled bed tem-perature in the ice ridge TR now appears explicitly. The boundary conditions remain mostlyunchanged from (4.14a) and (4.14b)2–(4.14d), only the far field boundary condition towards theice ridge (4.14b)1 is replaced byΘ→ 1−TrTR+1 + TRTRas Y → −∞. (C.19)C.2.1 Temperature solutions and migration velocity for α 1For α 1, we again observe a conductive boundary layer to form near the transition from freeslip to no slip. As before, we rescale (C.18) with (Y, Z) = α−1(Y˜ , Z˜), a = αa˜ and (V,W ) =α−β(V˜ , W˜ ), and obtain the rescaled version of (C.18):Vm∂Θ∂Y˜+ α−2Peδs∂Θ∂ts+ α−(1+β)Pe · (V˜ , W˜ ) ·(∂Θ∂Y˜,∂Θ∂Z˜+ α−11 + TRTR)−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= a˜ for 0 < Z˜, (C.20a)Vm∂Θ∂Y˜+ α−2Peδs∂Θ∂ts− κ(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= 0 for Z˜ < 0, (C.20b)with boundary conditionsΘ < 0 and −∂Θ∂Z˜∣∣∣∣++ κ∂Θ∂Z˜∣∣∣∣−= 0, for Y˜ < 0, Z˜ = 0, (C.21a)Θ = 0 and −∂Θ∂Z˜∣∣∣∣++ κ∂Θ∂Z˜∣∣∣∣−≤ 0, for Y˜ > 0, Z˜ = 0. (C.21b)153At leading order, these equations are identical to (4.21) and (4.22a)–(4.22b). The leadingorder outer temperature problem (‘outer’ to the leading order near-origin boundary layer justdiscussed) also remains unchanged from (4.24), so that (4.25) still holds.C.2.2 Temperature solutions for α 1 and Pe 1To reintroduce leading order advection, we put as before Λ = Pe−11+βα, see (4.26), and considerthe case of Λ strictly O(1). This changes (C.20) to:Vm∂Θ∂Y˜+ Λ−2Pe−1+β1+β δs∂Θ∂ts+ Λ−11+β · (V,W ) ·(∂Θ∂Y˜,∂Θ∂Z˜+ α−11 + TRTR)−(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= a˜ for 0 < Z˜, (C.22a)Vm∂Θ∂Y˜+ α−2Peδs∂Θ∂ts− κ(∂2Θ∂Y˜ 2+∂2Θ∂Z˜2)= 0 for Z < 0. (C.22b)Bearing in mind that β < 1, the time-dependent term is again asymptotically small, so that(C.22) again reduces to (4.27).To consider the outer temperature problem to (C.22) as in section 4.5.3, we again rescaleV = Pe−β/(1+β)Vˆ , W = Pe−2β/(1+β)Wˆ , Z = Pe−β/(1+β)Zˆ and Y ∼ O(1). For n = 1, omittingterms of order O(Pe−2/3) this givesVm∂Θ∂Y+ Λ−1Pe13 δs∂Θ∂ts+ Λ−1(Vˆ , Wˆ ) ·(∂Θ∂Y,∂Θ∂Zˆ)− Λ−1∂2Θ∂Zˆ2= a for 0 < Zˆ (C.23a)Vm∂Θ∂Y+ Λ−1Pe13 δs∂Θ∂ts− κ∂2Θ∂Zˆ2= 0 for Zˆ < 0 (C.23b)and for n = 3Vm∂Θ∂Y+ Λ−1Peβ1+β δs∂Θ∂ts+ Λ−1(Vˆ , Wˆ ) ·(∂Θ∂Y,∂Θ∂Zˆ)= a for 0 < Zˆ (C.24a)Vm∂Θ∂Y+ Λ−1Peβ1+β δs∂Θ∂ts= 0 for Zˆ < 0 (C.24b)to an error of order O(Pe(2β−1)/(1+β)). (C.23) and (C.24) differ from (4.29) and (4.30) by the154extra terms Λ−1Pe1/3δs∂Θ/∂ts and Λ−1Peβ/(1+β)δs∂Θ/∂ts, respectively. If this term is small,then the travelling wave model is valid and our original results hold. If this term becomes ofleading order, δs ∼ Pe−β/(1+β) for n = 3, then the boundary layer does not adjust sufficientlyquickly to neglect transient affects. For n = 3, (4.31) is therefore valid for δs  Pe−β/(1+β) only,but with −β/(1 + β) ≈ −0.2, this is not a very restrictive requirement.C.3 Change of model for sub-temperate slidingIn order to allow subzero sliding for y > |ym|, we now use the same relation as (3.7) on thefrozen side with basal shear stress τc > 0. This leads to (Schoof, 2010),Σ = τc u|u| , |u| > 0, u · n = 0|Σ| ≤ τc, |u| = 0, u · n = 0for |y| > ym, z = 0, (C.25)where we used the basal shear stress vector Σ = (Σ1,Σ2,Σ3) whose components are giventhrough Σi = τijnj − τklnknlni. Since the temperatures here are below zero, we assume thatτc is distinct from τb, the basal shear stress under the ice stream, (3.7), and we introduce theadditional scale [τc]BL = [τxz]BL, put τc = [τc]BLTc and obtain with the scales in (3.46)Txz = Tc U|U| , Tyz = TcV|U| , |U| > 0, W = 0√T 2xz + (s/δs)2 T 2yz ≤ Tc, |U| = 0, W = 0for Y < 0, Z = 0 (C.26)with |U| = (U2 + (s/δs)2V 2)1/2. We require that the switch from (C.26)1 to (C.26)2 occursin the boundary layer, so that the basal boundary conditons for the ridge remain unchanged.Rescaling (C.26) with the scales introduced at the beginning of Chapter 4 then gives (4.34).In addition to the ice flow, we have to consider the heat equation: sliding in the form (C.26)dissipates heat, and we have to drop the continuity constraint on the heat flux, boundarycondition (3.13b)2. The new boundary condition at the bed is therefore:T < Tm and kice∂T∂Z∣∣∣∣+− kbed∂T∂Z∣∣∣∣−= τij ·Dij for |y| > ym, z = 0, (C.27)155or in scaled formT < 0 and∂T∂Z∣∣∣∣+− κ∂T∂Z∣∣∣∣−= αTc|U| for Y < 0, Z = 0. (C.28)C.4 Local solution for free cold-temperate boundary(provided by Christian Schoof)Here we show that the cold-temperate free boundary in the ice must contact the bed tangentiallyat the origin in the case of a sharp transition from slip to no slip. Our approach is analogousto that in appendix A of Schoof (2012). Assume that the cold-temperate boundary is locally ata polar angle θ = γ to the horizontal axis. Close to the origin, heat transport is dominated byconduction, and we have−∇2T =a0(θ)rfor γ < θ < pi, (C.29)−∇2T = 0 forpi < θ < 2pi, (C.30)where a0 is a non-negative function, with boundary conditionsT =1r∂T∂θ= 0 on θ = γ, (C.31)[T ]+− =[1r∂T∂θ]+−= 0 on θ = pi (C.32)T = 0 on θ = pi (C.33)with [·]+− indicating the difference between limits taken as φ = pi is approached from above andbelow.As in Schoof (2012), we solve the problem using a complex variable method. Before we doso, we deal with the inhomogeneous term in Poisson’s equation for γ < θ < pi. In other words,we are looking for a particular integral that solves−∇2T0 =a0(θ)r.156The r-dependence of the inhomogeneous term makes this particularly easy; we need to look forT0 = rA(θ) with A satisfyingA′′ + A = −a0.For later convenience, we also demand that A(pi) = A′(pi) = 0. By the method of variation ofparameters, the solution isA = − sin(θ)∫ θpicos(θ′)a0(θ′) dθ′ + cos(θ)∫ θpisin(θ′)a0(θ′) dθ′.We define T = T0 + T˜ for γ < θ < pi, and T˜ = T for pi < θ < 2pi, so T˜ satisfies Laplace’sequation. To conform with the usual complex variable notation, we re-label the horizontal andvertical axes x and y, and put z = x+ iy. (This means that ‘Y ’ in the boundary layer model isreplaced by ‘x’, and ‘Z’ by ‘y’). Complex conjugates are denoted by a bar, so z = x− iy, andwe treat z and z as independent variables in the usual way (England, 1971). T˜ is harmonic,and hence T can be written asT =T0(r, θ) + φ(z) + φ(z) γ < arg(z) < pi,ψ(z) + ψ(z) pi < arg(z) < 2pi,where φ and ψ are holomorphic.Next, we implement the boundary conditions, beginning with (C.31). For later convenience,we note that T = 0 along θ = γ implies ∂T/∂r = 0 along the same boundary, and we find1r∂T∂θ= i[exp(iγ)φ′(z)− exp(−iγ)φ′(z)]+ A′(γ) = 0,∂T∂r= exp(iγ)φ′(z) + exp(−iγ)φ′(z) + A(γ) = 0,or, combining the two, we have Cauchy conditions in the formexp(iγ)φ′(z) =iA′(γ)− A(γ)2,where z is a point on the boundary, with arg(z) = γ.At the bed arg(z) = pi, the jump condition on T likewise requires [∂T/∂r]+− = 0, and we157have[1r∂T∂θ]+−= −i[φ′(x)+ − φ′(x)+] + i[ψ′(x)− − ψ′(x)−] = 0,[∂T∂r]+−= −[φ′(x)+ + φ′(x)+] + [ψ′(x)− + ψ′(x)−] = 0,where individual superscripts + and − now indicate limits taken as =(z) = 0 is approached fromabove and below, respectively. This implies that ψ′ can be extended to the wedge γ < arg(z) < piby identifying ψ′ = θ′. Then the boundary conditions on arg(z) = γ become in terms of ψ′ψ′(z) =exp(−iγ)[iA′(γ)− A(γ)]2on θ = γ. (C.34)Meanwhile, the Dirichlet condition (C.33) can be rewritten as ∂T∂r = 0 on the positive half ofthe real axis, orψ′(z) + ψ′(z) = 0 on θ = 2pi, (C.35)where the holomorphic function ψ′ is defined for γ < arg(z) < 2pi, and in both boundaryconditions (C.34) and (C.35), limits taken from within that domain are implied.Note that from (C.34) and (C.35), the real part of ψ′ is given everywhere along the boundaryof the domain, which allows us to solve for ψ′. To make use of this, we map conformally to thelower half plane by puttingζ = zpi/(2pi−γ), Φ(ζ) = ψ′(z),picking the branch of ζ that has a branch cut on the positive half of the real z-axis and mapspoints on the cut approached from below (i.e. with θ = 2pi) onto the real axis. Φ is defined for=(ζ) < 0. Let ξ be real. Then we can write boundary conditions on Φ asRe(Φ(ξ)−) =0 for ξ > 0[sin(γ)A′(γ)− cos(γ)A(γ)] /2 for ξ < 0on the real part of Φ, where the superscript − now indicates the boundary being approached158from below, and=(Φ(ξ)−) =cos(γ)A′(γ) + sin(γ)A(γ)2on ξ < 0. (C.36)The boundary conditions on the real part of Φ can reformulated as a Hilbert problem if weextend Φ analytically to the entire ζ-plane cut along the negative half of its real axis, puttingΦ(ζ) = −Φ(ζ) (C.37)for =(ζ) > 0. Then the boundary conditions on Re(Φ(ξ)−) simply state that Φ is analytic acrossthe positive half of the real axis, andΦ+(ξ)− Φ−(ξ) = − [sin(γ)A′(γ)− cos(γ)A(γ)]for ξ < 0, where the superscript + now indicates a limit taken from above the real axis. Thishas solutionΦ(ζ) = −12pii[sin(γ)A′(γ)− cos(γ)A(γ)] log(ζ) + P (ζ),where the logarithm is the usual branch with a cut along the negative half of the real axis, andP is an entire function, satisfying P (ζ) = −P (ζ) to ensure that (C.37) holds. Note that thisimplies that P (ξ) = −P (ξ) and P is imaginary on the real axis.Note that we have solved for Φ using only the fact that its real part is known on the ξ-axis.The fact we have an additional condition on its imaginary part of course arises from the factthat we are really solving a free boundary problem. In this case, solving that problem meansfinding the angle γ. To find γ, we require that the solution must also satisfy (C.36),=(Φ(ξ)−) =sin(γ)A′(γ)− cos(γ)A(γ)2pilog(|ξ|) + P (ξ)/i=cos(γ)A′(γ) + sin(γ)A(γ)2.To achieve this, we have to setP ≡ icos(γ)A′(γ) + sin(γ)A(γ)2159and sin(γ)A′(γ)− cos(γ)A(γ) = 0. Butsin(γ)A′(γ)− cos(γ)A(γ) =− sin(γ) cos(γ)∫ γpicos(θ′)a0(θ′) dθ − sin2(γ)∫ γpisin(θ′)a0(θ′) dθ′+ sin(γ) cos(γ)∫ γpicos(θ′)a0(θ′) dθ − cos2(γ)∫ γpisin(θ′)a0(θ′) dθ′=−∫ γpisin(θ′)a0(θ′) dθ′=0.We have 0 ≤ γ ≤ pi so sin(θ′) > 0, and a0 > 0 almost everywhere (the heat production rateis positive almost everywhere in the ice). The only way to make the last expression vanish istherefore to put γ = pi.160


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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


Related Items