A Numerical Model of Glaciohydraulic Supercooling Thermodynamics and Sediment Entrainment by Timothy Thomas Creyts B.Sc. (Honours); The Pennsylvania, State University. 1998 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF T H E R E Q U I R E M E N T S FOR T H E D E G R E E OF Doctor of Philosophy in The Faculty of Graduate Studies (Geophysics) The University Of British Columbia, July 26, 2007 © Timothy Thomas Creyts 2007 Abstract Beneath many glaciers and ice sheets, hydrology influences or controls a variety of basal pro-cesses. Glaciohydraulic supercooling is a process whereby water freezes englacially or sub-glacially because its internal temperature is below the bulk freezing temperature. Water super-cools when it is at its freezing point and flows from an area of higher pressure (lower ambient temperature) to an area of lower pressure (higher ambient temperature) without equilibrating its internal energy. The process is dependent on the configuration of the water flow path relative to the pressure gradient driving flow. I formulate the governing equations of mass, linear momentum, and internal energy for time-dependent clear water flow based on previous work (Clarke, 2003; Spring & Flutter, 1981, 19.82). Because field evidence and steady-state theory point to water distributing laterally across the bed, I modify this theory to account for an aperture that is much wider than deep, which I refer to as a sheet. Ice accretion terms are formulated with porosity because accreting ice has residual porosity. Ice intrusion into such a water sheet is not described in the literature, and I formulate intrusion based on previous work as well as ideas gained from subglacial cavity formation. In addition, I modify the clear water equations to include erosion and deposition of sediment along the glacier bed and incorporation of sediment into the accreted ice. Furthermore, water may leave the ice-bed interface and flow through the glacier pore space because subglacial water pressure is relatively high when supercooling occurs. To this end, I develop an englacial water-flow model that incorporates changes in ice porosity based on creep closure and ice melt or accretion. Simulations reveal behavior that cannot be inferred from simplified models. For example, while total ice accretion is comparable to field estimates, locations of simulated ice accretion along the ice-bed interface conflict with steady state models, which tend to overpredict accretion amounts. Simulations also indicate that much sediment deposition occurs prior to water being supercooled. Sediment deposition tends to smooth subglacial topography rather than enhance n Abstract it. Additional results and implications of numerical simulations are discussed. Table of Contents Abstract " Table of Contents i v List of Tables • • i x List of Figures x List of Symbols x i v Acknowledgements x y i i u Dedication X X V 1 1 Introduction 1 1.1 Glacier hydrology 1 1.2 Previous work ? 1.2.1 Storglaciaren 8 1.2.2 Matanuska Glacier • 10 1.2.3 Other glaciers 1 2 1.3 Importance and scope 13 1.3.1 Importance 13 1.3.2 Scope 14 1.4 Structure of the thesis 16 2 Subglacial flow equations 18 2.1 Conceptual model 20 2.1.1 Hydraulic geometry 20 iv Table of Contents 2.1.2 Clausius-Clapeyron equation 24 2.2 Balance laws 26 2.2.1 Mass balance 27 2.2.2 Momentum balance 32 2.2.3 Energy balance '. 35 2.3 The melt rate 36 2.3.1 Semiempirical correlations 37 2.4 Steady state accretion rate 40 2.5 Summary 41 3 Sheet closure mechanisms 42 3.1 Formulation 43 3.1.1 Regelation closure 44 3.1.2 Creep closure 47 3.2 Bed partitioning 49 3.2.1 Stress partitioning 50 3.3 Stress-velocity relation 51 3.4 Effective scales • • 53 3.4.1 Regelation area and volume scales 53 3.4.2 Creep length scale 55 3.4.3 Sediment grain size distribution 55 3.5 Example closure rates 56 3.5.1 Closure rates for beds with uniform grain size 56 3.5.2 Multiple grain sizes 61 3.6 Summary 65 4 Subglacial sediment transport 68 4.1 Balance laws for sediment transport • 70 4.1.1 Sediment mass balances 70 4.1.2 Fluid flow balance laws 73 4.2 Subglacial sediment characteristics 76 4.2.1 Mobile sediments 79 v Table of Contents 4.3 Equilibrium transport rules 81 4.3.1 , Sand formulae 83 4.3.2 Bed shear and sediment transport . 85 4.4 Conceptual model of sediment and ice accretion 87 4.5 Verification of the results 89 4.6 Summary 90 5 Englacial aquifer 92 5.1 Introduction 92 5.2 Balance laws . 93 5.2.1 Momentum balance 93 5.2.2 Mass balance 95 5.2.3 Energy balance 97 5.3 Pore geometry 99 5.3.1 Creep closure rate 100 5.3.2 Volumetric melt rate 101 5.4 Aquifer thickness evolution 102 5.4.1 Water table evolution 102 5.4.2 Surface boundary 103 5.4.3 Subglacial boundary 105 5.4.4 Downstream and upstream boundaries 106 5.5 Summary 106 6 Results and discussion 108 6.0.1 Longitudinal sections 108 6.0.2 Model details 109 6.1 Steady flow: simple hydraulic sheets 119 6.1.1 Constant recharge simulations 119 6.1.2 Energy Balance 138 6.1.3 Effect of porosity of the accreted ice 144 6.2 Comparison of channels and sheets 146 6.2.1 Results 147 Table of Content s 6.2.2 - Discussion . . • 149 6.3 Diurnal forcing 152 6.3.1 Results 153 6.3.2 Discussion 166 6.4 Sediment transport 168 6.4.1 Model details . 168 6.4.2 Initial and boundary conditions 170 6.4.3 Results 171 6.4.4 Discussion 186 6.5 Summary 188 7 Conclusions 189 7.1 Summary , 189 7.1.1 Model 189 7.-1.2 Simulations 190 7.2 Directions for future work . ._ 193 7.2.1 Technical directions 193 7.2.2 Topical directions 194 7.3 Outlook 195 Bibliography 196 v Appendices A Synthetic glacier section . 213 B Gra in parameterization for ice intrusion 216 C Solution of the equations governing ice intrusion 218 C . l Three equation model 218 C. 2 Two equation model 219 D Sediment mass balances 221 D. l Basal water mass balance 221 v i i Table of Contents D.I.I Channel 221 D.1.2 Sheet ' . . . 224 D. 2 Basal ice sediment concentration •. . . 225 D.2.1 Sheet • 226 D.2.2 Channel 227 E Semiempirical sediment transport relations 230 E. l Bed load . : 230 E. 2 Suspended load 234 F, Numerical formulation 238 F. l Sheet and channel numerical formulation 238 F . l . l Grid and discretization . 238 F.l.2 Boundary conditions 243 F.1.3 Initialization . ' .246 F.2 Sediment transport numerical formulation 246 F.2.1 Grid and discretization 246 F.2.2 Boundary conditions 249 F.2.3 Initialization 250 F.3 En'glacial aquifer numerical formulation . . . • 250 F.3.1 Vertical coordinate transformation 250 F.3.2 Discretization V •. . . . 251 ' F.3.3 Boundary conditions .• 256 F.3.4 Initialization 256 F. 4 Sheet closure ,. 256 F.4.1 Grain volume 260 G Additional discretization characteristics s 262 G. l Vertical coordinate transformation 262 G.2 Finite volume averages on non-rectilinear meshes 262 G.3 Interpolation of bed and surface topography 264 List of Tables 3.1 Closure parameters 58 6.1 Parameters synthetic glacier sections 112 6.2 Model parameters 114 6.3 Additional model parameters 115 6.4 Model parameters specific to sediment transport 169 A . l Parameters for the synthetic glacier section 214 F . l Subglacial water flow quantities on the centered and staggered grids 239 F.2 Sediment transport quantities on the centered and staggered grids , . 246 F.3 Englacial water flow quantities on the centered and staggered grids 252 F.4 Maximum interpolation errors for each of the largest grain sizes 258 ix List of Figures 1.1 Illustration of simple subglacial hydraulic systems 3 1.2 Downward view of a moulin 6 1.3 Synthetic geometry • 7 1.4 Schematic map of Storglaciaren 9 1.5 Longitudinal section of Storglaciaren 9 1.6 Basal ice at Matanuska Glacier 11 2.1 Subglacial flow path 22 2.2 Subglacial water cross sections 23 2.3 Phase diagram of water 25 3.1 Regelation closure 45 3.2 Creep closure • • 47 3.3 Example closure velocities for a single grain size 57 3.4 Total closure velocities for a single grain size with variable effective pressure. . . . 59 3.5 Total closure velocities for a single grain size with variable fractional area water. 60 3.6 Example closure velocities for a single grain size 61 3.7 Total closure velocities for a single grain size with variable effective pressure. . . . 62 3.8 Total closure velocities for a single grain size with variable effective pressure. . . . 63 3.9 Closure velocities for multiple grain sizes 64 3.10 Logarithmic grain size distributions 65 3.11 Stress partitioning for a case with ten grain sizes. . '. 66 3.12 Stress partitioning for ten grain sizes with grains exposed to their equator 67 4.1 Illustration of sediment transport and accretion 88 4.2 Hjulstrom diagram 90 x List of Figures 5.1 Englacial aquifer 102 6.1 Longitudinal Matanuska Glacier sections 110 6.2 Simulated synthetic sections I l l 6.3 Threshold output . 120 6.4 Threshold output: discharge and Reynolds numbers 122 6.5 Threshold output: temperature depression and ice accretion 123 6.6 Threshold output: melt rates 124 6.7 Above-threshold output 125 6.8 Above-threshold output: discharge 126 6.9 Above-threshold output: temperature depression and ice accretion 126 6.10 Above-threshold output: melt rates 127 6.11 Above-threshold output: heat balance terms 129 6.12 Above-threshold output: approximate heat balance terms 130 6.13 Below-threshold output 131 6.14 Below-threshold output: discharge and Reynolds numbers 132 6.15 Below-threshold output: A T b and melt rate. . 132 6.16 Flat-bedded output 133 6.17 Flat-bedded output: discharge and Reynolds numbers 134 6.18 Flat-bedded output: A T b and melt rate 135 6.19 Flat-bedded output: heat terms 136 6.20 Above-threshold output for fd = 0.20: discharge terms 140 6.21 Above-threshold output for fa = 0.12: discharge terms 141 6.22 Accretion characteristics for different friction factors 142 6.23 Energy balance comparison for different friction factors .' . 143 6.24 Accretion characteristics for different ice porosities 145 6.25 Effective pressures for subglacial sheets 147 6.26 Effective pressures for subglacial channels '. . . 148 6.27 Effective pressures for subglacial channels 149 6.28 Synthetic diurnal forcing 152 6.29 Daily cycles: Threshold case ub, ae 154 6.30 Daily cycles: Threshold case Re, m, closure rate 155 xi List of Figures 6.31 Daily cycles: Threshold case Hb, Zh:e, closure 156 6.32 Daily cycles: Above-threshold case ub, ae 158 6.33 Daily cycles: Above-threshold case A T , m, closure rate . . 159 6.34 Daily cycles: Above-threshold case Hb, Zb:e, closure . 160 6.35 Daily cycles: Above-threshold case Hb 161 6.36 Daily cycles: Below-threshold case. AT, m 162 6.37 Daily cycles: Below-threshold case. Lower discharge: Hb, Zb'e, closure 163 6.38 Daily cycles: Below-threshold case. Lower discharge: ub, ae 164 6.39 Daily cycles: Below-threshold case. Lower discharge: ae, closure rate 165 6.40 Daily cycles with erosion: Threshold case 172 6.41 Daily cycles with erosion: Threshold case, Hb, Closure, Azs 173 6.42 Daily cycles with erosion: Threshold case, Zb:e, A1 174 6.43 Daily cycles with erosion: Threshold case, Qmeit, hc, Qb, rh, \I/S 175 6.44 Daily cycles with erosion: Threshold case, A b , A b , Xb 176 6.45 Daily cycles with erosion: Above-threshold case 177 6.46 Daily cycles with erosion: Above-threshold case, Hb, Closure, Azs 178 6.47 Daily cycles with erosion: Above-threshold case, Zb:e, A1 . 179 6.48 Daily cycles with erosion: Above-threshold case, Qmelu hc, Qb,m,Vs 180 6.49 Daily cycles with erosion: Above-threshold case, A b , A b , Xb 181 6.50 Daily cycles with erosion: Below-threshold case 182 6.51 Daily cycles with erosion: Below-threshold case, Hb, Closure, Azs 183 6.52 Daily cycles with erosion: Below-threshold case, m, Closure rate, tys. . . . . . . . 184 6.53 Daily cycles with erosion: Below-threshold case, A b , A b , Xb 185 A. l Synthetic section parameters defined 214 B. l Grain geometry 215 E. l Shields diagram 231 F. l One dimensional grid for subglacial water flow 238 F.2 Ice geometry with upstream crevasse 243 F.3 Two dimensional grid for englacial water flow 251 List of Figures F.4 Grain size distribution used in the model 256 F.5 Comparison of calculated and interpolated closure schemes 257 F. 6 Average sheet thickness versus maximum sheet thickness 259 G. l Four-cornered two-dimensional grid. 262 List of Symbols Notation Superscripts are both geometrical and phasic. The geometrical superscripts b , e , and r are used for variables in the subglacial, englacial, and supraglacial water systems, respectively. In addition, an upstream crevasse is denoted by c . Phasic quantities, such as mass densities, will be referred to with the superscripts w , and s , for water, ice or glacier, and sediment, respectively. Throughout the thesis, if any term with a superscript is raised to a power, then the term is enclosed in parentheses before the index is placed on the term. Subscripts are less systematic and will be denned for quantities as they appear. Common subscripts include but are not limited to * and j as indices, a for areal, e for effective or equi-librium, T for thermal, H for hydraulic, and v for volumetric. Where closure of the ice into the subglacial water sheet is discussed, r indicates the process of regelation, and c indicates viscous creep of ice. Where sediments are discussed, the subscripts (, and s indicate bedload and suspended load, respectively; and & indicates a critical parameter relating to the threshold of sediment motion. Latin Letters A Flow law coefficient (compare with B below) s- 1 P a - n eq. (3.7), p. 48, Ac Cross-sectional area of a crevasse 9 eq. (F.10), p. 238 A Area of ice exposed to subglacial water and sed- 0 nr eq. (3.10), p. 49, iment ^ Total area of ith. grain size opposing ice intrusion m2 eq. (3.11), p. 50, As Area of sediment touching the glacier base m 2 eq. (3.10), p. 49, Aw Area of water touching the glacier base m 2 eq. (3.10), p. 49, Ae Effective area m2 eq. (3.3), p. 45, Aij Area of the ijth. grid cell m 2 eq. (G.2), p. 258, xiv List of Symbols B Flow law coefficient (= A~lln) Pa s1/" eq. (2.17), p. 29, a Grain Chezy coefficient m1'2 s- 1 eq. (E.3), p. 226, Ci Nondimensional coefficient [unitless] eq. (2.34), p. 36, c2 Nondimensional coefficient, (1 — f3CpPw) [unitless] eq. (2.41), p. 39, Sediment entrainment coefficient for basal ice [unitless] eq. (4.34), p. 88, Cx Arbitrary geometrical constant for a synthetic m l / 2 o x eq. (A.2), p. 208, overdeepening D Grain size diameter m eq. (4.18), p. 77, £>* Nondimensional particle parameter [unitless] eq. (4.21), p. 83, Do Reference grain size diameter m eq. (4.18), p. 77, Ds Characteristic grain diameter for suspended sed- m eq. (E.19), p. 231, iment Dx Representative grain diameter of the X t h per- m eq. (4.20), p. 83, centile of the grain size distribution E Enthalpy J eq. (2.6), p. 24, Ew Internal energy of water J eq. (5.16), p. 97, Ts Suspended load distribution correction parame- [unitless] eq. (4.28), p. 84, F ter Force on an obstacle N eq. (3.2), p. 45, F'1 Force of overlying ice N eq. (3.11), p. 50, Fs Force of sediment against ice N eq. (3.11), p. 50, Force of water against ice N eq. (3.11), p. 50, FT Fracture density per unit length m1 1 eq. (5.19), p. 98, G Error function for the ice intrusion solution varies eq. (C.l), p. 213, Hh Subglacial water sheet thickness m eq. (2.2), p. 21, Hb max Maximum value of Hb m eq. (F.45), p. 255 H1 Supraglacial water depth m eq. (5.29), p. 103, Supraglacial hydraulic conductivity tensor m s _ 1 eq. (5.30), p. 104, Hydraulic conductivity m s _ 1 eq. (5.1), p. 92, Thermal conductivity of water W m" 1 K - 1 eq. (2.33), p. 36, Ke:T Effective thermal conductivity W m"1 K - 1 eq. (3.2), p. 45, X V List of Symbols L Latent heat of water J k g - 1 eq. (2.29), p. 34, L% Adaptation length for sediment motion m eq. (4.4), p. 71, Mb:e Mass of accreted ice kg eq. (D.15), p. 221, M f Fluid mixture mass of water and sediment kg eq. (D.l), p. 216, M s Mass of sediment kg eq. (4.1), p. 70, M w Mass of water kg eq. (2.10), p. 27, M 0 Mass of the bed below zs kg eq. (4.5), p. 72, Nu Nusselt number [unitless] eq. (2.34), p. 36, N Number [unitless] eq. (3.11), p. 50, Ns Number of grains [unitless] eq. (4.18), p. 77, iV0s Reference number of grains of radius ro [unitless] eq. (3.23), p. 55, Number of sediments of the i grain size [unitless] eq. (3.21), p. 54, Pr Prandtl number [unitless] eq. (2.34), p. 36, P[ Melting perimeter m eq. (2.5), p. 23, Hydraulic perimeter m eq. (2.2), p. 21, Q b Subglacial water discharge m3 s _ 1 eq. (2.41), p. 39, Qmelt Flux of meltwater into a crevasse m3 s - 1 eq. (F.10), p. 238, n Bed to surface slope ratio [unitless] eq. (2.44), p. 40, Re Reynolds number [unitless] eq. (2.34), p. 36, i? Radius m eq. (5.24), p. 100, i? b Channel radius m eq. (2.2), p. 21, i*h Hydraulic mean radius m eq. (2.3), p. 22, 5 Instantaneous ice cross-sectional area, S1 — 5 W m 2 eq. (2.11), p. 27, 5 Englacial water cross section m 2 eq. (5.2), p. 93, Sh Subglacial channel cross section m 2 eq. (2.15), p. 28, S[ Channel cross-sectional area of ice m 2 eq. (2.11), p. 27, Sw Hydraulic cross section m z eq. (2.2), p. 21, T Transport stage parameter [unitless] eq. (4.20), p. 83, T Temperature K or °C eq. (2.6), p. 24, ipb Water temperature in the subglacial system T eq. (2.9), p. 26, Crevasse water temperature °C eq. (F.12), p. 239, xvi List of Symbols T - 1 mp Temperature of the melting point K or °C eq. (2.8), p. 24, rpw Temperature of water K eq. (5.17), p. 97, T0 Reference temperature K eq. (2.6), p. 24, V Volume m 3 eq. (2.6), p. 24, Ve Effective volume m3 eq. (3.2), p. 45, W Hydraulic/glacier unit width m eq. (2.2), p. 21, X General independent coordinate for the vertical m or s eq. (F.27), p. 246, coordinate transformation z Suspension parameter [unitless] eq. (E.12), p. 230, Z' Suspension parameter for sediment [unitless] eq. (4.29), p. 85, Zb:e Accreted ice thickness m eq. (2.19), p. 30, Threshold accreted ice thickness m eq. (D.17), p. 221, zc Water depth in an upstream crevasse m eq. (F.10), p. 238, Glacier thickness m eq. (2.18), p. 30, as Individual grain area m 2 eq. (3.11), p. 49, b Fracture width m eq. (5.19), p. 98, Specific heat of water J kg" 1 K " 1 eq. (2.29), p. 34, clh Representative water depth for suspended sedi- m eq. (4.26), p. 84, ment load dv Englacial pore closure rate kg m - 3 s _ 1 eq. (5.2), p. 93, fa Darcy-Weisbach friction coefficient for ice [unitless] eq. (2.25), p. 32, f'k Darcy-Weisbach friction coefficient for sediment [unitless] eq. (2.25), p. 32, fS J a Number of grains per unit size of a given texture m " 4 eq. (3.11), p. 49, per unit area of the bed fd Darcy-Weisbach friction coefficient [unitless] eq. (2.24), p. 32, fk Grain Darcy-Weisbach friction coefficient [unitless] eq. (4.31), p. 85, fd Form Darcy-Weisbach friction coefficient [unitless] eq. (4.31), p. 85, ft Fanning friction coefficient (/d/4)) [unitless] eq. (2.24), p. 32, xvii List of Symbols 9 Gravitational acceleration _ 9 m s eq. (2.21), p. 31, h Heat transfer coefficient W m" 2 K " 1 eq. (2.33), p. 36, i Subscript for zth grain size [unitless] eq. (3.11), p. 50, i Subscript for ith horizontal grid location [unitless] eq. (F.l), p. 234, 3 Subscript for jth vertical grid location [unitless] eq. (F.31), p. 247 k Permeability nr eq. (5.1), p. 92, ks Nikuradse sand bed roughness m eq. (E.17), p. 231 h Length scale m eq. (2.33), p. 36, le Effective length scale m eq. (3.8), p. 48, fh Melt rate per unit width kg m - 2 s _ 1 eq. (2.16), p. 29, m Fractal index (as an exponent) [unitless] eq. (3.23), p. 55, m Melt rate kg m _ 1 s _ 1 eq. (2.13), p. 28, mv Volumetric melt rate kg m - 3 s - 1 eq. (5.2), p. 93, n Flow law index (as an exponent) [unitless] eq. (2.17), p. 29, Porosity of ice [unitless] eq. (2.12), p. 28, < Mobile fraction of the bed [unitless] eq. (4.23), p. 83, nsp Porosity of the unlithifled bed [unitless] eq. (4.5), p. 72, na Fractional area of the bed not covered by sedi- [unitless] eq. (F.45), p. 255 ment Ox Index for a synthetic overdeepening [unitless] eq. (A.2), p. 208, P Pressure Pa eq. (2.6), p. 24, pb Subglacial water pressure Pa eq. (2.15), p. 28, pC Pressure in the upstream crevasse Pa eq. (F.12), p. 239 Pl Ice overburden pressure Pa eq. (2.17), p. 29, Water pressure Pa eq. (2.8), p. 24, q w Volumetric water flow rate per unit surface area m s_ 1 eq. (5.1), p. 92, q T Heat flux vector W m"2 eq. (2.33), p. 36, r Hydraulic potential of the glacier surface hydrol-ogy Pa eq. (5.28), p. 102, aj X 1200 3000 2500 2000 1500 1000 500 0 Distance upstream from terminus (m) Figure 1.5: Storglaciaren, Sweden, longitudinal section. E denotes the approximate equilibrium line location. R denotes the location of a riegel underlying the glacier. M O D and T O D denote the horizontal limits of the main and terminal overdeepenings, respectively. Vertical exaggeration is 1.5. Length is the same as the dashed horizontal line in Figure 1.4. Hydraulic studies at Storglaciaren led Hooke (1991) to suggest that overdeepenings play a substantial role in bedrock erosion. Part of his argument is that as water flows up the adverse 9 Chapter 1. Introduction slope, the discharge diminishes and sediment is deposited. The discharge lessens because higher water pressures force water to flow englacially. Sediment deposited when water leaves the subglacial pathways mantles the bedrock and prevents further erosion. For a sufficiently steep overdeepening, the water will move out of conduits and into a thin film along the ice-bed interface. Hooke envisioned the thin film as being mildly hydraulically conductive. He asserted that most of the drainage in the overdeepening must be englacial. Using additional borehole data, Fountain et al. (2005a) extended Hooke's (1991) idea and concluded that fractures may be the main hydraulic pathway through glaciers and, specifically, Storglaciaren. However, when Fountain et al. drilled boreholes and observed a hydraulic connection, they stopped drilling. As a result, it is difficult to ascertain what percentage, if any, of the hydraulic discharge occurs at the ice-bed interface. Fountain et al. did not provide estimates of discharge through the fracture system but proposed that the flow was significant. 1.2.2 Matanuska Glacier In a pair of papers assessing Matanuska Glacier, Alaska, a team spearheaded by Lawson and Alley described data related to subglacial supercooling (Lawson et al., 1998) and presented the basic mathematical formulation (Alley et al., 1998). At the base of Matanuska Glacier, a debris-rich package of basal ice is overlain by cleaner glacier ice (Fig. 1.6). The debris-laden ice exhibits a higher tritium content than the overlying englacial ice, which is tritium-poor (Lawson et al., 1998). The difference in tritium concentration signifies that the basal ice accreted after the early 1950s when large-scale production of tritium via free air nuclear weapons testing commenced. The tritium-poor englacial ice is older. Co-isotopic analyses of water and ice revealed that the accreted ice is isotopically lighter than the englacial ice and subglacial water discharge (Lawson et al., 1998). Isotopic analyses of frazil and anchor ice present in water issuing from the subglacial passageways have values equivalent to the accreted ice layers. The formation of the debris-rich basal layer directly from supercooled water is consistent with these data (Lawson et al., 1998). As a result, there is little doubt that the subglacial package of accreted ice and sediment is formed from actively supercooling water. Water flows through the overdeepening in conduits or a distributed drainage system. Near the terminus, water emerges at several localized vents. Radar profiles taken near these vents reveal what appear to be channels extending back toward the subglacial system (Arcone et al., 10 Chapter 1. Introduction Figure 1.6: Terminus of Matanuska Glacier. White to grey meteoric ice overlies dark-grey to black, sediment-rich basal ice that formed from glaciohydraulic supercooling subsequent to the most recent surface testing of nuclear weapons. Illumination is from the upper right. From Lawson et al. (1998). 1995; Lawson et al, 1998). Debris layers in the basal ice are laterally continuous, suggesting ac-cretion via a distributed system (Lawson et al, 1998). Frozen features resembling cross-sections of ice-choked subglacial conduits exist in the accreted ice but are rare. In addition, analyses of dye-tracing studies suggest that water distributes laterally beneath Matanuska Glacier (Lawson et al, 1998). Generally, the debris content ranges from 15-60% by volume but can range from practically none to entirely sediment (Lawson et al, 1998). Large clasts, such as pebbles and cobbles, found in the accreted ice are rounded and suggest a glaciofluvial origin, provided that those sediments have not been reworked by other subglacial processes. Fine sediment dominates the debris layers with fewer coarse-grained layers. The coarse layers can be clast-supported or ice-supported. Clast-supported layers may indicate regelation entrainment (e.g., Iverson & Semmens, 1995), but ice-supported layers form as water supercools and freezes (Lawson et al, 1998). Overall, material that is entrained via regelation accounts for a small percentage of the material in the basal layers. In addition, Lawson et al. discussed the similarity of parts of the basal ice facies dominated by silt. These layers resemble features that have formed as frazil ice floes entrain sediment. The presence of frazil ice is indicative of turbulent water flow (e.g., Daly, 1 1 Chapter 1. Introduction 1984) as is the presence of high concentrations of sediment in the flowing water. 1.2.3 Other glaciers Supercooled water and ice accretion are observed at other glaciers. Roberts et al. (2002) gave an overview of supercooling in Iceland. They separated the subglacial supercooling into two broad categories: supercooling under ablation conditions and supercooling under jokulhlaup conditions. Roberts et al. used ablation conditions to describe supercooling that occurs on a daily basis at Skaftafellsjokull under normal melt conditions. Data acquired at Skaftafellsjokull includes evidence of active frazil and anchor ice accretion. Here, the sediment-rich basal ice is several meters thick. In addition, water that flows up crevasses near the terminus supercools and accretes sediment. These conditions would be analogous to daily summertime conditions at Matanuska Glacier. Jokulhlaup conditions are strikingly different. Water under high pressure fractured the sur-face near the terminus of SkeiQararjokull during the 1996 jokulhlaup (Roberts et al, 2002). These high angle shear fractures acted as an overdeepened slope that Roberts et al. estimated was a factor of —44.5 times the surface slope, well beyond Rothlisberger & Lang's (1987) super-cooling criterion. Within the fractures, Roberts et al. identified intercalated platelets of ice and lozenge shaped ice crystals up to 40 mm in diameter that resembled deposits at Matanuska Glacier (Evenson et al, 1999). Recent observations at Taku Glacier, Alaska, included artesian water flowing from the ter-minus and freezing (M. Truffer, personal communication, 2003). There is also some accreted basal ice; however, detailed descriptions are not available. Observations are in concert with Taku Glacier's recent down valley advance and modification of its bed (Motyka & Echelmeyer, 2003; Nolan et al., 1995). Radar data point to the lower 1500 meters of the glacier lying in a sediment-floored overdeepening (Nolan et al., 1995). This area of the glacier is well below the equilibrium line, and the bed most likely receives abundant meltwater. Few data exist at other glaciers. Radar echo soundings of Muir Glacier, Alaska, reveal an overdeepening as well as accreted debris-rich basal ice (Alley et al., 2003a). However, there are no direct observations from the terminus. Observations from the terminus of Kviarjokull, Iceland indicate that subglacial supercooling is occurring and accreting debris (Spedding Sz Evans, 2002), but direct observations of supercooled water were not made. In addition, Spedding & Evans 12 Chapter 1. Introduction suspected that their observations were uncertain because co-isotopic studies of the accreted ice had not been completed at the time of publication. 1.3 Importance and scope 1.3.1 Importance Recently, Alley et al. (2003b) highlighted glaciohydraulic supercooling as being an important process through which glaciers modulate subglacial sediment erosion and deposition. Alley et al. argued that glaciers that rapidly erode and transport sediment will tend to develop an adverse slope because these glaciers cannot evacuate sediment as rapidly as they accumulate it. As a result, the bed steepens to its supercooling threshold which prevents the glacier from eroding its bed further. Hence, glaciohydraulic supercooling stabilizes the subglacial sediment erosion and deposition processes by preventing additional erosion. These arguments recapitulate and extend Hooke's (1991) ideas on the subject. However, if glaciohydraulic supercooling is a stabilizing feedback, why is it not more com-monly observed at glaciers with sedimentary bases? Questions such as this illustrate the lack of knowledge surrounding the phenomenon. Foremost, glaciohydraulic supercooling needs to be placed in the context of the larger subglacial and englacial hydrology. A general understanding of exchanges between different hydraulic network morphologies and how these influence erosion of the bed would be useful. In addition, Hooke et aUs (1988) contention that englacial pathways matter to the subglacial system appears to be corroborated by Fountain et a/.'s (2005a) work. Understanding why or how water leaves the subglacial system would aid in determining when supercooling is important to the subglacial hydrology. Fountain et al. made their observations upglacier from an adverse slope. If englacial water transport is fundamental to water flow through an overdeepening, then the water may not have the opportunity to be in contact with the bed. Obviously if water does not access the bed, then sediment at the ice-bed interface will not be mobilized. If this is the case, then the Alley et al. (2003b) hypothesis may need to be modified. Furthermore, it may be possible that the Fountain et al. (2005a) observations are enhanced by the presence of an overdeepening. Water may not be permitted to drain subglacially because the riegel acts as a dam. Water would then flow englacially. 13 Chapter 1: Introduction Ideally, issues such as these need to be resolved so that interpretations of glaciohydraulic supercooling can be applied to modern or past glaciers and ice sheets. In particular, the southern lobes of the Laurentide Ice Sheet occupied overdeepenings. For instance, the Lake Michigan Lobe may have had glaciohydraulic supercooling affecting the sediment budget, and the nearby Saginaw Lobe had active sediment transport in what appear to be, locally, subglacial sheets (Fisher et al., 2005) that could have been flowing up adverse slopes. A better understanding of glaciohydraulic supercooling may aid paleoclimatic and dynamic interpretations of ice sheets. For example large iceberg discharges, termed Heinrich Events, occurred quasi-periodically throughout the Pleistocene when northern outlets of the Laurentide Ice Sheet discharged icebergs to the North Atlantic (for a recent review, see Hemming, 2004). These icebergs had sediment-rich basal layers and rafted debris over much of the North Atlantic (Dowdeswell et al., 1995). Sediment package thickness decreases away from Hudson Strait, but lithologies point to an origin from regions that were ice-covered. The exact mechanisms of sub-glacial erosion and subsequent entrainment are not well understood. However, glaciohydraulic supercooling may be the phenomenon responsible for accreting sediment to the base of the ice sheet (Hemming, 2004; Hulbe et al, 2004). Water could have flowed up the sill of Hudson Strait, supercooled, and accreted water and sediment to the base of the ice sheet. The sill would act as the adverse slope necessary for glaciohydraulic supercooling. Similar situations may have occurred for other Pleistocene and contemporary ice sheets. 1.3.2 Scope Three areas of interest center on glaciohydraulic supercooling and overdeepenings: the shape of the hydraulic cross-section, exchanges between the englacial and subglacial system, and sediment transport within the overdeepenings. For the most part, information on these individual areas can be inferred from hydraulic, geomorphic, and geologic data. However, the rates, locations, and other analytic information throughout overdeepenings are difficult to quantify either from downhole englacial or near-terminus proglacial studies. The first area of interest is the cross-section that is favored by the hydraulic system. Data appear to indicate that a fast-flowing distributed sheet is preferred. The large body of theory on R-channels makes them, in a way, the default model of subglacial flow. Rothlisberger & Lang (1987) made no allowance for the absence of a conduit, while Hooke & Pohjola (1994) and Alley 14 Chapter 1. Introduction et al. (1998) favor distributed subglacial water systems. Furthermore, stratigraphic observations at the front of Matanuska Glacier seem to indicate that a distributed system, possibly a water sheet, is prevalent. There is also an informal "clogged artery" hypothesis. Because R-channels can act as hy-draulic arteries, the freezing point of the water may be well below that of the ambient ice. Freezing of the water then occurs rapidly, and the artery clogs itself and cannot reform faster than the water distributes into a sheet. Knowing the conditions for which the cross section is a conduit or a sheet would be useful, especially because both appear to be supported by different parts of the stratigraphic data (Lawson et al., 1998). The second area of interest is sediment transport. Both Hooke and Alley et al. (2003b) supposed that deposition occurs on the adverse slope. If this is true, the slope may change, possibly obtaining some kind of critical value; or the ice-bed profile may change, possibly assuming some kind of characteristic shape. Alternatively, there may be no feedback loop. While their arguments are persuasive, Alley et al. reason through how the system operates without presenting new field data or numerical results. I attempt to evaluate their proposed feedbacks associated with sediment transport at the bed. Glaciohydraulic supercooling can be an important component of subglacial sediment trans-fer. With its high rates of sediment accretion and glaciofluvial transport, Matanuska Glacier illustrates this significance; but questions still surround the process. For example, why is the process so well developed there relative to other locales? An understanding of the factors con-trolling subglacial sediment motion is incomplete. Clearly, surface and bed slopes are crucial (Rothlisberger & Lang, 1987), but how these interact with other hydraulic conditions to mobilize and accrete sediment is not clear. The third area of interest is the movement of water between the englacial system and the subglacial one. Principally, why does water leave a subglacial pathway to travel through an englacial one? Alternatively, the water might never enter a subglacial system, instead remain-ing englacial. Meltwater at Storglaciaren appears to favor englacial flow while meltwater at Matanuska Glacier appears to prefer to flow subglacially, but both locales show signs of glacio-hydraulic supercooling. The difference in the location of flow, englacial versus subglacial, must result from a physical mechanism that is active in one situation but not the other. Finding a pertinent trade-off would be extremely useful, especially if discernible characteristics could be 15 C h a p t e r 1. I n t r o d u c t i o n qualified for different types of glaciers. My approach is to develop and use a numerical model of en- and subglacial water flow. Much of the theory for the model is extracted from other glaciological studies or from other areas of the earth sciences. The notable exception is the closure of the sheet cross-section, which is entirely new. Attempts to link this same material either have not been tried or not been published. For example, the results from the simultaneous reconstructions of both englacial and subglacial water temperature from a hydrologic model are the first of their kind. On one hand, the numerical model and the physical model on which it is based are necessarily subject to assumptions to make the problem tractable. The diversity of material presented here rules against an exhaustive search of parameter space. I have attempted to minimize these effects or quantify them. On the other hand, the numerical model allows solutions over contiguous space and time that field measurements may not capture. The model allows detailed comparisons and contrasts of different variables and parameters both in space and time. In sum, there are notable deficiencies to numerical modeling, but there are also conspicuous advantages that I exploit. 1.4 Structure of the thesis This thesis is structured around the major components of the modeling effort. The aim is to elucidate factors controlling glaciohydraulic supercooling and to provide direction for future process-based studies. In general, shorter, simpler chapters are favored over longer, more elab-orate ones. • Chapter 1, the present chapter, focuses on introducing glacier hydrology and how glaciohydraulic supercooling fits into this larger picture. The process is identified and the scope of the work is established. • Chapter 2 introduces the fluid flow physics and the strategy for examining the sub-glacial hydrology. This physical model is modified from Clarke's (2003) effort. Rea-sonable approximations are used to determine pertinent length and time scales for glaciohydraulic supercooling. 16 Chapter 1. Introduction • Chapter 3 acquaints the reader with the way in which hydraulic sheets close and open in a vertical sense. These arguments provide a fresh look at regelation and creep closure from the work of Nye (1967) and Weertman (1972). In addition, the strong arguments put forth by Walder (1982) for the instability of sheets allow an easy access point for reassessing the stability of locally distributed sheets along the adverse slope. • The sediment transport relations are advanced in Chapter 4. Semi empirical sediment transport relations culled from the fluvial literature are recapitulated, modified, if necessary, and applied to the subglacial environment. • Chapter 5 captures the flow physics and their necessary approximations for englacial fluid flow in a Darcian model. This model describes the evolution of the pore space in the ice through melting or freezing along the perimeter of the pores. Formulation of the englacial model is included to provide a thorough elucidation of the water sources and sinks that may influence flow under and through overdeepenings. • Chapter 6 discusses the simulations and the implications of these simulations. I con-struct synthetic longitudinal sections that have characteristics of both Matanuska Glacier and Storglaciaren to simplify numerical results. Each simulation type pro-vides insight into the governing processes. Furthermore, different components of the previous chapters are turned on or off to produce stark contrasts and allow compari-son. Simulations focus on flow at the ice-bed interface. Englacial aquifer simulations are omitted for brevity and simplicity. • The conclusions to this thesis are found in Chapter 7. The summary of the major results can also be found here along with the important additions and modifications of our present understanding of glaciohydraulic supercooling. These are placed in the context of the larger body of glacier hydrology. Finally, ideas for future work and a synopsis of unfinished details are presented. Some of these chapters rely on useful auxiliary relations. These auxiliary relations appear in appendices at the end of the thesis. The appendices are presented in the order in which they are referenced in the chapters. 17 Chapter 2 Subglacial flow equations Physical models have been used as both diagnostic and prognostic tools in glacier hydrology. These models tend to be process-specific. Extensive examples of these types of models can be found in jokuhlhaup studies (e.g., Clarke, 2003, 1986; Fowler, 1999; Ng & Bjornsson, 2003; Nye, 1976). Fewer examples can be found in studies of cavities (e.g., Walder, 1986) or descriptions of morphological switching during surging (e.g., Kamb, 1987). A class of broader comprehensive numerical models is difficult to construct because of the feedbacks among ice dynamics, water flow, and basal conditions. In particular, capturing the essential physics of some processes, mor-phological switching (e.g., Kamb, 1987), for example, in numerical models is challenging and requires creative parameterizations (e.g., Flowers et al., 2004). Presently, the use of compre-hensive glaciohydrological models is limited to feedbacks between ice dynamics and subglacial hydrology (Arnold & Sharp, 2002; Johnson & Fastook, 2002) and feedbacks among different units of glacier hydrology (Flowers & Clarke, 2002a). The model I present here is not necessar-ily process specific, as it could easily be modified for other sub- and englacial uses, but I choose to present it with the specific intention of explaining glaciohydraulic supercooling. This model is not coupled to a model of ice flow and, therefore, attempts to capture feedbacks only within the glaciohydrologic system. Rothlisberger & Lang (1987) and Alley et al. (1998) described the physics of supercooling subglacial water. Aside from the hydraulic cross-section, both papers treated water flow similarly inasmuch as they make use of the same assumptions. For example, they both assumed that the water is always at the freezing point along the flow path. In addition, they assumed steady and uniform flow. Indeed, Alley et al. followed much of Rothlisberger & Lang's representation but added parameter analyses specific to glaciohydraulic supercooling. Any model of subglacial water flow that treats the thermodynamics appropriately can be used to simulate subglacial supercooling. Clarke (2003) provided one of these models. Clarke's work is based on the work of Spring & Flutter (1981, 1982) for flow through ice-walled conduits. 18 Chapter 2. Subglacial Row equations This model formulates the physics in terms of the continuum balance laws. Clarke modified the original theory by utilizing a common compressibility assumption to solve the two simultane-ous mass balance equations. In addition, he modified the conduit geometry, adding a simple sediment-floored conduit. Following Clarke's (2003) modification of the hydraulic cross-section, I introduce a dis-tributed water sheet as an appropriate hydraulic cross-section for water-flow through an overdeep-ening. There are several reasons for the introduction of this cross-section. Most importantly, accreted basal layers at Matanuska Glacier are broader than they are high (Lawson et al., 1998). Additionally, results from dye-tracing studies suggest that water distributes laterally at the ice-bed interface (Lawson et al., 1998). There is also evidence that water pressure is above flotation in parts of overdeepenings. Hooke (1991) suggested that water distributes laterally across the bed when overpressures occur. Furthermore, a turbulent sheet parallels the concept of broad-and-flat channels that were originally developed for water flow through overdeepenings (Hooke et al, 1990); however, in the case of sheets, it is not necessary to calculate the lateral extent of the channel. Alley et al. (1998) also envisioned a laterally distributed system but selected slower types of cross-sections (Fig. 1.1). I adopt a type of cross section that is consistent with turbulent conditions indicated by the presence of frazil ice and a considerable sediment load at the front of Matanuska Glacier (Lawson et al, 1998). Both fast and slow types of hydrology are not exclusive of one another as there is most likely longitudinal and lateral variation within the hydraulic system (e.g., Lappegard et al, 2006). Furthermore, water discharged through sheets as envisioned here can move fast or slow. Within this chapter, I review the formulation of the basic governing equations for the case of steady flows, which Rothlisberger & Lang (1987) and Alley et al. (1998) used in their anal-yses, and time dependent flows, which Clarke (2003) used. Equations from the Spring-Hutter formulation can be reduced to the steady-state representation, and this feature will be exploited in this chapter. In addition, I revisit the necessary fundamentals for simulation of water flow and glaciohydraulic supercooling. The Clausius-Clapeyron equation will be used throughout the thesis and requires illumination. The hydraulic geometry also requires parameterization, but these assumptions appear to be supported by observations reasonably well, suggesting that model results are pertinent to conditions in the subglacial system. 19 Chapter 2. Subglacial Bow equations 2.1 Conceptual model Water flows along the ice-bed interface (Fig. 1.3), in either a channel or sheet (Fig. 1.1). Flow in both R-channels and sheets can be formulated as a one-dimensional process using an along-path coordinate, s. The location of the water pathway is defined in x and z, the horizontal and vertical coordinate, respectively. The third coordinate, y, perpendicular to flow, is not simulated, and any variation in the flow path is assumed to be only in the x and z coordinates. The along-path coordinate is a line integral along the ice-bed interface in x and z. The change in any increment of the flow path is given by As= y/Ax2 + Az2. (2.1) The glacier overlies this water system, sealing it. The water layer has a substrate floor, which can be either rock or sediment. Data at both Matanuska Glacier (Lawson, 1979) and Storglaciaren (Hooke et al., 1997) suggest that much sediment exists at the base of these glaciers, such that I refer to the floor as being sediment. For the purposes of this study, the main difference between a rock and sediment floor is that a sediment bed may be mobile, whereas a rock bed would be static. I anticipate that both rock and sediment floors may be found under different areas of overdeepenings within the same glacier, but do not anticipate simulating this transition directly. The dynamical behavior of the water system and its interaction with the overlying ice and underlying sediment is described as each of the individual balance laws is introduced subsequently in this chapter. 2.1.1 Hydraulic geometry Generally, I consider geometries with cross sections that are circular, or Rothlisberger (R-) chan-nels (Rothlisberger, 1972), the semicircular R-channels envisioned by Clarke (2003), and water flow in sheets. Both R-channels and sheets can conduct water at a relatively high velocity. Other types of cross-sections have been described (e.g., Ng, 1998), but these representations are mathematically bold, and would require diligent modification of the Spring-Hutter formu-lation. Because circular channels have a well-described theory in the literature (e.g., Clarke, 2003; Fowler, 1999; Rothlisberger, 1972; Shoemaker, 1987; Shreve, 1972; Spring & Hutter, 1982; Weertman & Birchfield, 1983), these will not be discussed in as much detail as the subglacial sheets. Sheets are not nearly as well described in the literature (Walder, 1982; Weertman, 1972). 20 Chapter 2. Subglacial flow equations These sheets are a layer of water having a mean thickness and a width that is much greater than that thickness. Furthermore, a water sheet compatible with the work of Alley et al. (1998) or Hooke (1991) is not adequately characterized in relation to high velocity water flows. I resist the flow velocity characterizations of Hubbard & Nienow (1997) and envision sheets as being akin to broad-and-flat channels in flow characteristics (Fig. 1.1c) but more like films in cross section (Fig. l . ld). Within the glaciological community, the terms sheet and film are commonly applied for thin subglacial water layers. However, I draw the distinction that sheets are capable of conducting water rapidly. The term film will assume its conventional use to describe melt associated with regelation (e.g., Lliboutry, 1993; Nye, 1973a; Weertman, 1964). Rothlisberger & Lang (1987) performed their analysis for water flowing in conduits but acknowledged that any hydraulic cross section is possible. Alley et al. (1998) suggested a distributed water system, and as a geometrical simplification, performed analyses for a one-dimensional water layer. Alley et al. cited evidence from Matanuska Glacier that shows that the bulk of accreted ice is from water flow that is laterally distributed over the base (Lawson, 1979; Lawson et al, 1998). Because I am investigating channel and sheet cross-sections, both must be perpetuated throughout the theory. Water flows along the ice-sediment interface until exiting the glacier section downstream of an overdeepening. From input to output, an element of water is described in terms of the 1-D coordinate (Fig. 2.1). The base of the water lies at an elevation zs. The roof of the water is at an elevation zh. A glacier of thickness Z 1 overlies the water layer. Both R-channels and sheets can be defined in terms of their hydraulic perimeter P w , and cross-sectional area Sw, P S where Rb is the channel radius, Hb is the sheet thickness, and W is the sheet width. These cross-sectional geometries are summarized in Figure 2.2. In the along-path coordinate, these 21 w , 2irRh circular channel, (n + 2) Rb semicircular channel, (2.2a) 2Hb + 2W sheet, n(Rb)2 circular channel, = ^ ^(Rb)2 semicircular channel, (2.2b) HbW sheet, Chapter 2. Subglacial Bow equations Outlet Figure 2.1: Subglacial flow path. Not to scale. Front of the section presented in Figure 1.3. Water flows up the ice-sediment interface along path coordinate s. From left to right: Hh is the water layer height, which is equivalent to 2Rb. Z1 is the ice thickness. The along path coordinate s follows the water layer. The surface elevation is zx. The elevation of the glacier base is zh. The elevation of the base of the water layer is zs. The angle of the surface and the bed are a r and ah, respectively. cross-sections vary smoothly, growing larger or smaller depending on flow conditions (Nye, 1976; Spring & Flutter, 1982). In loose terms, the channels behave as tubes, while the sheets behave as ducts. Sheets have a cross sectional area; but for a sheet continuous in y, its thickness Hh determines the behavior of the cross-section. Equations governing flow become width-averaged quantities. This sheet thickness spans the average elevation of the subglacial sediments zs to the average elevation of the base of the glacier zh. These two horizons are locally parallel. The hydraulic radius, which is the ratio of the cross-sectional area to the wetted perimeter, Rh = Sw jPw, is commonly used in fluid flow engineering applications, Rh = 2 irRh circular channel, semicircular channel. (2.3) 2(TT + 2) For tubes with circular cross sections, the hydraulic radius is simply half the radius. For a semicircular conduit bounded on its base by sediment, the hydraulic radius is given as one quarter of the radius. Water sheets are relatively shallow. The width of the flow is then much 22 Chapter 2. Subglacial flow equations IT H b 1 W Figure 2.2: Subglacial water cross sections. White represents water flow area, light grey areas represent ice, and dark grey areas represent subglacial material, (a) R-channel surrounded by ice. (b) Semicircular R-channel. (c) Water sheet. Cross-sectional areas for both channels are equal. The cross-sectional area of the sheet is approximately half of the other two. Water flow is perpendicular to the plane of the page. greater than its height (W 3> Hh), and the hydraulic radius for sheets can be approximated by fib Rh = -y- (2.4) For the semicircular case, if the roughness of the bed is significantly less than the radius of the channel, then equation (2.3) is sufficient. For subglacial sheets, this approximation depends on how the water floods the roughness elements. For lack of a better parameterization, this hydraulic radius (eq. 2.4) will be used. The melting perimeter P1, will be different for each of the geometries, 2nRb circular channel, irRh semicircular channel, (2-5) 2Hh + W sheet. Only for the circular case will the melting perimeter be exactly the wetted perimeter. For the both the semicircular cross-section and the sheet, only the ceiling will melt. The bed is excluded from any melt relation. Clarke (2003) noted that the melt rate depends on the melting perimeter along the flow rather than simply the hydraulic perimeter. Assuming that the water layer is ice-roofed, the melting perimeter is the width of the layer. If the water layer is bound on either side by ice, then the small height approximation can be applied: melt on vertical ice walls can be ignored. However, if the sheet has sediment walls, then the melting perimeter is exactly the roof width. Either way, melt from the wall of the sheet is not immediately important. 23 Chapter 2: Subglacial Row equations Supercooled water need not form ice on the underside of the glacier. Examples from fluvial studies (e.g., Osterkamp & Gosink, 1983; Kerr et al, 2002) illustrate that much ice can form on rocks at the base of rivers. This fluvial process may be paralleled subglacially such that ice forms on the sediment below the water. Ice could encase some of the sediments, and then mobilize the sediments into the flow or float them to the ice-water interface where they freeze to the underside of the glacier. This idea merits attention but is not of primary importance. For simplicity, crystallization of ice can be assumed to occur along the melting perimeter. The melting and freezing perimeter are exactly the same for the remainder of this study but need not be. 2.1.2 Clausius-Clapeyron equation Figure 2.3 illustrates the phase diagram of water. Lines of coexistence separate water vapor, liquid water, and water ice. The unique feature of these lines is that the change in Gibbs free energy.required to cross from one phase to another is zero along them. By setting the change in Gibbs free energy to zero, the Clausius-Clapeyron equation can be derived for the change of phase (e.g., Krauskopf & Bird, 1995, p. 420), where the slope of the line dp/dT, is defined by the change in enthalpy AE, a change in volume AV, and a reference temperature for the phase change To. The inverse slope, dT/dp, is the quantity of interest. This slope is obtained by inverting equation (2.6) to obtain In moving from (2.6) to (2.7), the latent heat of fusion per unit mass L is substituted for the enthalpy; and the change in volume is replaced with the change in mass densities. The quantities p1 and p w are the usual mass densities of ice and water. Because the slope of the line of coexistence between water and ice is relatively smooth, the change in dT/dp can be held constant for values near To. The change in the zero point from the reference value To when the interfacial pressure pw, changes is given by (2.7) Tmp = 7b + Bp ,w (2.8) 24 Chapter 2. Subglacial How equations 10> IB Q. 0) 3 CO in „ „ 4 2 10 10' Solid (Ice Ih) Liquid (Water) Gas (Vapor) 240 250 260 270 280 290 300 Temperature (K) 310 320 330 340 Figure 2.3: Phase diagram of water. Lines indicate where both phases are stable. Al l three phases are stable at the triple point (~ 612 Pa, 273.16 K). The maximum pressure on this diagram (3.0 x 107 Pa) is equivalent to an ice thickness of 3300 m. Adapted from Wagner & Pruss (1993) and Wagner et al. (1994). where the coefficient (3 is equal to the inverse Clausius-Clapeyron slope dT/dp. The reference freezing point is 273.15 K in equation (2.7), but could also be 0 °C in (2.8) because this equation denotes only a relative change. From equation (2.7) the parameter (3 will be negative because ice is less dense than water at the freezing point. However, when these equations are applied to the subglacial environment, the density of water must include air within the water. The effect of air on the water density was noted by both Rothlisberger & Lang (1987) and Alley et al. (1998). The interfacial pressure was discussed by Spring & Hutter (1982). The two obvious choices are the overburden stress from the ice or the water pressure at the ice—water contact. Spring & Hutter used the water pressure, which makes sense given that the ice overburden may not govern the local variations in stress at the ice-water contact. I apply this interfacial pressure parameterization to equation (2.8) throughout the thesis. The Clausius-Clapeyron equation determines how the change in water pressure can modulate the melting point. Generally, (3 takes a value from —7.4 x 10~8 to —9.8 x 10 - 8 K P a - 1 depending on whether the water is pure or air-saturated, respectively. From the synthetic glacier section in Figure 1.3 with maximum ice thickness of about 47 m, assuming that the water pressure 25 Chapter 2. Subglacial Bow equations varies directly with ice overburden thickness, the maximum change of the melting point would be —0.03 or —0.04 °C for pure or air-saturated water, respectively. These example tempera-ture depressions scale linearly with ice thickness, but water pressure at the base may deviate significantly from ice overburden. One aspect of equation (2.7) that I have not discussed is the role of salts or other impurities. Salts typically depress the freezing point further. Equation (2.6) sets the change in Gibbs free energy to zero, making the change in entropy equal to the change in enthalpy divided by the reaction temperature. If salts are included, an additional term draws from the entropy term. Lawson (pers. comm., 2005) stated that all water quality measurements at Matanuska Glacier have pointed to the water being fresh. However, Bjornsson (2002) noted that a sulphurous odor preceded jokulhlaups at Skei3ararj6kull. The presence of chemical components suggests that supercooling under jokulhlaup conditions would necessitate the additional temperature depres-sion terms in the Clausius-Clapeyron equation. The effect of salts on the freezing temperature has been studied (e.g., Chen & Millero, 1977, 1986) but is not included because water quality data at most field sites are not readily available. 2.2 Balance laws The continuum balance laws provide the foundation for the model (Spring & Hutter, 1981, 1982). In total, there are four time-dependent equations: two mass balance equations, one momentum balance equation, and one energy balance equation. The four dependent variables are hydraulic cross-section, water pressure, water velocity, and water temperature. Each of these variables represents a local average in space and time. The balance equations are equivalent to the St. Venant equations that are commonly used in fluvial engineering. The main differences between the St. Venant equations and the Spring-Hutter formulation are that water flow is sealed and, to a degree, pressurized by an overlying glacier and that a phase change from water to ice is permissible (Spring & Hutter, 1981, 1982). Because the phase change is temperature dependent, an energy equation governing the temperature evolution appears that would not necessarily be present in the St. Venant equations. Furthermore, these equations have an accessible relationship to the previously published material of Rothlisberger & Lang (1987) and Alley et al. (1998). Rothlisberger & Lang pro-vided the basic relationship for melting or freezing of flowing water because of a change in the 26 Chapter 2. Subglacial Bow equations temperature-dependent melting point. In addition to assuming that the flow is steady and uni-form, they assumed that the water is always at the pressure melting point. These assumptions in a mathematical sense, using the linear subglacial water velocity uh in the s direction and the where equation (2.9a) denotes steady, uniform flow. This equation is somewhat misleading because it states that the source terms are equal to zero as well, which is the desired relationship. Equation (2.9b) indicates that the water is at the pressure melting point regardless of its location. The Spring-Hutter formulation has a reduced form that corresponds to equations (2.9a) and (2.9b). The relationship between both sets of equations will be presented in detail in subsequent parts of this chapter. Higher-order models of subglacial water flow have been proposed. For example, Bates et al. (2003) formulated flow based on the Reynolds averaged Navier-Stokes (RANS) equations. The main difference between a RANS representation and the Spring-Hutter formulation is that the dissipation terms from viscous stresses are parameterized as a shear stress along the hydraulic perimeter in the Spring-Hutter formulation. In general, the dissipation along the hydraulic perimeter is several orders of magnitude larger than the dissipation within the flow (Hinze, 1975). A RANS representation is also more complex in the case of subglacial flow because the exact hydraulic geometry is not known but would be required. Solution of a RANS representation is also computationally intensive. As a result, simulation of water flow via a RANS or other higher-order model representation is beyond the scope of the present work. 2.2.1 Mass balance Spring & Hutter derived the mass balance specifically for curvilinear tubes using the transport theorem. They assumed that the tube could vary in s. Spring & Hutter were chiefly interested in jokuhlhaups, and as a result, thought of the channel simply as melting open rather than freezing closed. For the phenomena described here, a negative melt rate is equivalent to an ice accretion rate. Because there are two phases, water and ice, there must also be two mass balance equations that include the processes of creep and melt closure to describe how cross-temperature of the subglacial water T b are (2.9b) (2.9a) 27 Chapter 2. Subglacial How equations section changes through time (Spring & Hutter, 1981, 1982). In a one dimensional system, the second dimension is approximated by having the cross section evolve. In the case of R-channels, the cross-sectional area evolves directly, while the other pertinent parameters can be calculated from that area. For example, hydraulic radius can be calculated from the area and the assumed cross-section. For the sheet, the second dimension is not critical; and as a result, it is not modeled. The mass balance relations then reduce to cross-sectional evolution and height evolution equations for their respective cross-sectional geometries. Equation (2.10) shows the general form for the mass of water in the hydraulic cross section, M w = f f pwdSds. (2.10) Following Spring & Hutter (1982), I define a difference in cross section, S = Si-Sw, (2.11) that represents the instantaneous misfit resulting from ice melted away per unit time. The mass of this misfit cross-section is M = JJ{(1- nj,) pl + n ^ w } dSds, (2.12) where nj, is the porosity of the ice. Equation (2.12) assumes that the surrounding ice is porous and that the mass is a mixture of water and ice. Moving equations (2.10) and (2.12) to a local balance per unit channel length along s via the material derivative yields dSw d nlpw {(1 ~ n') + { f + | (KS) } - - m - ^ - m . (2,3b) where m is the mass of melt per unit time per unit channel length, uj is the water velocity along s, and uls is the ice velocity along s. In these equations, the density of water and ice are assumed to be constant. Equations (2.13) show that the rate of change of mass in either cross section is equal to the amount of water produced from melting or freezing. Rewriting both of 28 Chapter 2. Subglacial flow equations these equations for dSw/dt yields 0SW ( l - nl) d + nlpw d = , p > \ . v m - j - « 5 w ) , (2.14a) dt (1 - n\) p'pw dsK 3 J' y } « - = ^ 7 + ( *)«*.• ( 2 1 4 b ) where equation (2.14a) is the mass balance for the water and equation (2.14b) is the mass balance for the ice. In moving from equation (2.13b) to (2.14b), the time derivative of the right hand side of equation (2.11) is used to substitute for the time rate of change of the area melted, dS/dt = dSl/dt — dSw/dt. Furthermore, the velocity of the ice along s is neglected. In consequence, [dS1 /dt) c ] o s e is the velocity of ice normal to water flow and represents ice closure. The form of the closure velocity [dS1 /dt) c l o s e will depend on whether the cross-section is a sheet or a channel. Solving these two equations for the hydraulic cross-section requires that a numerical solver set both results equal to the same value. This forced equality makes the problem numerically stiff. If instead, an approximation to one of the equations is made, then the equality can be relaxed. Clarke (2003, appendix A) detailed this relaxation. The salient feature of the relaxation is that the water is assumed to be compressible. One can then use the equation of state for water to substitute for the time rate of change of the water density. The net result is that the mass balance equation for water becomes an equation for the pressure of water pb, & 1_ id&, d(Sbub) _ (1-nL) p' + nLp" ] dt - ^ b S b \ d t + ds ( l - n L ) p i p w m J ' where uh replaces uJ, and 5 b replaces S l w. The notation Sb has sole reference to subglacial channels, and S'w will refer to the general forms of the cross section presented in equation (2.2b). Use of both Sb and ub makes notation simpler throughout the remainder of the thesis. Solving this equation in conjunction with (2.14b) and 7 b set as the compressibility of the water yields the exact solution produced by solving both equations (2.14a) and (2.14b). Relaxing 7 b physically approximates water being more compressible but enhances the convergence of the numerical solution. This relaxation has no noticeable effect on the overall mass balance. Clarke also noted that the mass balance of the ice could be used instead of the water balance. This approximation is equivalent to making the channel walls distensible, or stretchable. Compressibility assump-tions such as these are common in the fluid dynamics literature, and the compressibility has 29 Chapter 2. Subglacial Bow equations been used in hydrogeology to represent compressible aquifers (e.g., Bear, 1988, p. 204).3 For the subglacial sheet, the final equations are given by, 8Hb rh , fdHh\ dp* _ 1 jdH* d(H*ub) ( l - n ' j p ' + y - l ( 2 W h ) ~dt~~ ^ \ ~ d T + ~ds (l-n>)pipw mj • ^Abb> In these equations, fh is the mass of melt per second per square meter or m per unit width. The sheet closure velocity (dH/dt)close is the subject of Chapter 3 and will not be given additional attention here. Expanding equation (2.14b) using the closure relationship yields the evolution of the cross section for the R-channels, S S b - m - M W - ^ f l E ^ ! ) " , - ( M T ) dt p>(l-nj,) ~ " ' V nB where p' the is ice overburden pressure, n is the flow law exponent, and B is the flow law coefficient, or fluidity of ice. The term on the extreme right of (2.17) is the Glen's law formulation of creep closure that Spring & Hutter (1981) used for conduits. The mass balance system for the channels is given by equations (2.15) and (2.17). These equations can be contrasted with those for the sheet in equation set (2.16). I will present the momentum and energy equations twice also: once for the channels and once for the sheets. Because these equations have a melt rate term, they affect the overlying ice by melting or freezing along the melt perimeter. The thickness of the ice changes as ice is accreted or melted from the base. For the sheet, a simple mass balance term for overlying ice is given by ^ = - , * . . (2.18) at ( i - n j , y In this equation, the total derivative is equal to the negative melt rate term. Because I neglect ice flow, there is no convective velocity. While glaciers are rarely static, adding a convective term and the corresponding strain terms for the ice is beyond the scope of the present work. In addition, I use the notation Z1 for the thickness of the glacier itself (Fig. 2.1). When the porosity is reduced to zero, equation (2.18) reduces to a solid accretion rate. 3 I n Chapter 5, I show the details of the compressibility assumption in deriving the pressure evolution in the englacial water system. 30 Chapter 2. Subglacial flow equations An alternative to equation (2.18) is to define an accreted ice thickness Zb:e. The evolution of the accreted ice thickness is simply dzb-.e m dt sheet. (1 - nl) p1 Furthermore, the thickness of accreted ice for the channels is dZh[' = < 1 Rb + Zb:e Rb + Zb:e m + Za: b:e 2TT(1 - n')p' ' 2TT.R dsb * 1 ~bT close m _7r(l - ni)p" + nR5 \ — / cioseJ where the closure term on the right hand side accounts for the contraction of the accreted ice layer. Equation (2.19) does not account for this contraction because ice moves vertically into and out of the sheet without the radial divergence present in the case of channels. z»-b:e osb ~bT circular channel, semicircular channel, (2.19) (2.20) Porosity of the accreted ice Ice that accretes may be discontinuous and preserve a residual porosity. Direct observations of ice actively accreting subglacially are nonexistent. Instead, observations made at glacier termini are commonly used as proxies for the subglacial environment. Frazil ice near the terminus of Matanuska Glacier forms in water that ascends vents local to the terminus (Lawson et al., 1998). When water ascended a shear fracture extending from the bed to the surface of SkeiQararjokull frazil ice formed (Roberts et al., 2002). It is not immediately clear that water that ascends steep local features is representative of the subglacial environment. However, the formation of frazil is certainly a viable method of accreting ice and appears to be dominant in subglacial environments. Another method of accreting ice would be the formation of a mushy layer (e.g., Worster, 1997). Mushy layers typically form in low-velocity flows with the prominent analog being sea ice growing downward at the sea surface (e.g., Wettlaufer, 2001). Depending on conditions, the porosity in mushy layers can be high. In contrast, Gilpin's (1981) experiments of flow through closed conduits showed that accreting ice is not very porous. However, these experiments were run with high pressure gradients and high velocities. With a wide range of porosity to choose from depending on the choice of analog for the subglacial environment, it stands to reason that the accreting ice will be somewhat porous. Moreover, porosity of an accreted layer will be a function of flow conditions. 31 Chapter 2. Subglacial flow equations In the absence of an exact relationship for the evolution of the ice porosity, np becomes an adjustable parameter that affects the cross-section evolution, water pressure, and, for the sheet, ice thickness. The rate of change of the hydraulic cross-section is governed not only by melt but by the addition of water that is in the pore space. The rate of water addition from the pore space is the final term on the right hand side of equations (2.13). The rate of change of the thickness of the water sheet resulting from melting or freezing is rh/ [(1 — np)pl]. This relationship appears in equations (2.16a), and (2.18). 2.2.2 M o m e n t u m balance Steady flow While flow may be steady and uniform, the supply and flux of momentum are governed by the geometry of the ice—bed interface and the overlying ice. As is typical of steady conditions, momentum supply from gravitational head loss and a change in pressure gradient is balanced by shear along the perimeter (e.g., Bird et al., 1960, p. 182), — ^Ap b + pwgAzh^j Sw — PwAs TO, (2.21) where A p b is a change in water pressure, g is the gravitational acceleration, Azh is the elevation change of the ice—water interface, and ro is the shear along the hydraulic perimeter. Moving this macroscopic balance to a differential balance is straightforward, F w ds y ' The left hand side of equation (2.22a) is identifiable as the scaled negative gradient of the hydraulic potential, 0b = p h + p*gzh. (2.23) The Darcy-Weisbach formulation is an empirical relation for the shear stress along the hy-draulic perimeter,4 ro = lpwfd(uh)2, (2.24) 4 B i r d et al. (1960) gave their derivation in terms of the Fanning friction coefficient / / . However, Schlichting (1979, p. 597) gave a similar derivation using the Darcy-Weisbach friction coefficient fd- These two coefficients are related v ia / / = /d /4 . The difference of 1/4 results from a discrepancy in transferring the coefficient from Hagen-Poiseuille flow to turbulent pipe flow. The Schlichting form is more common in the fluvial geomorphology 32 Chapter 2. Subglacial How equations where fd is a friction coefficient. Clarke (2003) suggested that the friction coefficient was a combination of the ice and bed friction coefficients, fd circular conduit fd = < J d KpW } J d semicircular conduit (2-25) i (fd + fSd) ^eet where fd and / j are the friction coefficients for the ice and bed, respectively. In equation (2.25), I make the simple assumption that the roof and floor of the sheet account for an equal portion of the shear stress; but they need not be equal depending on roughness features of the ice or characteristics of the bed. For turbulent water flow, fd must be assigned as parameter. Values for fd for relatively smooth open channel flow range from 0.0012 to 0.06 (Schlichting, 1979, p. 617, Fig. 20.18). Clarke (2003) used hydraulic roughness parameters equivalent those for rough streams to simulate outburst floods. For laminar water flow, the friction factor is usually based on the Reynolds number. Because these relationships are based on data for smooth-walled conduits and have low accuracy (Bird et al, 1960, p. 187), the laminar relationships are not of interest here. A relationship for the velocity follows from equations (2.22b) and (2.24), b i 8 s b v u = — I sgnl i , , i 8 Hh\1 sgn( d(f)h\ d(j>b ds J ds d(f>b\ ds ) ds 2 channel, (2.26a) sheet. (2.26b) ,Pwfd 2 Stone & Clarke (1993) used a similar relationship for turbulent flow in a horizontal water layer beneath a glacier. Other shear stress parameterizations could be used to create equations (2.26). For example, the Manning-Strickler formula could be used in the place of the Darcy-Weisbach formula. Clarke (2003) used both but reasoned that the Manning-Strickler formula is more realistic for floods. literature; yet both are used and often carelessly given the notation / without proper identification. I thus adopt fd as the more common Darcy-Weisbach coefficient in place of the Fanning coefficient / / . 33 Chapter 2. Subglacial flow equations Time dependent flow Modification of the velocity evolution of Spring & Hutter (1981) to include the porosity of the ice gives, dub hdub {l-np)p[ + nppw mub 1 dpb dzb P w r 0 , . . ——— = — u — 1—^—rr—r*- - r « — r — channel, (2.27a) dt ds (1 - np) pl pwSb pw ds y ds Sbpv< ' v ' dub _ hdub (1 - np) pi + npp™ fhub 1 dpb dzb 2r 0 ~ d T - ~ U ^ i (l-nL) ^ ^ P^~ds~ 9 ds Hbp™ S h 6 e t ' { l M b ) where both equations include the porosity change discussed previously. These equations are derived from the material balance of the local momentum. Multiplying these equations by the water density and cross sectional area would yield the momentum balance in its original one-dimensional form. Beginning with the left hand side, the first two terms are the material derivative of the velocity. The next term, with m or rh, is the impulse lost (or gained) from melt (or accretion) along the hydraulic perimeter. The next two terms result from the change in the hydraulic potential along flow. The final term on the right is shear along the hydraulic perimeter. These equations are subject to the supplementary relations in equations (2.24) and (2.38b). Equivalent discharge The amount of meltwater that reaches the ice-bed interface will be independent of whether the water runs through sheets or conduits. Water discharge must then be equivalent for either a sheet or a channel. Exact equivalence of discharge might never occur subglacially. For example, in examining Figures 2.1 and 2.2, there are several competing features among the cross sections. In particular the momentum expended on shear along the hydraulic perimeter will be much greater for the sheet than for the conduits despite the expectation that hydraulic potential will be the same for all the cross-sections. For the steady case, the discharge Qb = Sbub can be set equal using equations (2.26). The radii of the channels can then be found in terms of the geometry of the sheet, 3 2 Rb = (# b )* (^ r Y circular channel, (2.28a) 3 Rb = (4TT + 8) s I — J W$ semicircular channel. (2.28b) 34 Chapter 2. Subglacial Bow equations Time dependent flow is not necessarily subject to these relationships because the flow will adjust according to the upstream boundary conditions. 2.2.3 Energy balance Time dependent flow Spring & Hutter (1981, 1982) derived the energy balance equation for one-dimensional pipe flow incorporating the ice-water phase change. Expanding the material derivative and rearranging the heat balance yields a partial differential equation for temperature, g T b _ g r b m ( ( 1 _ n i ) p . + n i p w ( ^ b ) 2 | p W T Q u h dt U ds p"c"Sb\+^> m p ( l - n l ) p i 2 j+Sbp"<%' (2.29a) dTb _ hdTb m f (1 - n') pi + (ub)2 \ 2 r0ub dt ~ ds P™c™Hb\h + C ^ l m v ( l - n j , ) ^ 2 J + ^ p w c w ' (2.29b) where the specific heat of water is , and equations (2.29a) and (2.29b) are for the channel and sheet, respectively. In these equations, the time derivative of the change in energy per unit volume is equal to the changes in latent heat L, sensible heat C™ATmp, and specific kinetic energy (ub)2/2 released through melting, m or rh, plus the shear heating along the hydraulic perimeter PwToub. The specific kinetic energy term results from the impulse term from equation (2.27) being dotted by the velocity to make an energy term. The primary difference between the energy equation for conduit flow and a distributed macroporous sheet is the hydraulic perimeter and cross-section of flow. Again, equations (2.3) provide supplementary relationships for one of the hydraulic cross sections. Steady flow: Rothlisberger heat transfer Because the water is always in contact with ice along its flow path, the only reasonable inter-pretation of the assumption in equation (2.9b) is that the temperature of the water must be equal to that of the overlying ice. Additionally, the overlying ice and underlying bed must be at the local melting point. As a result, the subglacial temperature must remain at the steady value of the pressure melting point. The assumption of the time rate of change being zero fixes the melt rate, and the energy balance becomes a statement of the melt rate. 35 Chapter 2. Subglacial How equations Because the water is at the melting point, the sensible heat term is unimportant in isolating the melt rate from (2.29). I also discard the energy imparted from the impulse term because it is approximately five orders of magnitude smaller than the latent heat term. Rewriting equation set (2.29) for the melt rate gives the equations, Shuh ( w w8Th , PwrQuh m = - ^j- { P W C ; - Q - - + — 5 5 — ) channel, (2.30a) Hbuh ( w w0Tb , 2r0ub ™=" ^ {pWc?lk7+ ^IW) sheet- ( 2 3 0 b ) Another modification to the melt rate is to use the chain rule in conjunction with the inverse Clausius-Clapeyron slope (eq. 2.8) for the change in temperature with pressure, ™ - a * (2 3i) By substituting equation (2.31) into equation (2.30) and parameterizing the shear stress by equation (2.22b) in equation (2.30), I obtain a melt rate that is dependent on the discharge through the water horizon, the water pressure, and the hydraulic gradient, m = " x ( ' " ^ ¥ + ¥ ) c h a n n e 1 , ( 2 ' 3 2 a ) + s , K e t . (2.32b) This melt rate falls directly out of the assumptions proposed by Rothlisberger & Lang (1987). As a result, I refer to this steady state melt rate as "Rothlisberger heat transfer" in a similar spirit to Clarke's (2003) reference to Nye's (1976) heat transfer assumption.5 2.3 The melt rate A heat flux formulated from Newton's law of cooling (e.g., Bird et al., 1960, p. 267) gives an estimate for the heat that can be conducted over a given length scale, which relates directly to the melt rate, qT = hAT, (2.33a) m = ^ . (2.33b) 5 Nye's heat transfer assumption is that dTb/dt = 0 whereas Rothlisberger's is dTb/dt = 0 (eq. 2.9b). 36 Chapter 2. Subglacial flow equations Here, h is the heat transfer coefficient, t\ is a length scale, A T is the temperature change across that length scale, and is the thermal conductivity of water. The heat flux qr is the magnitude of the heat flux vector qx normal to the phase change boundary, which is the melting perimeter in this case. The heat flux is in W m~ 2. For water flowing in either of the cross sections, there is an inner region of the flow where the temperature is relatively homogenous because the water is well-mixed. Near the boundary, heat is conducted across a boundary layer. The length scale t\ is the height of the conductive boundary layer Sjr, between the well-mixed inner region and the ice wall. The Nusselt number is the ratio of the total heat transfer to the conductive heat transfer. A Nusselt number greater than unity indicates that there is convective heat transfer. Multiplying the conductive heat flux in equation (2.33a) by the Nusselt number gives a heat flux that includes convection. After this multiplication, the length scale l\ ceases to be the conductive boundary layer and becomes the entire thermal boundary layer. 2.3.1 Semiempirical correlations There are numerous empirical relationships for the heat transfer coefficient from the engineering literature, in general, these fall into a broad class of correlations for the Nusselt Number, Nu = CiReroPr?, (2.34) where C\ is a nondimensional constant, and w and c are indices (e.g., Hutter, 1982; Schlichting, 1979). In this equation, the Reynolds number (Re) and Prandtl number (Pr) are given by the usual definitions, f channel, Re = < P. (2.35a) sheet, V LI W Pr=— = (2.35b) where fi is the viscosity of water, v is the dynamic viscosity of water, and is the thermal diffusivity of water. The Reynolds number is the ratio of the inertial forces to the viscous forces. The distinction between the thermal boundary layer and the inertial boundary layer is given by the,Prandtl number. For water at the freezing point, the Prandtl number is approximately 13; or the thermal boundary layer is 13 times the inertial boundary layer. 37 Chapter 2. Subglacial flow equations McAdams (1954) presented several formulae for flow in ducts, pipes, and sheets. From these relationships, the Dittus-Boelter heat transfer relationship for flow in turbulent pipes stands out because of its use in outburst flood studies (Mathews, 1973; Nye, 1976; Spring & Hutter, 1981, 1982). The Nusselt number becomes (McAdams, 1954, p. 219), Nu = 0.023 R e 4 / 5 P r 2 / 5 . (2.36) Alternative formulations exist for the Nusselt number, but these are less popular. For example, Isenko et al. (2005)6 used the relationship, Nu = 0.003 Re P r 2 / 5 , (2.37) to interpret experimental work. For the remainder of this analysis, I employ the more popular Dittus-Boelter relationship (eq. 2.36). However, any relationship in the form of (2.34) could be suitable. Because the height of the boundary layer is not obvious from the momentum equations (eq. 2.27), the Dittus-Boelter relationship substitutes the characteristic hydraulic cross-section length scale for the boundary layer thickness. For conduits, the characteristic length is 4Rh, which is the diameter of the circular channel or the radius of the semicircular channel. The same heat transfer in a macroporous sheet depends on the thickness of the water layer instead of the diameter. The mass melt rate per unit channel length can be formulated with appropriate substitutions of equations. The melt rate per unit area m is simply the equivalent melt rate for the channel (eq. 2.38a) divided by the melting perimeter, 0.023 Re^Pr^Jey rATmoP 1 m k T " ' m p channel, (2.38a) . 0.023 R e 4 / 5 P r 2 / 5 ^ A r m p m = H^L— sheet. (2.38b) Ultimately, the Dittus-Boelter relationship approximates heat transfer across the turbulent boundary layer. This relationship is found via equations (2.33a) and (2.36), 7 7 ^ — c h a n n e l , 0.023 R ^ 5 P r 2 / 5 ( 2 3 9 ) 0.023 R e 4 / 5 P r 2 / 5 S h e 6 t ' 6Isenko et al. (2005) gave their relationship in terms of a volumetric energy density per degree multiplied by the along path velocity. The Nusselt number presented here is equivalent to their formulation. 38 Chapter 2. Subglacial How equations A rough estimate of the thermal boundary layer for water flowing in sheets is 3 x 10 m— assuming Re = 2100, the transition from laminar to turbulent flows in plane-parallel apertures, Pr = 13, and Hh = 0.01 m. 7 For a circular channel a radius equivalent to a sheet with unit width (2.28), the thermal boundary layer is approximately 1 x 10~3 m. 8 While there is no guarantee that these examples are compatible with the balance equations presented thus far, this value does suggest limitations on the height and radius of the subglacial cross sections. Much of the heat transfer literature rests on a few correlations, one of these being the Dittus-Boelter relation. The Dittus-Boelter heat transfer relationship was originally designed for heat transfer within smooth-walled pipes. While the ice ceiling above the flow and the bed below are most likely not smooth, there are few published engineering Nusselt number correlations for rough boundaries. Other correlations exist for flat plates, concentric tubes, and tube bundles, but it is not clear that additional insight or accuracy can be gained by using these for rough-walled boundaries. Bird et al. (1960, p. 401) pointed out that these correlations can be applied but are less reliable in the transition region, defined approximately as Re = 2100-10000. Despite this limited range, the relationship proposed in equation (2.38b) appears to be the most relevant above the laminar flow regime. Within the laminar flow regime, Bird, et al. (1960, p. 399) suggested where an equivalent scaling to equation (2.39), assuming Hh/As as 0.01, gives the Nusselt number as approximately 6 for the laminar case and 27 using equation (2.36). In general, water flowing in subglacial sheets is not expected to be flowing below a Reynolds number of about 2100 for the model proposed here; however, I leave this option available should the flow ever drop below this value. 7The corresponding velocity from equation 2.35a is approximately 0.4 m s - 1 . 8The corresponding velocity from equation 2.35a is approximately 0.1 m s _ 1. channel sheet (2.40b) (2.40a) 39 Chapter 2. Subglacial flow equations 2.4 Steady state accretion rate Alley et al. (1998) arrived at an accretion rate via a heat balance, _rn = Q W ^ + { P W _ P L C 2 ) A B} ( 2 4 1 A ) C 2 = l - / i / c > w , (2.41b) where tana r is the glacier surface slope, and tana b is the energy slope of the water layer (Fig. 2.1).9 Furthermore, Alley et al. made small angle approximations such that a1 ~ tana r and a b ~ tanoA The discharge Qh in (2.41a) is constant such that both Hh and uh are both held constant. The bed slope and surface slope of the glacier are also held constant along the flow path, similar to Figure 2.1. The final cosa b term is a result of the flow path's inclination relative to the horizontal. Equation (2.32b) is essentially the relationship that Alley et al. (1998) outlined previously (eq. 2.41) without some of their additional assumptions. They reasoned that the water pressure gradient, dpw/ds is adequately approximated by the ice overburden pressure gradient dpl/ds. Furthermore, Alley et al. (1998) assumed a wedge-shaped ice configuration such as the area near the outlet in Figure 2.1. These assumptions lead to Hbub fdph . w . w dzb\ m = - — l — (p cpP + l)+p 9 ^ \ , (2.42a) m W • il^^L. ( p i ( S i n ar _ g i n ab^ (pwcw£ + !) + pw g i n ab^ (2A2b) After grouping the terms like Alley et al., an equation that parallels their work is quickly realized, m qHhWuh cos ah / h / w _ :. _ :. r.\ ,n A n . — r = r- (tano;b(p - C20) + C2P ( - tana )) . (2.43) p1 p'L \ / The assumption inherent in equation (2.43) is that cosa b ~ coso:r ~ 1. This assumption may not be valid depending on the slopes of the surface and bed but may be for relatively flat glaciers. Equations (2.41) and (2.43) illustrate that the Spring-Hutter formulation of the water flow equations can be reduced to a steady state representation. By setting equation (2.43) equal to zero, the crossing from melting to freezing can be found 9 I have modified Alley et al.'s (1998) notation to match the present formulation. In addition, the surface and energy slopes are taken to be positive numbers. 40 Chapter 2. Subglacial flow equations in terms of the surface and bed slopes, n _ s i n a b ( l + ^ P - V tano* sina r p w - ( l + Pc^pw) p1 tana r Rothlisberger & Lang (1987) gave this ratio as —2.02 for pure water and —1.30 for air-saturated water. Alley et al. (1998) gave the ratio as —1.7 for pure water and —1.2 for air-saturated water. Hooke (1991) also computed the ratio and found it to be —2.0 and —1.5 for pure and air saturated water, respectively. The large variability in these numbers results from small changes in the constants in equations (2.44) and (2.7). While the assumption that surface and bed topography are constant must be false, a simple ratio to predict glaciohydraulic supercooling from topography is extremely helpful. 2.5 Summary Throughout this chapter, I have focused on glaciohydraulic supercooling in a subglacial water system. Glaciohydraulic supercooling occurs because water flows upwards in a hydraulic sys-tem without warming sufficiently to maintain water. The change of phase is governed by the Clausius-Clapeyron relation (eq. 2.8). Three balance laws—mass (eqs. 2.15, 2.16, and 2.17), momentum (eqs. 2.26 and 2.27), and energy (eq. 2.29)—provide a mathematical representation through which water cross-section area or height, water pressure, water velocity, and water temperature can be simulated. Each of these equations has steady-state and time-dependent formulations. These formulations are subject to phenomenological conditions for heat transfer between the ice and water (eq. 2.32). In addition, the equations are subject to geometrical conditions on cross section (eqs. 2.2, 2.3, 2.4) and flow path (Fig. 2.1). For steady-state conditions there is a condition on the surface and base angles of the glacier, indicating the propensity towards supercooling (2.44). 41 Chapter 3 Sheet closure mechanisms The aim of this chapter is to provide a concise relationship for ice intruding into a water sheet via the processes of regelation and viscous creep. A water passageway can close or open depending on the difference between the water pressure and the stress in the ice acting toward the water. For example, if the water pressure is below that of the ice overburden, ice will intrude into the cross-section, contracting it. If the water pressure is greater than the ice overburden pressure, it may be possible for the water to push the ice, dilating the cross-section. For a conduit idealized as a tubular volume with a circular cross section, the relationship for closure is power-law creep (eq. 2.17). Creep closure of conduits in ice was observed and quantified by Glen (1952, 1954). To explain these observations, Nye (1953) constructed an extensive theory. Creep closure for conduits has proven effective in various forms and simulations (e.g., Nye, 1976; Spring & Hutter, 1981). In these situations, the difference between the ice overburden stress and the water pressure determines whether the cross-section is shrinking or growing. Typically, ice overburden squeezes the conduit shut. However, if water is being forced through at a pressure higher than the ice overburden pressure, then the ice walls may be pushed outward, opening the conduit. Examples of water dilating the cross-section occur during jokulhlaups when a large supply of meltwater may force its way along the ice—bed interface (e.g., Bjornsson, 2002). For sheets, the intrusion relationship is not well-defined. This ambiguity may be, in part, a result of sheets not appearing conspicuously at the front or surface of glaciers. As a water system nears the terminus of a glacier, water depressurizes in response to atmospheric pressure. Chan-nelization is then favored at the ice-bed interface (e.g., Fowler, 1987; Walder, 1982). Depressur-ization may explain the lack of observations of sheets emerging at the front of glaciers. Because observations of sheets are scarce, physical models help explain their development. Walder (1982) analyzed the dilation and contraction of water sheets and appraised their stability with respect to channelization. His exposition has not been extended but merits further attention because he was interested only in flat glacier beds. 42 Chapter 3. Sheet closure mechanisms Walder (1982) investigated the stability of a laterally continuous sheet of water at the glacier bed. He examined the two-dimensional, steady-state Navier-Stokes equations coupled to a simple diffusive thermodynamic model. In his model, ice sags into the water sheet, closing it. Walder approximated the rheology of ice by a linear viscosity rather than a nonlinear, power-law viscosity. He used linear perturbation analyses to determine what thickness of water would yield stable sheets under different conditions. Walder concluded that water sheets with flat beds were almost always unstable and would decay to channelized cross-sections. Sheets with rougher beds were stable to a thickness of about 4 mm. Unstable sheets would be discontinuous and give way to an arborescent drainage system. 3.1 Formula t ion Glacier sliding and ice intruding into a sheet are analogous processes whereby ice moves past obstacles at the ice-bed interface. The difference between sliding and intrusion is the direction of ice movement: horizontal for sliding, vertical for intrusion. This physical similarity opens the body of literature on sliding to closure of a water sheet. Glacier sliding is the result of two processes: viscous creep and regelation (e.g., Weertman, 1964). Physically, creep is the movement of lattice dislocations in ice crystals in response to applied stresses. Regelation is the movement of ice past an obstacle by melting on a high pressure side and then refreezing on a low pressure side. Water flows in a thin film from the high pressure side of an obstacle to the low pressure side. Weertman (1964) superposed the regelation velocity vr and the creep velocity vc to give the total ice velocity v, past an obstacle, v = v r + vc. (3-1) Equation (3.1) is the statement for glacier sliding speed, and this equation also represents the closure speed of ice into the sheet v = dHh/dt (see equation (2.16a)). Equation (3.1) holds for any individual obstacle that the ice is moving past. For ice at the melting point, regelation acts as a slip condition on the surface of the sediment. Creep acts on distances approximately the size of the obstacle that ice is moving past. Both regelation and creep act jointly over different length scales to move the ice along the bed. Weertman (1964) observed that, for very large obstacles, creep will be the dominant 43 Chapter 3. Sheet closure mechanisms mechanism. For very small obstacles, Weertman reasoned that the regelation velocity will be large. Near where the regelation velocity equals the creep velocity, both will be small. These obstacles will then control the speed at which the roof collapses. Nye (1969, 1970) suggested that a controlling obstacle size lies between 2 and 50 cm. In addition, Kamb (1970) concluded that a transition obstacle wavelength exists near 50 cm (Paterson, 1994, p. 142). A similar trade-off should occur for ice intrusion. 3.1.1 Regelation closure Regelation is a thermodynamic process where ice in contact with another material melts because of the local change in the pressure melting point. Closure of the water system via regelation is not understood nearly as well as the effect of creep closure. For a cross section idealized as a sheet less than the transitional obstacle size, regelation should be the sole mechanism acting to close a water layer. Nunn & Rowell (1967) performed laboratory experiments that have become classic examples of regelation. In these experiments a wire is loaded such that the pressure of the leading edge exceeds the local melting temperature of the ice (via eq. 2.8). Water produced flows in a thin film around the wire to the low pressure side, where it refreezes. Latent heat released from refreezing is conducted back through the wire, enhancing melt. Townsend &; Vickery (1967) performed similar experiments and obtained similar results by pulling a metal ball though a block of ice. Weertman (1972) and Nye & Frank (1973) gained insight from these simple experiments to explain a glacier sliding past obstacles at the ice-bed interface and, subsequently, to determine sliding speed over a rock bed. Nye (1967) constructed the theory for these experiments. He considered ice that is clean and separated from obstacles by a thin film of water and then solves for the force required for a steady regelation velocity. Nye noted if ice is below the melting point, regelation ceases unless there is an additional heat source. In this case, heat returning to the area of melting must first raise the temperature of the ice to the melting point before melting the ice. Because of this feature, I consider cases where the sensible heat is negligible. Other assumptions are inherent in regelation: ice is clean, there is a thin film of water separating the ice from the sediment grains, and transient effects are negligible. Also, the mechanical properties of ice, such as fabric orientation, should not change when ice melts and then reforms. 44 Chapter 3. Sheet closure mechanisms Ice flow past bedrock obstacles implies that ice flow is roughly parallel to the bed for sliding. For the closure of the ice—bed interface in the presence of a water sheet, the direction of ice motion must be into the sheet, or nearly vertical. Figure 3.1 depicts ice resting on a sediment grain at the ice-bed interface. On the leading edge of the grains, ice melts. The water generated then flows around the sediment grain toward the water sheet. "Regelation" in this sense is somewhat of a misnomer because water is not necessarily refreezing on the low pressure side. However, the theoretical basis for regelation is applicable provided that heat is conducted through the sediment grain to melt the overlying ice. o Figure 3.1: A simple picture of regelation closure. The water layer is the lightest gray. As the ice descends onto the grains, stress from the overlying ice is concentrated, and the ice melts. Grey arrows show the sense of ice motion. Black arrows indicate motion of water generated from ice melt. Nye (1967) gave the formulation for regelation in terms of a single cylinder and altered the shape factor to account for spheres. He then gave an approximate formulation for any shape, Vep'LVr F = (3.2) PKe:T ' where Ve is the effective volume of the shape, Ke:r is an effective thermal conductivity, vr is the speed at which the ice is regelating past the object, and 8 is the linearized pressure-melting coefficient (eq. 2.8). Rewriting (3.2) in terms of the normal stress available for regelation gives P1L Ve 8Ke:T A -Vr, (3.3) 45 Chapter 3. Sheet closure mechanisms where Ae is an effective area. The effective area is an area of ice—object contact. The negative sign in equation (3.3) indicates that the direction of closure is opposite to the positive z-direction, or downwards, when the ice overburden stress is greater than the water pressure. The definition of the effective area Ae in equation (3.3) involves both a sediment grain and an applied stress. Philip (1980) noted that the stress is the mean normal stress, and he defined the force as acting along the direction of ice flow. However, this is at odds with Nye (1967) who used the stress normal to the surface of the obstacle. In this case, Nye's formulation is preferable. Ae is then the surface area opposing the motion of the ice. If there is a known relationship between the surface and cross-sectional areas of sediments, then the cross-sectional area could be substituted. The effective area is not the area of all grains on the bed, but the effective cross sectional area of those sediment grains emerged into the ice and partially submerged in the water. Obstacles that are completely submerged by water do not have an effect. Similarly, Ve is the volume of an obstacle encased in ice at the ice-bed interface that the ice is regelating past. Philip (1980) investigated the effective thermal conductivity. He constructed an analytic model of ice regelating past arrays of equally spaced obstacles, cylinders or spheres, having an adjustable thermal conductivity. Philip solved for the effective thermal conductivity based on a given array geometry and an obstacle shape. If the thermal conductivity of an array is the same as that of ice, the effective thermal conductivity reverts to that of ice. As a result, obstacle spacing is irrelevant without varying the conductivities. For obstacles with thermal conductivities greater than that of ice, heat is preferentially moved through the cylinders. Con-versely, if the thermal conductivity of a cylinder is less than that of ice, a greater relative heat flux occurs through the ice. It is not clear that Philip's (1980) method of applying different conductivities provides adequate compatibility because sediment grains at the ice—bed interface are usually irregularly spaced. However, his study does give an impression of what may happen when the thermal conductivity is varied. In addition, Rempel et al. (2004) noted that varying the conductivity of the grains makes the problem unnecessarily difficult. Iverson & Semmens (1995) avoided the issue completely in their experimental study by setting the conductivities of ice and sediment equal. Because the thermal conductivity appears to have a small impact on the problem, I set the thermal conductivity equal to that of ice as has been done by Iverson & Semmens (1995). 46 Chapter 3. Sheet closure mechanisms 3.1.2 Creep closure Creep is the process whereby ice deforms under an applied stress. For sediment grains with mean dimension larger than the transitional obstacle size, the primary mechanism of ice intrusion into the macroporous water layer will be power-law creep. The most common form of the flow law for ice was given by Glen (1954) where a strain rate is proportional to a stress raised to a power. If ice rests on grains that concentrate stress then ice will preferentially creep around these grains. Figure 3.2 shows the conceptual framework for creep closure of a sheet with ice sagging between two large sediment grains. If the ice were lower in the figure, overburden stress would be dissipated because the ice would rest on more grains. This dissipation may reduce the rate of creep closure. Figure 3.2: A simple picture of regelation closure. Ice preferentially sags into the water layer with a larger spacing between particles. Same color scheme as Figure 3.1 Creep is driven by the stress deviator at the ice-bed interface. Equivalent shear r can be written in terms of the stress deviation from a hydrostatic state, where uc denotes a deviatoric stress available for creep. Weertman (1964) noted that there are two deviators. One deviator results from the compression on the obstacle. The other is a tension from the ice that passes the obstacle. This is clearly seen as (3.5) retains a 1/2. In addition, 2r 2 = cr2 + 2r 2 (3.4) 47 Chapter 3. Sheet closure mechanisms because shear from down glacier ice movement is neglected, T X Z in (3.4) can be neglected10, 2 r 2 = ( a c ) 2 , r=\j\ (°c)2. (3.5) Because the water sheet experiences only vertical motion, the bed must be in compression or tension. As a result, the 1/2 in (3.5) is not critical. Weertman (1964) also argued that the 1/2 is not important for glacier sliding where obstacles at the ice-bed interface are in both tension and compression simultaneously, and r = crc. (3.6) The strain rate in the ice can be approximated by the flow law (e.g., Paterson, 1994). Substituting the creep closure stress directly into the flow law yields the strain rate in terms of the creep closure stress, e2 = Asgn(ac)\oc\n, (3.7) where ez is the vertical strain into the sheet and A is a flow law coefficient. Multiplying (3.7) by a length scale produces a velocity, vc-Asgn(-ac)\oc\nle, (3.8) where le is the creep length scale. The effective length scale provides an approximation to the spatial integral of the strain rate without solving the integral directly. Solving the integral directly would necessitate solving for the stress distribution in the ice. An attempt to integrate the stress would not necessarily give a better approximation because an integration would be dependent on the boundary conditions, which are not straightforward for a bed of irregularly shaped objects of variable horizontal spacing. By using a simple length scale, equation (3.8) yields a tractable approximation of creep closure. The creep closure stress can be found by rearranging equation (3.8) f \ v c \ \ l , n ac =-sgn(vc) i-^-j . (3.9) 1 0 By neglecting shear, I preclude the circumstances that would allow a sedimentary bed to develop cavities. These cavities might take an elongate form (e.g., Weertman, 1986) that could then coalesce into laterally con-tinuous cavities (e.g., Kamb, 1987; Walder, 1986). This coalescence would lead to configurations similar to those discussed by Alley et al. (1998). 48 Chapter 3. Sheet closure mechanisms The natural scale for the effective length in equation (3.8) is the distance between the sediments at the bed because an ice roof will preferentially sag between supporting grains. Creep around the individual grains would be secondary to the sagging. This length scale will vary depending on the thickness of the water sheet. 3.2 B e d par t i t ioning The overlying ice rests on an unlithified sedimentary bed and water sheet. Within this water sheet, many smaller sediment grains are wholly submerged. Larger grains may penetrate the water to partially support the ice and stabilize the water layer. A combination of water pressure and sediment stresses jointly support the weight of the overlying ice. As a result, the total cross-sectional areas of water A w and sediment As sum to the total area of the ice, A[ = AW + AS. (3.10) In this case, A1 is a cross-sectional area of ice normal to the direction of closure. The ice seals the water sheet and partially rests on sediment grains. This equation is consistent with the ice-water configuration of a sheet from Chapter 2. Many different grain sizes may be present on the bed. Some grain sizes will be fully sub-merged in water and not touch the ice ceiling. The remaining grain sizes will support some of the weight of the glacier. The total contact area of sediments touching the ice A s is equal to the integral As= f f flasdaadAl (3.11) where is the number of grains per unit size (in area) of a given texture per unit area of the bed. The quantity f% is a number density. Each of a variety of grain textures has an individual fractional area as. One of the geometric necessities of equation (3.11) is that As is defined at the average elevation of A1, which is zh. The cross-sectional area normal to intrusion for one grain size is then dependent on the mean size of the grains, the level of submergence, and the number of that grain size on the bed. Implicit assumptions include the sediment grains being exposed to the same level and having a cross-sectional area and volume that can be approximated by a known function of the mean grain size. 49 Chapter 3. Sheet closure mechanisms 3.2.1 Stress partitioning If the ice is resting on more than one grain size across the water horizon then the ice overburden stress can be partitioned among the different grain sizes and the water pressure. A reasonable analog for stress partitioning emerges from work on sliding and cavitation at the ice-bed interface (e.g., Fowler, 1987; Lliboutry, 1979; Schoof, 2005). Fowler (1987) approximated a series of bumps at a glacier bed as a Fourier series. In his formulation, the largest wavelength bumps have the largest amplitude, and smaller wavelength bumps are superimposed on the larger bumps. When the glacier's sliding speed reaches a threshold, water cavities develop behind the smallest bumps. As the sliding speed of the glacier increases, larger cavities form behind the larger bumps. These larger cavities drown a fraction of the smaller bumps and cavities. In order to cope with the areal changes in normal stress and water pressure at the bed, Fowler used a stress renormalization technique introduced by Lliboutry (1979). A similar formulation can be used for ice intruding into a sheet of water. The downward force of the overlying ice F1 is partitioned among all grain sizes available at the bed and the water sheet, pi ^ F w + Fs^ AW =Awpw + A\oi + ... + AsNasN, N = A > W + ^ A ^ S , (3.12) i=l where Fs and F w are normal forces on the sediment and water, a1 = plgZl is the ice overburden stress, A\ is the total area opposing flow for the zth grain size, and erf is the stress on the ith grain size. Shear from ice movement is ignored as are any other stresses. Partitioning equation (3.12) into intermediate stresses Oi yields a renormalization, AW' =A\a\ + (A1 - A\) ax (3.13a) (A- A\) ax =A\a\ + (A - A\ - A%) o2 (3.13b) (A - ... - ASN^) aN-i =AsNo% + A w p w (3.13c) In these equations, larger i values indicate smaller grain sizes. Equations (3.13) are predicated on the largest grain sizes supporting a significant portion of the ice overburden stress. The 50 Chapter 3. Sheet closure mechanisms remaining stress is distributed among the smaller grains. The summation of the N equations (3.13) yields equation (3.12). Renormalization of the stresses is the crucial step that differenti-ates this work from that of previous workers who were interested in both regelation and creep (e.g., Weertman, 1964). Solving (3.13) for cr| yields a series of equations with a subtle relationship, * i = ^ t e - i - * o + 7h> (3-i4a) = % i * « i + "i> ( 3 - 1 4 b ) where ae:i is the relative effective stress for the zth grain size, and i = 1... N . The area AQ is equal to A 1 for the largest grain sizes, and er/v is equal to pw for the smallest grain sizes. One grain size: a simple case If there is only one grain size on the bed, then there is a simple stress division among the water and sediment. The stresses for regelation and creep are the deviatoric stress on the sediment grains and the effective stress, respectively. This stress must be obtained via a force balance on. the bed, whereby the mass of the overlying ice must be supported by either the water or the grains at the base of the glacier: JAX =aSAs+pWAW, as=^ae+pw. (3.15) In this case, the three stress components available for ice closure are the vertical stress in the ice a1, the water pressure pw, and the stress exerted by sediment grains on the ice as. 3.3 Stress-velocity relat ion There is one equation each for the total closure velocity (eq. 3.1), regelation (eq. 3.3), creep (eq. 3.9), and a constitutive relationship for stress at the bed (eq. 3.12). To solve these equations, three more relationships are necessary, o-r-.i =o-! - ri and (b) one when Hb < ri, where Vi is the individual grain radius. For the case of sediment grains lying on zs, the effective volume is then given by VKi = | [2n - Hb) ^ + riHb - (tfb)2} • (3.18) The effective cross section is the area of the disk, n{2Hbri- (# b ) 2 } for Hb > Ae..i={ I- • v ' J • (3.19) irrf for Hb 0.9v, C (creep regime) denotes vc > 0.9v, and M (mixed regime) denote that neither process is controlling. These and subsequent figures are constructed around the grain size(s) rather than realistic sheet water depths. Water depth of the sheet would not be as high as 0.4 meters but are probably less than one tenth that height (e.g., Lawson et al., 1998). There are several striking features in Figures 3.3 and 3.4. The first feature is that regelation is an important process for several parameter choices. In the case of Figure 3.3a, the effective stress driving closure is relatively low. As a result, the creep velocity is quite small, and the 1 1 Al l figures that report closure velocities use positive values. These are absolute values of the closure velocities. Al l velocities are downward toward the bed. No closure velocities for negative effective pressures are plotted. 57 Chapter 3. Sheet closure mechanisms Table 3.1: Closure parameters. Parameter Value Units Notes A 6.8 x 10~2 4 s- 1 P a " n (Paterson, 1994, p. 97) 4217.6 J kg - 1 K - 1 9 9.81 _2 m s Ke:T 3.3 W m -1 K - 1 L 3.336 x 105 J k g - 1 n 3.0 (Nye, 1953; Paterson, 1994) P{ 8.99 x 105 Pa Corresponds to 100 m of ice overburden P 7.440 x 10~8 K Pa" 1 From equations (2.7) and (2.8) Pl 916.7 kg m~3 Pw 1000.0 kg m- 3 regelation velocity is a greater part of the total closure velocity. However, throughout most of height of the water layer (as indicated by M), the creep velocity is greater than 10% of the total closure rate. In Figures 3.3b and 3.3c, the relative effective pressure is much higher. The result is that creep is important for most of the grain height. Near 0.4 m for both Figures 3.3b and 3.3c, there is a sharp upturn in the total velocity. This upturn is a result of the ratio Ae/Ve being large as the ice initially settles on a grain. Each creep velocity in Figure 3.3 curves slightly from 0.2 to 0.4 m. This curvature results from the cross-sectional area in equation (3.23) evolving until the equator is reached. Once ice intrudes below the grain equator of the grain, le is constant. As a result, from 0 to 0.2 m, the creep velocities are constant. Figure 3.4 shows the overall pattern for the total closure velocity. In the region where creep dominates, the velocities rise to the left as a result of the rate dependence on n. The bend in the black contours at 0.2 m in the creep region is a result of the change in le, as indicated above. The transition to a mixed regime and a regelation regime near the top of the figure results from the combination of both the ratio Ae/Ve and the ratio Aj/Af being large. From right to left, the transition from a creep regime to a mixed regime to a regelation regime results from a drop in the effective pressure that drives closure. 58 Chapter 3. Sheet closure mechanisms Normalized effective pressure (unitless) Figure 3.4: Total closure velocities for a grain size of rj = 0.2 m with variable effective pressure. Parameters are given in Table 3.1. The fraction of the bed area covered with water is 0.5. These velocities are total closure velocities (eq. 3.1). Black contours are spaced at intervals of 0.25 x 10~6 m s _ 1 . Numbers and thick white lines denote different closure regimes. Area R indicates regelation dominated closure (vr > 0.9u). Area C delimits creep dominated closure (vc > 0.9v). In area M, both creep and regelation are important for closure. Figure 3.5 shows the variability of the closure velocity for a change in the area of water. The practical course of creating a bed parameterization is to define the cross sectional area for a grain and then determine a number of grains based on a predetermined fractional area of water (Aw/A1). The smallest fractional area of water is 0.09 because this is closest hexagonal packing for circles on a plane. A creep regime dominates most of this figure. The transition from a creep to a mixed to a regelation regime toward the top of the figure is a result of the change in the ratio Ae/Ve. The gradual rise in the closure velocities toward the right hand side of the figure is the result of fewer grains and a larger grain spacing. The close contours in the upper right corner result from the product of Ae/Ve and Ai/A\ driving regelation. Figure 3.3b is a vertical slice of Figure 3.5 at a fractional area of 0.5. Because le is dependent on the number of grains (eq. 3.22), the fractional area water is not a direct proxy for the creep length scale. However, le does increase as the fractional area water 59 Chapter 3. Sheet closure mechanisms increases. Figure 3.6a shows the increase in closure rate as a function of le. This figure is a horizontal slice of Figure 3.5 at 0.15 m with the fractional area water converted to le. Figure 3.6a illustrates the creep velocity being linear in le despite the total closure velocity not being linear. According to equations (3.17), the creep velocity should vary linearly with le and the regela-tion velocity should vary linearly with Ai-\/A\. Figure 3.6a illustrates the linear dependence of the creep velocity on le. Figure 3.6b shows how that the regelation velocity is linear in Ai-\/A\. In this case, Ai-\/A\ = Al/As. The water depths correspond to Figure 3.5 and were chosen to show the relative importance of the creep and regelation terms. Fractional area water (unitless) Figure 3.5: Total closure velocities for a grain size of n = 0.2 m with variable fractional area water. Parameters are given in Table 3.1. These velocities are total closure velocities (eq. 3.1). Black contours are spaced at intervals of 0.25 x 10 - 6 m s . Letters and thick white lines denote different closure regimes. Area R indicates regelation dominated closure (vr > 0.9v). Area C delimits creep dominated closure (vc > 0.9v). In area M, both creep and regelation are important for closure. The regelation region saturates the colorbar and has a maximum value of 0.01 m s _ 1 . 60 Chapter 3. Sheet closure mechanisms 0 1 2 3 0 50 100 150 200 250 Effective length: lg (m) A r e a r a t i o : ^ . . / A * ( u n i t | e s s ) Total velocity Creep velocity Regelation velocity Figure 3.6: Example closure velocities for a grain size of ri = 0.2 m. Parameters are given in Table 3.1. The normalized effective pressure is 0.5. Figures (a) and (b) are horizontal slices of Figure 3.5. (a) Hb = 0.15, vc varies linearly with le; (b) min(A w /A') = 0.5, min(7e) = 0.35 m. Small grains As Figure 3.7 illustrates, closure for small grain sizes is dominated by regelation. The ratio Ae/Ve is the controlling term in the regelation velocity with velocities increasing on the right side of the figure as a result of the velocity dependence on Aj_i /Af . Figure 3.8 shows two cross sections through Figure 3.7, one horizontal and one vertical. The creep velocities throughout both figures are several orders of magnitude smaller than the regelation velocities. As effective stress increases, the creep velocity does increase but not to a magnitude where it becomes significant (Fig. 3.8b). The magnitude of closure velocities is approximately a factor of 100 greater for small grains than for large grains as noted in the colorbar (compare Fig. 3.4). 3.5.2 M u l t i p l e grain sizes Sedimentary beds beneath glaciers will most likely not consist of only one grain size. Larger grains will be interspersed among smaller grains. As ice descends, different grain sizes will be submerged to different levels in the water with the ice resting on the remainder of those grains. Figure 3.9 illustrates several scenarios for 5, 10, 25, and 50 different grain sizes lying on the 61 Chapter 3. Sheet closure mechanisms x 10 3 x10 0.1 0.2 0 3 0.4 0.5 0 6 0.7 0.8 0.9 Normalized effective pressure (unitless) Figure 3.7: Total closure velocities for a grain size of r, = 0.001 m with variable effective pressure. Parameters are given in Table 3.1. These velocities are total closure velocities (eq. 3.1). Black contours are spaced at intervals of 0.25 x 10 - 5 m s - 1 . All closure velocities in this figure are dominated by regelation (R), (vr > 0.9v). bed. The maximum grain size in each figure is 0.4 m, and the minimum grain size is 0.001 m. These figures are meant as illustrations of the processes. Grain size distributions for each of these figures are illustrated in Figure 3.10. The fractional areas for each grain size are the same within each distribution (eq. 3.23). Creep, mixed, and regelation regimes are determined by multiplying the individual creep and regelation velocities by the fraction of driving stress on each grain size (eq. 3.14b), summing the creep and regelation components for all grain sizes, and then comparing the magnitudes of those sums. As a result, the regimes indicated in Figure 3.9 are based on the relevant stress distributions. The velocities in Figures 3.9a-d are extremely smooth given the amount of variability mul-tiple grain sizes could introduce. In general, the closure velocities are less than 0.4 x 1 0 - 5 m s _ 1 for all distributions. Structure in the different closure regimes for 5 and 10 grain sizes (Figs. 3.9a,b) results from ice sequentially resting on smaller grain sizes. In Figures 3.9c and 3.9d, the creep regime is absent because each grain size has a relatively large value of Ae/Ve. The process of ice resting on the tops of grains is reasonably continuous for many grain sizes. 62 Chapter 3. Sheet closure mechanisms 0 0.5 1 1.5 2 0 0.5 1 Water depth (m) ^ Q-3 Normalized effective pressure (unitless) • • ^ ^ ^ Total velocity — — Creep velocity Regelation velocity Figure 3.8: Total closure velocities for a grain size of r-j = 0.001 m with variable fractional area water. Parameters are given in Table 3.1. These velocities are total closure velocities (eq. 3.1). (a) A vertical slice of Fig. 3.7 where ae/pl = 0.5 (b) A horizontal slice of Fig. 3.7 where Hh = 0.01 m For water depths of less than 0.1 m (e.g., Lawson et al., 1998; Seaberg et al., 1988) either a regelation or mixed regime is favored to close the sheet. Creep closure is more important for larger water layer thicknesses (e.g., Vogel et al, 2005). Figure 3.11 shows the fraction of stress available for each grain size for the case of ten grain sizes. These six figures represent the largest of those ten grain sizes. Between Figures 3.11a and 3.11b the introduction of the second grain size at 0.22 m divides the driving stress between the two grains. The fraction is not exactly half until the water depth is 0.18 m for high effective stresses. Similarly, introduction of the third grain size at 0.12 m water depth divides the the stresses among the first three grain sizes. Perhaps the most striking feature of these figures is that the stress used by each grain size is never equal across all grain sizes. Larger grains require larger stresses to achieve the same closure velocity associated with smaller grains. Grain burial Instead of using the whole grain parameterization for sediment shape on the bed, the half grain parameterization can be used (eqs. 3.20). This parameterization gives an estimate for grain 63 Chapter 3. Sheet closure mechanisms 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 —Tot/ i „Tot;„i „.Toti„.i ^ T o t / J a la a lo a la a la e e e e Figure 3.9: Closure velocities for multiple grain sizes. Figure labels are in the lower left corner. Parameters are given in Table 3.1. These velocities are total closure velocities (eq. 3.1). Black contours are spaced at intervals of 0.4 x 10 - 5 m s _ 1 . Letters and thick white lines denote different closure regimes, (a) Five different grain sizes, (b) Ten grain sizes, (c) 25 grain sizes. The closure velocity saturates the color scale in the upper right of the figure and has a maximum value of 9.9 x 10"4 m a - 1 , (d) 50 grain sizes. The closure velocity saturates the color scale in the upper right of the figure and has a maximum value of 2.0 x 10~4 m s _ 1 . burial. Relief introduced by the sediments on the bed is reduced. Closure velocity trends and magnitudes are similar to Figure 3.9 and are not plotted. Figure 3.12 displays the stress partitioning for the largest six grain sizes. The noticeable difference is that stress levels are significantly reduced as ice descends onto smaller grain sizes. The smaller grain sizes play less of a role in this case. Transitional obstacle size Weertman (1964) noted that only one of the two terms in either of equations (3.17) would be significant. This appraisal is based on the regelation velocity being proportional to an inverse length scale while the creep closure velocity being proportional to a length scale. For a very small length scale, regelation would be the dominant mechanism of closure. For a large length scale, creep closure should dominate the system. Thus, intersection of two length scales would produce a limiting length scale. The actual effective area, volume, and length scale will depend on the distribution, arrangement, and shapes of grain sizes at the ice-bed interface. 64 Chapter 3. Sheet closure mechanisms 10~3 10"2 10"1 Grain radius (m) Figure 3.10: Logarithmic grain size distributions for the 5, 10, 25, and 50 grain size cases. The total area of the bed is 100 m 2 with 9% of the bed unoccupied by sediment. The theory presented here does not strongly support or oppose Weertman's idea but offers an alternative view of the processes. When water depth is shallow and the effective stress is low, larger grains tend to support more of the available stress (Figs. 3.11 and 3.12). These grains lessen the available stress that the other grains could use to produce faster velocities. When water depth is shallow and effective stress is high, intermediate-sized grains carry nearly as much stress as larger grains. When the driving stress is low, perhaps one could say that the large grains sizes are the controlling obstacle size. 3.6 Summary In this chapter, a method of calculating the rate at which ice closes a subglacial sheet is in-troduced. The method is based upon steady-state formulations of two processes that occur subglacially: viscous creep and regelation. Creep is the deformation of ice from an applied stress, and regelation results from the change in the melting point from an applied stress. The rates of closure can be calculated numerically from equations (3.17) and (3.14). Parameteriza-tions are necessary for appropriate length, area, and volume scales (eqs. 3.18, 3.19, 3.20, and 3.22 ) that are based on a distribution of idealized grain sizes lying on the bed. Simulations show that instantaneous closure velocities range from about 1 x I O - 7 m s _ 1 to 1 x I O - 4 m s _ 1 but can be as high as 1 x 10 - 2 m s _ 1 . These example velocities do not account for the lessening of the 65 Chapter 3. Sheet closure mechanisms 1 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9 oTot/oj aToVo' aTot/aj e e e Figure 3.11: Stress partitioning for ten grain sizes. Fraction of stress available for closure that rests on each of the indicated grain sizes. These are the six largest grain sizes with velocities corresponding to Figure 3.9b. Stresses for r7...io plot below the threshold. Black contours are spaced at intervals of 0.05. (a) Fractional available stress for r\ = 0.2 m. (b) Fractional available stress for rd\b~ ~ds~ (4.11b) where A1 is the sediment concentration in the overlying ice and 5 b represents the cross section of the sediment-fluid mixture. If there is no supply of sediment from the bed or overlying ice, and the volumetric concentration of sediment in the flow is zero, then these equations revert to equations (2.17) and (2.15). Equation (4.11a) is the sum of the instantaneous misfit balances of the ice and sediment. Equation (4.11b) is the balance of the fluid in the cross section. These equations are derived in appendix D (eqs. D.8 and D.9) with an alternative formulation for the pressure evolution equation based on the mass balance of water (eq. D.12). Similarly, relations equivalent to equations (2.16) for the subglacial water sheet are dHh 1 dt (1 - A1) (1 - nl) p1 + nj,pw + A! ( l - ni) p s fdHb\ vl> 1 + 1 W + dp° ~dt 1 close ( l OP3' 7 b i ¥ b f dHh d (Hhuh) dt ds A b p s + (1 - A b)pw ( l - A 1 ) (1-nj,) p[ - Hb (ps m + d\h dt + 1 + nppv (1 - A1) (1 - ni) ? m + + 1 + nsvpw (1 - n p) p s (4.12a) , d A b ' ds (4.12b) Not surprisingly, these equations are similar to the equations for the channel. Their derivation parallels the derivation in Appendix D. The momentum equation is much more difficult to modify. The reason for this difficulty is the movement of sediment. Generally, sediment transport is correlated to the water flow 74 Chapter 4. Subglacial sediment transport using a semiempirical relationship. As a result, the sediment momentum balance becomes a semiempirical black box. There are a few manageable modifications, however. The simplest modification is to substitute the fluid density, p f = A V + ( l - A b ) p w , (4.13) in place of the water density in equations (2.27). Another modification, given equations (4.11b) and (4.12b), is the scalar multiplication of production and supply terms with the velocity uh. This modification creates additional impulse terms as material moves off the hydraulic bound-aries. The changes to equations (2.27) yield duh ~dt = — u duh ,dub ds 1 dph _ pf ds ® ds hdub p f ds ^ ds 1 + ( l _ A i ) ( i _ n i ) p i dzh PWT0 m + + 1 + I + ' + n > w + A'(l - n > s (1 - nj,) p1 + n > w + 2 ( 1 - ni) (ps - p1) A1 1 + 1 W ( l - A ' ) ( l - < ) p ' rh + (4.16a) 1 + 1 W n i p m (4.16b) ( 1 - A 1 ) (l-nl)p' where equation (4.16a) holds for fh > 0 and equation (4.16b) holds for rh < 0. Because I assume that the ice is static, there is no convective velocity, which is equivalent to the assumptions used in equation (2.18). If the melt rate is positive, then the supply term has a simple form (see eqs. 4.34) and equation (4.16a) simplifies to 1 dZl 8t ( l - A ' X l - n i y This equation has the satisfying effect of resembling equation (2.18). (4.17) 4.2 Subglacial sediment characteristics A necessary assumption is that the subglacial bed exposed to the overlying ice and water is unconsolidated sediment. While bedrock is undoubtedly exposed at the bed of many glaciers, in the ablation zones the bed is often sediment covered (e.g., Clark k. Walder, 1994; Lawson, 1981). In general, it is in the ablation zones where supercooling has been observed. Typically, subglacial sediments are tills. Sediment size distributions for tills formed under 76 Chapter 4. Subglacial sediment transport shear follow a fractal distribution, N S = N o { § ^ ) m ' (4-18) such that the dimension, or fractal index m in equation (4.18), is approximately 2.9 (Hooke k Iverson, 1995). Equation (4.18) relates the number of particles ./Vs, of a characteristic diameter D, to a known number NQ, for a reference size DQ. For a fractal dimension of exactly 3, the volume occupied by each of the grain sizes will be the same. For m < 3, the volume of material will preferentially be in the larger grain sizes. For m > 3, the converse is true and more of the total volume of the till will be in the smaller grain sizes. Because the sheet closure mechanism relies on the index being 2 for the surface sediments (Fig. 3.10), I assign the volumetric distribution consistent with this and choose m = 3 in equation (4.18). Khatwa et al. (1999) analyzed sediments in various diamictons from different glacial en-vironments. Their results show fractal distributions ranging from 2.55 to 2.97. Supraglacial deposits in this study had a fractal dimension near 2.83, while the subglacial deposits had a fractal dimension near 2.92. Khatwa et al. questioned whether the fractal nature of glacial sediments is really characteristic of the erosion and sedimentation process. In particular, all of the diamictons they analyzed had a reasonable fractal dimension, but not all of the diamictons showed a shear fabric. Benn k Gemmell (2002) demonstrated, through synthetic distributions of sediments, that the use of a fractal distribution is not an especially discerning analysis for tills. They showed superimposed Gaussian distributions with different mean values may resemble fractal distribu-tions with m ~ 3. Benn k Gemmell went on to suggest that tills with a fractal distribution of about three can be formed not only from shearing as Hooke k Iverson (1995) proposed but also from mixing of tills. A practical problem of a fractal distribution of clast sizes lies in the larger grain sizes. Most studies with field data (e.g., Hooke k Iverson, 1995; Fischer k Hubbard, 1999; Khatwa et al, 1999; Benn k Gemmell, 2002) truncate the fractal distribution at either 1 or 10 mm. This truncation makes it difficult to estimate the size fraction larger than fine pebbles. Estimating the larger fraction from proglacial deposits is difficult because the fines may be winnowed. Analyses of till recovered from under Ice Stream B suggest that the fraction of sediment larger than 4 mm is roughly 3% (Tulaczyk et al., 1998). A drilling program at Black Rapids Glacier produced several cores of variable quality. Pebbles and cobbles were retrieved from beneath the 77 Chapter 4. Subglacial sediment transport glacier with over 20% of the recovered sediments being larger than 4 mm (Truffer et al., 1999). Because of core quality problems, the Black Rapids data may provide only an upper limit to the percent of coarse material. While a fractal distribution may not be very indicative of the genesis of a till, it does provide a starting point for approximating the texture of subglacial sediments. Because of the generality of the fractal distribution, it serves as an adequate initial condition on the parent till below the ice. For grain sizes larger than four millimeters, the fractal distribution may not hold because these grains usually are not investigated. The work of Benn & Gemmell (2002) circumvents this inadequacy of the published data, allowing for the fractal distribution to be used. By using a fractal distribution, I assign an equal volume of solids in subglacial till to all clast sizes. Lawson (1979) performed sedimentological analyses of the accreted basal layers at Matanuska Glacier, Alaska. He found that the accreted ice had an average sediment content of 25% by volume that ranged in size from clay to pebble gravel with large fractions being medium silt to sand. Grain size distributions of accreted sediment commonly have multiple peaks in sand to silt size particles, 0.008 < D < 2.0 mm (Lawson, 1979, 1993). Often there is a spike in coarse silt, D ~ 0.05 mm, accounting for approximately 20% of the sediment in the basal layers. However, sediment distributions may be dominated by sand or even be up to 50% gravel and larger grains. Lawson also identified a melt out till directly below the accreted ice at the margin that has the same source as the sediment entrained in the ice. He described the melt out till as a "pebbly sandy silt" (Lawson, 1979, p. 39). On the basis of these observations, it seems plausible that the smaller grain sizes may be volumetrically more important. There are two possible source options for sediment accreted in the ice. Either the accreted ice reflects the exact composition of the subglacial layer; or the accreted ice reflects a preferential sorting of sediment. In the case of a preexisting melt-out till, the accreted sediment may match the source sediment exactly. A melt-out till is comprised of sediment that had already been accreted. In most other cases, the accreted sediment will not reflect the exact texture of the parent material. Subaerial fluvial systems readily sort sediment. Therefore, one expects flowing water at the base of a glacier to sort sediment as well. Lawson et al. (1998) also described individual processes such as the filtering of sediment by frazil ice floes near the terminus of Matanuska Glacier. This filtering mechanism would tend to sort sediment as it is frozen into the accreted ice. As a consequence, it is unlikely that the accreted sediment fully represents the 78 Chapter 4. Subglacial sediment transport source. There are no direct observations of sediments in overdeepenings. As a result, I have assumed that the subglacial material is a deformation till. Where an active subglacial hydraulic system exists, the sediments could possibly be sorted such that they resemble proglacial outwash. Law-son (1981) described a section of fluvially sorted material in front of Matanuska Glacier; however, the location of this section relative to the terminus of the glacier is unclear. In addition, data on subglacial till from Storglaciaren, Sweden, suggests that this overdeepening is occupied by deformation till (Hooke et al, 1997). 4.2.1 Mobile sediments Knowledge of the grain size distribution is necessary to state what is and is not mobile sub-glacially. From the foregoing discussion, two issues are critical: (1) grain size distributions in overdeepenings are not well characterized, and (2) the smaller grain sizes appear to dominate subglacial sediments. The result of the first issue is that the model presented here cannot ac-curately reproduce any data including the grain size distributions presented by Lawson (1979). This shortcoming is beyond control; and, if necessary, a community effort will be required to more accurately characterize subglacial beds in situ. The consequence of the second issue is that construction of a model can focus on the smaller grain sizes. Focusing on smaller grain sizes is a benefit because the number of processes occur-ring at the bed are limited. Armoring of a subglacial streambed, for example, does not merit consideration. Furthermore if smaller grain sizes are readily available at the bed in large quan-tity, this availability suggests that the source of much of the sediment in the flowing water is local. The amount of sediment in the flow is not limited by either source or upstream conditions. Sediment transport conditions are then approximately near equilibrium. Equilibrium transport occurs when flowing water is eroding and depositing equivalent amounts of bed material simultaneously. Laboratory experiments or field measurements of equilibrium conditions require that the flow be steady, the sediment supply not limit the sediment transport, and the transport conditions be statistically steady (Bagnold, 1966). Once these conditions have been met, relations are constructed to predict sediment transport based on flow and sediment characteristics. The equilibrium assumption pertains to flow conditions. For long time periods, flow conditions will not be steady, and it must be true that entrained sediment is preferentially 79 Chapter 4. Subglacial sediment transport eroded from or deposited on the bed. For example, daily cycles of water discharge do not meet equilibrium conditions. On a second-to-second time period, transport will be near equilibrium. Subtle changes in the equilibrium' load will lead to the daily variation. The assumption im-plicit in this statement is that the sediment concentration always responds faster than the flow conditions. For large hydraulic cross-sections, this assumption may not be valid. However, for relatively shallow subglacial water sheets, this assumption is acceptable. If sediment responds slower than the flow conditions or the amount of sediment available to be transported varies spatially, a correction to an equilibrium transport rule via an adaptation length (e.g., Armanini & di Silvio, 1988; Wu, 2004; Yalin, 1972) or an equivalent spatial lag in sediment transport (e.g., Einstein, 1968; Phillips & Sutherland, 1989) is possible (eq. 4.4). However, adaptation lengths range from a f^ew centimeters to a few meters for coarse sand (Phillips & Sutherland, 1989). The spatial effect introduced by not applying an adaptation length is quite small especially when the longitudinal glacier sections (Figs. 1.3 and 1.5) are hundreds to thousand of meters long. Larger grain sizes, such as gravel and coarser material, are undoubtedly mobilized by sub-glacial water flow. R-channels, in particular, may mobilize sediments as the cross section en-larges. If larger grain sizes are mobilized, then processes such as sorting and armoring may become important. Sorting and armoring are observed in esker cross-sections (e.g., Shreve, 1985) that indicate these processes occur subglacially. To simulate armoring, the bed must be divided into vertical levels. Each of these levels.contains grain size information. This type of bed discretization is referred to informally as a Hirano bed (e.g., Parker & Sutherland, 1990). Sediment transport rules for a partially mobile sediment are required when multiple grain sizes may be mobile along the bed (Wilcock & McArdell, 1997, 1993; Wilcock, 1997). These processes are not beyond the scope of this thesis. However, the data available to constrain a model of subglacial water flow with multiple grain sizes using a Hirano bed is limited. The prudent course in this situation is to use a relation for total sediment transport before returning to finer de-tails such as armoring and grain distributions. The consequence is that larger grains, primarily gravels, will not be a primary concern here. In addition, the height of subglacial cross-sections along overdeepenings is expected to be small (Hooke, 1991; Rothlisberger & Lang, 1987) resulting in a size limitation to the material transported. In this instance, if the shear stress does cross the threshold necessary for entrain-80 Chapter 4. Subglacial sediment transport merit pf larger grains in the water layer, these grains will not be mobile. The grains could be rotated within the flowing water such that they touch the ice ceiling and are frozen in place along the melting perimeter. Another mechanism that may limit sediment available to be transported is the closure re-lationship presented in Chapter 3. If ice rests on larger grains, these will be locally fixed. If the shear stress exceeds the critical shear necessary for entrainment, the grain will not move. Finer sediments around the fixed grain will be winnowed. Ice could then preferentially move around the larger grain, possibly entraining it (Iverson, 2000; Iverson & Semmens, 1995). Large grains with the fines winnowed do appear in the accreted basal layers at Matanuska Glacier, Alaska (Lawson et al, 1998). Because of these complications the closure mechanism presented in Chapter 3 is decoupled in terms of grain size distribution from the model presented here. 4.3 E q u i l i b r i u m transport rules Sediment transport is commonly split into two categories: suspended load and bed load. This division is common within the fluvial literature and follows analyses of transport equations by Gomez & Church (1989) and Garcia & Parker (1991). Lawson (1993) and Loso et al. (2004) also make this distinction for sediments with a subglacial source. Suspended load consists of sediment grain sizes that are small enough to be suspended by the turbulence within the flow. Bed load is the sediment that travels in contact with or near the bed via rolling, sliding, or saltation. For most turbulent water flow, silt- and clay-sized particles (D < 0.06 mm), will be in the suspended load, while gravel and coarser sediments (D > 2 mm) will be in the bed load. The break between bed and suspended load often falls in sand-sized particles (0.06 < D < 2 mm); and the break is dependent on both the sediments and water flow (e.g., Church, 2006). The distinction is not necessarily sharp as there is an unbroken transition from bed to suspended load. Another sediment categorization is bed material load and wash load. Bed material load is that part of the sediment load that interacts with the bed. The sediment sizes in the bed material load can be found in the bed material. Wash load consists of the finer portion of the suspended load that rarely interact with the bed. These sediments can be held in suspension for quite a long distance and may not represent the local bed material. Komar (1989) noted that little attention has been given in the literature to wash load; and the amount of wash load 81 Chapter 4. Subglacial sediment transport in the flow is governed by its availability (Bagnold, 1966; Komar, 1989). While the sediment mass balance equations (eqs. 4.10) do not explicitly capture this feature of sediment transport, an approximate calculation of the wash load/bed material load fraction is possible based on the material transported. From the sediment size considerations, it is clear that two transport relations are necessary: one for bed load and one for suspended load. Ideally, all grain sizes below 4 mm are of interest. However, these gravel-sized sediments are still relatively large. Avoiding gravel, the maximum size would be 2 mm. Typically, sediment enters suspension around 0.5 mm (e.g., Garcia & Parker, 1991). The two relations that best account for suspended sediment are those proposed by van Rijn (1984b) and Smith & McLean (1977) (Garcia & Parker, 1991). Because there is a complementary bed load formula for van Rijn's suspended load formulation, I use both in tandem to estimate sediment transport. These relationships provide an appropriately scaled estimate of sediment transport that is preferable to simpler formulae that have less accuracy, van Rijn's (1984a; 1984b) formulae will handle sand-sized suspended load (0.1 < D < 0.5 mm) and sand-sized bed load (0.2 < D < 2.0 mm). These sand-bed formulae will determine the equilibrium sediment flux cfe in equation (4.4). At first glance, determining q~l from fluvial relationships seems tenuous. These flux rela-tionships were developed for rivers and canals that are deep relative to subglacial water flows. Nevertheless, these relationships will most likely yield the best estimates for subglacial flows. Sand-bed flume experiments often have water depths of three to 15 centimeters (e.g., Abbott & Francis, 1977; Lee & Hsu, 1994) that are within the range of water depths expected along overdeepenings. Furthermore, field data suggest that sand is commonly mobilized subglacially. The bulk of accreted sediment at Matanuska Glacier is in the silt/sand regime (Lawson, 1979; Lawson et al., 1998), and glaciers with larger subglacial passageways, Skei5ararj6kull for exam-ple, often discharge a significant amount of sand-sized sediment (e.g., Russell et al., 2006). In what follows, I briefly review van Rijn's formulae. For clarity, only an overview appears here. The essential parts of the formulation can be found in Appendix E. Additional information can be found in the original papers as well as the review of suspended load transport algorithms by Garcia & Parker (1991). 82 Chapter 4. Subglacial sediment transport 4.3.1 Sand formulae van Rijn published a series of papers that explored sand movement as bed load (van Rijn, 1984a), suspended load (van Rijn, 1984b), and in creating bedforms and modifying alluvial roughness (van Rijn, 1984c). The papers build a framework for sand motion in fluvial environments. The two later papers build on the bed load paper. As a result, I briefly review the bed load before examining the suspended load. Bed load van Rijn developed a physical model of saltation at the bed based on first principles. He then compared model simulations against experimental data of saltation along the bed. His model reproduced observed jump lengths and jump heights of sand grains reasonably well. The jump height varies from a few millimeters to about a centimeter. The jump length varies from a few millimeters to a few centimeters (Abbott & Francis, 1977; Lee & Hsu, 1994; van Rijn, 1984a). From this physical model, van Rijn developed a bed load transport relationship. The relationship is valid for grain sizes in the range (0.2 < D < 2.0 mm). The bed load transport is qe-.b = ubSb\b.b, (4.19) where ub is the mean sediment velocity in the bed load layer, 6b is the thickness of the bed load layer, and \ b is the volumetric sediment concentration. The bed load velocity is ub = 1.5T0f5 | - lj °/7Ta5A>o, (4-21) where £>* is a nondimensional particle parameter. It is a ratio of the gravity forces on a particle in fluid relative to the viscous forces and can be derived from the grain Reynolds and Froude numbers. From the discussion above, 5b should be approximately one centimeter for sand. 83 Chapter 4. Subglacial sediment transport Finally, the volumetric sediment concentration in the bed load layer is A b : 6 = 0.18A b : 0 -£, (4.22) where A b 0 is a reference concentration for the bed, A b 0 = n s a ( l - n p . (4.23) van Rijn considered the entire bed mobile and gave A b 0 as 0.65. Because the closure relationship presented in Chapter 3 requires larger, immobile clasts, I introduce the mobile fraction of the bed nfj. This term states what areal fraction of the bed material is mobile. van Rijn (1984a) supported equations (4.20), (4.21), and (4.22) with examples from exper-imental data. Substituting equations (4.20), (4.21), and (4.22) yields the final form of the bed load flux, qe:b = 0 . 0 8 1 A b : 0 ^ | - l)g} (4.24) van Rijn substituted A b . 0 = 0.65 directly into equation (4.24), using equation (4.23) instead gives g e : b = 0 . 0 5 3 < | ^ { ( ^ - l ) 5 } (4.25) as the final form of the bed load flux; If TT.^=1, as would be the case for a fully mobile bed, then the coefficient is simply 0.053, which is what appears in van Rijn (1984a). Suspended load The formula for the suspended load is qe:s = u b c i b A b s , (4-26) where ub is the velocity of the subglacial water flow, db is a representative water depth, and A b is the volumetric suspended sediment concentration. The depth of flow in equation (4.26) is db=t ^•Rb semicircular channel, 4 (4.27) Hb sheet, where the depth of a rectangle with an equivalent cross-sectional area approximates the depth of a semicircular channel, van Rijn (1984b) noted that his relationship is valid for grain sizes in the range (0.1 < D < 0.5 mm). Grain sizes smaller than 0.1 mm, silt and finer particles, 84 Chapter 4. Subglacial sediment transport undoubtedly enter suspension, but van Rijn's algorithm does not account for these grain sizes because he used data for predominantly clean sands. The volumetric concentration is, A e b : s = ^ A b : 0 (4.28) where Fs corrects for the distribution of suspended sediment not being vertically uniform. The reference suspended sediment concentration A b 0 is defined either at half the bedform height or at the Nikuradse equivalent roughness height. For a flat bed, the bed load layer height 6h is the suspended sediment reference height. Together equations (4.26) and (4.28) are the equation van Rijn (1984b) wrote for the suspended load flux. (4.29) where 8S is the suspension reference level and Z' is a nondimensional ratio of the momentum from a falling particle to the momentum available to keep that particle in suspension. The reference concentration is A b 0 = 0 . 0 1 5 ^ ^ (4.30) where D50, D*, and T were defined for bed load. Al l equations presented for both bed load and suspended load are subject to the supplementary relationships in Appendix E. 4.3.2 Bed shear and sediment transport The shear stress imparted on the fraction of mobile sediments at the bed determines whether those sediments are entrained in the flow. The shear stress available to sediments is commonly broken up into a grain stress necessary to move sediment, and a form stress (e.g., Yang, 1996, p. 70). Because fluid velocity and density are the same for both stresses, the friction coefficients are additive such that fd = f'd + fd, (4-31) where f'd and f'd are the grain and form roughness coefficients, respectively. The boundary shear stress from equations (2.24) and (2.25) becomes TO = 4 +fa + rd). (4.32) 85 Chapter 4. Subglacial sediment transport where TQ and TQ are the grain and form shear along the bed. Partitioning the shear stress along the bed is necessary to compute the sediment load. This partitioning is important for estimating the grain Chezy coefficient in equation (E.3). The grain Chezy coefficient is related to the grain Darcy-Weisbach friction coefficient via As a result, estimates of f'd are necessary. In alluvial rivers, form resistance generally accounts for bedforms such as ripples, bars, and dunes. In rough streams, large clast sizes account for a significant portion of the form shear, with relatively little grain shear available to move smaller grain sizes (Church et al, 1998). In a subglacial water system, not only are fluvial bedforms possible, but the assumed fractal nature of the grain distribution suggests that large grains bear much of the total stress available for the bed. As a result, a reasonable assumption is that the grain stress is approximately 15-20% of the total stress available to the bed. It is possible to have slightly lower or higher values depending on conditions. These percentages are only suggestive and need to be refined via field studies. The fraction of the bed that is mobile is limited to the grain sizes less than about 2 mm in accordance with van Rijn (1984a). Furthermore, the closure relationship proposed for the bed has 18 grain sizes (Fig. F.4). Because the bed is assumed to be fractal with m — 3 (eq. 3.23), each of the grain sizes occupies an equivalent volume of the bed. As a result, half of the bed material is assumed to be mobile. By making this assumption, I implicitly assume that the suspended load grain sizes below $ = 4 behave as suspended sand (van Rijn, 1984b). In general, this is not true, but it does allow an expeditious procedure for simulating smaller grain sizes. Small grain sizes will be mobile, and treating them as suspended sand provides a sensible mobility criterion. Depending on how the bed elevation evolves, the grain size distribution could also evolve. There is little data against which to check such an evolution. In consequence, I make no allowance for the erosion or deposition of fine-grained material relative to the coarse-grained material. Thus, dynamically armoring the bed via sand mining—removal of fines from space between coarse clasts—is not possible. Measured friction coefficients are lower for mobile beds. Vanoni & Nomicos (n.d.) performed flume experiments on both mobile and immobile sand beds. Their intention was to examine the influence of roughness elements on the flow. They found that, relative to equivalent clearwater (4.33) 86 Chapter 4. Subglacial sediment transport flows, bed mobility lowered the friction coefficient by 5 to 28%. Using similar experiments, Gao & Abrahams (2004) found that the shear along the bed was lowered by up to 10%. While modulation of the friction factor with sediment concentration may occur, it is difficult to examine variations when the effective magnitudes of the friction coefficient are not known for subglacial water flow. In consequence, no modification is made to the friction coefficient to account for sediment transport. 4.4 Conceptual model of sediment and ice accretion Sediment accretion is not well understood. Using a sediment concentration of subglacial water of 5 g L - 1 (Lawson, 1993), then the equivalent volume of accreted sediment should be approxi-mately 2%. Because sediment concentrations in the ice are approximately 25% by volume but can range from practically 0% to over 95%, Lawson et al. (1998) invoke a filtering mechanism to increase the concentration of sediment accreted. Figure 4.1c illustrates subglacial sediment transport and this filtering mechanism. Using the ideas of Alley et al. (2003b) and Hooke (1991) as a template, Figure 4.1 illus-trates most of the ideas captured by the balance laws for sediment transport and equilibrium transport rules. In Figure 4.1a, sediment is eroded from the bed and transported in the fluid. As the fluid begins to ascend the overdeepening, water pressure increases (Hooke & Pohjola, 1994). The increased pressure forces water to distribute laterally along the ice-bed interface. Alternatively, water may escape through an englacial aquifer. The reduced flow rates along the bed imply a lessened velocity. As a result, sediment transport is diminished (Fig. 4.1b). The fluid continues to ascend the adverse slope, and the water temperature drops below the bulk freezing temperature and supercools. Frazil ice forms in the flow and can flocculate. The floes float to the top of the water layer. Figure 4.1c shows schematically how these floes can trap sediment. The flocculated layer is exceptionally porous and allows water flow to pump sediment through those pores. In Figure 4.Id 1 2, the freezing front descends, and sediment is trapped in the ice. Obviously, all of these processes are continuous, or nearly so, but the four subfigures 1 2 I n Figure 4.Id, I have extended the fractures to symbolize the residual porosity originally introduced in equation (2.11). Obviously, the fractures wil l not extend to the accreted ice without incorporating a rule for fracture generation in a moving glacier. 87 Chapter 4. Subglacial sediment transport Figure 4.1: Illustration of sediment transport and accretion. Lines through the overlying glacier represent porosity in the form of fractures (Chap. 5). Filled circles represent sediment grains, and thin white rectangles represent frazil ice. (a) Sediment is eroded from the bed and transported in the flow, (b) Water begins to flow up the adverse slope and loses sediment, (c) Frazil ice forms in the flow and flocculates along the top of the water layer. Sediment is filtered through a porous lattice of frazil crystals, (d) Ice seals the porous lattice, and sediment is accreted. illustrate snapshots. Furthermore, Figures 4.1c and d would be continuous in a vertical sense, with the freezing front descending as frazil floes are added to the base of the glacier. Processes occurring in the flocculated layer are responsible for sediment accretion, and to model sediment accretion properly a frazil ice component is necessary. However, the nucle-ation, enlargement, flocculation, and eventual solidification of a frazil ice layer is a complex and demanding topic (e.g., Ashton, 1982; Daly, 1984; Martin, 1981; Wettlaufer, 2001). Because incorporation of a frazil ice component into the model is beyond the scope of this thesis, \&1 and & in equations (4.10) become adjustable parameters when the melt rate is negative and ice is Chapter 4. Subglacial sediment transport accreting. The supply to or from the ice is for m > 0 channel (4.34a) C 3 A b m for m < 0 A y - for rh > 0 sheet (4.34b) C 3 A b m for m < 0 where C 3 is an arbitrary constant. An assumption inherent in this mass supply term is that there is no sediment in the water in the pore space of the ice. Ideally, C 3 would be related to the degree of supercooling, the flow Reynolds number (or alternatively, the turbulence intensity which has been linked to frazil ice formation (Daly, 1984)), as well as the concentration of sediment in the water. These heuristic suggestions would make C 3 nonlinear. However, because there is no data to support these suggestions, implementation would be tenuous. As a result, I retain C 3 as an assigned parameter. 4.5 Verif icat ion of the results The amount of sediment accretion can be tuned via the free parameter C 3 . Tuning the sediment accretion can yield results that may not be related via a physical mechanism to field observations of sediment accretion. Furthermore, the equilibrium transport relations of van Rijn (1984a,b) have prediction limitations. When compared to data, van Rijn's bed load algorithm successfully predicts an amount that is 0.5 to 2 times the actual sediment transported for 77% of 600 data sets (van Rijn, 1984a). Similarly, the suspended load algorithm predicts a value that is between 0.5 and 2 times the suspended sediment concentration for 76% of 800 data sets (van Rijn, 1984b). Amounts of subglacial sediment accreted can vary substantially because of the variability associated with the equilibrium transport relations and balance laws presented here. A fully independent assessment of subglacial sediment transported and accreted may not be possible; however, some checks are available. One of these checks is a comparison of sediment eroded, transported, and deposited to the Hjulstrom diagram (Fig. 4.2). The Hjulstrom dia-gram illustrates velocities necessary to erode, transport, and deposit grains of a given diameter. While not a stringent test, the Hjulstrom diagram provides a viable way of checking the sedi-ment concentration versus the flow conditions. Because the wash load can be calculated from 89 Chapter 4. Subglacial sediment transport the balance laws, the Hjulstrom diagram may be most stringent for smaller grain sizes in the transportation regime. 10 D (mm) Figure 4.2: Hjulstrom diagram. (After Graf, 1971, p. 88). Hjulstrom (1935) derived these grain size-velocity relationships empirically for beds with uniform grain size. The double line separating the transportation and erosion regimes shows that this boundary is not necessarily sharp. Dashed vertical lines indicate the range of validity for equilibrium bed load and suspended load relationships (van Rijn, 1984a,b) 4.6 Summary This chapter presents a model of subglacial sediment erosion and deposition. Using the mass balance for a continuum, rates of sediment transport per unit flow width are derived (eqs. 4.10) as well as the elevation evolution for the subglacial bed (eq. 4.7). If the subglacial fluid density changes significantly because of sediment transport (eq. 4.13), then the fluid flow balance laws take a different form for both the channel and sheet (eqs. 4.11, 4.12, 4.14, 4.15). Suspended load and bed load are common divisions of types of sediment transport. For each type of trans-port there is one set of equations (4.10), resulting in two equations for sediment transport in subglacial channels and two equations for sediment transport in subglacial sheets. Because sed-iment textures at the base of the glacier are difficult to quantify and fines seem to dominate, only equilibrium transport relationships for sand are employed (eqs. 4.24 and 4.26). These rela-tionships are subject to a number of auxiliary equations (Appendix E). Furthermore, sediment 90 Chapter 4. Subglacial sediment transport transport in the sheet is not explicitly coupled to the closure relationship proposed in Chapter 3. In addition, sediment accretion to the base of the glacier is parameterized based on the amount of sediment in the water (eqs. 4.34). Finally, verifying the sediment load within the subglacial water system is not straightforward. One simple check is to compare results to the Hjulstrom curve. 91 Chapter 5 Englacial aquifer 5.1 Introduct ion Subglacial passageways will not always be available at the ice-bed interface. In these instances, water will move through either an englacial aquifer in the overlying ice or a subglacial aquifer below the sheet. An englacial aquifer poses unique problems in that the porosity and perme-ability may evolve as water freezes or melts depending on local conditions. Flow then responds not only to the change in head driving the system but also to the change in permeability with time. Glaciers with an active hydrology have several different types of aquifers based on the form of the porosity. In the glacier surface snow, there is a multigranular porosity where meltwater percolates through the snow (e.g., Illangasekare et al, 1990; Pfeffer et al., 1990). Veins at ice-grain boundaries are another type of porosity (Nye & Frank, 1973). Veins are found throughout most glaciers and ice sheets but contribute relatively little to the overall hydrology because they have extremely small apertures. Fractures comprise a third type of porosity. Recently, Fountain et al. (2005a) suggested that fractures not only dominate englacial flow but may also be the dominant hydraulic pathway through glaciers. Fractures form in glaciers for any number of reasons (e.g., Ffambrey &; Dowdeswell, 1997; Ffambrey et al, 2005; Lawson, 1996), but they are commonly generated as ice flows over steep bedrock. For example, crevasses develop at Storglaciaren as the ice descends bedrock riegels (Fig. 1.5) (e.g., Fountain et al., 2005a; Stenborg, 1973). These crevasses do not completely close, creating an ubiquitous fracture network. These fractures are expected to govern water flow within an englacial system. There are only a few existing models of englacial water flow. Nye (1976) discussed flow in veins. He provided a physical and thermodynamic formulation of water flow through veins. Hutter (1982) developed a model of ice flow and water transport within the cold and temperate portions of an ice sheet. His formulation of the temperate ice touched on the mass and mo-92 Chapter 5. Englacial aquifer mentum balances of the water and ice as well as the jump conditions at the phase boundary. Fowler (1984) added to Flutter's work by including a formal heat balance for the temperate ice. Fowler's work is particularly interesting because of this feature. Flowers & Clarke (2002a) created a vertically integrated model of englacial water flow through fractures. I expand upon Fowler's work with the intention of simulating flow in fractures. Water temperature, pore space, and permeability evolution are of particular interest. Because there may be supercooled water at the base of a glacier, this water may be able to move upward, freezing to block pore space. Alternatively, warmer englacial water may move downward, curbing subglacial supercooling. These two scenarios illustrate how the englacial hydrology may affect the subglacial hydrology or vice versa. Much of this chapter focuses on the derivation of appropriate relationships to expand previous work to be able to investigate hydrologic scenarios pertinent to glaciohydraulic supercooling. 5.2 Balance laws Chapter 2 featured steady-state and time-dependent variations of three balance laws: mass, momentum, and thermal energy. For an englacial system, the balance laws closely parallel a subglacial water system. However, the balance of momentum is the most important because it leads to Darcy's law. 5.2.1 Momentum balance Darcy's law governs the mass and momentum balance for flow in a porous medium, Kw q w = - ^ r V c T , (5-ia) K^J^PZI ( 5 . L B ) h where q w is the specific discharge of water vector, KJJ is the hydraulic conductivity, k is the permeability, and p is the static viscosity of water. The hydraulic potential assumes its usual form, c/>w = pwgz + pw (eq. 2.23). A melt rate term must be included in the momentum equation prior to applying the simplify-ing assumptions inherent in Darcy's law. Where the pore space is changing with time, equation (5.1a) may not provide a proper balance. The goal of expanding the momentum equation is to 93 Chapter 5. Englacial aquifer alter Darcy's law for a change in pore space via melting or freezing. Alternatively, inclusion of a melt rate may not warrant a change in Darcy's law. Either conclusion is useful. The material derivative of linear momentum per unit volume is given by, p W i ? + p W u W ' = ~{mv ~ ^ ~ V r ~ ^ T' (5'2) where u w is the velocity vector of the water through a cross section S, mv is a volumetric melt is a pore closure rate via viscous deformation of ice, ro is the shear along the hydraulic perimeter P w . The right hand side represents an impulse from melt —mvuw, an impulse from moving ice dvuw, a change in potential V(/>w, and a loss in momentum from the water flowing within the matrix PWTO/S. The hydraulic perimeter PW is related to the surface area of the pore space. The melt rate mv is positive when matrix ice melts or negative when pore water freezes. Because flow is predominantly slow and steady and du™/dt is usually assumed small, the left hand side of (5.2) can be neglected. Recognizing that Darcy's law is a balance between hydraulic potential and dissipation from shear, the shear term in equation (5.2) can be approximated using equation (5.1a), - ( m v - d v ) n w - V r - ^ d w = 0. (5.3) H Rewriting this equation for the specific discharge gives owa L J Equation (5.4) is further simplified substituting the pore water velocity with q w / nL . A simplified expression for specific discharge with the momentum source terms is q w = ? — T T T V 0 W (5-5) w f 1 + * 5 ( m , 7 4 , ) l I P w 5 < J v2"p If the melt rate is zero, equation (5.5) reverts to Darcy's law (eq. 5.1a). In practice, the impulse term should be negligible such that the product K^j {mv — dy) /{p^grip) The equation of state for ice is similar and varies according to p ^ P o e x p j y ^ - P o ) } . (5T2a) where p0 is the reference ice density at the reference ice pressure pQ, and 71 is the bulk com-pressibility of ice. These reference values are the density of ice at atmospheric pressure and atmospheric pressure, respectively. The pressure in this case is the stress exerted on the ice frame. Using a pressure in the equation of state necessitates that the stress state be isotropic. If the ice were deforming viscously, strictly speaking, the stress tensor would not be isotropic, but the isotropic case would provide a reasonable approximation. In addition, the compressibility of the porous ice matrix is not that of pure ice. This compressibility should be lower than that of pure ice but should converge to the value of ice as the porosity approaches zero. Expanding equations (5.8) and (5.9c) and solving them for (mv — dv) gives two equivalent equations. Setting these equations equal and substituting the"derivatives of the equations of state for the change in water and ice density gives an evolution equation for the pressures, + ( 1 - H * H % = ^ V • ( p W q W ) - <'W - & lH> ( 5 J 3 ) where 7g is the effective compressibility of the ice matrix. While the values of the pressure in the ice and water may be different, if there are no changes in the aquifer, then locally 96 Chapter 5. Englacial aquifer dpw/dt = dp1 /dt; and equation (5.13) can be solved for the water pressure, -^VWr-Zd-^ t^ ^^ '-^ }^- <514) Bear (1988, p. 204) derived a similar form of this equation without source terms. Equation (5.14) is also reminiscent of the equation for water flow in a deforming medium (e.g., Domenico & Schwartz, 1997, p. 463), but there are two large differences. In water flow in a deforming medium, the matrix has a velocity; and the material derivative is not equal to a source term. If a flowing ice matrix is also incorporated into the present formulation, then the deforming matrix and the pore space change via melt/freezing would both be necessary and the equations would be slightly different. With both the equation for the evolution of the porosity and the evolution of pressure in hand, one final assumption can be made: that the local change in density is negligible. This assumption simplifies equations (5.9c) and (5.14) to £4 = m v - d dt pl dt nj,7« + (1 - njjTj ( ' " 'V" f ) dt The mass balance for the water is not used for the evolution of the porosity for the simple reason that the ice matrix defines the porosity. Inclusion of the divergence of the flux helps to determine the pressure and would, essentially, be double counted if used in the porosity evolution equation. While the change in the porosity is expected to be small, one weakness of this formulation is that the mass balance of the ice should be fully coupled to ice flow, especially for the viscous deformation term. This coupling is well beyond the scope of the present work, but equation (5.15b) is a reasonable approximation because the hydraulic system responds much more quickly than the ice flows. 5.2.3 Energy balance The thermal energy balance can be broken into parts in a fashion similar to the mass balance equation. In this case, the mixture of water and ice is the sum of the constituent balances. With the porosities expected to be low, the mixture temperature T should be close to the ice temper-ature Ts. Below the water table, water content at Falljdkull, Iceland, for example, is commonly 2-4% (Murray et al., 2000), indicating that the bulk temperature is controlled by the temperate 97 Chapter 5. Englacial aquifer ice. Because neither the mixture temperature nor the ice temperature will vary significantly near the melting point, the total heat of the water must be determined independently of the mixture heat. The water sheet below the ice and the boundary condition at the surface of the glacier will affect the temperature of the water. The temperature evolution of the ice fraction may be important, but the suspicion is that the ice temperature change is not important. Sensible heats near the freezing point are typically orders of magnitude smaller than latent heat terms. In order to simplify the problem, I assume that the ice is isothermal and at the melting point, as is commonly observed at many temperate glaciers. I thus ignore the thermal energy balance of the ice entirely. Other formulations for polythermal glaciers exist (e.g., Fowler, 1984; Hutter, 1982) but are beyond the scope of the present work. Using elements of the derivations of Fowler (1984) and Spring & Hutter (1982), the total internal energy of the water Ew in the pore space is = f pwnlewdV, (5.16) Jv E iv where the internal energy density can be approximated as ew = cpATmp+L. This approximation neglects the kinetic energy as well as any variations in the kinetic energy. For Darcian flow, ignoring the kinetic energy is justified. Moving equation (5.16) from a total internal energy to a local energy balance yields, pWCPNP^T =pWc>P^f + PWCP^ • VRW» (5-17A) w = - (mv - dv) (L + c;ATmp) + - t £ - ( g w ) 2 ; ( 5 1 7 b ) npJ\H where qw is the magnitude of the water flux vector. Both the momentum balance (eq. 5.1) and the mass balance (eq. 5.7) have been used to simplify these local balances. The only source term is the final term on the right hand side of equation (5.17b), which represents viscous dissipation. Because of Darcy's law, the viscous dissipation can also be written as K^(V4>w)2/(nppwg). Other source terms, such as diffusion of heat through the water via Fourier's law is omitted be-cause porosity is low. Rewriting equations (5.17) as an englacial temperature evolution equation yields, 9TW _ q w V T W mv-d g (q"\2 ( , 1 R ) This equation is qualitatively similar to equation (2.29). 98 Chapter 5. Englacial aquifer 5.3 Pore geometry Numerous attempts have been made to relate hydraulic conductivity, permeability, and porosity (e.g., Bear, 1988, chap. 5). Specifically, researchers have attempted to link physical models of different pore geometries to Darcy's law. The goal is to be able to say how the hydraulic conductivity varies with a change in the pore volume. A relationship between porosity and hydraulic conductivity is necessary to close the mass balance and momentum equations. Flowers & Clarke (2002a) chose the intersecting fracture model of Snow (1968). This model has several advantages. Snow's model is relatively simple in its treatment of both the geometry and the flow physics. In addition, because englacial hydrology is expected to be a fracture network, Snow's model is immediately relevant. The hydraulic conductivity of a network of intersecting fractures is where Fe is a number of fractures per unit length, and b is an average fracture width (Snow, 1968; Flowers & Clarke, 2002a). The porosity is the average fracture width multiplied by the number density: These relationships are strictly for laminar flow between two approximately parallel-sided frac-ture walls. Snow (1968) originally applied this relationship to horizontal packer tests of crys-talline rocks with low permeability and porosity. He suggested that the maximum fracture size be less than two millimeters. Furthermore, equation (5.19) assumes that all of the fractures are hydraulically conductive. For Storglaciaren, Sweden, approximately half the fractures are hydraulically conductive with remaining fractures storing water (Fountain et al., 2005b). In consequence, equation (5.19) represents an upper limit of englacial hydraulic conductivity. It is possible to reduce the hydraulic conductivity by multiplying equation (5.19) by the fraction of active fractures. Of the three quantities in equation (5.20) that describe the hydraulic properties of the ice, Fe is the most problematic. Recently, Hambrey conducted a study of fractures at Trapridge Glacier, Yukon Territory, Canada. He found that fracture densities ranged from 0.6-2.6 m - 2 (~0.8-1.6 m _ 1 ) (pers. comm., 2006). Fracture densities for the terminus of Matanuska Glacier most likely fall at the lower end of these values, and fracture densities for the overdeepenings of Storglaciaren P^gFeb* p 12 (5.19) < = Feb. (5.20) 99 Chapter 5. Englacial aquifer can be significantly less. Between the other two remaining parameters np and b, the porosity of the ice below the water table is better constrained and can range from approximately 0.01 to 0.10 (Pettersson et al, 2004, and references therein). Fracture widths are difficult to ascertain because hot water drilling often enlarges the fracture before it is measured (e.g., Fountain et al., 2005b). However, these are constrained to a maximum of 2 mm as noted above. Rather than evolve the number density of fractures, which would require additional con-ditions for fracture opening and closing, it is preferable to continue simply with the thermal evolution of crack width. From (5.20), the time evolution of the porosity becomes, d4 = Fef, (5.21) dt *dt' v ; where the additional bdFt/dt is intentionally absent from the right hand side. Without a change in the number density of fractures with time, this number density becomes a constant; and the change in porosity is directly linked to the change in fracture width. The evolution of the fracture width can be calculated from equations (5.15b) and (5.20). 5.3.1 Creep closure rate Nye (1953) explored the viscous creep of isotropic ice using Glen's power law for ice rheology (Glen, 1954). Nye used these relations to show how cylindrical holes, spherical holes, and conduits in ice close. The extension to an isotropically closing fracture in ice is not necessarily straightforward given that the divergence terms in a circular, spherical, or cylindrical coordinate system lend themselves immediately to solution via separation of variables (Nye, 1953). Assuming incompressibility and that tractions are negligible, the strain rate tensor for a crack becomes du/dx = -du/dz — 0. This relationship implies that the closure rate is independent of crack width and depends only on the pressure difference between the ice and the water. However, experience gained from conduit closure suggests that the hydraulic cross section does matter. For example, Ng (1998, chap. 3) derived a closure relationship for elliptical channels, but a complementary derivation is beyond the scope of this work. Rather, it will be more useful to parameterize both Ng and Nye's work with respect to fractures in a porous medium. Considering Nye's (1953) relationship for contraction of a tunnel in ice and rewriting it for the change in the diameter of that tunnel yields a one-dimensional closure rate, i = - » w - , . - ) ( ^ ) " . < - ) 100 Chapter 5. Englacial aquifer which is similar to the closure term in equation (2.17). Given that a circular aperture should close faster than a long fracture because of the circular convergence of the ice, this relationship would provide a reasonable upper bound for the closure of a fracture. Conversely, if the fractures are opening via deformation of the ice, then this relationship will overestimate the change. Arguably, a shape factor could be inserted into (5.22) to account for the change in cross section. The pressure of the ice in equation (5.22) is simply the ice overburden pressure, pl = pwg (z1 — z). From equations (5.21),and (5.22) the closure rate dv becomes, where p' corrects the units for equation (5.15a). Equation (5.23) arbitrarily takes closure to be positive when fracture apertures are decreasing (i.e., closing). 5.3.2 Volumetric melt rate Heat transfer in laminar flow regimes has been well-studied. Analytic solutions for both the flow distribution and the temperature distribution have been derived for Hagen-Poiseuille flow through tubes and laminar flow between flat plates. Pigford (1955, eq. 32) gave the heat transfer coefficient for nonisothermal conditions for a tube, where R is the radius of the tube, s is its length, and S is its cross-sectional area. If the transverse length of a fracture can be taken as equivalent to its height, then an equivalent relation for plane-parallel fractures is Equation (5.24) is the heat transfer coefficient for only one crack. Multiplying this equation by yields the heat transfer density. Substituting the heat transfer density into equation (2.33b) yields the volumetric melt rate, The temperature difference is simply that of the interstitial water and the melting point, with the melting point governed by the Clausius-Clapeyron equation (eq. 2.8). The melting point is dependent on the interfacial pressure between the water and the ice. (5.23) (5.24) (5.25) 101 Chapter 5. Englacial aquifer 5.4 Aquifer thickness evolution At the top of the englacial aquifer, water may leave or enter the glacier, depending on the availability of water and the saturation of the ice. Water can also enter the aquifer from the sheet or conduit at the base of the glacier. These two sources of water as well as the flow of water within the englacial aquifer may cause the total thickness of the aquifer to change height. Figure (5.1) shows the conceptual configuration of the englacial aquifer where the thickness of the aquifer is defined as Ze = ze — z°. This figure complements figure (2.1) in defining the different geometrical quantities. Glacier Outlet Sediment Figure 5.1: Schematic englacial aquifer. Front of the section presented in Figure 1.3. Dotted line indicates the height of the water table. From left to right: Zx is the glacier thickness. Z e is the englacial aquifer thickness, z1 is the height of the glacier surface. ze is the height of the englacial water table. z b is the height of the base of the glacier. zs is the height of the base of the water layer. Hr is the thickness of water above zT. If Hr > 0, then ze = zT. In this illustration, the front of the glacier has a seepage face. 5 .4 .1 Water table evolution The rate of change of englacial water table elevation ze, is a kinematic boundary problem. The formulation is well known (e.g., Fowler, 1997, pp. 79-80), 102 Chapter 5. Englacial aquifer where all values are valid on the boundary. Equation (5.26) holds for cases of the change in the water table elevation is below the surface of the glacier. When the englacial water table reaches or exceeds the surface of the glacier, there is no condition on equation (5.26) to either keep the water table at the surface or restore the water table to the surface. As a result, there will need to be a condition that restores the water table to the surface of the glacier when it is exceeded. A penalty function forces the englacial system to monotonically approach the surface of the glacier once the water table has exceeded an arbitrary elevation, where r r is an arbitrary time constant. Flowers & Clarke (2002a) used a similar formulation to describe the exchange of water between the en- and supraglacial water systems instead of a form of Darcy's law (eq. 5.28). Penalty functions are a mode of mitigating inconsistencies between a physical formulation and necessary simplifications. In this instance, the physical formulation allows the water table to exceed the glacier surface, which is not realistic. Other methods of solving the confined kinematic boundary problem exist, and equation (5.27) is simply numerically expedient. 5.4.2 Surface boundary The boundary at the glacier surface is phreatic with the fluid pressure set to atmospheric pressure, if there is no supraglacial flow at the surface of the glacier. If there is water flow at the surface then conditions are assumed for the flow of water into the ice. No attempt is made to model the capillary fringe above the phreatic surface, nor is any attempt made to capture the full physics of fluid flow in the snowpack, firn, or unsaturated zone in general. There are two options for the surface boundary: the water table is at the surface of the glacier or the water table is below the surface of the glacier. If the water table is at the surface of the glacier, then the pressure condition on the aquifer surface is simply atmospheric if there is no surface water at the surface. If there is water on the glacier surface, then the pressure condition is taken as the height of the surface water layer. These pressure conditions affect the surface formulation of Darcy's law, a" - K » ^ (5 28^ W - - ^ ^ , ( - ) 103 Chapter 5. Englacial aquifer where d^/dz is the hydraulic potential difference across the surface of the water table. This hydraulic potential does not have a pressure component when when there is not a supraglacial water system. In some cases, there may be a perched aquifer on the surface of the glacier. For example, meltwater ponds form at the surface when water flows over the surface of the glacier but cannot drain through moulins, fractures, or other drainage features. If the flux of water across the surface of a glacier is greater than the flux from the surface to the englacial system, then there may be a perched aquifer feeding the englacial aquifer. In terms of the Darcian model presented here, the flux of water to the englacial system may be approximated as, where q\ is the vertical water flux at the glacier surface and Hr is the thickness of the supraglacial water system. If there is no supraglacial system, there is no flux from the surface. Rather than model the unsaturated zone explicitly, I assume that water is instantaneously added to the water table with an amount governed by Darcy's law. Because the unsaturated zone has a pressure equal to atmospheric pressure, the pressure difference is approximated by the hydrostatic pres-sure at the base of the supraglacial water layer taken over the height of the unsaturated zone (upper condition in (5.29)). If there is no surface water layer, then this term decays to zero. If there is no unsaturated zone such that the englacial water height is equal to the surface elevation of the ice, ze = z\ then q\ takes the form of equation (5.28). In groundwater applications, the hydraulic conductivity of the unsaturated zone is gener-ally lower than that of the saturated medium. Because equation (5.19) relates flow through a parallel fracture, it is not immediately clear what part of the physics should be altered. To com-pensate, the unsaturated zone hydraulic conductivity will be the equivalent saturated hydraulic conductivity multiplied by an arbitrary retarding coefficient. If water emerges at the surface of the glacier, then the water may flow supraglacially. Typi-cally, water will form streams at the surface of the glacier in a way similar to subaerial streams. Because these streams are not of particular interest here, water flow on the glacier surface is for HT > 0, for HT = 0, (5.29) 104 Chapter 5. Englacial aquifer handled as a diffusive flux in the same fashion as Flowers & Clarke (2002a), a m 5L_=_V-(u r fr)+£, (5.30a) u r = - K ^ V [Zx + Hr) , (5.30b) where u r is the surface velocity, and K.rH is the surface hydraulic conductivity tensor. Surface energy balance terms and precipitation terms can be added to this balance for a more inclusive surface hydrology (Flowers & Clarke, 2002a). These terms are avoided here. Streams preferentially move water directly downhill on glaciers, melting out supraglacial stream beds. As a consequence, the hydraulic conductivity tensor in equation (5.30b) is rarely isotropic. However, because this relationship is used only to accommodate water flow, I assume that the conductivity is isotropic in the same manner as Flowers & Clarke (2002a). 5.4.3 Subglacial boundary The lower boundary is the water layer. Flux into the aquifer is governed by Darcy's law applied at the base of the glacier. This water flux not only provides water to the englacial aquifer, but drains water from the subglacial sheet. Obviously, water can flow from the englacial aquifer to the water layer or in the reverse direction, depending on the hydraulic potential gradient. There is no necessary flow boundary condition at the lower bound except the pressure condition, which is simply the water pressure in the sheet. The mass of water in the subglacial sheet will be governed by equation (5.31a) while the flux into the base of the aquifer is governed by its negative in equation (5.31b), (5.31a) dt pwg dz v ( 6 ' 3 1 b ) pwg dz The pressure difference governing the flux is simply the difference in pressure between the sheet and the aquifer. Similarly, temperature can be imported or exported depending on the water flux. Both the pressure and temperature are Dirichlet conditions in that they are specified at the boundary of the aquifer as the conditions in the sheet. These boundary conditions are not dependent on a transition from laminar to turbulent flow when the water leaves the aquifer and enters the sheet or the converse when water leaves the sheet and enters the aquifer. 105 Chapter 5. Englacial aquifer Melting or accretion of ice by the subglacial water will either expand or contract the aquifer at its base. This accretion is described by the melt rate equation for the subglacial sheet, dt ~ (l-ni)pi- ( 5 - 3 2 ) This equation is the same as equation (2.18) presented for the glacier thickness; however, here, the equation has a different function as a change in the height of the base of the glacier. The thickness of the englacial aquifer, dt ~ dt dt' is simply a derivative of the definition of the aquifer thickness. Either the ordinary or partial derivatives could be used in (5.33) because the convective velocities are zero for these terms because the glacier is assumed to be static. 5 . 4 . 4 Downstream and upstream boundaries At the downstream end of the glacier, ice and water meet. If the water table is high enough englacially, then the glacier may have a seepage face. The condition on the seepage face is that the pressure is atmospheric. As a result, the hydraulic potential is simply the height at the outlet, with no pressure component. Computationally, this outlet is a simple Dirichlet condition on the seepage face. The upstream boundary is also prescribed as a Dirichlet condition. This condition can vary with time, and as a result, requires a function to describe its evolution. An upstream condition is somewhat artificial because water would normally enter the glacier from the surface. The usual upstream condition would be a flux condition. For example, a zero flux at a glacier head wall, or an equivalent condition would be more physical. However, because I am only interested in sections of glaciers at this point, head will be the prescribed quantity. 5.5 Summary Throughout this chapter, I have focused on the hydrologic and thermal evolution of the englacial aquifer. Together the three balance laws—momentum (eqs. 5.1), mass (eqs. 5.15), and energy (eq. 5.18)—provide a mathematical representation through which water flux, aquifer porosity, water pressure, and water temperature can be simulated. These formulations are subject to 106 Chapter 5. Englacial aquifer auxiliary conditions for creep closure of the pore space (eq. 5.23) and heat transfer between the ice and water (eq. 5.25). Boundary conditions at the surface of the glacier are based on Darcy's law (eq. 5.28) and on the conservation of water along a kinematic boundary (eq. 5.26). If the water reaches the surface of the glacier, a penalty function (eq. 5.27) limits the height of the aquifer to the glacier surface. However, water can leave the glacier to flow over the surface of the glacier via an enhanced diffusion scheme (eq. 5.30). At the base of the glacier, water can enter or exit the glacier from the sheet (eqs. 5.31). If ice is accreting, the englacial aquifer is also changing (eq. 5.32). Upstream and downstream boundaries have conditions on hydraulic potential. Together, all of these equations form an internally consistent scheme for simulating englacial water flow in temperate ice. 107 Chapter 6 Results and discussion In this chapter, I illustrate some of the behaviors of the physical model using numerical simu-lations. This departs from previous chapters where the focus was the construction of a physical model. Details of the numerical formulation can be found in Appendix F. Because much of the theory in the previous chapters has either not been assembled in the present format or has not been developed for sub- and englacial water flow, I separate these into clear water flows (Chapter 2) and sediment laden flows (Chapter 4), but I omit water flow through an englacial aquifer (Chapter 5). The relative importance of the terms in several of the governing equations need to be illus-trated. For example, the relative importance of the energy balance terms in equations (2.29) governing the thermal evolution of the water are important. Furthermore, for parameters that require careful selection, an exhaustive search is not possible. Fortunately, an exhaustive search is also not warranted. Glaciers such as those discussed in Chapter 1 have a finite range of physical characteristics. Numerical simulations presented here capture relevant processes and rates for these physical characteristics. 6.0.1 Longitudinal sections The subglacial water flow model is not designed to give prognostications for any particular glacier. Rather, the model is intended to diagnose and elucidate controlling processes. Ideally, simulations should capture relevant features of the physical model. Because the ice-surface slope and the bed slope are of first-order importance, suitable longitudinal sections are necessary. Longitudinal sections from the terminal overdeepening at Matanuska Glacier are shown in Figure 6.1. Lawson et al. (1998) published these sections as approximate hydraulic sections. Each has a water vent emerging at its downstream end. These sections display much information, both in the surface and bed structure and also in the supercooling ratios plotted below the sections. The bed slope in these sections is often —3 to —5 times the surface slope, where 108 Chapter 6. Results and discussion negative values imply that the bed slope is adverse to flow (Figs. 6.1b,d). Therefore, these bed slopes are well-above the steady-state supercooling threshold. Surprisingly, at 50 m from the terminus, Figure 6.1c displays a ratio that is below the steady-state subglacial lake criterion of — 11.0 (Clarke et al., 2005, using values from Table 6.2). Obviously, simulating these sections is useful if local process-based results are desired for Matanuska Glacier. However, the level of information in the Matanuska sections is unwarranted for elucidating the behavior of the governing equations. As an alternative, any number of parameters can be chosen for synthetic sections, and these sections have the benefit of being spatially smooth. Smoothness in space aids in the convergence of the numerical solver. Four relevant sections are shown in Figure 6.2. These sections are meant to mimic some of the features of Matanuska Glacier and other glaciers while retaining an uncomplicated structure. Relevant parameters for each of these sections appear in Table 6.1. Of these, the surface-to-bedslope ratio is the most telling. The section in Figure 6.2a has a threshold bed and surface slope ratio, —1.70, for the parameters listed in Table 6.2. The section in Figure 6.2b has a surface to bed slope ratio of —3.00 and a maximum ice thickness of 35.99 m. Both sections in 6.2c and 6.2d have surface to bed slope ratios that are below the steady-state threshold necessary for supercooling. These non-supercooled cases will be used to illustrate several key points. None of the synthetic sections is meant to illustrate entire glaciers. Rather, these are meant to illustrate processes occurring near a glacier terminus. 6.0.2 Model details Parameters and constants that are applicable to all model runs are presented in Table 6.2. Many of these values are well-established in the literature and need no further clarification. However, several require additional clarification. The natural variability of the pressure melting coefficient B, lies between —7.4 x IO""8 and —9.8 x 10~8 K P a - 1 depending on whether the water is pure or air-saturated, respectively. Alley et al. (1998) investigated effects of applying different values and found that the accretion rate is sensitive to 3 when water flux is high. More ice accretes when water is air-saturated. I thus choose a value of 3 that corresponds to pure water to yield conservative estimates of accretion. This value is consistent with the densities in Table 6.2 and equation (2.7). 109 Chapter 6. Results and discussion Distance upstream from terminus (m) T 1 1 r Distance upstream from terminus (m) Figure 6.1: Interpreted longitudinal sections from Lawson et al. (1998). (a) Figure 9a from Lawson et al. (1998). Height=0 m is the minimum ice elevation in this figure. The aspect ratio is 1 : 1. (b) Supercooling ratio derived from bed and surface data in Fig. 6.1a. Thick dashed line is the steady-state supercooling threshold (TZ = —1.70, eq. (2.44)). Values below this line indicate locations of supercooling for steady-state hydraulic conditions. Thin dot-dashed line is TZ = 0. (c) Figure 9b from Lawson et al. (1998). Height=0 m is the minimum ice elevation in this figure. The aspect ratio is 1 : 1, but the x-axis is not the same scale as (a), (d) Same as (b) but data from Fig. 6.1b. Minimum value at 50 m is TZ = —15.6. 110 Chapter 6. Results and discussion 250 200 150 100 50 0 Distance upstream from terminus (m) Figure 6.2: Simulated synthetic glacier sections. Subglacial water flows from left to right along the line at the interface of the ice and bed. (a) Threshold-overdeepened case with TZ = —1.70 (Figs. 1.3 and A. l ) . (b) Overdeepened case with TZ = —3, which is above the steady-state supercooling threshold., (c) Overdeepened case with TZ = —0.85, which is below the steady-state supercooling threshold, (d) Flat-bedded case. Dashed lines in (b), (c), and (d) are the surface and bed in (a). The aspect ratio is 1 : 1. Ill Chapter 6. Results and discussion Table 6.1: Parameters used to create the synthetic glacier sections in Figure 6.2 Threshold Above-threshold Below-threshold Flat-bedded Figure 6.2a 6.2b 6.2c 6.2d Start of overdeepening: 150.00 150.00 150.00 n/a Start of wedge: XW 107.59 107.59 122.84 n/a Elevation at outlet: 21.85 21.85 14.82 0.00 Upstream ice thickness: A 46.86 35.99 46.86 46.86 Bed to surface slope ratio: K -1.70 -3.00 -0.85 0.00 Surface slope: tancnr 0.10 0.06 0.13 0.19 Bed slope: tana b 0.17 0.17 0.11 0.00 Values of x, z, and Z have units of meters. The other parameters are unitless. For all sections, the coefficient Cx = 0.002 ml~2°x, ice thickness at the outlet Z0 = 0.00 m, the transition coefficient ox = 1, the location of the downstream water outlet XQ = 0.00 m, the location of the upstream water boundary x\ = 250.00 m, and the elevation of the upstream water boundary zf = 0.00 m. See Appendix A for a more detailed explanation. 112 Chapter 6. Results and discussion Another parameter that deserves attention is the compressibility of the subglacial hydraulic system jb. Its physical limit is the compressibility of the mixture of water, ice, and sediment in the flowing water and should be close to the compressibility of water, 4.5 x 1 0 - 1 0 P a - 1 . Clarke (2003) noted that this parameter is essentially numerical and governs the rate of convergence of the computational scheme. However, 7 b does not govern accuracy. I have tested numerous values for this parameter over several orders of magnitude using the sheet version of the code, with no noticeable difference in the results. As a result, I assign 7 b to the value of 1.0 x 10~7 Pa" 1 (Clarke, 2003). The flow-law exponent n also deserves some discussion. There is evidence that n is commonly less than 2 at low driving stresses (e.g., Duval et al, 2000). Low effective stresses are often found subglacially. For most examples contained within this thesis, closure is not dominated by creep. Rather, sheet closure is predominantly regelation, and creep is less important (Fig. F.5). In the case of negative effective pressures, n = 3 will overestimate the rate of opening for faster changes in sheet thickness. In this case, n = 3 provides conservative estimates of pressure response from ice closure. Thus, a higher value of n is not problematic in terms of the model. Furthermore, smaller values of n can be substituted easily into the formulation presented in Chapter 3. The grid spacing Ax is relatively arbitrary. Essentially, Ax needs to be small enough to capture the relevant processes. However, if Ax is too small, then the matrix inversions necessary to advance in time for the numerical differentiation formula become prohibitively expensive. In the limit of an extremely tight spatial discretization, time stepping is not possible because the matrices have high condition numbers. Knowledge of a horizontal length scale of any of the processes discussed in previous chapters would limit Ax. The most conspicuous of these is the ice-bed configuration of the longitudinal section to be simulated. For the synthetic sections, I take Ax as 5 meters unless otherwise noted. I have tested other discretizations and found the results presented here to be largely independent of the grid spacing In addition to parameter choices listed in 6.2, other parameters are also necessary. These are given in Table 6.3. The crevasse-cross sectional area is relatively small. This low value provides little buffering in terms of the response time of the subglacial system to the recharge rate. In terms of a simple geometry, the crevasse cross-sectional area corresponds to circular moulin of approximately 3.6 m diameter. The Darcy-Weisbach friction coefficient is high relative to subaerial streams. However, Clarke (2003) showed that appropriate flood discharge hydrographs 113 Chapter 6. Results and discussion Table 6.2: Model parameters Parameter Value Units Notes Physical constants: A 6.8 x I O - 2 4 Pa"" s- 1 Flow law coefficient (Paterson, 1994, p. 97) 4217.6 J kg- 1 K - 1 Specific heat of water at constant pressure 9 9.81 m s z Gravitational acceleration K% 0.5610 W m" 1 K - 1 Thermal conductivity of water L 3.336 x 105 J kg" 1 Latent heat of water n 3.0 [unitless] Flow exponent (Nye, 1953; Paterson, 1994) Pw 1000.0 kg m~3 Mass density of water at 0 °C Pl 916.7 kg m- 3 Mass density of ice P 1.781 x 10~3 N s m~2 Viscosity of water Derived constants: B 5.28 x 107 Pa s1/™ Flow law coefficient (= A1,n) Pr 13.39 [unitless] Prandtl number P -7.44 x I O- 8 K Pa" 1 Pressure melting coefficient (eqs. 2.7, 2.8) n -1.70 [unitless] Threshold bed to surface slope ratio (eq. 2. Numerical parameters: 7 b 1 x IO" 7 Pa" 1 Subglacial water compressibility Ax 2.5-5 m Longitudinal grid spacing At 30 s Maximum forward time step 0 m Minimum accreted ice thickness 114 Chapter 6. Results and discussion for Steele Glacier, Yukon Territory, Canada can be reconstructed using values of fa in the range 0 . 1 2 - 0 . 2 0 . I take 0.16 as the mean value. The crevasse recharge rate is taken to be 0.1 m 3 s - 1 . This value corresponds to a maximum summertime discharge for Matanuska Glacier, as described above. Finally, I explicitly state that the width of the sheet is 1 m. Thus, all simulations are per unit glacier width. Values of width-averaged parameters discussed in Chapter 2 are returned to their usual dimensions using this value. Table 6.3: Additional model parameters Parameter Value Units Crevasse cross-sectional area Ac 10 m 2 Darcy-Weisbach friction coefficient (eq. 2.24) fa 0.16 [unitless] Crevasse recharge rate Qmeit 0.1 m 3 s _ 1 Unit glacier width W 1 m Boundary conditions Meltwater from the glacier surface can flow into the subglacial water system through fractures, moulins, or crevasses. Because there is no surface component to the model, a feeder crevasse is constructed at the upstream end of the longitudinal section. The feeder crevasse acts as a reservoir. This reservoir is a simplification of upstream hydrology and is meant to loosely approximate upstream conditions. The crevasse has a cross sectional area normal to z, and a recharge rate feeds into the crevasse. The pressure difference between the crevasse and the subglacial water system determines the flux of water into the subglacial system. Upstream crevasse water temperature is assumed to be at the equilibrium melting point. If instead sensible heat were imported from the surface, the subglacial hydraulic system would melt open. Because the upstream crevasse is a simplification of upstream hydrology, I make the expedient assumption that there is no sensible heat; and water enters the subglacial system at the melting temperature. Details of the crevasse formulation appear in Appendix F.1.2 and are illustrated in Figure F.2. While the crevasse condition is a useful addition to the model, the selection of the crevasse cross-sectional area and the recharge rate are arbitrary. A large upstream crevasse buffers the evolution of the hydrostatic pressure at its base. Large crevasses represent an upstream 115 Chapter 6. Results and discussion hydraulic system with much storage while small crevasses represent the opposite. Given these representations, the cross sectional area of the upstream crevasse need not be physically exact but only justifiable for the simulation. Surface melt water flows through a glacier for only a fraction of a year. Opening and closing of the subglacial system are not treated here because they are extremely difficult to simulate. To simplify, I assume that the subglacial water system is capable of accepting water at the time that water flows into the crevasse. Furthermore, time integration is limited to a melt season. Duration and intensity of the melt season depends on climatic and geographical parameters which are not treated here. For simplicity, I make the expedient assumption that a melt season is approximately 100 days long (e.g., Pearce et al., 2003). There are two limits for steady state solutions for the height of water in the crevasse: (1) ei-ther the water height is at the glacier surface, or (2) the water is at the height of maximum bed elevation downstream of the crevasse. For long time intervals, the recharge rate will determine if and when the crevasse water elevation approaches either of these limits. If the crevasse water level is below the maximum downstream bed elevation, then water will flow in reverse from the subglacial system to the crevasse. In this case, recharge will also play an important role in returning the crevasse water elevation to a permissible level. In sum, the crevasse water height must be between these two limiting elevations in order for water to flow through the subglacial water system. For longer times, the numerical model seeks an equilibrium among crevasse water pressure, storage, and subglacial discharge. If the recharge rate is high relative to discharge to the subglacial system, then the crevasse height can rise to the glacier surface. Discharge through the subglacial system increases until a critical subglacial cross-section is reached. The crevasse then drains to a steady level. This decay often takes tens to hundreds of hours depending on the crevasse size and water recharge rate. If the converse is true, and the recharge rate is low, then the crevasse water elevation can approach its minimum value. Generally, long integration times are necessary for the water system to creep closed in this case. Depending on the simulation, recharge into the crevasse can take a value of 0.001-0.3 m 3 s - 1 m _ 1 based on field data. Discharge data from Matanuska Glacier and Storglaciaren are relatively important because glaciohydraulic supercooling at these locales is relatively well-described and the width-averaged water fluxes constrain upstream model conditions. At Mat-116 Chapter 6. Results and discussion anuska Glacier, summer discharge measured approximately 200 m downstream of the terminus can exceed 80 m 3 s _ 1 (Lawson et al., 1998). The width at the terminus is approximately 5 km, giving a width averaged discharge of about 0.02 m 3 s _ 1 m _ 1 . However, a large portion of the terminus is not active. Therefore, Alley et al. (1998) suggested that summertime discharge may meet or exceed 0.1 m 3 s _ 1 m - 1 . In addition the subglacial system may leak water to an englacial or. subglacial aquifer (Alley et al., 1998; Hooke et al, 1988), which then discharges to the hy-drologic system downstream of the overdeepening. Consequently, the width-averaged discharge at the ice-bed interface is most likely lower than 0.1 m 3 s _ 1 m _ 1 . During winter, proglacial water discharge continues at a rate of approximately 1 m 3 s _ 1 , either from water stored sub- or englacially. Accretion occurs but is slow relative to summertime conditions. At Storglaciaren, two streams, Nordjokk and Sydjokk, drain the glacier. These two streams have maximum summer discharges of approximately 2.5 and 1.5 m 3 s - 1 , respectively (estimated from Hooke et al, 1988; Hooke & Pohjola, 1994; Seaberg et al, 1988). The width of the glacier at the riegel is approximately 850 m (Fig. 1.4). The width-averaged discharge is therefore approximately 0.005 m 3 s _ 1 m _ 1 , which is considerably lower than the equivalent value for Matanuska Glacier. This is not surprising given that Storglaciaren is in a drier climate and is a small fraction of the size of Matanuska Glacier. Other glaciers have similar width-averaged discharges. Width-averaged discharge on Columbia Glacier, Alaska, approach 0.1 m 3 s _ 1 m _ 1 , but the water pathways are not well known (R.A. Walters et al. as cited in Alley et al, 1998). Bering Glacier, one of the largest glaciers in the Chugach Mountains, Alaska, has width averaged discharges that can exceed 0.3 m 3 s _ 1 m _ 1 in the glacier's non-surge phase (Fleisher et al, 1998). Findelengletscher, Switzerland has width averaged discharges of approximately 0.001-0.01 m 3 s _ 1 m - 1 and the water there travels through a subglacial cavity system (Iken & Bindschadler, 1986). Haut Glacier d'Arolla, Switzerland has a peak summer discharge of approximately 7 m 3 s - 1 m _ 1 and a width near the terminus of 1 km (Swift et al, 2005) leading to a maximum width averaged discharge of approximately 0.007 m 3 s - 1 m _ 1 . Bench Glacier is a relatively small glacier (~7.5 km2) in the Chugach Mountains. Maximum summer discharge is approximately 7 m 3 s _ 1 and the width is approximately 0.9 km (Riihimaki et al, 2005). These numbers yield a maximum width-averaged discharge of 0.008 m 3 s - 1 m - 1 . There are no overdeepenings at Bench Glacier, and probably none at Haut Glacier d'Arolla and Findelengletscher; however, these glaciers are nearly the same size as Storglaciaren 117 Chapter 6. Results and discussion and provide roughly equivalent width-averaged discharges in the range 0.001-0.01 m 3 s _ 1 m _ 1 . Columbia and Bering Glaciers are closer to the same size and have similar climatic conditions to Matanuska Glacier. Because the recharge rate is arbitrary, selection of the direction is largely by trial and error within the bounds of the physical system. However, Alley et al. (1998) showed that accretion is favored for higher discharge systems. In response to their estimates, I take 0.1 m 3 s _ 1 m _ 1 as a baseline discharge in an attempt to reproduce their accretion rates. At the outlet end of the longitudinal section, water pressure must be prescribed. Two simple options are either atmospheric pressure or ice overburden pressure at the terminus. Numerous field studies of supercooled vents demonstrate that water flow is artesian (e.g., Lawson et al., 1998; Pearce et al., 2003; Tweed et al, 2005). Consequently, a value of pressure slightly above atmospheric or ice overburden can also be used. However, the degree of overpressure is difficult to estimate a priori. Unless otherwise noted, pressure is prescribed as atmospheric. Initial conditions Nonzero initial values are necessary for pb, uh, T b , and, depending on cross-section, either Hh or 5 b . Initial conditions are discussed in Appendix F.l.3. A reasonable initial assumption for the water pressure is the ice overburden pressure. The water velocity can be assumed to be the velocity of steady state via equation (2.26). Initial temperature can be set via the Clausius-Clapeyron relationship (eq. 2.8). Values of ph, uh, and T b are readily purged in numerical simulations provided that initial values are reasonably close to an expected solution. The values of Hh or 5 b are more difficult to initialize because there are no accepted values of either for overdeepenings. Moreover, these values should vary through the overdeepening and depend on the mass balance of the water layer. In previous studies, Clarke (2003) and Ng (1998) chose Sb to match data from jokulhlaups. However, matching data from outburst floods is not a goal of this study. Alley et al. (1998) posited that distributed water systems in overdeepenings were cavities and used a value of 0.1 m taken from field studies (Walder & Hallet, 1979). For the sheet model, Hb values of 0.1 m are possible but high, and I thus choose 0.1 m as a default initial maximum water depth. This value translates into an average water depth of 0.087 m (see Appendix F.4.1). 118 Chapter 6. Results and discussion 6.1 Steady flow: simple hydraulic sheets The basic model presented in Chapter 2 with auxiliary relationships presented in Chapter 3 governs clear water flow. In this case, the subglacial hydraulic system is opened via a high melt rate. Closure occurs through accretion via freezing or a combination of creep and regelation. One important scenario for the clear water model is a steady forcing where the recharge to the upstream crevasse is kept constant. This boundary condition represents the simplest forcing function for the subglacial water system. Responses among terms within the governing equa-tions and different geometries are most easily illustrated for steady forcing because the subglacial water system is forced to find a configuration that matches input conditions, fllustrated config-urations may not be attainable in reality because of diurnal, seasonal, and other variations. I discuss the effects of these fluctuations later in this chapter. This boundary condition is meant to illustrate governing processes rather than give site- or scenario-specific information. Furthermore, discharge is included in the steady state accretion rate (eq. 2.42a). However, there is not a steady state as is implied in this equation because the melt rate must either be positive or negative for flowing water. Only for the trivial case of no water discharge is there an obvious steady state. It is thus useful to explore the behavior of the accretion rate under steady forcing. Unless otherwise noted, these simulations are specifically for hydraulically conductive sheets. I contrast conduits later in this chapter. 6.1.1 Constant recharge simulations Overdeepened case: Threshold section Clear water results for sheet thickness Hb, water pressure p h , water velocity ub, and water tem-perature T b are plotted in Figure 6.3 from the model presented in Chapter 2. In Figure 6.3a, from approximately 250-175 m from the terminus, the water depth decreases. From approxi-mately 175-60 m, the water sheet bulges; and from 60-0 m, the water sheet thins below its initial value. In Figure 6.3b, the water pressure is smooth throughout the section, with a relatively smooth bend in the contours at approximately 25 d and farther than 150 m from the terminus. Water pressure decreases smoothly from the inlet pressure to the outlet. The water velocity in Figure 6.3c, shows variations similar in overall form to Figure 6.3a, but opposite in sign. Values 119 Chapter 6. Results and discussion for the velocity are high—approximately 1 m s - 1 (Fig 6.3c). Where the water sheet is relatively thick, velocity is slow relative to areas where the water sheet is thin and vice versa. The water temperature illustrated is Figure 6.3d and bears some resemblance to Figure 6.3b; but there is additional structure in the plot for water temperature that is not evident in the water pressure. 100 £> 75 50 25 100 75 50 25 100 75 50 25 100 75 50 25 m\\\\V^yj)V -• i d I LI i 0.12 0.1 .§. 0.08 1c 0.06 40 30 «T i 20 E • l0 - °a 1.5 _ 1 0.5 % 0 -0.01 p -0.02 r -0.03 w E 250 200 150 100 50 Distance upstream from terminus (m) Figure 6.3: Threshold section output for a 100-day simulation under constant recharge condi-tions. Water flows from left to right along the line at the interface of the ice and bed. Time is along the vertical axis, (a) Water depth. Contour interval is 0.01 m. (b) Water pressure. Con-tour interval is 3 m. (c) Water velocity. Contour interval is 0.12 m s _ 1 . (d) Water temperature. Contour interval is 0.003 °C. From the four basic output variables a working hypothesis appears. Water sheet depth and water velocity are related to each other via a simple continuity relationship such that Qh = WHhuh. Thus, as the water sheet thickness increases, the water velocity must decrease to maintain a constant discharge. Furthermore, water pressure and water temperature appear to be related to one another via the Clausius-Clapeyron relationship. However, additional effects appear in the results for the temperature that are not clearly evident in the water pressure results. 120 Chapter 6. Results and discussion Discharge results shown in Figure 6.4 point to the monotonous nature of inlet and along-path discharge. In Figure 6.4a, the water height of the crevasse, a surrogate for inlet pressure head, (solid black line) is plotted against the left axis, and discharge from the crevasse to the subglacial sheet (solid gray line) is plotted against the right axis. The pressure head has a strong transient response for approximately the first 10 d before relaxing to a constant pressure head of 34.3 m. This value is well below the maximum allowable value of 46.86 m which is the upstream ice thickness. Furthermore, the inlet pressure head is less than the notation height of 42.96 m. Discharge to the subglacial sheet is constant after approximately three hours of integration. The along-path discharge is plotted in Figure 6.4b with time on the vertical axis. Approx-imately 0.2% of the array values fall out of the range 0.1 ± 0.0001 m 3 s _ 1 , and most of this variation occurs within the first 10 d. Once inlet pressure variations subside, the variations in water discharge also subside. The Reynolds number varies directly with discharge for sheet water flow (eq. 2.35a). Because discharge does not vary significantly, Figure 6.4 only shows"small variations in Reynolds number. The consequence is that the Reynolds number does not appear to be very telling in this case. For example, there are no large variations from the non-overdeepened portion of the section to the overdeepened portion. The important result is that water discharge is fully turbulent because the Reynolds number exceeds 2100. Key results from the discharge information are (1) sheet water discharge attempts to mimic the crevasse recharge rate, (2) the discharge is roughly constant over the section such that QB = WHbub where W = 1 m, and (3) Reynolds numbers are fully turbulent. In Figure 6.5a, the difference in temperature over the course of the model simulation ranges from about —0.001 to approximately 0.010 °C. The temperature depression shows that most of the flow path is above the melting point except for a small fraction that lies up the overdeepening. The area bounded by the zero degree isotherm is of particular interest. This swath starts approximately 20 d from the onset of the simulation. Eventually, it steps back to approximately 50 m from the terminus over the course of the simulation. Simultaneously, water upon exit from the terminus is not supercooled at 100 d. A similar feature in the same area also exists in the velocity simulation (Fig. 6.3c). This swath is more clearly evident in longer simulations (not shown), and the swath steps upstream to the onset of supercooling at about 107 m, which is the start of the wedge (Table 6.1). 121 Chapter 6. Results and discussion 250 200 150 100 50 0 Distance upstream from terminus (m) Figure 6.4: Discharge and Reynolds numbers in Figure 6.3. (a) Left axis: Elevation of water in the crevasse (black line). Black dotted line is the maximum water height and corresponds to the ice surface elevation. Black dashed line is the flotation limit. Right axis: Water discharge from the crevasse into the subglacial sheet (gray line), (b) Subglacial water discharge. Contour interval is 1 x 1 0 - 6 m 3 s _ 1 . Approximately 0.2 % of the array is not in the range 0.1 ±0.0001 m 3 s _ 1 . The colorbar has only one tick owing to its small range (0.099995, 0.100003). (c) Reynolds numbers corresponding to (a). The precision is indicative of the lack of variability and does not reflect additional knowledge of subglacial water flow. Contour interval is 0.25. The net accretion from this area of supercooling is shown in Figure 6.5b. A total of ap-proximately 0.017 m ice accretes near the terminus. The maximum accretion is nearly 0.015 m at the terminus at approximately 75 d. Because the supercooled swath is stepping upstream (Fig. 6.5a), accretion is not monotonic or sustained throughout the simulation; and some ac-creted ice is melted in the final days of the simulation. However, this simulation is a threshold simulation. The expectation is that there will be no or very little supercooling, which is evident in the amount of accretion. Relevant melt rates are shown in Figure 6.6. Figure 6.6a shows the model melt rate computed using the Dittus-Boelter semiempirical relationship (eq. 2.38b). Rates vary from about —0.6 x 1 0 - 5 to slightly over 8 x I O - 5 kg m - 2 s - 1 . The swath of negative melt rates corresponds directly 122 Chapter 6. Results and discussion » 50-P 25-200 150 100 _i i 50 0 Distance upstream from terminus (m) Figure 6.5: Temperature depression and ice accretion for the threshold case in Figure 6.3. (a) Temperature relative to freezing. Positive values indicate melt. Contour interval is 0.001 °C. (b) Amount of ice accreted. At t = 100 d, the maximum ice accreted is not at the terminus, but is approximately 5 m from the terminus. Contours indicate 0.001, 0.005, 0.010, and 0.015 m. to the swath of temperature depression (Fig. 6.5a). This feature is a result of the melt rate being directly dependent on the temperature depression. There is a secondary dependence on (uh)4'5 It is possible to compute alternative melt rates based on the output variables. Specifically, the approximate melt rates based on the steady state assumptions are useful for comparison (eqs. 2.9) (e.g., Rothlisberger & Lang, 1987). Figure 6.6b shows a melt rate using ice overburden pressure as an approximation from equation (2.42b). In this case the melt rate varies from about —0.4 x 1 0 - 5 to slightly over 21 x I O - 5 kg m - 2 s _ 1 . The pattern results from the small variation in the overall glacier ice thickness during the model run. High melt rates occur where ice is thickest. Tight contours from 150 to 100 m occur because of the onset of the overdeepening. Generally, the approximate melt rate is a factor of ±2-3 greater than the model melt rate. Furthermore the pattern is starkly different. Instead of approximating the water pressure as the ice overburden pressure in the steady state melt rate, it is possible to use the water pressure directly. Figure 6.6c shows the approximate melt rate using this substitution. Instead of being a better fit to the model melt rate, the result is worse. Average misfit melt rates are closer to ±5-20 times the model rate along the adverse slope. However, locally the misfit can be much greater. and (Hb)-V5. 123 Chapter 6. Results and discussion Distance upstream from terminus (m) Figure 6.6: Melt rates for the threshold case in Figure 6.3. (a) Melt rate using equation (2.38b). Contour interval is 0.5 x 10~ 5 kg m - 2 s _ 1 . (b) Approximate melt rate using ice overburden pressure from equation (2.42b). Contour interval is 1 x 10~ 5 kg m ~ 2 s _ 1 . (c) Approximate melt rate using water pressure from equation (2.32b) Contour interval is 5 x I O - 5 kg m ~ 2 s _ 1 . Overdeepened case: Above-threshold section The above-threshold section fully activates the glaciohydraulic supercooling phenomenon. In Figure 6.7a, the water depth shallows along the adverse slope from about 100 m to the terminus. The darkened curved contours are the area of supercooling; and it steps upstream towards the onset of the wedge at about 107 m. The minimum water depth is below 0.07 m. Beneath the flat-bedded portion, the sheet expands to above 0.18 m at 150 m and 100 d. The pressure evolution displayed in Figure 6.7b has a strong upstream contour from ap-proximately 0-50 d. This feature results from the input pressure taking a much longer time to approach a steady value (Fig. 6.8, left axis). The input pressure head evolves to a value of approximately 32 m, which is below ice flotation. However, the pressure head only reaches a steady value after 200 d (not shown). As a result, the pressure evolution in Figure 6.7b has not equilibrated after 100 d. For the first four days, the crevasse water elevation is at the glacier surface and above flotation. These first four days are clear in all of the plots for the above-threshold overdeepened section. After approximately 4 d, the water system has reorganized, and the upstream crevasse 124 Chapter 6. Results and discussion 200 150 100 50 0 Distance upstream from terminus (m) Figure 6.7: Above-threshold section output for a 100-day simulation under constant recharge conditions. Water flows from left to right along the line at the interface of the ice and bed. Time is along the vertical axis, (a) Water depth. Contour interval is 0.01 m. (b) Water pressure. Contour interval is 3 m. (c) Water velocity. Contour interval is 0.1 m s _ 1 . (d) Water temperature. Contour interval is 0.002 °C. pressure evolution follows a damped decay pattern. For the threshold case, the damping takes 30 d (Fig. 6.4a, left axis). For the above-threshold case, the damping to a steady value takes more than 200 d. In both instances, the discharge into the subglacial system from the crevasse equilibrates relatively rapidly to the recharge rate (Figs. 6.4a, 6.8, left axes). Because water discharge is nearly constant across the overdeepening (Fig. 6.8, right axis), the velocity increases where the water depth decreases. This behavior is illustrated in Figure 6.7c. Water velocities up the adverse slope exceed 1.4 m s _ 1 where the water depth is shallow. Where the sheet bulges, water velocities are lower. From Figure 6.7d, water temperature again varies in a fashion very similar to water pressure. However, water temperature relative to the melting point appears to be more closely related to the water depth and velocity (Fig. 6.9a). The maximum amount of ice accreted is slightly above 0.039 m at 100 d and 70 m (Fig 6.9b). This ice thickness is over double the maximum 125 Chapter 6. Results and discussion _ 40 £ 38 | 36 1 34 2 32 LLI 30 10.2 Maximum 10 20 30 40 50 60 Time (d) 70 80 90 100 Figure 6.8: Boundary conditions for the above-threshold case (Fig. 6.7). Left axis: Elevation of water in the crevasse. Black dotted line is the maximum water height and corresponds to the ice surface elevation. Black dashed line is the notation limit. Right axis: Water discharge. Gray line is the input into the subglacial water system. Light gray dotted line contained within the gray line is the effluent subglacial water discharge. ice accreted ice thickness for the threshold case. Furthermore, the area of accreted ice extends along most of the adverse slope by 100 d. There is a drop in the accreted ice thickness near the terminus with time because the water temperature drops below freezing, rises above freezing, and then drops below freezing again. This feature is shown in the location of the zero degree isotherms in Figure 6.9a with a swath that is slightly above zero separating the two zero degree isotherms that are farthest right. 0.006 -j. 10.004 t . 0.002 % 0.000 200 150 100 50 Distance upstream from terminus (m) Figure 6.9: Temperature depression and ice accretion for the above-threshold case in Figure 6.7. (a) Temperature relative to freezing. Positive values indicate melt. Contour interval is 0.001 °C. (b) Amount of ice accreted. At t = 100 d, the maximum ice accreted is not at the terminus, but is approximately 75 m from the terminus. Contours are 0.001, 0.01, 0.02, and 0.03 m. 126 Chapter 6. Results and discussion In Figure 6.10a, the modeled melt rate follows the temperature depression inasmuch as the zero degree isotherms are also the zero-values for the melt rate. The patterns for both the temperature depression and melt rate are similar. The steady state melt rate shown in Figure 6.10b appears to be governed more by the ice overburden pressure than the actual flow. The approximate melt rates are approximately ±2-3 times the modeled melt rates along the adverse slope. 250 200 150 100 50 0 Distance upstream from terminus (m) Figure 6.10: Melt rates for the above-threshold case presented in Figure 6.7. (a) Melt rate using equation (2.38b). Contour interval is 0.5 x 10 - 5 kg m - 2 s _ 1 . (b) Approximate melt rate using ice overburden pressure from equation (2.42b). Contour interval is 1 x I O - 5 kg m - 2 s _ 1 . (c) Approximate melt rate using water pressure from equation (2.32b) Contour interval is 2.5 x 10~5 kg m~ 2 s _ 1 . Results are plotted in Figure 6.10c for the case when water pressure is used in the approx-imate melt rate instead of ice overburden pressure. This plot resembles many of the others inasmuch as there is a swath. However, the melt rates do not come close to giving the proper rates or areas of melt. Except for the obvious statement that water freezes along the adverse slope, it is clear that this approximation does not capture the areas or rates very well relative to the model melt rate. Reconstructing the time-dependent local balance of temperature is useful because it is possi-ble to discern the governing terms in the heat balance. From equation (2.29b), the local balance 127 Chapter 6. Results and discussion is + c w A r (1 - nj,) p» + n > w ( t i b ) 2 ~ V rap 7 " T T — : — Sensible heatv V ' Specific kinetic energy Hb p w c w (6.1) Viscous dissipation The latent heat, sensible heat, and kinetic energy terms all have the same melt rate prefactor outside the left, vertical brace. The actual values for each of these five terms is not necessarily important, but the relative magnitudes and patterns are important. These are plotted in order of appearance in Figure 6.11. Immediately apparent from Figure 6.11 is that the divergence and viscous dissipation terms (Fig. 6.11a and e, respectively) are clearly dominant and nearly balance each other with ranges up to ±50 x 10~5 °C s _ 1 . In Figure 6.11b, the latent heat term is the next largest term. The specific kinetic energy term is surprisingly large in Figure 6.lid. The sensible heat in Figure 6.11c is negligible because it is approximately five orders of magnitude smaller than the largest terms. Both the viscous dissipation and specific kinetic energy terms depend on the water velocity. From equation (2.24), the viscous dissipation goes as (ub)3 while the kinetic energy term goes as (ub)2. In consequence, it is not surprising that these two terms bear a resemblance to Figure 6.7c. Furthermore, because the velocity is closely coupled to the change in water thickness via Qb = WHbuh when discharge is constant, it is not surprising that these terms are coupled to the ice accretion rate. In addition, the latent heat term is coupled directly to the ice accretion because this term is directly linked to the phase change. In Figure 6.7a, the heat divergence term is more difficult to interpret as-is because it re-lates how much heat is available to move through the subglacial water system. Substituting the Clausius-Clapeyron relationship (eq. 2.8) gives —uhdTb/ds ~ —uh3dp°/ds. The result is that an approximate heat divergence follows from the subglacial water pressure gradient. This approximation appears in Figure 6.12a, and the residual appears in Figure 6.12b. The residual has a similar pattern to the approximate melt rate in Figure 6.10 which is not a coincidence given that the same approximation appears in both figures. The obvious conclusion is that 128 Chapter 6. Results and discussion 100 g 75 | 50 F 25 100 g 75 * 50 i= 25 100 g 75 ® 50 i= 25 10 r ? in o 5

0.1 m) occurs at the upstream end of the adverse slope. At the upper end of the adverse slope, sedimentation occurs, but is slow relative to the upstream end, especially after day 40, when little sedimentation occurs. Figure 6.42a illustrates that accretion is low for the threshold case. Total accretion is slightly above 0.01m near terminus. This result is consistent with all other simulations of threshold sec-tions, where the accreted ice thickness is near zero. Furthermore, accretion begins downstream of the onset of the adverse slope. After 25 d, accretion begins at 55 m from the terminus, which is downstream of the onset of the supercooled slope at about 107 m. The total volumetric concentration of sediment in the ice varies near the terminus but is low. Relatively high concentrations occur for low accreted ice thicknesses. These concentrations occur prior to 37 d. After 37 d, concentration drops off in the accreted ice. The daily cycles illustrated in Figures 6.43a-c detail the relationship of crevasse recharge, crevasse water elevation, and subglacial discharge for days 80-85. Despite large variations in the recharge, there is little variation in the discharge to the subglacial system from the upstream crevasse. Downstream discharge mimics upstream discharge in the same manner as 172 Chapter 6. Results and discussion 250 200 150 100 50 0 Distance upstream from terminus (m) Figure 6.41: Daily cycles with erosion: Threshold case, (a) Thickness of the subglacial water sheet. Contour interval is 0.01 m. (b) Total ice closure. Contour interval is 3 x 1 0 - 6 m. (c) Change in bed elevation. Positive values indicate sedimentation. Contour interval is 0.01 in. all previous simulations. In addition, the crevasse water elevation lags the recharge rate by about three hours. During nighttime, the crevasse water elevation drops nearly linearly because the subglacial discharge is approximately constant. The daily variations in melt rate illustrated in Figure 6.43c show that accretion occurs from about midnight to noon. Melt occurs during the remainder of the day. Along the flat-bedded portion, melt occurs at all times but is highest from 1600h to 1700h. Accretion occurs only on the adverse slope. Because water discharge is constant, the system is able to sustain accretion at night. This result contrasts with the clear water case where discharge was zero and accretion was negligible at night. Figure 6.43e shows the total sediment supply, both bed load and suspended load, from the bed. Where the supply is negative, sedimentation occurs. During nighttime, there is little sedimentation. During the day, there is much sedimentation with most of it occurring in the early part of the afternoon along the flat-bedded portion of the section. Sedimentation drops off sharply from 140 to 75 m from the terminus as the water loses its sediment load. Because the water is losing its sediment load, there is little sediment that can be accreted. Thus, in Figure 6.42b, the drop in sediment concentration occurs because accretion happens at night while sediment transport occurs during the day. 173 Chapter 6. Results and discussion 75 50 25 75 50 25 a • 1 _ — • 1 1 1 • I -• — J 0.01 40.005 X 2 - k 250 200 150 100 50 0 Distance upstream from terminus (m) Figure 6.42: Daily cycles with erosion: Threshold case, (a) Accreted ice thickness. Except for the first contour at 0.001 m, the contour interval is 0.01 m. (b) Concentration of sediment in the basal ice. Contour interval is 2 x 10~5. Blocky contours upstream prior to day 25 are a numerical artifact and result from low accretion ice thickness. Prior to day 37, sediment transport and accretion occur at nearly the same time (not shown). During the early part of the simulation, accretion still occurs during the shoulder times of morning and evening, but there is enough sediment in the flow that elevated concentrations are captured by the accreted ice. Figures 6.44a-c show the total sediment load in the subglacial system, the suspended load, and the bed load, respectively. Sediment load is largely dependent on the daily discharge being relatively high. Even though the discharge varies by only about 0.001 m 3 s _ 1 from daytime to nighttime, this is enough to stimulate sediment transport. In addition, concentration in the water exceeds concentration in the ice by two orders of magnitude. Both suspended and bed load sediment concentrations are low. These are approximately factors of 5-10 lower than concentrations observed at Matanuska Glacier (Lawson et al., 1998). Water velocities range from 0.3 to 0.5 m s - 1 . In terms of the Hjulstrom diagram (Fig. 4.2), these velocities should imply erosion. However, the Hjulstrom diagram was developed for rivers with higher hydraulic radii than simulated subglacial sheets. For comparison, Reynolds numbers have a daytime peak at about Re=2700 during days 80 to 85. At night, Reynolds numbers drop to about 1600. In consequence, the subglacial system is tending toward a low-flow state that is not accurately captured by the Hjulstrom diagram. In sum, the hydraulic system becomes choked with sediment over the course of the simu-lation. Silting of the hydraulic system prevents significant drainage of water. Because of the 174 Chapter 6. Results and discussion Q m e l t Q Distance upstream Distance upstream . 3 -1 . (m3 s~1) from terminus (m) from terminus (m) (m s ) Figure 6.43: Daily cycles with erosion: Threshold case, (a) Recharge into the subglacial water system, (b) Water elevation in the upstream crevasse. Dark dotted line is the glacier surface elevation. Dark dashed line is the notation elevation. Light dashed line is the elevation of the glacier terminus, (c) Water discharge. Dark line denotes water discharge into the subglacial water system. Light dashed line overlying the dark line denotes water discharge at the terminus, (d) Melt rate. Contour interval is 2 x 10 - 6 kg m - 2 s - 1 . (e) Sediment supply from the bed. Positive values indicate erosion. Contour interval is 3.6 x 10 - 6 kg m - 2 s - 1 . blockage, the numerical solver eventually fails because there is no optional routing for the dis-charge. Sediment choking is compatible with field observations. Furthermore, as the subglacial system chokes with sediment, accretion switches from occurring during morning and evening to occurring at nighttime. As a result, ice that accretes at nighttime has low sediment concentra-tions. Overall, ice accretion is low when the hydraulic system chokes. 175 Chapter 6. Results and discussion x s b (gL - 1 ) * * (gm- 1 s - 1 ) 1 2 3 4 0.1 0.2 0.3 0.4 0.3 0.6 0.9 Distance upstream Distance upstream Distance upstream from terminus (m) from terminus (m) from terminus (m) Figure 6.44: Daily cycles with erosion: Threshold case, (a) Total concentration of sediment in subglacial water. Contour interval is 5 x 10~5. (b) Concentration of suspended sediment in subglacial water. Contour interval is 0.1 g L - 1 . (c) Bed load transport rate in subglacial water (converted from volume concentration and velocity). Contour interval is 0.2 g m _ 1 s _ 1 . 176 Chapter 6. Results and discussion Overdeepened case: Above-threshold section The crevasse water elevation reaches the surface of the glacier after 21 d (Fig. 6.45a). During the remaining days of the simulation, crevasse water elevation reaches the surface of the glacier during the daytime. In contrast to the threshold section, the above-threshold section is able to drain completely at night. The crevasse is able to drain in the above-threshold case because the ice thickness is about 11 m thinner than the threshold case. After 21 d, the maximum daily discharge to the subglacial hydrology from the crevasse decreases from 0.05 to about 0.008 m 3 s _ 1 . Minimum nighttime discharge is negligible. There is not a particular discharge that the above-threshold case strives to maintain; however, the subglacial discharge decreases significantly. 0 25 50 75 100 Time (d) Figure 6.45: Daily cycles with erosion: Above-threshold case, (a) Water elevation in the upstream crevasse. Uppermost dotted line is the elevation of the glacier surface. Dark dashed line indicates the equivalent notation elevation. Lowermost, gray dashed line is the elevation of the terminus, (b) Discharge from the crevasse to the subglacial water system. Minimum water depths are below 0.02 m but above 0.01 m after 80 d (Fig. 6.46). The water depth responds nearly in unison over the entire length. Downstream water depth evolution precedes the upstream end by a day or two. Much of the evolution of the water sheet occurs within the first 21 d when the sheet drops over 0.03 m. Figure 6.46b illustrates that ice intrusion into the sheet is negligible. Maximum intrusion occurs upstream and is about 2 x 10~5 m. In terms of water depth, ice intrusion is almost irrelevant. However, sedimentation is not negligible. Figure 6.46c shows that over the course of 100 d, maximum sedimentation of about 0.1 m occurs at the start of the adverse slope. Sedimentation 177 Chapter 6. Results and discussion also occurs along the adverse slope for most of the first 50 d. During the final 50 d, the erosion steps upstream from the terminus to about 100 m from the terminus by the end of the simulation. 250 200 150 100 50 Distance upstream from terminus (m) Figure 6.46: Daily cycles with erosion: Above-threshold case, (a) Thickness of the subglacial water sheet. Contour interval is 0.01 m. (b) Total ice closure. Contour interval is 3 x 10 - 6 m. (c) Change in bed elevation. Positive values indicate sedimentation. Contour interval is 0.01 m. Ice accretion is also important in reducing water depth over the course of 100 d. Figure 6.47 shows that maximum accretion is at the terminus and is over 0.06 m at 100 d. This amount of accretion is roughly double the amount for daily water cycles with clear water flows. The spatial lag between the onset of the adverse slope and the onset of accretion is also smaller for the case with sediment motion. The volume concentration of sediment in the ice varies but is greatest near the terminus. Maximum values are about 1.3 x 10~4. After about 60 d, the concentration of sediment in the ice begins to decrease. Meltwater recharge to the upstream crevasse, crevasse water elevation, and discharge in the subglacial water system are illustrated in Figures 6.48a-c, respectively. Response of the crevasse water elevation follows the onset of daytime conditions at 0900h. Maximum crevasse water elevation is reached by noon, and crevasse water elevation begins to drop at 2000h. The system has a time lag such that the minimum water elevation is reached at about 0700h. The daily cycle then repeats. 178 Chapter 6. Results and discussion Distance upstream from terminus (m) Figure 6.47: Daily cycles with erosion: Above-threshold case, (a) Accreted ice thickness. Except for the first contour at 0.001 m, the contour interval is 0.01 m. (b) Volume concentration of sediment in the basal ice. Contour interval is 2 x 10~5. Over the last five days of the simulation (Fig. 6.48d), the melt rate is largely negative, indicating of accretion. During the morning, daytime, and evening, melt rates are strongly positive along the flat-bedded portion of the section. From 0700h to 0900h, during the end of the night phase, melt rates are negative along the flat-bedded portion. However, the lack of discharge precludes accretion. Downstream, accretion occurs from about 100 m to the terminus and is strongest during the morning and evening. This result is consistent with results from clear water diurnal forcing for the above-threshold section. Figure 6.48e reveals that sediment supply from the bed is negative upstream for most of the day. Along the final 75-100 m from the terminus, erosion occurs during the daytime. Little sediment is available to the hydraulic system during the morning, evening, or night portions of the diurnal cycle. Small amounts of sediment are present throughout the night, and these are deposited nearly uniformly along the section. Sediment concentrations illustrated in Figures 6.49a-c, show that sediment transport de-creases over the final five days during the daytime. Total sediment concentrations are highest upstream. As water flows downstream, sediment is deposited along flat-bedded portion of the section. Suspended sediment load and bed load are lower than field studies, but are higher than the concentration of sediment accreted in the ice. This result corresponds to sediment transport occurring largely within the daytime portion of the diurnal cycle. Very little sediment is mobilized during the morning or evening when accretion rates are higher. Maximum daily discharge, therefore, crosses a threshold at approximately 60 d when less 179 Chapter 6. Results and discussion Figure 6.48: Daily cycles with erosion: Above-threshold case, (a) Recharge into the subglacial water system, (b) Water elevation in the upstream crevasse. Dark dotted line is the glacier surface elevation. Dark dashed line is the flotation elevation. Light dashed line is the elevation of the glacier terminus, (c) Water discharge. Dark line denotes water discharge into the subglacial water system. Light dashed line overlying the dark line denotes water discharge at the terminus, (d) Melt rate. Contour interval is 2 x 10~6 kg m - 2 s _ 1 . (e) Sediment supply from the bed. Positive values indicate erosion. Contour interval is 3.6 x 10~6 kg m - 2 s _ 1 . sediment is mobilized. As a result, less sediment is captured by the accreted ice, and the accreted ice becomes cleaner over the final 40 d of the simulation. Maximum daily water velocities over the final five days of the simulation reach 0.5 m s _ 1 . In terms of the Hjulstrom diagram (Fig. 4.2), these velocities should produce erosion albeit by a nominal amount given a D50 of 0.001 m. Maximum daytime Reynolds numbers for the final five days peak at about 5000, which is a consequence of shallow water depths. As a result, the Hjulstrom diagram is not necessarily a valuable tool for estimating whether subglacial water sheets will erode or deposit sediment because the diagram is intended for use with higher hydraulic radii and higher Reynolds numbers. Still, the Hjulstrom diagram does indicate that 180 Chapter 6. Results and discussion >. b x10" 4 Xb (gL"1) X* (gm- 1 s- 1 ) 0) E P 1 1 i 1 — 1 1 1 1 J J ! ! !_] "TT7 /— T n A 1// / ^ l-L—f my / / r n / m , / / / / ^\ ..) c a b 250 200 150 100 50 Distance upstream from terminus (m) 200 150 100 50 Distance upstream from terminus (m) 200 150 100 50 0" Distance upstream from terminus (m) 100 99.5 99 98.5 98 _ 97.5 | 97 96.5 96 95.5 95 Figure 6.49: Daily cycles with erosion: Above-threshold case, (a) Total volume concentration of sediment in subglacial water. Contour interval is 3 x 10 . (b) Concentration of suspended sedi-ment in subglacial water. Contour interval is 0.1 g L - 1 . (c) Bed load transport rate in subglacial water (converted from volume concentration and velocity). Contour interval is 0.2 g m _ 1 s _ 1 . the erosion-transportation-sedimentation thresholds are bunched together, which is what the above-threshold simulations indicate. In sum, the Hjulstrom diagram is loosely analogous to what occurs in the subglacial alluvial system but is not particularly revealing. The above-threshold system also chokes with sediment. However, this choking takes longer to accomplish because the system moves less sediment than the threshold section by about a factor of two. The net result is that the hydraulic system is able to evacuate water more effectively despite the additional effect of accretion during along the adverse slope. 181 Chapter 6. Results and discussion Overdeepened case: Below-threshold section Crevasse water elevation for the below-threshold section rises to the surface of the glacier by simulation day 9 (Fig. 6.50a). After 21 d, the crevasse water elevation rises from the minimum value of the terminus elevation. This result parallels the threshold section where the crevasse also did not drain to its minimum value. As shown in Figure 6.50b, maximum daytime discharge for the below-threshold section decreases over the course of the simulation. Nighttime discharge rises above zero beginning on day 21. This result is consistent with the upstream crevasse retaining additional water at night. By the end of the simulation, subglacial discharge varies between 0.003 and 0.005 m 3 s _ 1 . This variation is greater than that of the threshold section and adds to the numerical stability of the below-threshold section for the course of the entire 100 days. 0 25 50 75 100 Time (d) Figure 6.50: Daily cycles with erosion: Below-threshold case, (a) Water elevation in the upstream crevasse. Uppermost dotted line is the elevation of the glacier surface. Dark dashed line indicates the equivalent flotation elevation. Lowermost, gray dashed line is the elevation of the terminus, (b) Discharge from the crevasse to the subglacial water system. Most of the evolution of the subglacial sheet is complete by day 21, as illustrated in Fig-ure 6.51a. For the final 79 days of the simulation, the sheet shallows by less than 0.01 m. In this instance, ice intrusion is negligible because the total closure over 100 days is tens of microns (Fig. 6.51b). The shallowing of the water depth primarily results from deposition of sediments along the section (Fig. 6.51c). Maximum sedimentation is about 0.15 m and occurs at the upstream end of the adverse slope. Along the adverse slope, there is a subdued pattern where additional 182 Chapter 6. Results and discussion sedimentation occurs approximately 75 m from the terminus. Distance upstream from terminus (m) Figure 6.51: Daily cycles with erosion: Below-threshold case, (a) Thickness of the subglacial water sheet. Contour interval is 0.01 m. (b) Total ice closure. Contour interval is 5 x 10 - 6 m. (c) Change in bed elevation. Positive values indicate sedimentation. Contour interval is 0.01 m. Figure 6.52 illustrates the components of the mass balance of the sheet for the final five days of the simulation. The melt rate is always moderately positive in both space and time. Only during the late night part of the diurnal cycle does the melt rate approach zero. This observation draws sharp contrast to the clear water example where accretion was possible during the morning and evening parts of the diurnal cycle. Melt during the day opens the water system at about double that of the nighttime closure rate (Fig. 6.52b). Sediment supply from the bed is roughly twice the magnitude of the melt during the daytime (Fig. 6.52c). There is no sedimentation at night when water velocity is too low to move sediment. Sedimentation concentrates during the morning portion of the daily cycle along the flat-bedded portion of the section. During nighttime, competition between the closure and melt rates acts to close the sheet. However, onset of melt and cessation of closure during the morning open the sheet. Accommodation space created by melt is quickly filled in by the movement of sediment during the morning. As Figure 6.53a illustrates, the below threshold case carries much sediment. Maximum daily total volumetric sediment concentration is approximately 0.001 near the inlet. During daytime 183 Chapter 6. Results and discussion Melt rate x 10"5 (kg m"2 s"1)Closure rate x 10"8 (m s"1) M / S x 10"5 (kg m"2 s"1) 0.5 1 1.5 -8 -6 -4 -2 -4 -3 -2 250 200 150 100 50 0 200 150 100 50 0 200 150 100 50 0 Distance upstream Distance upstream Distance upstream from terminus (m) from terminus (m) from terminus (m) Figure 6.52: Daily cycles with erosion: Below-threshold case, (a) Melt rate. Contour interval is 1.16 x 10~6 kg m~ 2 s - 1 . (b) Closure rate. Contour interval is 1 x 10~1 0 m s - 1 . (c) Sediment supply from the bed. Contour interval is 0.75 x 10~5 kg i n - 2 s - 1 . conditions, the concentration decreases along the section. This drop increases along the adverse slope. At night, sediment transport is negligible. The evening transition from sediment-laden to low-sediment conditions occurs over a few hours and is not as sharp as the morning transition from low to high sediment load. Relative to the threshold section, water flow in the below-threshold section carries signifi-cantly more sediment up the adverse slope. Along the flat-bedded portion of the section, bed load is approximately three times greater in the case of the below-threshold section (Fig. 6.53b). Suspended load along the flat-bedded portion of the below-threshold section is also about three times greater (Fig. 6.53c). Despite the fact that the below-threshold section flushes more sediment through the system, the water system chokes with sediment. The added sediment reduces the water depth consider-ably during the first part of the simulation. Nighttime melt opens the channel enough to create accommodation space which is then filled with sediment during daytime conditions. 184 Chapter 6. Results and discussion xbxi ( A 5 ) 1Z yx0 - xw J If the upstream ice thickness ZQ is known and ox = 1, then xw takes the value xw = xo + ^ ^ + ^ Y x l - x o ( n + l) + nxu)2 + — {Z0-Zl)j . (A.6) 213 Appendix A. Synthetic glacier section 40 20 '3 -20 1 : 1 Ice ; Substrate ; 250 200 150 100 50 Distance upstream from terminus (m) Figure A . l : A synthetic glacier section. Subglacial water flows from left to right along the line at the interface of the ice and substrate. Parameters are given in Table A . l . The aspect ratio is 1 : 1. Table A . l : Parameters used to create Figure A . l Parameter Value Units x 0 0.00 m Xyj 107.59 m 150.00 m XI 250.00 m 21.85 m zl 0.00 m 4 0.00 m A 46.86 m Cx 0.002 m l - 2 o * Ox 1 [unitless] TZ -1.70 [unitless] tana r 0.10 [unitless] tana^ 0.17 [unitless] 214 Appendix B Grain parameterization for ice intrusion For spheres that are intersected by a plane above their equator, the volume above that plane is a spherical section less the volume of the inscribed cone below it. The cone is inverted with its vertex at the center of the sphere and its base along Hh. The spherical section plus the cone is below the plane, and the spherical cap is above the plane. The plane intersects the sphere at a height Hh which is greater than the sphere radius v\. Figure B . l : Figure illustrating the relationship of rj, Hb, and ip For Hb > ri, the volume is that of the spherical cap. The cap is the volume of the spherical segment through a polar angle ip minus the inscribed cone below the cap, p2ir pn /.cos-1((iyb-r i)/r i) 7 T / u \ / u / u\ 2 \ Ve:i = / / / R2siniPdipdRde--{Hh-ri) 12riHh-(Hh) ), = I (2n - Hh) (2r2 + nHh - ( i / b ) 2 ) , (B.l) where Ve-i is the effective volume, R is a radius, and tp and 9 are angles. The effective area is simply the area of the base of the cone, (B.2) 215 Appendix B. Grain parameterization for ice intrusion For spheres intersected by a plane below their equator, the relationship takes a similar form. The spherical cap is the bulk of the sphere when Hb < rv In this case, the section of interest is the spherical section plus the inscribed cone. The cone has its vertex at the center of the sphere and its base along Hb. A similar integration that yielded equation (B.l) yields, rl-K rn /.T-cos-'((r 1-i/ b)/r 1) . / , 2 \ Ve..i = B2 sini>dipdRdO+ -(ri-Hb) ( (Hb) -2riHb) Je=oJR=oJip=o 3 V / \ \ ; J = | (2n - Hb) (2r2 + nHb - ( t f b ) ^ . (B.3) It is not surprising that equations (B.l) and (B.3) are exactly equal. Once the ice has passed the equator, the area normal to the force is simply the area of the equator. The effective area is then that of the maximum cross section of the sphere, ae-.i = 7T7- 2. (B.4) If the geometry is a hemisphere resting on the same datum as Hb, then the spherical cap takes a similar form to equation (B.l), rln pn pcos-l(Hh/n) / / , \ 2 \ Ve:i = R2 sin ip dip dRd8 - -Hb [rj - (Hb) = l(2rl~3Hbr2 + (Hb)3y (B.5) Without an extra length of radius, the effective area is similar to equation (B.2), ae:i = TT (r2 - (Hbfy (B.6) 216 Appendix C Solution of the equations governing ice intrusion Because the closure relationship is dependent on the number of grain sizes and the closure velocity, the solution must be simultaneous for all of the variables. Solving for the bed closure velocity and the relevant stress terms requires a nonlinear, iterative solver. In this case, I employ a Newton scheme. The Newton solver requires error estimates G, on the functions as well as the Jacobian terms. C . l Three equation model One method is direct solution of the 2./V + 1 equations (3.17), (3.12), and either (3.13) or (3.14b) for the 2N+1 variables: of, C T J , and v. The first two conditions and their accompanying Jacobian terms (dGi/d{ospOj, v}) are straightforward, G1(V, 0) into equation (D.8) yields a coefficient of 1/ {(1 — A')(l — nj^p1}, which is similar to the equivalent coefficient in equation (2.17). On the surface, equations (D.9) and (2.15) look quite different. However, both equations are essentially the same with the exception of the space created by the supply of sediment from the bed. The small corrections for sediments in the flow could largely be ignored without sacrificing accuracy in the pressure solution. These corrections include: eliminate the derivatives of the sediment concentration in the flow, and eliminate the supply of sediment as ice melts, ff the bed were aggressively aggrading or eroding, it would be difficult to eliminate either of the bed sediment supply terms. However, if the bed was fairly stable, then these terms could also be eliminated. 222 Appendix D. Sediment mass balances A n alternative formulation The balance of mass for the water, M w = f [ (1 - Xb)pwdSds, Js Js* (D.10) can be drawn separately from equation (D.l). Carrying out the same steps above, the rate of change of the cross section is dSh d (Sbub) 1 b „w dt + ds 1 + (1 - A b) pw d>P dt ,dA b ~ds~ 1 - AO (1 - nl) p< (1 - n*) ff (D.ll) Rewriting this equation in the spirit of equation (D.9) gives dpb dt 7 b 5 b f dSh d (Sbub) b „w dt ds + 1 + n 'p w (1 _ A') (1 - nj,) pj ( 1 - A b ) p ' m + S°p d\b hdXb + ub dt ds „ S - W (D.12) This equation bears more similarity to the clear water equation (2.15) than equation (D.9) does. As an addition simplification, the derivatives of A b could be eliminated without affecting the rate of change of the pressure. D.l.2 Sheet Writing the equations for the sheet is a matter of converting the relevant equations from the foregoing discussion. The evolution of the height of the sheet is dHb 1 dt (1 - AO (1 - nj,) pi + n ip w + A' ( l - nj,) ff fdHb\ * s 1 + 1 W npp ( l - A 1 ) ( l - n ' V j + + V dt yclose ( l - n s ) p s ' (D.13) 223 Appendix D. Sediment mass balances where the closure relationship is derived in Chapter 3. The pressure equation can also be appropriated from the discussion above, dph ~dt 7b#b f OHb d (Hbuh) 1 dt ds \bf? + {l - \h)pv #b ( p s - P w ) dA b hdXh dt + 1 + < P W (1 _ A i ) ( i _ n i ) pi rh + + 1 + < p w (1 - n») p* as (D.14a) dph ~dt 7 b#b f a # b a (# b u b ) a i ds (1 - A b )p w aA b b ! H + U d\h ds + 1 +

where equation (D.14a) comes from the total mass balance and equation (D.14b) comes from the mass balance of the water. These equations are subject to the same discussion of approximation that applied to equations (D.9) and (D.12) as far as eliminating terms. D.2 Basa l ice sediment concentration Chapter 4 requires equations that govern the thickness of the accreted ice and the sediment concentration in the ice. Sediment can accrete to the base of the glacier when water freezes there. Alternatively, sediment can be lost from the glacier when ice melts to enter the water flow. The result of this jump condition is that there is an asymmetrical condition on the concentration of sediment in the ice. The asymmetry appears in the mass balance for the accreted ice layer. I make the expedient assumption that the basal ice layer has a constant vertical sediment concentration. Field evidence shows that this is not the case with various grain size distributions and stratigraphic relationships (e.g., Lawson, 1979; Lawson et al, 1998; Roberts et al, 2002). However, to accurately reproduce these features, it would be necessary to record the sediment concentrations more accurately via a tracing algorithm (e.g., Staniforth & Cote, 1991) or a stretched grid (e.g., Hildes, 2001). Both of these methods are beyond the scope of this thesis for the simple reason that it is imperative that the scientific questions are well-posed, both in specific goal and mathematical formulation. For tracing methods to be accurate, the scalar to be traced would need to be smooth, which is not the case with the melt condition. Stretched grids have the limitation that numerical diffusion can occur and they do not necessarily conserve mass. As a result, I avoid both in favor of an average concentration in the accreted basal layer. 224 Appendix D. Sediment mass balances Another assumption inherent in the derivations presented here is that the along-path ice velocity is neglected. This assumption is consistent with the mass balance equations presented in Chapter 2. D.2.1 Sheet For a subglacial water sheet, the mass of an accreted layer at the glacier base is M b:e= j f fZ {(l-n^il-X^^ + il-r^Xy + r^p^dzdyds, (D.15) J s Jy Jz^ where z b : e is the elevation of the englacial transition between accreted ice and overlying meteoric ice. In equation (D.15), nlp(x, y, z) and X'(x, y, z, t) are spatially varying quantities.15 Through-out the remainder of this exposition, these quantities are local averages, but I avoid the tedium of creating extra notation for these averages. Melting conditions: mass loss When the melt rate is positive and melting is occurring, then there is no change in the sediment concentration in the basal ice. The local form of the rate of change of mass and the change in accreted ice thickness are, respectively, AnW\fdz»« dzb ( . . . i (Bzb'e dzb\ { (1 - ni) (1 - A1) p' + ( l - nL) XY + npp™} | — - — j 1 + n 'p-( l - A O C l - n i y dZb-' dt (1 - raj,) (1 - As) ^ + (1 - nL) Aips + npP« 1 + ( l - A i ) ( l - n i ) p ' m + ¥y (D.16a) m - r - ^ j , (D.16b) where Zb:e = zb:e - zb. Physically, Zb:e is the thickness of the basal ice layer. Equation (D.16b) is equivalent to equation (4.16). If the right hand side of equation (D.16b) is negative and Zb:e is less than an arbitrary height, then dZbe 1 ^ dt = ^ 0b : e - ^ b : e ) , (D.17) where ZQ:G is a minimum thickness and r b : e is a time constant. Equation (D.17) is a penalty function that converges to ZB:E. 1 5 A n argument also exists that p\ p w , and p s are spatially varying quantities. However, at no other point i n the thesis do I suggest this. The densities are taken as constant. 225 Appendix D. Sediment mass balances Freezing condit ions: mass gain When ice is accreting, both the concentration of sediment in the ice and the thickness of the basal ice can change simultaneously. In consequence, two simultaneous equations are necessary: one for the basal ice thickness, and one for the englacial sediment concentration. The local balance from equation (D.15) is r ~\ f)7b:h {(1 - nj,) (1 - A1) j + (1 - nj,) AV S + n ^ } - ^ - + ( l - nj,) (ps - p1) Zb ,d\l ~dt I + < P W m + ¥}. (D.18) ( l _ A i ) ( l - n j , ) p i Because there are two time derivatives in equation (D.18), the local form of the mass of the sediment in the ice is also necessary, d(psX[Zh:e) dt dA1 dt Zb:e A dZhe & dt + "p7 (D.19) Substituting equation (D.19) into equation (D.18) gives dZb:( dt (1 - nj,) p' + nj,p« (1-4) 1 1 + < P W ( l _ A i ) ( l _ n i ) p i (D.20) D.2.2 Channel M b:e J JQ {{l-nL)(l-\i)pl + (l-rti\V + rtppv'}rdrddds, (D.21) For a subglacial semicircular channel, the mass of an accreted layer at the glacier base is CK rRb+Zb:e 'BP where 0 is the azimuthal angle, Zhe is the thickness of the accreted layer, and r dr d6 ds is the differential volume of accreted ice. The quantities np and A1 represent local averages throughout the remainder of this derivation. 226 Appendix D. Sediment mass balances M e l t i n g conditions: mass loss Assuming that the concentration of sediment does not change with time in the accreted layer as ice melts, the local rate of change of mass becomes { (1 " 4) C1 "7 A ' ) P1 + (! - nv) X>PS + n>W} {H\ [(Rh + Zh:eT ~ {RhT] } = 1 + < P W { (1 - nj,) (1 - A') p1 + ( l - nj,) \{ps + nj,pw} ^nZbe ( l - A o a - n i y m + V'1}, (D.22a) dRb dt + 1 + close TT [Rh + Z B : E ) dZb:i dt < p w ( l - A O ( l - n ' ) p ' m + & } , (D.22b) where dRb/dt\c\ose is the closure of the channel resulting from creep. This term is not the full change in radius from the mass balance of the channel because the melt terms are on the right hand side of these equations. Separating for the rate of accreted layer thickness change gives dZbe 1 ( • dt + vr {Rb + Zbe) {(1 - nj,) (1 - A') p1 + ( l - nj,) A'p s + nj,pw} 1 Zb:e dRb I tf1 1 + nip™ ( l - - A i ) ( l - n j , ) p i m > — Rb + Z^e dt (D.23) close Freezing conditions: mass gain The local form of the mass balance with a change in the concentration of sediment in the ice becomes | ( { ( l - n i ) ( l - A V + ( l - n i ) A y + ^ ^ ^ = 1 + < PW ( 1 - A ' ) ( 1 - npp1 dRb (D.24a) {(1 - nj,) (1 - A') p' + (1 - nj,) A y + nj,pw} {*Z b » + TT (Rb + Z B : E ) ^ } + \\2RbZb* + ( > f } {(1 - nj,) (ps - P') f } = -close 1 + ( l_ A i ) ( l -n l , )p> (P. 24b) 227 Appendix D. Sediment mass balances The local balance is of sediment in the ice is .2 / , \ 2V £ { ^ [ ( t f + z - ) ' - ( r f . y -7T 2RbZb:e + (zb«)' OX1 7T : 2Z b:e a #b <9i close + 2 ( V + Z b : e ) 9 Z h : ' dt (D.25a) (D.25b) From the local balance of sediment in the ice, the rate of change of sediment concentration is dX[ dt 2RhZb:e + (Zh:e) I TT/9: 2 1 TT rtS 7b:e dt dZb:' close + ( * b + ^ * ) 2 f r ! ] } -(D.26) Substituting equation (D.26) into equation (D.24b) and rearranging gives dZh:' dt TT {Rh + Zbe] {(1 - raj,) p{ + raj,pw} 1 + ( l - ^ l 1 - ^ ) - 1 ra'pw (1 _ Ai)(l _ n i ) p i Zb:e dRb m pb + Zb-e dt (D.27) close This equation is similar to equation (D.20) with the main difference being that creep closure is included in equation (D.27). 228 Appendix E Semiempirical sediment transport relations van Rijn (1984a,b) developed procedures to apply sediment transport relations. The aim of this appendix is to present the transport relations rather then derive them. Because of this aim, the procedures are presented first; and the supporting equations follow. This format is operationally succinct. The derivations use theoretical and empirical considerations. I refer the reader to the original papers for these derivations. E . l B e d load Procedure 1. Compute the particle parameter D* using equation (E.l). 2. Compute the critical bed-shear velocity according to the Shields criterion (eq. E.7, Fig. E . l , eq. E.5). 3. Compute Chezy coefficient related to grains C using equation (E.8). 4. Compute the effective bed shear velocity using equation (E.3). 5. Compute the transport stage parameter T (E.2). 6. Compute the equilibrium bed load transport qe:b using equation (E.9). Formulation van Rijn advocated two nondimensional numbers to describe sediment motion. The particle parameter describes the nondimensional grain size and is derived from the particle Froude and 229 Appendix E. Semiempirical sediment transport relations Reynolds numbers, where D is a characteristic grain size and v is the kinematic viscosity of water. If this rela-tionship is applied to the entire alluvial bed, D can be replaced with /J50 the median of the grain size distribution. This substitution implies that D50 is an accurate representation of the sediment distribution. The second nondimensional number is transport stage parameter, which characterizes an excess shear stress available to move sediment, j _ (u*) ~ iu*,cr) ^ 2) (u*,cr) The transport stage parameter is a normalized difference of shear stresses, where the bed shear velocity accounts for the lessening of shear stress away from the bed. In (E.2), / 1 \ ,b 9'' u * = u " I ~QJ I > (E.3) is the bed shear velocity related to the grains where C is the grain Chezy coefficient. The maximum value of the bed shear velocity related to grains is the bed shear velocity, = ( E- 4) where ro is the shear along the hydraulic perimenter (eq. 2.24). In equation (E.l), the critical bed shear velocity from the Shields criterion ulcr 6CT = y —5— , (E-5) where 0cr is a critical mobility parameter. The Shields threshold criterion is commonly used to express the critical shear stress necessary for motion, Tcr = (ps-Pv,)gDOcr. (E.6) Shear stresses below this critical value will not cause motion of the sediments. While arguments against this criterion exist (e.g., Yang, 1996, p. 22-23), no relationship is in more common usage 230 Appendix E. Semiempirical sediment transport relations for the initial motion of sediments. An analytic form of the critical mobility parameter is D* < 4.5, 4.5 < £>* < 10, 10 < D* < 18, (E.7) 18 < £>* < 144, 144 < £>*, where I have modified the original form given by van Rijn (1984a) to make the curve smoother. A value of 0.6 for D* > 144 indicates tightly packed, uniform beds. The gray area in Figure E . l represents a lower value because of mixed-sized beds. Values higher than 0.06 are possible provided that the bed is structured in a way that reinforces the grains at the bed (for a review, see Church 2006). 10° o, b 10" 1 0 _ : 0 1 2 3 10 10 10 10 D, Figure E . l : Shields diagram. An illustration of equation (E.7). The gray area represents the range discussed by Church (2006) for a mixture of different sediment sizes. The Chezy coefficient related to the grains C in the bed shear relationship is where D90 is the size of the 90th percentile of the grain size distribution. The equivalent sand roughness is 3D90 as reviewed by van Rijn (1984a). The grain C is a form of the Chezy coefficient parameterized to the original Nikuradse experiments of water flowing over a sand bed (e.g., Henderson, 1966, p. 94). The grain Chezy coefficient can be converted to a grain Darcy-Weisbach friction factor (eq. 2.24) via (fd)~1^2 — C'/\/8g. However, the form of equation (4.33) 231 0.24 D " 1 for 0.14 D~0M for 0.04 D-010 for 0.013 L>°-29 for 0.055 for No motion Appendix E. Semiempirical sediment transport relations is the one that is most compatible with the equations presented thus far, and I opt for that formulation rather than equation (E.8). The bed load flux (eq. 4.19) is equal to the product of the concentration of bed load material (eq. 4.22), bed load layer thickness (eq. 4.21), and bed load velocity (eq. 4.20). In terms of the transport stage parameter and the grain size characteristics, the bed load flux is where semiempirical relations are used for the concentration and bed load height relations. This equation is the relation given in eq. (4.24). The bed load flux is in units of cubic meters per second per meter channel width. (E.9) 232 Appendix E. Semiempirical sediment transport relations E.2 Suspended load van Rijn (1984b) developed a method of calculating the suspended sediment load for particles in the size range 0.1-0.5 mm. The procedure presented here is taken from van Rijn (1984b). Procedure 1. Compute the particle parameter £>* using equation (E.l). 2. Compute the critical bed-shear velocity u*tCr using equations (E.5) and (E.7). 3. Compute the transport stage parameter T (eq. E.2). 4. Compute the reference level 5S using equation (E.17). 5. Compute the reference concentration A b 0 using equation (E.16). 6. Compute the particle size of suspended sediment using equation (E.l9) or another method. 7. Compute the fall velocity of suspended sediment ws using equation (E. l l ) . 8. Compute the 3s factor (eq. E.l3). 9. Compute the overall bed shear velocity u* using equation (E.4). 10. Compute the ip factor using equation (E.15). 11. Compute the suspension parameters Z and Z' using equations (E.12) and (E.14), respectively. 12. Compute the correction factor Ts using equation (E.23). 13. Compute the volumetric concentration of suspended sediment A b using equation (E.22). 14. Compute the suspended load transport qs using equation (E.21) 233 Appendix E. Semiempirical sediment transport relations Formulation The initiation of suspended load is commonly given in terms of the fall velocity of a sediment particle as well as the shear velocity of the water. Based on experimental results, van Rijn (1984b) formulated his own criterion for the initiation of suspended load, for 1 < -D* < 10, 0Aws for 10 < £>*, (E.10) where «*,ers is the critical shear velocity necessary for entrainment into suspension, and ws is the fall velocity of a grain. For large grain sizes, u*, c r s takes the value 0Aws. The sediment fall velocity is ' i p" p™ - 1 J g& "TSV 0.01 1 + i )gD3 I.I for D < 0.1 mm, for O . K D < 1.0 mm, for 1.0 < D. (E. l l ) A suspension parameter Z relates the downward fall of sediments to upward turbulent motions, Z=.-^—, (E.12) where KS is von Karman's constant and 8s describes how the sediments interact with the individ-ual turbulent eddies. Values of 3s greater than unity indicate that the sediments are propelled to the outside of the individual eddies. The justification is that sediments at the outside of the eddies mix more readily in the flow (van Rijn, 1984b). The 0s parameter is „8\ 2 (3s = 1 + 2 ( — u* We for O . K — < 1. it* (E.13) If the fall velocity is greater than the bed shear velocity in equation (E.13), then the particle cannot be in suspension. Particles in the flow will occupy space, damp turbulence, and reduce the particle fall velocity. As a result, a simple correction to the suspension parameter is introduced, Z' = Z + (p, (E.14) 234 Appendix E. Semiempirical sediment transport relations where

5o) 0 - 3 ( l - exp {-0.5T}) (25 - T ) , (E.18) where db is the water layer thickness. If the bed is flat, the bed load layer thickness (eq. 4.21) can be used. The mean size of the suspended load will most likely differ from the mean grain size of the bed. To account for this difference, Ds is a characteristic suspended grain size, Ds = D50 [1 + 0.011 (pa - 1) ( T - 25)]. (E.19) where os is the geometric standard deviation of the sediment distribution. The standard devi-ation is os = 0.5(LWr>50 + D16/D50), (E.20) where Dsn and D\Q are the 84th and 16th percentile of the grain size distribution. Calculating crs from equation (E.20) is excessive for this study. Arbitrarily choosing a pertinent grain size from the fall velocities or the grain size range 0.1-0.5 mm, the which is the range of validity, might yield more insight. 235 Appendix E. Semiempirical sediment transport relations The suspended load can then be calculated as the product of the mean flow velocity u , the water depth d b , and the suspended load concentration A b 0 , qe..s.= uhdh\he,s. (E.21) The volumetric sediment concentration is AL = ^ : 0 , (E.22) where the correction factor Ts is (E.23) 236 Appendix F Numerical formulation Balance law numerical formulation Chapters 2, 4, and 5 are formulated in a similar fashion using the balance laws. Because all of the equations are implemented using MATLAB, the formulation and averaging schemes are an important component of the numerical implementation. These discretizations are made for use with the numerical method of lines (Schiesser, 1991). Use of MATLAB implies that the solution in time is relatively simple to implement using the ODE15S package. This package is specifically for stiff numerical systems and uses a numerical differentiation formula (NDF) that is similar in form to the more common backward differenti-ation formulae (BDF) (Shampine 8z Reichelt 1997; see also Ascher & Petzold 1998, Chap. 5). Both NDF and BDF are multistep methods. In addition, ODE15S uses an adaptive step to advance in time. F . l Sheet and channel numerical formulation Because the sheet and channel are one dimensional in that the quantities of interest vary in s, both have similar discretizations and boundary conditions. I employ a staggered grid and assign each of the variables to either the cell centers or cell edges (Patankar, 1980). F . l . l Gr id and discretization Implementation of the governing equations (eqs. 2.15-2.17, 2.27, and 2.29) requires an along-path coordinate system. Generally, the path is regular in x, but not necessarily regular in s. Figure F. la illustrates a case where As is regularly spaced in x. However, a change in z for any of those points without a change in x would make the grid irregular (Fig. F.lb). The irregular grid appears in all the s-derivatives. 237 Appendix F. Numerical formulation a b Figure F . l : One dimensional grid for subglacial water flow. Centered values reside at points denoted by integer values such as i — 1, i, and i + Staggered values inhabit points at integer values plus or minus one half such as i - TJ and i + \- (a) Illustration of a regular grid in x and s. (b) A similar illustration where the grid in s is irregular. Note that the grid in x is exactly the same for both cases. For a general variable that lies on the centered grid (Table F. l) , first derivatives in s lie on the staggered grid such that dQ ffl dx Ail Ax ds dx ds Ax As Xi+l - Xi J V Si+1 - Si (F.l) For a general variable u> that lies on the staggered grid the first derivative is du = dudx^AuAx M+i-^i-A /Sj+i-Sz-A ds dx ds Ax As l \xi+i - xt_i J \ si+i_ - s{_x J ' While these equations appear to be intuitive from the perspective of the chain rule, the As/Ax terms need to be calculated from the along-path coordinates. If s(x) only, then this term can be set prior to finding the numerical solution. However, if s(x,t), as in the case of an eroding bed, then As/Ax must be calculated at every time step. Often there will be quantities that require values on both the centered and staggered grid. For example, the melt rate is needed on both grids. In such cases, quantities on the centered grid can be averaged according to Vi =H*W . - s,_ i) sheet, (F.3b) Vi =Sh (s{+i — s i - i ) channel, (F.3c) 238 Appendix F. Numerical formulation Table F . l : Subglacial water flow quantities on the centered and staggered grids. Centered Staggered zb(x) Glacier base elevation zs(x) Glacier sediment bed elevation zr(x) Glacier surface elevation Z1(x) Ice thickness p1(x, t) Ice overburden pressure pb(s, t) Subglacial water pressure 4>h(s,t) Hydraulic potential Tb(s, t) Water temperature np(s,t) Ice porosity *4^-(s, £) Gradient of the base -r^j-(s,t) Potential gradient ^^-(s,t) Temperature gradient ub(s,t) Water velocity Channel-specific quantities Rb(s,t) Channel radius Sb(s,t) Cross sectional area m(s, t) Melt rate P w Hydraulic perimeter Pl Melting perimeter Sheet-specific quantities Hb(s,t) Mean sheet thickness rh(s, t) Melt rate per unit width d ^ s ^ ^ ( s ' Flux divergence \^^-(s,t)j Closure rate \ / rinse 239 Appendix F. Numerical formulation where the overline denotes an average and Vi is the volume of the ith cell. Similarly, quantities on the staggered grid are averaged to the centered grid according to LUj V, 1 LO-, 1+V 1 LO • 1 Vi+i2 =Hb+1W(si+l-Si) sheet, channel. (F.4a) (F.4b) (F.4c) Using equations (F.4) and (F.4), it is now possible to discretize the time derivatives for the time derivatives. The discretized forms of the mass balance equations (eqs. 2.15-2.17) are dHh m d t * ( l - 4 : i ) ^ dSb mi dt + dlP dt closeti 25 bsgn(pj-p! nB Px (l -4 i ) dt i sheet, (F.5a) channel, (F.5b) 7 b ^ b I dt i + si+h ~ si-h ( l - n i : i ) p'pw ) sheet, (F.5c) dph dt i rbSb dt + ( l - ^ i ) ^ / m7- > channel. The melt rate and melt rate per unit width become 0.023 ( R ^ ) 4 / 5 Pr2?5K%ATmp..iPi rrii = 0.023 ( R ^ ) 4 / 5 Pr2/5K^ATt HJL mp:i channel, sheet, (F.5d) (F.6a) (F.6b) which are not interesting except that the Reynolds number needs to be averaged such that, channel, •h:i+i P vb (F.7) sheet. P In this equation, assignment of the Reynolds number to the staggered grid limits the number of averaging steps necessary to obtain the melt rate on both the staggered and centered grids. 240 Appendix F. Numerical formulation Physically, the melt rate must match on the staggered and centered grids; however, numer-ically, the melt rate is more important on the centered grid than on the staggered grid. As a result, the melt rate is defined on the centered grid. In the momentum equation, the melt rate terms are small and can be approximated by the melt rate averages. Using the average melt rate ensures that they are smooth on the staggered grid. The melt rate terms do not appear in the steady state momentum equations (eqs. 2.26) because the assumption that the melt rate is small is valid under most flow conditions. When discretized, the momentum equations (2.27) become ,2 / , \ 2 duh ~dt S i + l - Si Pwsu 1 *i+\ ~dt Si+l - Si , 2 "2> -9-'z+1 ,b P * i r 0 . i + i Si+l , 2 sUpw 2rn channel, (F.8a) Si+l - Si -9 ~b _ ~b ^i+1 H Si+l ~ Si sheet, (F.8b) where the first term on the right hand side is upwinded. Similarly, discretizing equations (2.29) yields 2 ~dti 1 2 Si- Si-1 mi nvi nvj c b P cp D i L + . C p A T m p ; j — ( l - n j ^ + n j ^ f o ) ' PX + 5J5 p w c w channel, (F.9a) ~dt, T t - l t , fhi "2 Sn Up Hb L + c* AT r rap:i (l - nj,,) p[ + np:ip™ (ub) ( l - < i ) Px + 2 (r0nb) i sheet. (F.9b) Overall, these equations strike the proper balance for the discretization of internal energy. How-ever, one of the discouraging aspects of this energy discretization is that the dTh/ds terms lie on the staggered grid and are simply moved to the centered grid. This feature guarantees 241 Appendix F. Numerical formulation upwinding but introduces a small amount of error to the second-order accurate scheme in space. Any other scheme would require Ti-i, Ti, T j + i , which arbitrarily bridges the solution from the 2 — 1 to i + 1 nodes. Obviously, when water flow is opposite to the upwinding, the scheme presented in equations (F.9) needs to be modified. F.1.2 Boundary conditions Upstream Two conditions are necessary at the upstream boundary. These conditions are commonly em-ployed as a condition on water pressure and water discharge. Without a surface melt and infiltration condition applied to the model presented here, it is useful to borrow ideas from outburst floods. Ice-dammed lakes release through a conduit where the temperature and configuration of the lake are either known or assumed at the upstream end of the channel (e.g., Clarke, 2003). Instead of a lake, it is possible to create a similar situation by imposing a crevasse at the upstream end of the subglacial water system such that, dZ° Qme\t — Q b m i n \ ~dt = ~AC ' ( F - 1 0 ) where Z c = zc — zs is water depth in the crevasse, zc is the water elevation in the crevasse relative to a fixed datum, Q m e i t represents recharge from a supraglacial water system, Qb is a water flux to the subglacial water system, A° is the representative cross-sectional area normal to z, and dzs/dt is the rate of change of the elevation of the bed below the crevasse (see Fig. F.2). Equation (F.10) neglects any crevasse opening or closing mechanisms via melting/freezing or mechanical opening/closing. Because this equation does not directly account for the boundaries at the surface or base of the crevasse, penalty functions similar to equation (5.27) prevent the crevasse water height from being above the surface of the glacier or below the bed elevation, dZc dt 7 C j . „ c ( Qmelt Q \ Z +r ^ J yc I _c I Qmelt Q \ £ + T \ A° ) (F. l l ) < 0. In each of the cases, r c is a characteristic time scale for the crevasse. These relations can be implemented such that the penalty function is used when the water depth approaches the limits in (F . l l ) . Alternatively, the rate of change of the crevasse water elevation can be prescribed directly without the use of a recharge rate. The advantage of this method is that the crevasse water 242 Appendix F. Numerical formulation height is known for all time. The disadvantage is that the water elevation change is not based upon the outflow from the crevasse, which would be expected. I avoid this representation. Figure F.2: Ice geometry with upstream crevasse. Water flows into the crevasse at a predeter-mined rate, Q m e i t (gray arrow) and flows out of the crevasse into the subglacial water system at a rate QB. From the water depth in the crevasse, the upstream conditions of pressure and temperature are dpc w dZc IF = P 9 ^ dTc _ dpc ~dt ~dl' (F.12a) (F.12b) Equation (F.12a) states that the upstream pressure condition is hydrostatic, and equation (F.12b) states that the upstream temperature evolution is based on the Clausius-Clapeyron equation. The discharge from the crevasse into the subglacial system is defined as Q\ =S\U\ channel, 2 2 OX =H$WU\ sheet, (F.13a) (F.13b) where the implied boundary conditions16 are that dSb/ds = 0 and dHb/ds = 0 for the channel and sheet, respectively. By setting the two mass balance equations (eqs. 2.14 and eqs. 2.16) 1 6 In these equations, numbered subscripts indicate the relative position of the ith location on the grid. For example QI is exactly Qh_ 3 . Near the outflow boundary, i may take a value based on TV. 1 2-i 243 Appendix F. Numerical formulation equal, the boundary conditions on dub/ds become = - ^ { ( ^ - ^ ) m 2 - 2 5 2 b s g n ( p i - ^ ) ( ^ ) l C h a n n e l< ( R 1 4 a ) i n\ i vp 1 pwJ v 9t yclose:j ds | ^2 ds These relationships substitute directly into equations (F.8) for the divergence terms. Downstream Only a pressure boundary condition is necessary downstream. This condition can be either a Dirichlet or Neumann condition. For the discretization employed here, this condition is applied as a Dirichlet condition meaning that the outflow pressure must be known. Commonly, water exits the glacier at atmospheric pressure, which is a natural choice for the pressure. In some situations, water exits the subglacial water system into a lake. In addition, at some of the hydraulic vents at the margin of Matanuska Glacier, water flow is artesian. Both the lake case and the artesian case can be prescribed as Dirichlet boundary conditions. The other boundary condition option is to prescribe the derivative as a Neumann boundary condition. If Neumann boundary conditions are applied, it is possible that open channel flow might occur. However, the equations presented in Chapter 2 are not compatible with open channel flow. Thus, not only does a Dirichlet boundary condition approximate the bound, but this condition also prevents the equations from behaving nonphysically. Open channel flow is observed at many glaciers near the terminus, but its subglacial extent is not clear. Both the ice mass balance (eqs. F.5a and F.5a) and the thermal energy balance (eqs. F.9) require a melt rate at the exit. A Reynolds number appears in this melt rate and requires an appropriate velocity at the final centered cell, b J y 2 UN - Ob N channel, (F.15a) UN = -ik1 sheet> (R15b) N where iV denotes the final grid-centered cell, and Q b is the average discharge per unit width. The discharge between cells is conserved, and as a result, dividing the discharge at the final staggered point by the final cross-section is a reasonable approximation to the velocity there. With the melt rate term defined and because the thermal energy balance is already upwinded 244 Appendix F. Numerical formulation with respect to its divergence term, this equation can be applied at the outflow boundary without complication. Similarly, the ice mass balance can be applied directly to the final centered grid cell. F.1.3 Initialization Because a numerical model integrates over time, the initial conditions are paramount. Initial conditions that are ill-posed will not be corrected by numerical integration. Of the fields that need to be specified, the cross-sectional area 5 b , and the water layer thickness Hb, are the most difficult to ascertain. Essentially, reasonable guesses for these terms are necessary. Both Clarke (2003) and Ng (1998) chose Sh to match known examples. Clarke (2003) also tested initial values of Sh and concluded that the overall sensitivity was limited. It stands to reason that Hb must also be chosen arbitrarily. For the velocity, pressure, and temperature equations, natural initial conditions stem from applying the steady state velocities (eqs. 2.26), setting the water pressure equal to ice overburden pressure, and then applying the Clausius-Clapeyron relation (2.8), respectively. F.2 Sediment transport numerical formulation Governing equations for sediment in the subglacial water system are given by equations (4.11), (4.12), (4.14), and (4.15). There are also the balances of sediment in the water (eqs. 4.10) and the rate of change of the bed elevation (4.7). Other equations that are necessary are the supply from the ice (4.34) as well as the semiempirical sediment transport relations (eqs. 4.19-4.30). Many of these equations have a numerical formulation that parallels the formulation for water flow without sediment transport. As a result, I only touch on the portions of the numerical formulation that are either exclusive to sediment transport or require additional clarification. F.2.1 Gr id and discretization The grid is the same as that of the subglacial water flow (Fig. F. l ) . Additional quantities that fall on this grid are listed in Table F.2. In the four sets of governing equations (eqs. 4.11, 4.12, 4.14, and 4.15), there are necessary modifications of the discretization presented for flow equations presented above because of the inclusion of the sediment concentration A b and the source terms Vf/S, and $ s . In the mass 245 Appendix F. Numerical formulation Table F.2: Sediment transport quantities on the centered and staggered grids. Centered Staggered das ds L . There are no additional discretizations necessary for the balance of internal energy that are not discussed above. The rate of change of the sediment concentration requires two modifications. The definition of sediment flux appears in equation (4.3a) as ) (1 - nj,) p! + nj,pw + A 2 ( l - nj,) ff 1 p w (1 - A*) (1 - nj,) pi + nj,pw + AJ, ( l - nj,) (f 1 + 1 W nvp (1-A ' 2 ) ( l - n p ) p\ 2 dub +25 2 bsgn(p 1-p 2 D) 1 pw nB + "^2bp f3 £ - 1 p w 1 ~a7 - S b u b ^ ff (1 - A*) (1 - nj,) pi + nlpw + A 2 (1 - nj,) ff (F.23a) + + 1 2 1 ^ " (1 - A') (1 - nj,) p ; + nj,pw + A 2 (1 - nj,) ff + I — - 1 1 + 1 W < P m3 2 fdHb\ (ps \(d? „ b ) ) b 9 A b , (1 - A 2) (1 - n p) p1 (F.23b) where equations (F.23a) and (F.23b) are for the channel and sheet, respectively. The sediment divergence terms dXb/ds, dqs/ds, and dqs/ds can be estimated by assuming that the water flows in with an equilibrium concentration of sediment. Water velocities are estimated by assuming that the discharge does not vary significantly over the first cell, such that Ui — Q 3 2_ 52 b 2 channel, sheet. (F.24a) (F.24b) Because sediment concentration is at equilibrium, the first cell cannot evolve its bed elevation via equations (4.7) and (4.4) because * s = 0. However, the spatial derivative of the bed figures 248 Appendix F. Numerical formulation into equation (4.14). As a result, on the first cell dzs/dsi = 0, which is equivalent to saying that the erosion or aggradation rate in the first cell is equal to that of the second cell. F.2.3 Initialization The initialization of the subglacial velocity via equation (2.26) allows calculation of both bed and suspended load via the semiempirical sediment transport relationships. The other important initialization is the concentration of sediment in the ice A1, and it is not clear what this should be set as. In consequence, I take it as a free parameter. F.3 Englac ia l aquifer numerical formulation Equations (5.1), (5.15), and (5.18) form the system of governing equations. These governing equations are subject to the supplementary relations in equations (5.19), (5.23), and (5.25). The conditions at the surface of the aquifer and base of the aquifer given in equations, (5.28-5.32, 5.26, and 5.27). The values of pressure and temperature at the top of the glacier are assumed to be atmospheric and 0°C, respectively. Pressure and temperature values at the base of the aquifer are given by equation (eqs. 2.16b and 2.29b). Note that at this point it must be clear that the aquifer is only over the sheet. This prevents a correction for the radial divergence of water over a channel. F.3.1 Vertical coordinate transformation Because the thickness of the englacial aquifer can evolve through time, and because the aquifer is assumed to be contiguous, a vertical coordinate transformation is useful. This coordinate transformation, and its resulting change in derivatives, is presented in Appendix G. Quantities that vary in the vertical dimension transform via equation (G.la). The irregularly spaced vertical derivatives in z become regular vertical derivatives in £ w for numerical operations. Terms containing £ w map directly back to z via r = (F.25) where £ w takes values in the range (0, 1) dependent on the assignment of z. The vertical derivative of (F.25) is d£w 1 ^ - = —• (F.26) dz Ze v ' 249 Appendix F. Numerical formulation Equations (F.25) and (G.la) can then be used to create the other derivatives according to (F.27) dX Ze \dX K dX where X represents either x or t depending on the situation. When the independent variable is time t or longitudinal position x, the derivatives gain an extra term based on this equation. Because z° and Ze do not vary vertically, the derivatives in equation (F.27) are relatively simple to calculate. Simple calculations such as these are some of the advantages of a vertical coordinate transformation. Each of the four governing equations for flow within the aquifer assumes a form in £ w . Darcy's law (eq. 5.1a) changes slightly with the change of variable, Qx = 1 K%d(j)v Zepwg dt, ' dw 1 dcf>v jy-w pwg l dx Equations (5.15) and (5.18) become z- di dzb ^dZe h £ dx dx dnp 1 (F.28a) (F.28b) (F.29a) 1 (F.29b) (F.29c) di Ze\dt " * dt The divergence of specific discharge V - q w that appears in equations (5.15b) and its transformed equivalent (eq. F.29b) expands to dq] 1 dql d±_+ dZ^ dx dx + 1 dq\ (F.30) dx Z e di \dx ' s dx ) ' Ze di ' The other equations necessary for numerical implementation of the aquifer (eqs. 5.23, 5.25, 5.28, 5.30, 5.31, 5.26, 5.27) are either applied along the boundaries or do not have a dependence on z. As a result, these equations are not transformed via equation (F.25) or its derivatives. F.3.2 Discretization In general, the englacial aquifer grid is positioned with respect to the subglacial sheet. The staggered grid is similarly determined. With two grids, distinction must be made between the 250 Appendix F. Numerical formulation values that reside on the centered grid and those that reside on the staggered grid. Figure F.3 illustrates this grid. In addition, this figure illustrates the difference between the z- and £w-grids. The advantage of the £w-grid is that it is regular and uniform. j+l j+'A j N i-'A i+'A W*M-a-l-B-x-:»«-i3-|--:- • • • • • • I T ^ y -i-1 i i+1 s Figure F.3: Two dimensional grid for englacial water flow. Centered values reside at points denoted with integer values such as i — 1, i, and i + 1 for x and j — 1, j, and j + 1 for £ w . Staggered values inhabit points at integer values plus or minus one half such as i — \ and i + \ for x and j — \ and j + \ for £ w . (a) Illustration of a regular grid in x and z. (b) A similar illustration where the grid in £ w is uniform. For a general variable i}, that lies on the centered grid in x-z coordinates, the vertical, horizontal, and time derivatives can be approximated in x-£ coordinates as, 1 i\j+l m _d£ld£ _ Ail d£ dz 9£ dz ~ A£ dz di} _dn di} d£ dx dx d£, dx Ax „pi+l,j ' X i + i - a ; . Zf A ^ i + i j \ x i + i - X i ' S j di} di} di} d£w (F.31a) (F.31b) (F.31c) (F.31d) dt dt ' d£ dt dt Zf A£ v h J In equations (F.31b) and (F.31d), equation (F.27) has been substituted for the relevant deriva-tive in x or t. In equation (F.31b), the-di~l/d£ term lies on the vertical points of the staggered grid on + \ ) . As a result, this term needs to be averaged to the (i + gridpoints in equation 251 Appendix F. Numerical formulation Table F.3: Englacial water flow quantities on the centered and staggered grids. Centered Staggered nlp(x,£,t) Ice porosity Vertical porosity gradient pw(x,i,t) Water pressure dg-(x,i,t),^(x,i,t) Pressure gradient V q w Flux divergence q^(x,i,t),qj(x,i,t) Specific discharge qw Magnitude of qw T™(x,ti,t) Water temperature Temperature gradient r(x,i,t) Hydraulic potential Potential gradient Gradient of glacier base Aquifer thickness gradient Ft, b Fracture parameters K% Hydraulic conductivity dy, rriy Source terms (F.31c) using a four point average of differences, ~ A O i f An An + + Ail. + An (F.32) where the relevant vertical derivatives in £ are given in equation (F.31a). The four differences require six field points in this equation. Derivatives in t in equation (F.31d) require a simpler stencil for the vertical averages, An (F.33) 1 / — 1 \ If A £ w is constant, then equation (F.33) reduces to Ail, _ ~~ ^ij-l A£wid ~ 2A£W ' Both equations (F.32) and (F.34) utilize the fact that £ w is a uniform grid because additional weightings are not necessary for the averages. At the boundary, either forward or backward differences are necessary. (F.34) 252 Appendix F. Numerical formulation Discretizations for the water fluxes given by equations (F.28) follow as, q •, i • = 1 *H:i,j+\ *,3 (F.35a) 1 d(pv H+l i+\,3 Xi+l Xi 7 ' +2 ' J Xi+l - Xi >. (F.35b) The d(j)w/d^i+i • term can be constructed using equation (F.32). In addition, the hydraulic conductivity in equations (F.35) is defined on the centered grid because equations (5.19) and (5.20) are dependent on np, Fe, and b that reside on the centered grid (Table F.3). Typically, harmonic averages are used for flux perpendicular to layering and arithmetic averages are used for flux parallel to layering. For the anisotropic, nonlayered case, geometric averaging is appro-priate. Without prior knowledge of the internal structure of the glacier, I assign the averages as geometric, KH:ij+± -ynH:i,j*H:i,j+l \ +1 V. ,3 + h KH:i+y ~yKH:i,jKH Vi,jVi+i,j i+h,3 where V is the individual volume of the cell. These volumes are defined by, Vi+l]+h =W (Xi+i - Xi) - tJ) , (F.36a) (F.36b) (F.37a) (F.37b) where W indicates that the volumes have unit width. Combinations of subscripts are quickly realized using equations (F.37) as a template. The mass balances and the thermal balance (eqs. F.29) are all discretized using a second-order finite volume technique. The evolution of porosity becomes, dn: p dt ij m V.l,J dv.ij dnp 1 p> d £ i d Z ? \ d t t where the closure rate (eq. 5.23) has an obvious discretization dt (F.38) dv.ij = A p r i j s S n (Pij - Pij) \Pi,3 'Pi i,3 1 nB (F.39) 253 Appendix F. Numerical formulation In this equation, there are no averages or gradients that require a special formulation. Addi-tionally, the melt rate from equation (5.25) becomes, = 1.27 np:i.jbi.j (F.40) where qfj requires elucidation. Because q w lies-on the staggered grid, qf is an averaged magni-tude such that % 2 { f e + i j + ^ - i , ) 2 + f e , + | + ^ - i ) 2 } (F.41) The final term that is necessary in equation (F.38) is drip/d£ij, which is simply an average of the differences via equation (F.33). The second mass balance equation and the thermal balance are discretized in a similar fashion to equation (F.38), 1 dp™ dt ij dt ij n. dp :p:ij<*w + (1-nr.ij)<* 1 fdz°_ ~dt n. m •v.i,j P cp 'lp:i,3 1 dTw (F.42a) dTw dx ij d v ^ ( L + c;ATmp:h3)+ 9 dt i *j dt ., (F.42b) In equation (F.42a), averaging of dpw/d£,ij to the centered grid is accomplished via equation (F.33). In equation (F.42b), <9Tw/<9£jj is similarly averaged using equation (F.33), and qj- is obtained via equation (F.41). The difficulty in discretizing equation (5.18) is the first term, — ( q w • VT^/n-p, which enlarges because of the coordinate transformation in equation (F.29c). The result of this difficulty is that there are six averages in equation (F.42b) with one of these being repeated. Luckily, these are simple averages where q™ij, dTw/dxij, dzb/dxi, and dZe/dxi are averaged from the horizontal staggered grid to the centered grid. Similarly, dTw/d£ij and q™t j are simple vertical averages from the vertical staggered grid to the centered grid. Implementation of equations (F.38) and (F.42) and their supporting equations uses compu-tational matrix operators as well as MATLAB'S sparse matrix functions. The discretization is formally second-order accurate, but averages may slightly reduce this level of accuracy. 254 Appendix F. Numerical formulation F.3.3 Boundary conditions The upstream and downstream Dirichlet conditions on temperature and pressure outside the aquifer are necessary. Downstream, the pressure is atmospheric, and the temperature is taken to be zero degrees. Upstream, the hydrostatic pressure in the crevasse is the boundary condition. Glaciers undoubtedly exchange water between a more rapidly filling and draining system of crevasses and a slower englacial hydrology, and other models incorporate this feature (e.g., Flowers, 2000). F.3.4 Initialization From the initial water height, the initial temperature profile can be found by assuming that the water is at the equilibrium zero. This assumption requires that the pressure in the water be known. The obvious choice for the pressure is simply setting the englacial water pressure equal to local hydrostatic pressure, pe(x,z,0) = pwg(ze-z). (F.43) For a regular but nonuniform grid in z, an average pressure over the nonuniform portion of the grid is necessary. The technique used to obtain this average pressure is shown in equation (G.5). Temperature is then found from pressures using the Clausius-Clapeyron relationship. Initial porosities and fracture densities must be set arbitrarily. F.4 Sheet closure The discussion of subglacial sediments in Chapter 4 notes that large grains are poorly represented by most sampling techniques. However, in Chapter 3, it is clear that larger grain sizes are important for sheet closure. A reasonable way through this conundrum is to extend the grain size distribution using the fractal distribution (eq. 3.23 or 4.18) along the = 0 (D — 0.001 m, coarse sand) for a total of 17 grain sizes. Extending the scale to $ = —9 255 Appendix F. Numerical formulation (D = 0.512 m, boulder) gives a reasonable distribution for the closure relation (Fig. F.4). While not a perfect match, any reported data for larger grain sizes would fall along or could be easily converted to the $ scale. 101 103 105 107 109 1011 Number of gra ins Figure F.4: Grain size distribution used in the model. Vertical axis on the left size is linear and corresponds to the black line. Right vertical axis is logarithmic and corresponds to the gray line. Numbers along the gray line are the corresponding $ values. Both lines contain the same information. The number of grains corresponds to 100 m 2 but will change depending on the grid cell size. The relations among the points and curves will remain the same, however. When the value of water pressure is lower than the value of ice overburden stress, sheet closure occurs according to equations (3.17). Solution of these equations requires using Newton's method for N +1 variables, where N is the number of grain sizes on the bed exposed to the ice. Performing these iterations for all spatial grid cells becomes computationally expensive within a few time steps. Equations (3.17) describe a steady-state process, and the erosion rules presented in Chapter 4 do not change the distribution of sediments governing sheet closure. These two details can be exploited to make a lookup table whereby the steady state values are computed for a range of reasonable overburden stresses and water pressures prior to integrating the water flow model through time. The closure rates should be insensitive to low erosion rates because sheet closure depends more strongly on larger grains. Erosion will affect the grains smaller than <3> = — 1. Furthermore, the effective length scales (eqs. 3.18-3.22) are insensitive to a change in the grid size provided that the grain size distributions are linked directly to A1 and that partial grains 256 Appendix F. Numerical formulation are used in the grain size distributions. Physically, partial grains are those sediment grains that intersect more than one grid cell or intersect the unit width of flow in the y-direction. Closure velocities for the synthetic longitudinal section are presented in Figure F.5 (see also Fig. A. l ) where the half grain case is used (eqs. 3.20). The numerical solution in Figure F.5a is calculated every 0.002 m from 0.001 m to 0.255 m. Driving stress is calculated at 25 equally spaced points between zero and ten meters of ice equivalent driving stress (equivalent to spacing 0.42 m or 3750 Pa). From 10 meters of ice equivalent driving stress to the stress, there are 75 equally spaced points (equivalent to 0.56 m of ice equivalent driving stress or 5000 Pa). The driving stress grid is intentionally tighter at lower driving stresses because water pressures in subglacial sheets should be near the ice overburden stress. The maximum ice thickness for the synthetic section is 46.86 meters (Table A. l ) . Because ice may accrete throughout a model run, the maximum driving stress in Figures F.5 is 51.3 meters or 110% of the maximum ice thickness in the synthetic section. x 10 5 25 45 5 25 45 5 25 45 o,/(p'g) (m) oe/(p'g) (m) o./fo'g) (m) Figure F.5: Comparison of calculated and interpolated closure schemes. Contours for Figures (a) and (b) are spaced at 1.0 x 10~6 m s _ 1 . (a) Calculated closure velocities. Regelation, mixed, and creep regimes are delineated by white contours and denoted by white R, M , and C, respectively, (b) Closure velocities interpolated from Figure (a). Closure regimes are the same as in (a), (c) Normalized residual between the calculated and interpolated values. Contour intervals are 0.025. 257 Appendix F. Numerical formulation Between Figures F.5a and F.5b, there is little noticeable difference. Figure F.5b shows the values interpolated halfway between the numerical solution points discussed above. Figure F.5c shows the normalized residual. Horizontal lines are contours at 0.00 normalized error. The importance of this figure is not in the values, but in the locations of the maximum errors. The contours capture the seven largest grain sizes. The maximum errors for each of these grain sizes at the radius height is recapitulated in Table F.4. The error results from the change in closure velocity as the number of grains touching the overlying ice changes. For r\ the maximum error occurs at a low driving stress. Obviously, if the water depth meets or exceeds the maximum grain size, there will be a numerical instability as the ice floats off its bed. Thus, error from the largest grain size is not of principal concern because that grain size is problematic for physical reasons ab initio. For all other grain sizes, maximum error occurs at the maximum driving stress. For water depths at or below r% the absolute value of the maximum error is less than 0.03 such that a lookup table appears to be a reasonable approximation to the closure relationship. Furthermore, the errors introduced by interpolation via a lookup table are much smaller than the probable errors introduced by assuming a grain size distribution. Note that when the value of water pressure is greater than the value of the ice overburden stress, the closure relation reverts to creep closure (eq. 3.8). Because creep is not dependent on a Newton-Raphson numerical iteration, creep velocities are computed whenever necessary in the appropriate time step. Table F.4: Maximum interpolation errors for each of the largest grain sizes. Grain Radius (m) Error ae/{p[g) (m)' r i 0.256 1.220 0.21 r2 0.128 0.002 0.21 rz 0.064 0.001 51.3 r 4 0.032 0.0004 51.3 r;> 0.016 0.0003 51.3 re 0.008 0.0001 51.3 r7 0.004 -0.0002 51.3 Negative error values indicate that the interpolated solution is greater than the numerical solution. 258 Appendix F. Numerical formulation F.4.1 Gra in volume Because the largest grains take a portion of the volume of the sheet the maximum thickness of the sheet .ffm a x , is not the same as the average thickness Hh. However, the closure relationship acts on the maximum thickness rather than the average thickness. Because the distribution of sediments that define the closure relationship is known a priori, the relationship between the maximum and average thickness is also known a priori. For the half grain case with a distribution given by Figure F.4, the average thickness is / 18 Hb = -L \ n^AH^ + (A - n-A) ( tfmax - ]T K„ i=l (F.45) where ri% is the fractional area of the bed that is not covered by sediment, 18 is the number of grain sizes, and Ve:i is the effective volume in the ice. The fractional area of the bed not covered with sediment is equivalent to the fractional area of the bed covered by water in Chapter 3. When the effective volume is zero, the volume of a grain in the flow is exactly half the volume of the sphere. The volume of the water sheet is then solely dependent on n™, which forces the thickness of the water sheet to always be finite. This relationship is plotted in Figure F.6. max 0.25 Figure F.6: Average sheet thickness versus maximum sheet thickness. The black line gives the relationship in equation (F.45) where Ve:i is given by equation (3.20a). The dashed gray line gives Hh = i ? m a x . The average sheet thickness is less than the maximum sheet thickness. 259 Appendix F. Numerical formulation The initialization of the sheet thickness and the closure relationship are done simultaneously so that there are no inconsistencies in the volume of either. For example, the amount of obstacle that is within the water layer must not cause the volume of water to be less or more than the average value of Hh indicates. 260 Appendix G Additional discretization characteristics G . l Ver t i ca l coordinate transformation For coordinates that change during the course of a simulation, mapping a grid that changes with time to a set, regular grid is particularly useful. The idea is that a grid may be irregular, but the these irregularly spaced points map to a regular grid. The regular grid remains constant through time while the irregular coordinates evolve with time. Haltiner & Williams (1980, p. 14-15) formulated grid mappings in general terms for partial and ordinary derivatives. For a function f(X,z), where X represents time and/or additional space coordinates and z is its vertical dependence, it is useful to have a set of mathematical rules for transforming the derivatives from z to £ or vice versa. I recapitulate Haltiner & Williams's (1980) as well as parts of Marshall's (1996) derivation for completeness of the change of derivatives, where df(x,z) =df(x,t) df(x,Q dt_ dX dX dt dX' y ' ' df(X,z) =df{X,£) dt dz dt, dz' d2f(x,z) = d jdf(x,t) dt] dz2 dz\ dt dzy _o2f(x,t) (dtV df(x,t) dt d_ (dt dt2 \dzJ dt dz dt \dz d2f(X,t) (dtV • df(X,Q d2t ( c . ~ dt2 \d~z) + dt dz2' [^ j In these equations, / and / have the same value but are be labeled as different variables because of the change in vertical coordinate. (G.lb) G .2 F in i t e volume averages on non-rectilinear meshes For any two-dimensional mesh, the four corners around point (i, j) are defined as — — + (i-\,j + \), and (i + \, j + (Fig. G . l a). Each of these points has a corresponding 261 Appendix G. Additional discretization characteristics x and z that define mesh location. For a mesh with only a vertical transformation, x • 1 • 1 and x- i •, i have the same numerical value; and x4,i - i and x-,i „•, 1 have the same numerical value (Fig. G.l). The numerical values of z, i „• i , _•, i , z,-, i _ _ i , and z4,i , i can all be v 7 1 2'J 2 2'J '2 2 " 2 ' 2 " ~ r 2 different. i-'Aj+'/i i+'AJ+'A i-'Aj-'A i+'AJ-'A i-'AJ+'A i-'Aj-'A i-l i+'AJ-'A i+1 Figure G . l : Four-cornered two-dimensional grid, (a) Illustration of an irregular grid in x and z. (b) A similar illustration where the grid in x is uniform. The finite volume method relies on the smooth scalar quantity f(x, z) being represented by a volume-average of that quantity placed at the centroid of the cell. Mathematically, this integration looks like a two-dimensional average, I fXi+l Tapper f(xij,zij) = — / f(x,z)dzdx, ^ x , _ l S lower (G.2) where the - 5 J + 5 « + 2 > J - 2 + 2 ' J 2 1 2 >•? 2 ' " 2 J " 2 2 J + 2 ~^ 2 2 2 '-^ 2 * 2 2 (G.5) Equation (G.5) can be expanded to include six vertices of an irregular solid volume. In its present form, equation (G.5) is also the weighting per unit glacier width. G.3 Interpolation of bed and surface topography For surface and bed topography data from glaciers without regular sampling, it is important to interpolate the irregular sampling to a regular grid. For the purposes of presentation, the irregular grid has M points and will be denoted with M when the grid is indicated. Similarly, the regular grid N, has N points. The M points of the irregular grid have an ordinate, z\, and an abscissa, X j . From these values, an along path coordinate is, in a simple Pythagorean sense, AsM:i = {(4l : i+l - 4 l : i ) 2 + (xU:i+l ~ Xu-.i^Y , (G.6) where the beginning of the grid is set at sua = 0. The resulting values of SM-.I are M - l SM:i+l = «M:1 + A S M : i - (G^) i=l The differenced regular grid between cell centers is A SM:M ~ sl / p i 0 \ A s n = M - l • ( G ' 8 ) As opposed to the irregular grid, all points on the regular grid have a common spacing, SN:i = SM:1 + (« ~ 1)ASN, . (G.9) but both grids start at SM-.I = SN-.I = 0. 263 Appendix G. Additional discretization characteristics The staggered grid is computed via s m + i = s N : i + ^ for i = 1 : JV - 1. (G.10) If any of the interpolated values for S N fall on the SM values, then X N ; J and z^-.i are assigned the values of XM-A and ZMA, respectively. If a new grid point, s^-.j falls between two old grid points, SM-.J and S M - . J + I , then a linear interpolation scheme is used to find the new values for x^-j and ~s XN:j = {SN:j ~ S M A ) 1" Xy[:i, ( G . l l a J SMA+l ~ SMA 4 - 7 = (sN:j - SM:i) ~ 1" 4l : r (G.llb) SM:i+l - SM:i Similarly, the surface topography of the glacier is interpolated with a linear scheme such that -ZJYI is interpolated to z^. The thickness of the glacier is thus Z\ = z^:i — z^:i. Terms on the staggered grid x N . J + i and are created in much the same way as s N : J - + i in equation (G.10). From these relations, dz/dsN.+ i can be created. 264