Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Hydraulic properties of subglacial sediment determined from the mechanical response of water-filled boreholes Waddington, Brian Stewart 1993

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

Item Metadata


831-ubc_1993_fall_waddington_brian.pdf [ 5.12MB ]
JSON: 831-1.0052960.json
JSON-LD: 831-1.0052960-ld.json
RDF/XML (Pretty): 831-1.0052960-rdf.xml
RDF/JSON: 831-1.0052960-rdf.json
Turtle: 831-1.0052960-turtle.txt
N-Triples: 831-1.0052960-rdf-ntriples.txt
Original Record: 831-1.0052960-source.json
Full Text

Full Text

HYDRAULIC PROPERTIES OF SUBGLACIAL SEDIMENTDETERMINED FROM THE MECHANICAL RESPONSE OFWATER-FILLED BOREHOLESByBrian Stewart WaddingtonB.A.Sc. (Mechanical Engineering) University of British ColumbiaA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF GEOPHYSICS AND ASTRONOMYWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAJune 1993© Brian Stewart Waddington, 1993In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)Department of Geopky5iG5 and Asir°notttYThe University of British ColumbiaVancouver, CanadaDate  JUne ae., 19F3DE-6 (2/88)AbstractThe hydraulic properties of basal sediments must be known if the surge mechanism ofsoft-bedded glaciers is to be understood. Freezing of water-filled boreholes drives wa-ter into the subglacial bed and the associated pressure effects give information aboutsubglacial hydraulic properties. A numerical model describing the mechanical responseof an unconnected borehole and the bed beneath it to this freezing forcing was devel-oped, using a non-linear transient viscoelastic ice rheology and an approximate modelof top-down freezing. The resulting system of equations was solved using the methodof lines. Results agree well with analytic solutions, if parameters are correctly chosen.Forward modelling of pressure records from three 1992 boreholes and three from otheryears indicates that the till underlying Trapridge Glacier has a hydraulic conductivity of1.3 x 10-9-1.35 x 10' m s-1, in agreement with previous indirect estimates. The modelwas also used to investigate the response of a borehole to sudden pressure changes. Theresponse is very fast compared to pressure sensor sampling rates: thus the "true" basalsignal is essentially unaffected by the presence of the borehole, except during the initialfreeze-in. The decay events seen in Trapridge Glacier pressure records are caused by theresponse of the basal drainage system, not the kprehole.Table of ContentsAbstractList of Tables^ viList of Figures^ viiiList of Symbols^ xiiAcknowledgements^ xvii1 Introduction 12 Ice Deformation Model Development 72.1 Introduction ^ 72.2 Governing Equations ^ 82.2.1^Rheology of Shyam Sunder and Wu ^ 92.2.2^Balance Equations ^ 152.2.3^Definition of Strain and Strain Rate^ 152.2.4^Compatibility of Strains ^ 152.3 Application of Limiting Assumptions 162.3.1^Isotropic Ice ^ 162.3.2^Very Slow Flow 232.3.3^Radial Symmetry ^ 232.3.4^Small Vertical Gradients (Except Hydrostatic Pressure) ^ 26ili32.3.5^Infinitesimal Strains ^2.4^Initial and Boundary Conditions 2.4.1^Initial Conditions ^2.4.2^Boundary Conditions Bed Equations, Forcings, and Model Synthesis3.1^Basal Groundwater Flow ^3.2^Borehole Equations 3.3^Forcings ^3.3.1^Freezing ^3.3.2^Sudden Pressurization ^3.3.3^Pressure Step Change 3.4^Summary of Coupled Initial Value Problems ^3.5^Method of Lines Solution ^3.5.1^Spatial Coordinate Transformation ^3.5.2^Finite Differencing ^3.6^Computer Programming of Solution ^313434353738424646474848495055584 Model Testing 604.1 Introduction ^ 604.2 Ice Deformation Model ^ 614.2.1^Glen-Law Viscous Ice 614.2.2^Linear Elastic Ice ^ 714.3 Basal Groundwater Flow Model ^ 764.4 Coupled Model ^ 774.4.1^Sealed Borehole with Glen-Law Ice ^ 804.4.2^Effect of Viscoelastic and Transient Viscoelastic Ice ^ 82iv4.4.3 Ice—Water—Bed Model ^  854.5 Conclusions ^  875 Model Results^ 885.1 Introduction  885.2 Freezing Forcing L(t) Obtained from Thermistor Records ^ 895.3 Acceptable Parameter Ranges ^  965.3.1 Geometry ^  965.3.2 Ice Properties  ^995.3.3 Bed Properties ^  1025.4 Ice Parameters from Blind Hole Fits ^  1025.5 Bed Parameters from Unconnected Hole Fits ^  1105.5.1 Effect of Ice Properties ^  1145.5.2 Effect of Freezing Forcing  1175.5.3 Other Boreholes ^  1205.5.4 Summary of Bed Properties ^  1205.6 Response to Sudden Pressure Forcing  1246 Summary^ 128References^ 131Appendices^ 135A Summary of Model Equations^ 135List of Tables2.1 Notation and physical definitions of Shyam Sunder and Wu rheology. . . . 104.1 Standard rheological parameters used for model testing ^ 624.2 Optimum ddriv3 error control parameters ^ 634.3 Dependence of Glen-Law strain on coordinate transformation . 704.4 Optimum parameters for Glen-Law ice ^ 714.5 Dependence of linear elastic strain on coordinate transformation ^ 764.6 Optimum parameters for coupled problem ^ 875.1 Freeze-in times for 92T01 and 92T02 935.2 Transient parameters from Shyam Sunder and Wu [1990] ^ 1005.3 Grain-size-independent transient parameters ^ 1015.4 Computed transient parameters ^ 1015.5 Summary of input parameter ranges 1045.6 Bed parameters from fit of 921124 freeze-in ^ 1145.7 Bed parameters from fit of 921124 freeze-in, with various ice properties 1155.8 Bed parameters from fit of 921124 freeze-in, with various freezing forcings 1185.9 Bed parameters from fits of various freeze-ins ^ 121viList of Figures1.1 Types of Trapridge Glacier boreholes and freeze-in pressure records . . .1.2 Response to freezing forcing  ^43.1 Bed geometry and coordinate system ^  393.2 Ice finite-difference grid  ^553.3 Bed finite-difference grid  ^564.1 Effect of maximum radius on strain in Glen-Law ice ^  664.2 Glen-Law stress distributions for various coordinate transformations .^684.3 Effect of number of nodes on strain in Glen-Law ice ^  694.4 Effect of maximum radius on strain in linear elastic ice  734.5 Elastic stress distributions for various coordinate transformations .^744.6 Effect of number of nodes on strain in elastic ice ^  754.7 Aquifer response for various numbers of nodes  784.8 Aquifer response at various times ^  794.9 Sudden pressurization of a sealed borehole surrounded by Glen-Law ice^834.10 Effect of ice rheology on pressure decay in a sealed borehole ^ 844.11 Pressure decay for various ice/bed combinations ^  865.1 92H24 freeze-in temperature record (92T01)  915.2 921126 freeze-in temperature record (92T02) ^  925.3 Multi-year pressure record for 90P02 ^  945.4 921124 length versus time: observed and fitted  ^955.5 Parabolic borehole ^  97vii5.6 Freeze-in pressure record for 92P02 (in borehole 921126) ^ 1035.7 Initial part of 921126 freeze-in (from 92P02) ^  1055.8 Freeze-in pressure records considered in this thesis  1065.9 Best fit to 921126 freeze-in pressure record ^  1085.10 Best fit to 921145 freeze-in pressure record  1095.11 Best fit to 921126 freeze-in pressure record, with slower-decaying L(t)^1105.12 Freeze-in pressure record of 921124 (from 91P12) ^  1115.13 Best fits to 921124 freeze-in pressure record, with various r, ^ 1135.14 Comparison of freeze-in with softest, stiffest, and rigid ice  1155.15 Best fit to 921124 freeze-in pressure record with softest and rigid ice . .^1165.16 Best fit to 921124 freeze-in pressure record, with weaker freezing forcing . 1185.17 Best fits to 921124 freeze-in, with rb = 0.025m and various rc ^ 1195.18 Best fit to 891150 freeze-in pressure record ^  1215.19 Best fit to 901154 freeze-in pressure record  1225.20 Best fit to 911146 freeze-in pressure record ^  1235.21 Pressure record of 90P07 ^  1245.22 Sudden pressurization at t = 1 day of hole 921124 ^ 1265.23 Sudden pressurization at t = 1 day of hole 921124 impermeable bed. . . ^ 127List of Symbolsa()^a5^coefficients of L(t) expression (m)A Glen Law flow law parameter (Pa's1/3)A^kinematic hardening parametergrain-size independent kinematic hardening parameter.b1 • • • b5^exponents of L(t) expression (s-1)scalar drag stress associated with isotropic hardeningB0^initial scalar drag stressgrain-size independent initial scalar drag stressc1^c6^anisotropic parametersc* normallizing coefficient, c* = c2c3 c3c1 c1c2ice grain size (mm)hydraulic diffusivity (m2s1)Dijki^orthotropic elastic stiffness tensor (Pa)Young's modulus (Pa)60^unit vectorfunction (of which finite difference is taken)f.^body force vector (N)gravitational acceleration, g = 9.806 m s'Gipet^fourth-order stress transformation tensorhydraulic head (m)ho^initial hydraulic head (m)Hijkl^fourth-order strain-rate transformation tensor1Xtemperature-independent isotropic hardening parameterH'^grain-size independent temperature-independent isotropic hardeningparameterhydraulic conductivity (m s")length of borehole (m)mass (kg)mo^initial mass of water in borehole (kg)maq^mass of water expelled from borehole into bed (kg)Maq,ramp^mass of water expelled into bed from 0 < t < tramp (kg)Mf^mass of water frozen to top of borehole (kg)Mf,ramp^mass of water frozen to top of borehole from 0 < t < tramp (kg)Mramp^mass of water in borehole at t = tramp (kg)mw^mass of water in borehole (kg)number of ice nodesMa^number of bed nodesporosity of basal sedimentflow-law exponenthydrostatic pressure (Pa)Po^initial hydrostatic pressure (Pa)Pb borehole water pressure (Pa)Pb^water pressure in basal cavity (Pa)p f or ce^forcing pressure (Pa)qi specific discharge (ms')creep activation energy of ice (J mol-1)Q H^creep activation energy of ice, for T> —10 ° C (J morl)QL^creep activation energy of ice, for T < —10°C (J mol-1)Q.^mass flow rate of water (kg3s-1)Qw^volume flow rate of water (m3s-1)radial spatial coordinate (m)rb^borehole radius (m)initial radius of borehole (m)radius of basal cavity (m)?max^radius of outer ice boundary node (m)transformed radial coordinate in iceRa^transformed radial coordinate in bedRgaa^gas constant, Rgas = 8.314 J mo1-1K-1elastic back stress tensor (Pa)applied loading tensor, pseudo-deviatoric stress tensor (Pa)pseudo-deviatoric reduced stress tensor (Pa)back stress measure tensor (for isotropic materials this reduces to thedeviatoric stress tensor) (Pa)time (s)tramp^ramp-up time for pressure forcing (s)temperature (K)T263^constant temperature, T263 = 263.12 Kui displacement vector (m)vi^velocity vector (m s-1)V^viscous stress factor (Pa)Vo^temperature-invariant viscous stress factor (Pa)VO,H^temperature-invariant viscous stress factor applicable for temperaturesxiT> —10°C (Pa)V^volume (m3)Vb^volume of borehole (m3)volume of basal cavity (m3)Vf^volume of ice frozen to top of borehole (m3)Vw^total volume of borehole and basal cavity (m3)weighting factor of transformed spatial derivativesxi^spatial coordinates (m)X^power of radial coordinate transformationvertical spatial coordinate, z = x3 (m)compressibility of basal sediment matrix (Pa-1)compressibility of water (Pa-1)7^engineering shear strain (s-1)Kronecker deltatotal strain tensor (s-1)Eli^elastic strain tensor (s-1)steady-state viscous strain tensor (s-1)Et..^transient viscous strain tensor (s-1)constant which scales the equivalent stresscylindrical angle coordinate and spherical angle coordinate (rad)steady-state viscous softness parameter (Pa-1)At^transient viscous softness parameter (Pa-1)density (kg m-3)pw^density of water (kg m-3)Po^density of water at pressure Po (kg m-3)xiilor^density of newly frozen ice (kg m-3)cri;^stress tensor (Pa)stress difference tensor, reduced stress tensor (Pa)creq^equivalent uniaxial stress (Pa)Cr cl,eq^equivalent stress difference (Pa)Crij^deviatoric stress tensor (Pa)E'2 second invariant of deviatoric stress tensor (Pa)spherical angle coordinate (rad)(del dt)0^reference strain rate (s-1)(c/e/dt)d,eq^equivalent uniaxial strain rate (s-1)d dt^material time derivativea/Ot^spatial time derivativeAcknowledgementsMy deepest thanks are due to my supervisor, G.K.C. Clarke; his inspiration, guidance,and encouragement were essential to the completion of this work. I also thank all thoseinvolved in the Trapridge Glacier field program, especially Dr. E.W. Blake, U.H. Fischer,S.J. Marshall, Dr. T. Murray, and Dr. D.B. Stone, for their part in collecting subglacialpressure data, and for their insightful suggestions. I also thank my family and friends fortheir support during the inevitable ups and downs of this project. This work has beensupported by grants from the Natural Sciences and Engineering Research Council andthe University of British Columbia.xivChapter 1IntroductionThe question of how glacier surges work is a major problem in glaciology. It is now wellestablished that subglacial hydrology plays a pivotal role — surges result from acceleratedbasal sliding induced by sustained high water pressure [Raymond, 1987]. How this highwater pressure comes about and how it is sustained are still unclear. The surge ofVariegated Glacier has been modelled in terms of a linked cavity drainage system betweenthe ice and an impermeable rock bed [Kamb, 1987], but not all surge-type glaciers overliea hard bed. Some surge-type glaciers are underlain by a bed of soft, unconsolidated,permeable sediment [Clarke et a/., 1984], as is the fast-flowing Ice Stream B in WesternAntarctica [Engelhardt et al., 1990].Previous investigations have established that Trapridge Glacier, a small sub-polarsurge-type glacier currently in its quiescent phase, is underlain by permeable, deformingsediment [Clarke et al., 1984; Blake, 1992; Blake eta., 1992]. Stone [1993] suggests thatTrapridge Glacier rests on a layer of relatively low-permeability till, perhaps one metrethick, which overlies till and alluvial deposits of much higher permeability. Smart andClarke [submitted] propose that in the winter, and at other times of low surface-water in-put, water generated at the bed by frictional and geothermal heat percolates through thelow-permeability layer, into the higher-permeability layers beneath (from which it easilyescapes by groundwater flow). There also exists during the summer a high-permeability(although spatially discontinuous) drainage system immediately under ice, above the low1Chapter 1. Introduction^ 2permeability-layer [Stone, 1993]. This drainage system has been investigated using bore-hole response tests carried out in connected boreholes [Stone, 1993; Stone and Clarke,1992]. (In a connected borehole the water level drops when the drill reaches the bed,indicating that the drill has intersected some sort of subglacial drainage system. Typicalfreeze-in pressure records are rather different for the different types of holes (Fig. 1.1).Connected holes immediately assume the drainage system pressure (which often displaysdiurnal oscillations as shown in Fig. 1.1). Unconnected holes experience a pressure risestarting 6-12 h after drilling, when the hole freezes over. Blind holes experience a muchmore rapid pressure rise with repeated sudden drops, believed to be caused by hydraulicfracturing of the ice.)Knowledge of the hydraulic properties of the low-permeability layer is important tounderstanding the surge mechanism of Trapridge Glacier because this layer controls waterescape from the bed during much of the year. (Interestingly, a number of observed surges,such as that of Variegated Glacier, have started during the winter when water input islow [Kamb et al., 1985]. This is precisely the time of year when the low-permeabilitylayer controls the Trapridge drainage.) The only existing in-situ estimate of the hydraulicconductivity of this layer is an upper bound of 3.6 x 10-7 m s-1 obtained by Smart andClarke [submitted] from consideration of the minimum observable water level drop inan unconnected borehole. Clarke [1987] has estimated the hydraulic conductivity tobe 1.27 x 10-8 m s-1, based on the granulometry of samples taken from the bed, andassumptions about porosity. Tavi Murray (unpublished data, 1992) obtained a valueof 2.2 x 10-8 m s-1 from consolidation tests of forefield till, which is presumed to besimilar to the basal material. None of these estimates is completely satisfactory and thisfact motivates my attempt to establish a more direct, in-situ measurement of subglacialhydraulic properties.In this thesis I examine the freeze-in of unconnected boreholes and use this analysis— I100500< 2000ccD 100cncncc 300200100a10^1^(77// cChapter 1. Introduction^ 30^ 1TIME (DAYS)Fig. 1.1: Types of Trapridge Glacier boreholes and freeze-in pressure records. Boreholesare classified as: (a) connected if the water level drops when the drill reaches the bed,indicating that the hole intersects a highly-transmissive basal drainage system, (b) un-connected if the hole reaches the bed but the water level does not immediately drop,indicating that the hole bottoms in less permeable material, and (c) blind if the holedoes not reach the bed. Note the differing time and pressure scales.Chapter 1. Introduction^ 4aFig. 1.2: Response to freezing forcing. These cartoons illustrate the two ways that excesswater displaced by freezing can be accommodated: (a) by ice deformation, and (b) bywater flow into the bed. In blind holes only ice deformation operates. In connectedholes water flow into the bed dominates. It is unclear which mechanism dominates inunconnected estimate the hydraulic conductivity of this layer. When unconnected boreholes freezeshut, water in the borehole is compressed and its pressure rises. The magnitude of therise depends on the hydraulic properties of the bed, as well as the rate of freezing of theborehole, and the rate of ice deformation around the pressurized hole (Fig. 1.2).I construct a numerical model of the coupled ice—borehole—bed system which describesthe mechanical aspects of a freezing borehole, including the ice deformation and ground-water flow of water into the bed, and an approximate freezing forcing obtained fromtemperature measurements made as boreholes freeze. The inclusion of a proper thermalmodel, such as that of Jarvis and Clarke [1974] or Humphrey and Echelmeyer [1990], in-stead of the approximate freezing would improve the model, but was precluded by limitedtime. I make many other simplifications to render the problem more tractable, includingusing only a single layer of deforming ice, neglecting non-hydrostatic background stress,Chapter I. Introduction^ 5and describing water flow into the bed in terms of Darcian flow in a uniform half-space(neglecting bed deformation and bed heterogeneity). A complexity that I retain is theuse of a transient viscoelastic rheology [Shyam Sunder and Wu, [198913] to describe theice; the inherently transient nature of the problem precludes the use of a simpler Glen-Law description. Compressible water is used in the borehole so that blind holes can bemodelled. The resulting system of non-linear partial differential and algebraic equationsis solved using the method of lines [e.g., Schiesser, 1991].I test the model to check its accuracy in some test cases, and then use it to estimateice and bed properties, by iterative forward modelling until observed freeze-in pressurerecords are satisfactorily fit. The sensitivity of the results to uncertainties in inputs suchas freezing rate, ice properties, and basal-cavity size (excavated by the hot water drillat the base of the borehole) was investigated. Hydraulic conductivity proved to be mostsensitive to the freezing rate and the basal-cavity size, and essentially insensitive to theice properties. I estimate the hydraulic conductivity of the low-permeability layer to be1.3 x 10-9 m s' to 1.3 x 10-9 m s-1. This range is consistent with the previous estimates.No satisfactory estimate of bed compressibility was obtained because model results forthis parameter were so highly variable.Boreholes in Trapridge Glacier freeze from the top down, since the coldest ice isnear the surface [Clarke and Blake, 1991]. Pressure sensors located 0.5 m above thebed continue to indicate water pressure for 1-2 yr after insertion, before they becomeisolated by freezing. Hence above every pressure sensor there is a water column whichvaries in length over time. The question arises of whether the response of this boreholecould influence the pressure observed in unconnected holes, essentially filtering the "true"basal pressure signal. This question must be resolved before the pressure signals inunconnected holes (and those holes that, during winter, become unconnected) can beconfidently interpreted.Chapter 1. Introduction^ 6To settle this issue I use the model to investigate the effect of the borehole responseon the pressure observed in unconnected boreholes after the initial freeze-in period. Themodelling results indicate that the borehole has a negligible effect on the observed pres-sure (except during the initial period of rapid freezing).Chapter 2Ice Deformation Model Development2.1 IntroductionTo predict the pressure response of an unconnected borehole as it freezes shut after hot-water drilling, the mechanical response of the ice—bed—borehole system must be modelledtogether with the freezing forcing itself. The model must describe the deformation ofice in response to elevated pressure, as well as the flow of water from the bottom ofthe borehole into the bed. In this chapter the ice deformation model is developed. InChapter 3 the bed and borehole equations are developed and the numerical solution tothe coupled problem is described.The problem in this thesis is clearly transient, so an appropriate transient ice rheologymust be used. This rheology must describe the response to time-varying pressure at theborehole wall, including pressure which first rises, then falls. Therefore the strain recoverybehaviour of ice [Duval, 1978] must be properly described. The simple Andrade-Lawexponential form [e.g. Paterson, 1981, p. 26] adequately describes the behaviour of iceunder a suddenly applied steady load (i.e., in a creep test), but does not model strainrecovery. The viscoelastic fluid relation of Morland and Spring [1981] and the viscoelasticsolid relation of Spring and Morland [1982] model only the response to constant load orconstant strain rate. A number of formulations seek to model full response, includingstrain recovery [Sinha, 1978; LeGac and Duval, 1980; Ashby and Duval, 1985; ShyamSunder and Wu, 1989a; Shyam Sunder and Wu, 1989b]. Of these the physically based7Chapter 2. Ice Deformation Model Development^ 8phenomenological model of Shyam Sunder and Wu [1989a, 1989b] appears to best predicttest results [Shyam Sunder and Wu, 1990]. For the present work I accept the model ofShyam Sunder and Wu [1989b].It must be noted that none of the foregoing transient rheologies describes ice deforma-tion at large strains (tertiary creep), where recrystallization is important. Instead theydescribe the behaviour from zero strain up to the time when the minimum strain rateis attained (at between 0.5% and 2% [Budd and Jacka, 1989]). Thus the model that Idevelop is valid only for small strains.The freezing-borehole problem is inherently three-dimensional, but I simplify it tomake it tractable. I assume the background ice stress to be purely hydrostatic. I considerthe borehole to be a long cylindrical hole, so that end effects can be ignored and icedeformation assumed radial. Furthermore variations in ice properties (e.g., those causedby temperature change) and stress state with depth (the hydrostatic pressure in theborehole changes more quickly with depth than the background hydrostatic pressure inthe ice due to the different densities of ice and water) are ignored, so that the deformationcan be considered constant with depth. This allows the ice deformation problem to bereduced to a single layer of radially symmetric flow, a one-dimensional problem.The balance of this chapter outlines the derivation of the ice deformation model. Thedevelopment is lengthy and complex, an inescapable consequence of the complexity of thephysical model. Readers who have little appetite for mathematical details are encouragedto proceed directly to Chapter 4.2.2 Governing EquationsThe governing equations fall into three basic categories: balance equations (mass andmomentum), equations defining strain, and the theological equations. These are the^Chapter 2. Ice Deformation Model Development^ 9starting point for the derivation of the ice deformation model.2.2.1 Rheology of Shyam Sunder and WuThe Shyam Sunder and Wu (SS86W) rheology can be described as a transient, nonlinear,viscoelastic rheology with internal state variables. Shyam Sunder and Wu [1989b] presenttheir rheology, in rate form, for a viscously incompressible orthotropic viscoelastic ice.Since I will treat isotropic ice only, some conversion work is necessary.This rheology was developed with the following assumptions: orthotropic polycrys-talline ice [Shyam Sunder and Wu, 1989b, p. 224], only primary and secondary creepare described (damage mechanics, no recrystallization) [Shyam Sunder and Wu, 1989b,p. 224], and strains and rotations are small [Shyam Sunder and Wu, 1989b, p. 225]. Theexclusion of recrystallization effects, and hence the inability to describe tertiary creep,make this rheology unsuitable for modelling large-strain deformations. Because I amprimarily interested in small strains, this limitation is acceptable. No existing rheologiescan model the full range of primary, secondary, and tertiary creep.The SSSEW rheology is summarized on p. 231 [Shyam Sunder and Wu, 1989b], andrestated, in my notation (which differs slightly from theirs), below. Notation and physicaldefinitions are given in Table 2.1.Strain Rate Decompositionclei;^de7i^de:idt^dt^dt^dtEquations of State (in rate form)dui;^delidt^Mid dtdafij =- dui;^dRi;dt^dl^dt(2.1)(2.2)(2.3)Chapter 2. Ice Deformation Model Development^ 10SS&W thesis parameter descriptionA A "kinematic hardening parameter" (my terminology)al ... a6 c1...c6 anisotropic parametersB B scalar drag stress associated with isotropic hardeningD Dijkl orthotropic elastic stiffness tensorE E Young's modulusG Giiki fourth-order stress transformation tensorH Hijki fourth order strain-rate transformation tensork .i-/ temperature-independent "isotropic hardening parameter"(my terminology)R R i ; elastic back stress tensorS * S.i'i applied loading tensor, pseudo-deviatoric stress tensorSa* S'ciii pseudo-deviatoric reduced stress tensor (for isotropic materialsthis reduces to the deviatoric stress tensor)SR* Sliij back stress measure tensor0 77 constant which scales the equivalent stressE Ei; total strain tensorEeEt<17ct,iielastic strain tensortransient viscous strain tensorc, c:i steady-state viscous strain tensorA A, "steady-state viscous softness parameter" (my terminology)Ad At "transient viscous softness parameter" (my terminology)cr a-ii stress tensorad ,ij stress difference tensor, reduced stress tensorTable 2.1: Notation and physical definitions of Shyam Sunder and Wu rheology.Chapter 2. Ice Deformation Model Development^ 11whereEvolution EquationsdR•;^deki^ — AE 'kJ^•dt^3^3 dt(2.4)deYildt(del dt)0de:j I dt(de I dt)0dBdtAtSliHE ( de)crd,eq kdt ) t.egThe equations of state (2.2) and (2.3) may also be expressed in absolute terms== Dijki67;cri;^ (2.8)ach =^Ria-^ (2.9)It will prove convenient to rearrange the strain decomposition. Letting the viscous strainbegivesand^v ^t= fij ^dey-^de5^def.1 1 +t ^dt^dt==E75 "Equations (2.11) and (2.12) replace (2.1).The preceding equations and those appearing throughout this thesis are in tensornotation. Unfortunately many of the same variables are defined in Shyam Sunder andWu [198913] in "engineering" vector notation. In this convention the symmetric stressand strain tensors are represented as six-element vectorsjr. = Cryy) Crzz, Crxy, Cryz, CrzxlE^ T[eyy, Eyy, CZZ 7xy, 7yz 7zxi(2.13)(2.14)Chapter 2. Ice Deformation Model Development^ 12where -y = 2E are the engineering shear strains, equal to twice the conventional shearstrains. Time derivatives of shear strains are similarly defined. It is also apparent,although not explicitly stated, that the various pseudo-deviatoric stress vectors, such as^:5* =^S:z, 2 S:y , 2 Sy'.. , 2S:r. I T.^ (2.15)^are similarly defined. Likewise^and -cid are defined similarly to^e^es, et, §*, St7i,and ,§1 are defined similarly to e. In the engineering notation used by SS&W fourth-ranktransformation tensors are replaced by 6 x 6 symmetric matrices.The term Hoc/ in the elastic back stress equation (2.2) is defined as 3 (c1 1- 4) c;^--3c1c2c3^--3C1C2C3c*2C*2 C*23 (c1 c3) 4^--3c1c2c3c*2^c*23 (ci^c3)c*20H 1 (2.16)SYM.(2c4)1(24)1(2c6)[Shyam Sunder and Wu, 1989b, p. 229], where c* c2c3 + c3c1 + c1c2. The term (de/dt)0may be set equal to unity without loss of generality [Shyam Sunder and Wu, 1989b,p. 228].The following terms of the viscous strain-rate equation (2.5) are defined as [ShyamSunder and Wu, 1989b, p 228]= (-1 )N N —1v^eq (2.17)2^3 (idCreq TG (2.18)Chapter 2. Ice Deformation Model Development=1^ 1(ci +c3 )^— -3 C1^— -3 C1Ikci^c3)^--3c1 013(2.19)1G =5 (C1 + C3)(2.20)2c4SYM. 2c52c6where= Vo exp (Q/NRg„,T) (2.21)77 =^c2) (2.22)A flow law exponent value of N = 3 is generally accepted [e.g., Shyam Sunder and Wu,198913] and is used in this thesis. With this definition of n, creq cryy when the stateof stress is uniaxial (i.e., a = (0, cryy, 0, 0, 0, 0)). This will result in creq being the usualequivalent uniaxial stress. Equation (2.21) is the commonly used Arrhenius Equation forthe temperature dependence of the viscous stress factor V (which is loosely equivalent tothe flow law parameter A in the Glen flow law), and is generally accepted to be valid fortemperatures below —10°C [e.g., Shyam Sunder and Wu, 19894 At higher temperaturethe viscous stress factor changes more rapidly than (2.21) predicts, due to the presenceof water at grain boundaries. Hence a different expression for V must be used above—10° C. I follow the approach of Paterson [1981, p. 39], and use the Arrhenius relation,but with a higher activation energy Q H (letting Q in (2.21) become QL). The high- andlow-temperature expressions must, however, give the same value at V at —10° C, whichwould not be true if the Arrhenius relation with different Q but the same Vo were used.This difficulty is overcome by computing a new Vo,H from the V(-10 ° C) found usingChapter 2. Ice Deformation Model Development^ 14the low-temperature expression. HenceV(-10 ° C) V0,11 exP (QH/NR9as71263) — Vo exP (QL/NRga5T283)where T263 = 263.12 K. ThusVO,H = VO exP [QL/NRgasT263] exp [—QH/NRgasT2631= V0 eXp [(QL QH)/NRgasT263] •^ (2.23)Therefore, for T > —10° C,V(T) = VO,H eXP [QH/NR9a87]Vo exP [(QL — QH) /NR9a.T263] exP [QH/NRga,T] .^(2.24)The following terms of the transient viscous strain-rate equation (2.6) are defined as= 3 (  1  )N,TN-1d,eq3crd2,eq = —CrdGe/d^ (2.26)71^:511:1 = Grrd (2.27)= GR^ (2.28)[Shyam Sunder and Wu, 1989b, p. 230]. The (c/e/dt)d,eq term in the scalar drag stressevolution equation (2.7) is defined as(c/c/dt)d,eq(de/dt)0[Shyam Sunder and Wu, 1989b, p. 231].(Crcl,eq) NBV(2.29)At (2.25)Chapter 2. Ice Deformation Model Development^ 152.2.2 Balance EquationsThe relevant balance equations here are those of mass and momentum. Energy balance isnot important to the mechanical model. In Cartesian coordinates we have mass balance[e.g., Mase, 1970, p. 126]op a--(97 5Ti(pvi). 0and momentum balance [e.g., Mase, 1970, p. 128]&xi;^dv;Ph+ xi = P dt(2.30)(2.31)2.2.3 Definition of Strain and Strain RateThe rheology of Shyam Sunder and Wu [1989b] was derived under the assumption ofsmall strains [Shyam Sunder and Wu, 1989b, p. 225]. Hence it is acceptable to use theinfinitesimal definition of strain here, instead of the more complex one for finite strains.The definition of infinitesimal strain is [e.g., Mase, 1970, p. 128]1 (auj Dui=" 2 Oxi axiThe strain rate is (for any state of strain) [e.g., Mase, 1970, p. 113]deij^1 (Ovidi^2Oxi+ Oxi2.2.4 Compatibility of Strains(2.32)(2.33)Arbitrarily chosen strains will not in general yield a single-valued, continuous displace-ment field (because the definition of strain provides six equations to determine the threecomponents of displacement). A real displacement field is guaranteed to exist if thecomponents of strain satisfy the compatibility equations of strain. These are [e.g., Mase,Chapter 2. Ice Deformation Model Development^ 161970, p. 92J:02611 , 402 f 22^2  u612n8X1^(▪ 9X1 UXiCIX202E22  02633^02623 1-^204 04^10M20X3.92633^026112azaX?• 834 OX3OX1a^8623 8631^0E12^azeii(2.34)(2.35)(2.36)(2.37)(2.38)(2.39).^az, ( oxi ax2 ax3^8x28x3a (^+8623 _ 0631 0612^82622 _ .^ax2 8x1 ax2 8z3^axoxi(ax3 ax, ox2 ax349^8623 , ,8631^.8612-t- r,a8.2:83:2 •. ^Similar expressions for the compatibility of strain rates also exist, but are not requiredbecause compatibility of strain rates is assured by the compatibility of strains.2.3 Application of Limiting AssumptionsThe following assumptions will be applied in this section to simplify the problem equa-tions: isotropic ice, very slow flow, radial symmetry, small vertical gradients, and in-finitesimal strains.2.3.1 Isotropic IceIt will be assumed that the ice surrounding the borehole is initially composed of randomlyoriented crystals, and hence is macroscopically isotropic. While this is unlikely to be truein glacier ice, the assumption of isotropy leads to enormous simplification of the rheologyand no information exists that would justify a more complicated assumption.For an isotropic material, (2.8) describing the relationship between stress and elasticand2/3 —1/3 —1/32/3 —1/3^02/311= (2.43)SYM.Chapter 2. Ice Deformation Model Development^ 17strain reduces to Hooke's Law for isotropic materials [e.g., Mase, 1970, p. 143]= 2;Leli AEL6ii^ (2.40)where Sij is the Kronecker delta. For subsequent work I require elastic strain as a functionof stress [e.g., Muse, 1970, p. 143]1 ^Ae7 = 2/1cYij2p(3A -I- 2) akkk. (2.41)For an isotropic material the coefficients c1...c6 in (2.16), (2.20), and (2.22) are allequal to unity [Shyam Sunder and Wu, 1989b, p. 228]. Substituting these into (2.22)yields n = 2. Furthermore the transformation matrices G and H becomeSYM.222G = (2.42)2/3 —1/3 —1/32/3 —1/3^02/31/2These transformation matrices can now be used to evaluate the pseudo-deviatoric stressvectors .C"`, §:ri, and fill".Chapter 2. Ice Deformation Model Development^ 18Substituting for G in (2.19) gives2/3 —1/3 —1/32/3 —1/3^02/32SYM.^222^1^1axx 30Jl - azz+ 0. ___- 3 ^ 3 yy 3 zz1^1^2- 5rCrxx - ayy 5crzz2o-xy2tryz2o-zz=Crxx0" VYazzCrzyCfyzCrZX(2.44)The right hand side of (2.44) can be proven to be equal to the deviatoric stress vector,as follows: The definition of the deviatoric stress is, in tensor notation,1Crij = C•ij -01k8ij• 3and expands to the individual terms1 ,cru. =^— —^+ 0'22 + ass)31G-22 = 1722^(^a22 033)31 ,C733 = 0.33 - -3 Vrii (r22 + 0.33 )V12 = (712U23 = 0-23Cr31 = U31-(2.45)Chapter 2. Ice Deformation Model Development^ 19These can be translated into engineering vector notation, using the definition (2.15)1 (ri„40-•^■-• yy^t'zz)—1 (cr.. -I- 2o ^crzz)- (crz. — crw 2crzz)2crzi,2o.yz2cr z.Cr = (2.46)Note that (2.46) is identical to (2.44), the expression just derived for ,§*• Therefore it isapparent that for isotropic materials kr = CI, and that the matrix G is equivalent to adeviatoric operator (because G = 1-6) . Thus===or, in tensor notation,CrijA"ci^dS CrijSR*(2.47)(2.48)(2.49)Now the definition of the equivalent stress can be simplified. Substituting for 77 andGO in (2.18) using 77^2, (2.19), and (2.47),2eq — 3 —T—cr CTChapter 2. Ice Deformation Model Development^ 203 r„.2" eq — —2 La.., airy, zz, Crzy, (ryz, o](2Crxx — Cryy Crzz )—1 (a^2o-yy — Cize3 xx—1 (a — cryy 2U)3 xx 1^1^1 1=er_ :x — —2o-xxoyy — —2o-zzcr.. io-z„ayy ti — —2 Crvy trzy1^1— —2 crzzo-zz — —3 cryyc•zz crz2z 2azy 2cryz 2crzz_1 (0.2. _^4_ 0.2 ) 4_ _1 (0.2 _ 20.yy0.zz 0.z2z)2^ 2 YY1^,+— (T" — 2crzzaz. cr:x)^3o-yz 3azx2^zz—1 (Cr2 — 2CrxxCryy cr2^—1 (cr2 — 2cryycrzz o)2 k " 2 11111+— (0-2 — 20-zza= cr:z) + 3 [crzy Cryz crzz]2^zz1—2 [(azx — ayy)2 + (Cfyy — Cree) 2 + ((Tye — crxx) 2} + 3 [o-zy o-yz crzx] (2.50)The right hand side of (2.50) reduces to the conventional definition of equivalent stress.The conventional definition of equivalent stress is [e.g., Mase, 1970, p. 182]creq2 = 3E' —^-2 — 2 ta (2.51)where Ei2 is the second invariant of the deviatoric stress tensor [e.g., Fung, 1969, P. 81]which expands to1= —6I—^-r v0-22)^722 — u33)) ,^k0.33 ,^kai2-1- 0.23 + 0'30-Hence12^ ,^ , f^N^.3 ICreq — —2 [kali — Cr22)^k(022 — Cr33) --r- “733 —^2 } -I- k0-12+ Cr23 a31) (2.52)_Chapter 2. Ice Deformation Model Development^ 21Note that (2.52) is identical to (2.50) just derived. Hence (2.51) may be used for thedefinition of creq• Therefore, from (2.5), (2.17), (2.47), and (2.51), the steady-state viscouscomponent of the strain rate is given bydeyi dt(del dt)0jLJjJ3 ( \N N-11eq 3^)N ^(N —1)12—2 V(2.53)This is identical in form to the commonly used Glen Flow Law for "steady-state" viscousdeformation of ice [see e.g., Paterson, 1981, p. 31]. Similarly, it can be shown that2^3 1 cCr d,eq^crip'ij• (2.54)The evolution equation (2.4) for R must be examined. The vector equivalent ofhod de! I dt is Hdetl dt. Thus in vector form (2.4) becomes, with 77 = 2,= —2 AEHd;t(2.55)dt^3^dtExpanding H det/ dt using (2.43) for the isotropic definition of H gives2/3 —1/3 —1/32/3 —1/3^02/ /dtdetzz/dt2 detxy/dt2 detyz/dt2 dezz / dtH t =dtZetdt(2.57)Chapter 2. Ice Deformation Model Development^ 22d-;t—dt23- (2 dexxl dt — dcyyt Idt — dezzldt)(—dexzldt + 2 deyyt /dt — detzz/dt)(—det../dt — /dt + 2 detzz/dt)(Idyl dtdeytzldtde/dt(2.56) Because de!a/dt is proportional to 'crej, for which the first invariant vanishes, the firstinvariant of de:i/dt also vanishes. Hence dext.Idt deyyt Idt dezzldt = 0 and=dtor, in tensor notation,dd.! dtdeyyt dtdezzldt2 dexyl dt2 d4zIdt2 detz.ldt(2.58)(2.59)de de^de •Hi •'kl^= 8i1c87^=dt dt^dt •Therefore (2.4) becomesdRii 2 AEd4idt^3^dt •In summary, the rheology for isotropic ice is described by the equations==^dh EZi^ (2.60)de' V de2i^,J^Bj dt^dt^dt(2.61)1 Ae • =^Cli3^2/.z (3A + 2f.t)akkSia^ (2.62)— v —^ (2.63)Chapter 2. Ice Deformation Model Development^ 23dRi; 2 AE detii(2.64)dt^3^dt()—(16 (2.65)dt 0=^d (dE)^ (2.66)-11t oE 1 crd,eq)N(2.67)cra,eq BV )(1-3--117 )N a N-1d,eq^ (2.68)where( de)= 1s-1^ (2.69)3,/=^ (2.70)creq^ i- a Jai.;2,....2^3/ di d (2.71)' d,eq = — Cri• Cri .2 3 3A, =^( 1 )N crit-i (2.72)2 Vl^q3 ( 1 N N-1At =^.17)^dcr ,eq •^ (2.73)2 ./32.3.2 Very Slow FlowThis assumption simplifies the momentum balance equation by eliminating the inertialterm, thus reducing it to a force balance equation. Slow flow implies that accelerationsare very small, and p dvildt < pfj,acrii/Oxi. Thus momentum balance (2.31) becomesacri •Pf, + 3^ = 0-^ (2.74)uf Xi2.3.3 Radial Symmetrydtde!idtdBdtHere I convert from Cartesian to cylindrical coordinates, and apply the assumption ofradial symmetry.Chapter 2. Ice Deformation Model Development^ 24First consider the definition of strain (2.32). The gradient of a vector field in termsof cylindrical coordinates is [Prager, 1961, P. 36]Ouioxi(r) (r) °Ur^(r) (9)^(r) (z) aUz- ei ei^- F ei ej + ei ei or(0) (r) 1 CU?^)^(0) 0)1 (au+ e e 0i `•^— — us + e• e • — — + ur ao r 00^r+e(z)e(r) OUr ecoe(e)Oue e(.)e(z)auz3 az^3 az^3 az(e) (z)1 auzei ej —r 00(2.75)With radial symmetry alae 0 and ue = 0 so that the gradient of the displacementfield becomesauOxi eloe(r)2t-`1. + e(r)e(:9) (0) + e(r)e(z)auz3 Or^3^2 3 Or^+er er (0) + er^.e(8) + e• e3 (0)r^3+elz)er -b--+ elz)er (0) + e1 e 5(z) aau z .z (2.76)Substituting this definition into the definition of strain (2.32) yieldsOurerr =^ (2.77)Orfee = (2.78)auzerr Oz^ (2.79)ere = 0 (2.80)ee. = 0^ (2.81)1 (auz au?)=^+ (2.82)erz 2 ar^azSimilarly the strain-rate components becomede,.,.^Dv,.(2.83)at^Orde00vr^ = (2.84)dtde„^avz(2.85)(2.86)(2.87)(2.88)dt^azdere =0dtdee. = 0dtderz^1 (avz avr)dt^2 ar + az )Chapter 2. Ice Deformation Model Development^ 25Consider the mass balance equation (2.30). Expanding it by carrying out the spatialdifferentiation givesop^Ova^Op+19:9X; + vi azi = 0.Noting that Op/at vi 8p/Oxj = dpldt, the material time derivative of density [ .g.,Prager, 1961, P. 74], we haveavj dpP axj^dt^°.(2.90)The divergence of a vector field expressed in cylindrical coordinates is [Prager, 1961,p. 36]avj av, 1 ( ay,^avzax,^Or + -7: a° +vrj^az •With radial symmetry vo = 0 and alae 0 so that mass balance (2.90) becomesav, vr avz\ dpP^+ 7 +^+ Ti °•t (2.91)(2.92)Next consider the force balance equation (2.74). The divergence of a tensor fieldexpressed in cylindrical coordinates is [Prager, 1961, p. 37]acrii^e(r) 00r,. + 1 acrer + (Tr,^coo + aux,-^ ,-ax,^3^Or^r ae^az+e ^ + 59Ocrre +1 ace°( + aro + as,^au+  ze^3 ^ar^r ae^az+e^Or(aurz +1 (threz 4_ crTz) + OCTzz)^3^a ^r 09^az • (2.93)(2.89)Chapter 2. Ice Deformation Model Development^ 26With radial symmetry alae^0. Because dere/dt = cleoz/dt = 0 (and therefore E,.9 =fez = 0), and the ice rheology is isotropic, are = crez = 0. Hence the divergence of thestress field reduces toao-i;ax,(r) Carr 1 ,^aazr)= e. ^ + ((Tr^cree)+3^ar^r az+e(z) Car.. + crr. acrzzOr^r^az ) •Thus the force balance condition yields two non-trivial equations(2.94)acrrr^1Pfr ^ +(crrr aeo) acrzr) = 0^(2.95)ar^r azoar..^CrTZ^aCr2.2 pfz+ ^ +^+ = 0.^(2.96)Or^r^azGravity is the only body force acting on the ice. Therefore the force balance equa-tions (2.95) and (2.96) become(aarr + 1 (arr ow) (9;zr)Or^r(aarz arz acrzz)pg + Or^r^az= 0^ (2.97)= 0.^ (2.98)2.3.4 Small Vertical Gradients (Except Hydrostatic Pressure)With the exception of the elastic response to pressure change, the ice rheology is inde-pendent of hydrostatic pressure. For linear elasticity, hydrostatic pressure does not affectthe solution, except to cause a small change in ice density, which is second order, andhence insignificant. We can therefore neglect the effects of vertical gradients in crx,. andCTZT •Applying this simplification the mass balance relation (2.92), becomes(2.99)Chapter 2. Ice Deformation Model Development^ 27and the force balance relations, (2.97) and (2.98), becomeaurr , 1 fOr^r-r- kurr — Cr t 9 0) = 0aa. pg + ^ = 0.azSince crzz = —p 'crzz, and a'azziaz = 0, (2.101) becomesOPPg = 0,oz( 2.100)(2.101)(2. 1 02)the standard equation for hydrostatic pressure. Since the hydrostatic pressure has aninsignificant effect on the flow solution, this equation has no influence, and I shall hence-forth treat the state of stress in the ice as being independent of depth.With (9/az = 0, the strain expressions (2.77)—(2.82) simplify to(2.103)(2.104)(2.105)(2.106)(2.107)(2.108)and the strain-rate expressions (2.83)—(2.88) becomeOurErr^OrcooE rr =0Er =fez = 0EZT - 0derrdt(1E09 yrdt=0000.dtclErodtdEezdl(kyrdt=(2.109)(2.110)(2.111)(2.112)(2.113)(2.114)Chapter 2. Ice Deformation Model Development^ 28Hence, there exists a state of plane strain in the r-9 plane. It is also apparent thatÔ, and are principal axes of strain (those axes for which all non-diagonal terms ofthe strain tensor vanish). For isotropic ice these are also the principal axes of stress.Therefore in this coordinate system the shear terms of all tensor quantities (o-ii, 14,eii,^4i, 4, and ey;) and all their time derivatives and deviators vanish. Only thediagonal terms remain. The notation for the diagonal terms of the tensors can now besimplified, so that o„ becomes or, and so on. In the simplified notation the ice rheologybecomesE,. = E +^ (2.115)€8 = fee + Eto' (2.116)0 = cez + evz^ (2.117)^de:).^de:.^de:= (2.118)t^dt + dt^dee'^dee de:= (2.119)^dt^dt + dtdel^det^des=^z + z (2.120)t^dt^dt1^AE,e. = —21iar 2p, (3A + 2/t) (ar ± ao + az )^(2.121)1=^A(2.122)Ee^—a, ^211,^2p (3A + 2/4) (ur ± CO ± Cr z)1=^A(2.123)ez^— Cr z ^2p^2p, (3A + 2) (a'. + a° ± gz)d (2.124)CI r = ar — Rrd = - (2.125)^go^as  Red= uz -^ (2.126)a z RzChapter 2.^Ice Deformation Model Development2dt 3dR9 2dt 3dRz 2dt 3de'(2.127)dtdct(2.128)dtdEt(2.129)dtde) 0 (2.130)de,'9^(1) ^ (2.131)Asicrs^0dtdesz^Azicrz (7ide) 0 (2.132)dtderAt' ()^ (2.133)ad cAdt^dtdee = Azi4 ((A) (2.134)dt^dtdetz=^'acz‘ (1) (2.135)Atdt29wheredB^1 y cr ddt BV^eq(2. 136 )cr„ =(2.137)(2.138)(2.139)(2.140)(2.141)(2.142)1s'3 ( N N_i2\V) "-23- (iur2^'1792^icrz2)1ar — — ar ao ± az )31— —3 (ar^+ az)1o-z — — (q,.^aer + az)3(ToCrzChapter 2. Ice Deformation Model Development^ 303(  1^ )N d N-1= (2.143)At^creg2 BV^adeg .^(3IF^ 2 + isre + ,012) (2.144)2 \i d ad..^— —1 (ard + 4 + azd) (2.145)ar^r 3/ad^1= erect _ _3 (0.7.d + crcj + azd) (2.146)e/ d^ad=^— _1 (ard + 4 + 0.zd) (2.147)cr. .^3It is also apparent that the mass balance equation is superfluous for small strains.If the strain-rate definitions (2.109)—(2.114) are substituted into the mass balance equa-tion (2.99), we have(de.^deo dez)^dpP = 0. (2.148)dt^dt^dt )  The term in brackets, dekkldt, is the cubical dilation rate, that is the rate of change ofvolume of an element. Because the mass of an element (m = pV) is constant^d ^dVdt = ° = P dt^dtanddp^p dVdt^V dt.Substituting (2.149) into (2.148) gives(2. 149)( de,.^deo^dez^p dV^0.^P dt^dt^dt ) V dtSimplifying and rearranging,( de,.^dco^dez)^1 dVdt + di + dt ) V dtwhich is just the definition of the cubical dilation rate. Hence it is apparent that themass balance equation is unneccessary in this problem, as the ice density is of no interest.The strain-rate definitions also become superfluous and will not be considered further.Chapter 2. Ice Deformation Model Development^ 31The compatibility of strain equation can now be developed. One possibility would beto apply the foregoing limiting assumptions to the strain compatibility conditions (2.34)—(2.39). A simpler approach, which yields the same result, is to start with the straindefinitions in cylindrical coordinates. From the definition of tangential strain (2.104) itfollows thatV.,. =Teo.^ (2.150)Substituting this into the definition of radial strain (2.103), and differentiating, gives,er^— krce)ar569rOrOr569r^— e,. = 0. (2.151)An equivalent expression for strain rates can be similarly derived, but is unnecessary,because (2.151) sufficiently constrains the strains.2.3.5 Infinitesimal StrainsThe assumption of infinitesimal strains has already been invoked a number of times inthe preceding derivation. Here I use it once more to set the advective part of spatialtime derivatives to zero, and allow the preceding equations to apply at spatially fixednodes. As a result of this assumption all the material time derivatives may be replacedby spatial time derivatives. Hence the ice rheology becomes(2.152)(2.153)(2.154)Chapter 2. Ice Deformation Model Development^ 32aft,:^ae^SE.=^r + rat^at^at(34 = afte .94at at + &af: _ aft 06:at — at + &OR,.—=^at^ToTMe^2^ete—AEat^3^at^aRz^204—AE^at a3^at(2.155)(2.156)(2.157)(2.158)(2.159)(2.160)(2.161)(2.162)(2.163)(2.164)(2.165)(2.166)(2.167)(2.168)(2.169)(2.170)(2.171)(2.172)c^ A1= —cr, ^ fry + cre + crz)^r 2p^2p (3A + 2p)e^1=Aee —^2p^2p (3A + 2p) (ar^+ crz)1 A= 2ttaz 2(3A + 2) (cr, cre + az)"r =crea^era —or^crz — Rz1 s-1(de 0 =/c•r/ao===At =Erpar 1Dr + -r^- o-e) 0Definition of Strain(2. 1 86)(2.187)Chapter 2. Ice Deformation Model Development^ 33aB —= E()N d N1at^Cr^.BV "whereo-eqd = \/_3 (42 + /act + ,(712)2'4a.r = 4 — 1 (4 + 4 + 4)3 k'40'0 = 4— _1 (4 + 4 + 4)3 \„,s1• - ' z = crd _ _11^1 (crd 4_ 4 + 4).^3^T .(2.173)(2.174)(2.175)(2.176)(2.177)(2.178)(2.179)(2.180)(2.181)(2.182)(2.183)(2.184)The the other relevant equations are summarized below:Compatibility of Strainsaceco+ r^- c,. = 0^ (2.185)DrMomentum BalanceChapter 2. Ice Deformation Model Development^ 34Ures = —r(2. 1 88 )ez = 0.^ (2.189)The above equations are in no particular order. They can be grouped as follows toform the field (differential) equations of an initial value problem, with algebraic constraintequations. The differential equations are (2.155)—(2.157), (2.164)—(2.166), and (2.173).The constraint equations are force balance (2.186) and plane strain (2.189). All the otherequations are used to calculate the intermediate variables needed for the differential andconstraint equations.2.4 Initial and Boundary ConditionsHere I will state the initial conditions and boundary conditions, which together with thefield equations just developed, fully specify the initial value problem in the ice. Thesefollow directly from the assumptions of the model, as outlined in section Initial ConditionsThe initial state of the ice is that of hydrostatic equilibrium with no viscous strain. Hence(r, 0) = (78 (r, 0) = o^(r, 0) = —Po (2.190)(r, 0) =^(r, 0) = Evz (r, 0) = 0 (2.191)R,. (r, 0) = Re (r, 0) = Rz (r, 0) = 0 (2.192)B (r,O) = Bo (2.193)where po is the background hydrostatic pressure, and B0 is the initial value of the scalardrag stress, which depends on the grain size of the ice [Shyam Sunder and Wu, 1989a,p. 56].Chapter 2. Ice Deformation Model Development^ 352.4.2 Boundary ConditionsThe condition on the inside boundary is fairly obvious. The radial compressive stressmust be equal to the water pressure in the borehole, and the displacement at the boreholewall must equal the change in borehole radius. Thuscr,. rb , t = —Pb^ (2.194)u,. (rb,t) = Arb.^ (2.195)A boundary condition in terms of strain can be obtained from this by substituting (2.195)into the definition of tangential strain (2.188), yieldingE9 (re,, t) =Arb(2.196)rb At the outside boundary there are numerous possibilities. Although the outsideboundary in the problem as formulated so far is at infinite radius, in the numericalproblem this boundary will have to be placed at some large, but finite radius. The choiceof outside boundary conditions will be made with this fact in mind.Some possibilities, all of which are true at infinite radius, are as follows: zero strain,specified hydrostatic stress (all stress components specified), or specified radial compres-sive stress only (equal to hydrostatic pressure). The problem with zero strain is thatif the boundary is moved to finite radius there will be no flow unless the ice is elasti-cally compressible. It would be useful if the model could be run with elasticity turnedoff, to facilitate comparison with analytic solutions. The second and third options aresimilar. The specified-radial-stress option does have the advantage that some analyticsolutions are available for thick-walled cylinders subject to internal and external pressure,for which it is the appropriate boundary condition [Kjartanson et. al, 1988, p. 257; Jaegerand Cook, 1969, p. 128]. Therefore I will apply the specified-radial-stress condition tothe outside boundary.rb—Po.(rb, t)^—PbArbEs (rb, t) =cr,. (oo , t) =(2.197)(2.198)(2.199)Chapter 2. Ice Deformation Model Development^ 36Thus, to summarize the ice boundary conditions:This completes the description of the ice deformation initial value problem. In Chap-ter 3 the bed and borehole equations are developed, and synthesized with the ice defor-mation model developed in this chapter, to form the complete coupled model.Chapter 3Bed Equations, Forcings, and Model SynthesisThe flow of water in a thin sheet aquifer when flow is induced by borehole responsetests has been modelled by Stone and Clarke [1992] and Stone [1993] using a method-of-lines solution on a radially symmetric grid. Smart and Clarke [submitted] considerthe same problem that I do here; namely, to determine the hydraulic conductivity of thelow-permeability layer beneath Trapridge Glacier by considering the behaviour of waterin an unconnected borehole. They find an upper bound to the hydraulic conductivityby considering the bed to be a uniform half-space bounded by impermeable ice, exceptat the bottom of the borehole where a hemispherical water-filled cavity exists. Theyconsider the steady-state analytic solution for Darcian groundwater flow.In this model I also approximate the bed as a homogeneous half-space, with hydraulicproperties equal to those of the low-permeability till layer. This is a reasonable approx-imation as long as drainage channels are not too close to the hole, and the radius ofthe hole-bottom cavity is small relative to the thickness of the low-permeability layer. Ifurther assume that the flow is described by the standard diffusion equation for transientgroundwater flow in a saturated porous medium [e.g., Freeze and Cherry, 1979, p. 65].Many assumptions underlie this equation, including that Darcy's Law applies, and thatthe only motion of the soil matrix is the small vertical movement caused by elastic de-formation of the matrix as effective pressure is changed. The ice and bed are coupled byconservation of mass of compressible water in the borehole.The borehole freezing thermal problem has been considered thoroughly in connection37Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 38with the interpretation of thermistor measurements and hot water drilling in cold ice[Jarvis and Clarke, 1974; Humphrey and Echelmeyer, 1990; Humphrey, 1991]. Modelssuch as the finite-element model of Humphrey and Echelmeyer [1990] describe the freezingof ice to the wall of a borehole, and hence predict its radius as a function of time,background ice temperature, and drilling method. This is the best way to obtain thefreezing forcing. Unfortunately time constraints prevented me from implementing sucha solution, so instead I assume that all freezing is confined to the top of the borehole,so that the hole gets shorter, but its diameter remains constant. The borehole length asa function of time is described as a series of decaying exponentials, having coefficientsthat can be determined from temperature measurements made as a borehole freezes. Ialso allow for arbitrary sudden pressurization of the borehole, to permit the modelling ofsudden pressurization and decay events.In this chapter I develop the bed model, borehole model, and description of freezingand other forcings, then synthesize these with the ice deformation model to form thecomplete coupled model, and formulate the problem for numerical solution using themethod of lines.3.1 Basal Groundwater FlowConsider a homogeneous, isotropic, permeable half-space bounded on top by impermeableice, and by water in a hemispherical cavity at the bottom of a borehole. It is postulatedthat such a cavity is formed during the hot water drilling process. Let the coordinatesystem be spherical, with the origin at the centre of radius of the hemispherical cavity(Fig. 3.1). While it is known that the till under Trapridge Glacier is deforming [Blakeet al., 1992], it will be assumed that this deformation is slow enough that the standardChapter 3. Bed Equations, Forcings, and Model Synthesis^ 39Fig. 3.1: Bed geometry and coordinate system.Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 40equations for groundwater flow in a porous medium still apply. Taking the till defor-mation into account would introduce much more complexity than is justified in such arough model.The standard water mass balance equation in a saturated porous medium [e.g., Bearand Verruzjt, 1987, p. 63] gives^-v^Pwg (a +^7113)-07Oh^(3.1)where is specific discharge (or volume flux) of water, pw is density of water, a is soilmatrix compressibility, n is porosity, and 13 is the compressibility of water. Stone [1993,pp. 35-36] explicitly states the assumptions required by this expression. The generalform of Darcy's Law (which applies to low-Reynolds-Number flow through a stationarymatrix) [e.g., Bear and Verruijt, 1987, p. 38] is writtenOh^= —Kij(zi) Q^ (3.2)where Kij is the second rank tensor of hydraulic conductivities. I will further assumethat the basal material is hydraulically homogeneous and isotropic, in which case Darcy'sLaw simplifies toOhqi = —Kor, in vector notation,^—KVh^ (3.4)where hydraulic conductivity K is a scalar constant.Substituting Darcy's Law into the mass balance equation (3.1) gives the well-knowngoverning equation for transient flow in a compressible porous mediumOhICC1211 = pwg (a + ni3)^ (3.5)Ox; ( 3 . 3 )Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 41This equation can be expressed in spherical coordinates and simplified to take advantageof symmetry. First express Viz in spherical coordinates [e.g., Spiegel, 1968, p. 126]v2h 1 a (7. 2 ah\ +  1  a 1. 8ah\ +  1 02h^r2 Or k Or^r2^e ao^ao^r 2 2 9 802.^(3.6)With spherical symmetry moo 0 and moo 0, hence (3.6) simplifies to give1 0 al^V2h =^(r2-±.) .^ (3.7)r2 Or^OrThus (3.5) becomesah1 a,^(2r- —) —^Ohr2 ar^ar^Pwg(a+n13)Wiwhich can be rearranged to giveOh ^K ^1 a ( 2ah\at - pw g (a + nfl) r2 Or kr Or ) •The initial conditions and boundary conditions for the basal groundwater flow prob-lem follow directly from the assumptions that the bed is initially in a no-flow (uniformhydraulic head) state, and that the background head in the bed changes much moreslowly than the head in the borehole. The head at the cavity wall must also be equal tothe borehole head. Henceh (r, 0) = po I pwg, r, < < 00h (oo,t) = polpwgh (re, t) = pb(t)Ipwg, t > 0where re is cavity radius (the interface with the borehole), polpwg = ho is backgroundhead in the bed, and pb(t) is pressure in the cavity at the borehole bottom.Therefore, to summarize, the full initial value problem describing the response of auniform permeable half-space, initially at constant head, to an arbitrary pressure forcing(3.8)(3.9)Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 42at a hemispherical interface is given byOh ^K et =(2ahpwg(a+ nt3) r2 -6; r^r < 00 (3.13)subject to h (r, 0) = pol pwg^ (3.14)h(oo,t) = pol pwg (3.15)h(r,,t) = pb(t)I pwg.^ (3.16)3.2 Borehole EquationsThe borehole is filled with compressible water, which acts to couple the ice and bed. Thecoupling is provided by the water mass balance equation for the borehole. The mass ofwater in the borehole plus the mass of water lost by flow into the bed or freezing (seesection 3.3) must be equal to the initial mass of water in the borehole. Thus(3.17)where mw is the mass of water in the borehole at time t, maq is the mass of water whichhas flowed into the bed between time 0 and t, mf is the mass of water which has frozenbetween time 0 and t (see section 3.3), and mo is the mass of water in the borehole attime 0:The mass of water in the borehole at any instant is derived as follows. The compress-ibility of water must be taken into account; hence the density of water may in generalvary in space as well as time. Thusmw = f w pw dVvwhere Vw is the sum of the volume of the borehole VB and the volume of the basal cavityVc. The depth of the cavity is small, so that the density of water in it can reasonablyChapter 3. Bed Equations, Forcings, and Model Synthesis^ 43be considered constant. Density in the borehole cannot, however, be considered constantsince the pressure change with depth is of similar magnitude to the pressure change withtime.First, consider water mass of the cavitym(t) = pc(t)Vc= pc(t)17rr!^(3.18)where pc is cavity water density, given by the assumed equation of state for waterPc(t) = Po exP {/(3 (Pb (t) 1)0)]Thus2m(t) = 1rrpo exp (Pb (t) — po)]Next, consider the water mass of the borehole itself(3.19)(3.20)mb(t) = iv pw (z ,t) dV^ (3.21)wheredV = rr(OdzPw(z,t) = po exP (P(z, t) — Po)1.^(3.22)Here p(z,t) isp(z, t) = pb(t) — foz pw(z, t)g dz.^ (3.23)Because pw changes very slowly with pressure, and hence with z, and the pressuresconsidered are never especially large, pw can be considered a constant Po in (3.23).Hencep(z,0 = pb(t)— pogz.^ (3.24)Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 44Substituting (3.4) into (3.22)pw(z, t) = po exp [fi (Pb(t) — pogz — po)}Po exP [fl (Pb(t) — Po)] exP [—fipogz] .^(3.25)Now substituting (3.25) into (3.21) and simplifyingmb(t) =rL(t)jo1 Po exP (Pb(t) — Po)] exP [—PPogz]rrl(t) dz1(t)Po exP b6 (Pb(t) — po)] irrl(t) I^exP [—PPogz] dzJo(t)PoexP [P(Pb(t) — 790)] [_opog exP (—Opogz)]1 rr(t) exp [p(pb(t)— p0)][1 — exp (—PpogL(t))1 .gL(t) (3.26)Note that the height of the water-filled borehole L is time dependent. This allows forthe approximate freezing forcing developed in Section 3.3. Borehole radius rb(t) can befound from the tangential strain in the ice at the borehole wall. The ice strain boundarycondition at the borehole wall (2.198) leads tofe(rb,t) _Arb(t) rb(t)rb(t) — rb,o(t)?b(t)^•For infinitesimal displacementsrb(t)— rb,o(t)^rb(t)— rb,o(t)rb(t)^rb,0andrb(t) = Tb,0 (ce(rb,t) +1) .^ (3.27)Combining (3.18) and (3.26) the total mass of water in the borehole and cavity is,mw = mb inChapter 3. Bed Equations, Forcings, and Model Synthesis^ 45mw = R-7-1(t) exp V) (pb(t) — AO] [1 — exp (—,3 pog LW)]13 g2 ,^r„ ^„+ —Irr-po exP E(Pb (t) — P0).1 •3 '(3.28)Evaluating this at time t = 0 gives mo, henceMO = IrT/'° exp [,8 (po — RI)] [1 — exP (-0^-3-PogLo)1+ 1'7'W° exP [fl (Po — Po)]„ 2 313g2=- rrb,o,0 r1 _ exp FfipogL0)} + i2irr:poe0)3g^12 2--= 7rrb'0 [1 — exp (-13pogLo)] + '71-7. cP0-3i fig(3 .29)Now let us consider the flow of water from the borehole into the bed. First, Darcy'sLaw must be converted to spherical coordinates and simplifying assumptions applied. Inspherical coordinates the gradient can be written [e.g., Spiegel, 1968, p. 125]v _ e(r) (9 + e(9)-1 a + „(0)  1 ^aOr^r ao e r sine aoand for spherical symmetry 0/89 = 0 and a 1 ao = 0, so thatV' = e(r) aOrSubstituting (3.31) into Darcy's Law (3.4) yields the radial specific discharge,Ohqr = —A --.Or (3.31)(3.32)Equation (3.32) gives the volume flow rate per unit of area (the flux) across a hemispher-ical surface centred on the origin. The total volume flow rate through such a surfaceis simply the flux multiplied by 277-,2, the surface area of a hemisphere. Therefore thevolume flow rate into the bed (the flow rate at the cavity surface) is(3.30)(3.33)Qw (r,) = —271-7.1C :-711r,Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 46and the mass flow rate isQ, = —27rperif^ (3.34)rcSubstituting (3.19) for pc givesQm = —21-po exp (pb(t) — po)lr^ (3.35):K —reThe total mass of water lost by the borehole to the bed is the mass flow rate integratedover time. For compatibility with the rest of the problem, this can be•expressed as adifferential equationthrtaq^Oh= —2rpo exP [fi (Pb(t) — P0)]r2K —at ar (3.36)subject to maq(0) = 0.3.3 Forcings I shall consider three distinct types of forcing: volume forcing caused by' freezing to thetop of the borehole, sudden pressurization, and step pressure change.3.3.1 FreezingFreezing of borehole water causes a volume forcing because water expands as it freezes.Freezing is modelled in this thesis approximately by top-down freezing, at an arbitraryrate, with no freezing to the borehole walls. Under this simplifying assumption freezingforcing enters the model by specifying L(t), the height of the water-filled borehole asa function of time. The most important place L(t) appears is in the term M1 in theborehole mass-conservation equation (3.17).Let the height of the water-filled borehole as a function of time be given by somearbitrary function having coefficients that can be assigned to approximate the desiredChapter 3. Bed Equations, Forcings, and Model Synthesis^ 47freezing rate. As a functional form, I have chosen a series of decaying exponentials whichdisplays the appropriate behaviour of initially fast decay followed by slower decayL(t) = a() aieb't a2eNt a3eb3t a4eNt a5e 5 t^(3.37)where ai and bi can be arbitrarily chosen. The mass of water Mf which has frozen to iceis then the volume of ice produced multiplied by the density of ice. The volume of newlyfrozen ice V1 isMt) = irrl(t)—dL dt.dt (3.38)Because rb changes much more slowly than L, rb(t) may be replaced by rb,o, a constantV1(t) = irr [L(0) — L(t)1^ (3.39)thusmf = rrt,o [L(0) — L(t)]^ (3.40)3.3.2 Sudden PressurizationThe sudden-pressurization forcing is included to allow modelling of sudden pressurizationand decay events observed in Trapridge Glacier pressure records. In principle a suddenpressure change could be modelled by appropriate choice of initial conditions. The bore-hole pressure would be something other than Po, and stresses and strains in the ice wouldbe those corresponding to the analytic elastic solution. The difficulty with this approachis that, in general, the analytic elastic stresses and strains will be inconsistent with thesystem of model equations (they will not be possible solutions of the model equations)due to errors introduced into the numerical solution by finite differencing. Many nu-merical integration schemes fail when started with inconsistent initial conditions. Analternate approach which avoids this difficulty is to start with normal equilibrium initialChapter 3. Bed Equations, Forcings, and Model Synthesis^ 48conditions, and then impose a very rapid, but smooth, change in the borehole pressure.I have chosen a shifted half cosine for the ramp-up function, due to its smoothness,(t = PO Pf°rce {1 -1- cos {7 ^t^I)]2^tramp (3.41)where pf, is imposed pressure change (so the pressure during the ramp-up period isPo + pforce) and tramp is ramp-up time (half the cosine period).3.3.3 Pressure Step ChangeIt is desirable to be able to simulate a creep test, in which the borehole pressure is sud-denly changed, then remains constant. This allows comparison with analytic solutions(even though the infinitesimal-strain assumption of this model will cause inaccurate re-sults at large strains). Creep tests are simulated by using (3.41) to ramp up the pressure,and then holding borehole pressure constantpb(t) = Po pforce for t > tromp.^ (3.42)3.4 Summary of Coupled Initial Value ProblemsI now have all necessary equations to describe the coupled initial value problems: icedeformation, groundwater flow in the bed, and coupling through the water mass balancein the borehole. The initial value problem in the ice is specified by (2.152)-(2.189), withinitial conditions (2.190)-(2.193) and boundary conditions (2.197)-(2.199). The initialvalue problem in the bed is given by (3.13), with initial condition (3.14) and boundaryconditions (3.15)-(3.16). The borehole coupling equations, which include the forcings,are (3.17), (3.28), (3.29), (3.35), and (3.35).Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 493.5 Method of Lines SolutionIn the preceding sections the equations of the coupled ice and bed initial value problemshave been derived. The ice and borehole equations are non-linear, and there are bothpartial differential and algebraic equations to be simultaneously solved. There is nohope of an analytical solution, so numerical methods must be turned to. Few numericalmethods can handle such a problem. The finite element method has been used to solveproblems involving the ice rheology of Shya.m and Wu, at least in one dimension [ShyamSunder and Wu, 19894 The numerical method of lines (discussed at length in [Schiesser,1991]) is capable of solving non-linear initial value problems. It can also handle algebraicconstraints (if an appropriate numerical integration routine is selected), and coupledinitial value problems also present no special difficulty. I have chosen to use the methodof lines, primarily as a result of my familiarity with the method and its comparativesimplicity.Briefly, the method of lines involves discretizing the region into a grid, replacingcontinuous functions of the spatial variables with discrete functions at the grid nodes,replacing spatial derivatives in these equations with finite-difference approximations, andthen solving the resulting system of ordinary differential equations, and possibly algebraicequations too, by numerical integration.In this section a coordinate transformation will first be applied to the spatial variable,in the ice and bed, to concentrate the nodes where the field variables are changing rapidly.Then the grids will be established, and partial differential equations over the regionsreplaced by systems of ordinary differential equations. The equations that apply at theinterior and boundary nodes will be stated. Finally the resulting system of equations forthe coupled problem will be presented in a form suitable for programming.Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 503.5.1 Spatial Coordinate TransformationThe method of lines requires that the ice and bed regions be discretized, so that finite-difference spatial derivatives can be taken. In both the ice and bed the field variables(e.g., stress and strain in ice, and head in the bed) change much more quickly at smallradius, close to the borehole. It is readily apparent that at infinite distance from thehole, all these variables will become spatially constant. Intuitively one would expect amore accurate approximation to the true result if the nodes are more numerous wherevariables change rapidly, namely at small radius. Some possible ways to achieve a closernode spacing at small radius include: (1) increasing the total number of nodes withoutaltering their distribution, (2) applying a spatial coordinate transformation and spacingnodes evenly in the transformed variable to concentrate nodes at small radius, or (3)an adaptive grid could be used to generate new nodes when and where they are needed[Schiesser, 1991].The first alternative has two problems: (1) Increasing the number of nodes increasesthe computing resources required. (2) More importantly, increasing the number of nodes,and hence the number of ordinary differential equations, increases the problem's stiffness[Schiesser, 1991, p. 22]. (A stiff problem is one in which there are widely differingtime scales, requiring the integrator to use very small time steps to achieve reasonableaccuracy.) Hence the problem may become too stiff to be solved. The adaptive gridmethod is very good for problems in which the solution shape is not known a priori, orwhere the solution shape changes with time as in the present problem. Unfortunately ithas the disadvantage of being complex to program. As a result of the drawbacks of theother two alternatives given above, I chose to transform the spatial coordinate r, in theice and bed.One approach would be to transform the spatial variable such that the field variablesChapter 3. Bed Equations, Forcings, and Model Synthesis^ 51would change by comparable amounts between each pair of adjacent nodes. In this casethe resulting solutions for the field variables would approximate straight lines in thetransformed space. Although I do not in general know what the ice and bed solutionshapes will be, the analytic solutions for a number of simple cases can be used as guides:1. For steady-state radial groundwater flow in a whole-space or half-space, h is pro-portional to lir (see [e.g., Carslaw and Jaeger, 1959, p. 247] for the solution to theequivalent heat conduction problem).2. For steady-state or transient flow of linear (elastic, viscous, or viscoelastic) icesurrounding a cylindrical hole, stress and strain are proportional to 1/r2 [Nye,1953]3. For steady-state viscous flow of Glen-Law ice surrounding a cylindrical hole, stressis proportional to 1/7-2/3 and strain is proportional to 1/r2 [Nye, 1953].Thus it appears reasonable to choose a transformed coordinate R = 1/r in the bed.The best choice for the ice is less obvious. Therefore a range of transformations will beallowed: R = ?X where x is a finite nonzero real number, and R = In r. The logarithmictransformation is included because the steady-state solution to the diffusion equationwith cylindrical symmetry is logarithmic [e.g., Carslaw and Jaeger, 1959, p. 189], andone might expect some diffusive behaviour in this problem.Here I apply the spatial coordinate transformations to the ice deformation equations.The only expressions which explicitly involve r are initial conditions (2.190)—(2.193),boundary conditions (2.197)—(2.199), compatibility of strains (2.185)e,€9 +^— Er = 0orand force balance (2.186)ourr ^— 179 = .orChapter 3. Bed Equations, Forcings, and Model Synthesis^ 52The initial and boundary conditions are dealt with by direct substitution. The only termin the compatibility of strains and force balance equations that needs to be transformedis ra/ar. The power-law transformed expression is derived as follows:R = rxso thatr = R-i.^ (3.43)Hencea^aRT— =ar^r Or aRa ( x) a7*-5-; r^aRrxrx-i aORxrxaORa= x R aRThe equivalent logarithmically transformed expression is derived as follows:(3.44)R = ln rso thatHencer = eR.a^aR ar Or ORa^ar-ar (ln r) aR=r aR(3.45)a(3.46)0RChapter 3. Bed Equations, Forcings, and Model Synthesis^ 53A convenient way to express the two transformation forms is to leta , a'where W is a weighting function given by(3.47)W = XR if R = rx^ (3.48)W = 1 if R = ln r. (3.49)Thus the transformed ice deformation problem is obtained by substituting (3.47) into(2.185) and (2.186), with (3.48) or (3.49), and substituting the definition of r into theinitial and boundary conditions. Hence compatibility of strains (2.185) becomesand force balance (2.186) becomesfe — — = 0WORacr+ cr — cre = O.aR(3.50)(3.51)The R = 1/r transformation of the bed problem is similar In the bed equations rappears only in initial condition (3.14), boundary conditions (3.15) and (3.16), and thefield equation (3.13)Oh _^10 ( ,ah)at pwg(cx+ ni3) r2 Or e. OrThe initial and boundary conditions are dealt with by direct substitution. In the fieldequation the only term to transform is1 a ( 20h)r2 -87r^Or •First note thatChapter 3. Bed Equations, Forcings, and Model Synthesis^ 54andOhOh OROr = OR Or_2 0hTOR'(3.52)Therefore1 a ( oh)-Or)1 0 r 21 _2ahm7-;ir 8 ( Oh\1 aR a ph)r2 8r OR 1:1)--1—r_2) 02hr2 OR2_4 02hr 01124 02hR OR2. (3.53)The transformed problem in the bed is therefore obtained by substituting (3.53) into(3.13) and substituting the definition of r into the initial and boundary conditions. Thefield equation (3.13) for head becomesOh _^-02hat — pwg (a + nfl)R4 OR2.^ (3.54)The borehole equations also require transformation. Equation (3.35) for the flow rateof water into the bed must be changed as it includes a spatial derivative term. Applyingthe R = 1/r transformation and using (3.52) to convert (3.35) yieldsQ. = —271-po exP (pb(t) — po)] r:K^ah)aRah= —27rpo exp [fi (pb(t) — Po)ir:K aRrcrc(3.55)ah= 2rpo exP (Pb(t) — po)) K oRChapter 3. Bed Equations, Forcings, and Model Synthesis^ 55R =Rb^BOREHOLE rARH^WALL^N1-1^2^M-1^MFig. 3.2: Ice finite-difference grid.3.5.2 Finite DifferencingThe next step in the method-of-lines formulation is to establish grids in the ice andbed, on which the partial differential equations in It and t can be replaced by systemsof ordinary differential equations in t only. The differentials with respect to R will bereplaced by finite-difference approximations.First consider the ice. The innermost node is at the It value corresponding to theborehole radius rb. I will let the outermost node be at Rmaz, the R value correspondingto some TmaI, to be be chosen later. This approach preserves flexibility, and retains theability to model the response of thick-walled cylinders. The remaining nodes will bespaced evenly (in R) between the inner boundary node and the outer boundary node,allowing simple definitions of finite-difference operators. As shown in Figure 3.2, thereare M nodes, identified by subscript i, with i = 1 at R Rb and i = M at R Rmaz. Atall nodes the differential equations (2.173), (2.155)—(2.157), and (2.164)—(2.166) apply,as well as the algebraic plane strain relation (3.50). At the interior nodes (i = 2, M — 1),force balance (3.51) also applies. At the borehole-wall boundary node (i = 1) the normal-stress condition (2.197) applies, unless a pressure forcing is applied in which caseR Rma),.—Pbor ^ (3.56)Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 56CAVITYWALLRa= RcrARal Ra=0N^•^I2^Ma-1 MaFig. 3.3: Bed finite-difference grid.where Pb is given by (3.41) for sudden pressurization, or by (2.5) for constant boreholepressure. At the outermost node (i = M) the radial-stress boundary condition (2.199)applies.The bed is discretized in a similar way to the ice, except that the outer boundary nodeis Ra = 0, corresponding to 7.. = co (Fig. 3.3). Henceforth the subscript, a will be usedto identify quantities pertaining to the basal grid, as opposed to the ice grid. As shownin Figure 3.3 there are M. nodes, identified by subscript j, with j = 1 at R. = 1/re andj = Ma at Ra = 0. At interior nodes (j = 2,M. —1) the groundwater flow relation (3.51)applies. At the cavity wall node (i = 1) the borehole-pressure boundary condition (3.16)applies. At the outermost node (i = M) the background head boundary condition (3.15)applies.Finite-difference approximations are now needed for all spatial derivatives in the iceand bed equations listed above, as well as the water mass flow expression (3.16) whichappears in the borehole equation. Examination of these equations reveals that finite-difference approximations are required for the first derivative at all ice nodes (innermost,interior, and outermost) and at the innermost basal node. Finite-difference approxima-tions of the second derivative are required at bed interior nodes. There are numerouspossible choices of finite-difference operators; I use second-order (3 point) operators forChapter 3. Bed Equations, Forcings, and Model Synthesis^ 57both first and second derivatives. Although accuracy could be improved by taking higher-order approximations, the stability of the resulting method of lines solution would bereduced [Schiesser, p. 110].A general expression for the finite-difference approximations of first and higher deriva-tives, along with tabulated coefficients, is given by Abramowitz and Stegun [1972, p. 914]die f(x)dxkk! m!hk E Aif('i) (3.57) where k is the order of differentiation, m is the order of the finite-difference approxima-tion, h is the step size (or node spacing), j is the node at which the derivative is beingapproximated, and Ai are the (tabulated) coefficients. From these, I obtain the followingexpressions for second-order finite differences:1. First derivative, centred difference (applies at interior nodes)- 2 A x2. First derivative, forward difference (applies at inner boundary node)dfdx(3.58)dfdx-fi+2 + 4fi+i - 3fi2 Ax(3.59)3. First derivative, backward difference (applies at outer boundary node)dfdx3f= --- 4A-1 -I- A-2 2 Ax(3.60)4. Second derivative, central difference (applies at interior nodes). (Second derivativesare never taken at boundary nodes.)d2 fdx2f=1-1. - (Ax)2(3.61)This completes the method-of-lines formulation of the coupled ice-bed-borehole prob-lem. The resulting system of equations is summarized in Appendix A.Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 583.6 Computer Programming of SolutionThe final step in the problem solution is the solution of the system of equations of Ap-pendix A, using a numerical integration routine. I have chosen ddriv3 of Kahaner andSutherland (public domain software, available from Los Alamos National Laboratory),an implicit integrator which uses the Backward Differentiation Formula method, an ex-tension of the Euler method [Schiesser, 1991, p. 65]. ddriv3 is similar to the stiffly stableintegrator sdriv2 described in Kahaner et al. [1989, pp. 295-297, 341-346], except thatit is double precision and capable of handling systems of differential and algebraic equa-tions.The solution has been implemented in the program bfi.I.4 (Borehole Flow Infinitesi-mal, 14th version), written in Sun Fortran77. An earlier version of the program, bfil3,was used for the tests of Chapter 4. It contained bugs which only manifested themselveswhen the background pressure Po 0 0. The tests of Chapter 4 are valid because all thetests were done with Po = 0.The program allows the user to set many parameters at run time including:1. Forcing type can be step pressurization (creep test), sudden pressurization followedby free relaxation, sudden pressurization followed by freezing forcing, or freezingforcing alone.2. Problem geometry including borehole radius, borehole height as a function of time,and basal-cavity radius.3. Ice properties including elastic constants, steady-state viscous parameters, andtransient viscous parameters.4. Ice temperature at each node.Chapter 3. Bed Equations, Forcings, and Model Synthesis^ 595. Basal hydraulic properties including hydraulic conductivity and compressibility+ 743.6. Water compressibility.7. Grid parameters including number of ice and bed nodes, and spatial transformationin the ice.8. Initial conditions including background pressure and initial ice scalar back stress.9. Pressure forcing including forcing pressure and ramp-up time.10. Output control parameters including frequency of output, optional output of icestresses, basal hydraulic heads, and debug output.11. ddriv3 control parametersParameters are communicated to the program via an input file.The program can be run with the transient viscous part of the rheology turned off(leaving a power-law viscoelastic rheology), or the elastic part turned off (leaving a tran-sient viscous theology), or with both the elastic and transient viscous parts turned off(leaving a Glen-Law viscous theology). An effectively linear elastic rheology can beachieved by suitable choice of viscous parameters (very stiff). The ice-borehole systemcan be modelled independently of the basal system by setting hydraulic conductivity andcompressibility a ± nig to zero. The bed-borehole system (i.e., ice assumed rigid) can bemodelled by suitable choice of ice parameters (very stiff). The outputs of the programare: a summary of the input parameters ( . sum file), time series of the field variables( . out file), and optional debugging information ( .dbg file), all stress components at allice nodes at all time steps ( str file), and the head at all aquifer nodes at all time steps( .hed file).Chapter 4Model Testing4.1 IntroductionAny computer model must be tested to ensure that its results are correct. This isdone by comparing model output to available known solutions (usually analytic), andby qualitative assessment when no know solutions exist. The performance of a method-of-lines solution is limited by both the accuracy and stability of the solution. In generalaccuracy is increased by decreasing the grid spacing (increasing the number of nodes),increasing the finite-difference approximation order, increasing the integration-methodorder, and decreasing the temporal step size [Schiesser, 1991, pp. 19-31]. In this problemthe finite-difference approximation order has already been fixed, so that only the gridspacing, integration-method order, and temporal step size remain to be varied. In ddriv3,as in many other integrators, the temporal step size is adjusted automatically to keep theestimated error at each time step within a user-specified tolerance. Hence in this casethe error tolerance replaces the temporal step size as the user-specified parameter. Stiffproblems may also require the integrator to take very small time steps to meet the errortolerance. Stability can also be a problem, causing the model solution to diverge from thetrue solution, or causing*temporal numerical oscillation [Schiesser, 1991, pp. 122-129].The model testing reported on here has two purposes: to determine the control param-eter values which give optimum stability and accuracy, and to determine if the optimumsolutions are acceptably close to known solutions. These questions are first addressed60Chapter 4. Model Testing^ 61for the ice part of the model, then for the bed part, and last, in a qualitative way, forthe coupled model. The ice and bed models are tested and results compared to ana-lytic solutions. It is shown that although numerical oscillation of the ice solution can bea problem and stiffness limits accuracy for certain choices of ice properties, acceptablesolutions can in general be obtained with the correct choice of input parameters.4.2 Ice Deformation ModelThe ice model was tested by running creep tests and systematically varying input pa-rameters. These parameters include: ddriv3 error control parameters, spatial coordinatetransformation, ratio of maximum radius to borehole radius rmazIrb, and number of nodesM. The effect on the Glen-Law viscous and linear elastic solutions was investigated. Thevalues of ice parameters used for all tests in this chapter are listed in Table Glen-Law Viscous IceThese tests were run by modelling one-year creep tests, with transient and elastic partsof the rheology disabled. First I ran a number of tests to determine the optimum ddriv3error control parameters. The ddriv3 integrator can be run using several different one-step error control strategies: absolute error control, relative error control, or relativeerror control with user-specified definitions of "problem zero". The last case is identicalto relative error control, except that any dependent variable having a value less thanits problem zero is considered equal to zero when its relative error is computed. Thishas the effect of relaxing the error tolerance on variables having small absolute values.This option is used to reduce the computing time that would be wasted in accuratelycomputing numbers which are essentially zero. The integrator allows the problem zeros tobe set separately for each variable. I have taken advantage of this capability because theChapter 4. Model Testing^ 62Parameter Value SourceVo 6590 Pa s113 Shyam Sunder and Wu [1989a]estimated from data of Sinha [1978]A 0.0142 Shyam _Sunder and Wu [1989a]fitted from data ofJacka [1984]H 0.02 Shyam Sunder and Wu [1989a]fitted from data of Jacka [1984]0.286 Shyam Sunder and Wu [1989a]fitted from data of Jacka [1984]Rga, 8.314 J mo1-1K-1 Shyam Sunder and Wu [1989a]QL 67000J mori Paterson [1981, p. 34]QH 139000J mol-1 Paterson [1981, P. 34]it 3.3005 x 109Pa calculated from seismic velocitiesof Paterson [1981, p. 78]A 6.3608 x 109Pa calculated from seismic velocitiesof Paterson [1981, p. 78]ice temperature 0 °Cwater temperature 0 °CPw 999.9 kg m-3 Streeter and Wylie [1979, p. 54]13 4.4 x 10-1°Pa-1 Streeter and Wylie [1979, p. 54]Table 4.1: Standard rheological parameters used for model testing.Chapter 4. Model Testing^ 63Parameter bf i 13 Name ValueMaximum integration order mxord 5Error control method specifier ierror 4Relative error tolerance eps _ 1 x 10- 2Problem zero for m,,,q, B, R, and a ewt (1 ,3-6 ,10 ,11) 1 x 10-6Problem zero for it ewt (2) 1 x 10-10Problem zero for et' ewt (7-9) 1 x 10-40Table 4.2: Optimum ddriv3 error control parameters.dependent variables have widely differing characteristic magnitudes (for instance stresseshave typical magnitude 100-104, whereas viscous strains have typical magnitude i0-10-s.) It was found that the error tolerance (eps) and problem zeros (ewts) were limitedby the step size required to keep the error within the tolerance. ddriv3 terminates withan error message whenever the step size has to be reduced more than 50 times in theattempt to take one time step, or if the time step size falls to less than the roundofferror in time. All step size errors observed during the model testing occurred during thepressure ramp-up. It was found that the solution changed by less than one part in 107when eps was changed 1 x 10' to 1 x 10. The strictest parameters that could be rununder most conditions are summarized in Table 4.2.It was found that the Glen-Law ice deformation problem is very stiff; many combina-tions of parameters resulted in problems so stiff that a step size error resulted, even withthe loose error tolerance of eps = 1 x In contrast, when elasticity was included alltrial runs with eps = 1 x 10-5 succeeded.The ratio of maximum node radius to borehole radius (r,„,,z/rb) has a large effecton how well the model solution approximates the true solution for infinite ice. This isChapter 4. Model Testing^ 64because bfil3 really solves the problem of flow in a thick-walled cylinder, which hasa solution different to that of infinite ice, especially for small (rmazIrb). Hence even ifthe numerical solution were perfect it could not be expected to do any better than theanalytic solution to the thick walled cylinder problem.In this section I focus on tangential strain at the borehole wall (comb) as a measure ofsolution accuracy. It is through tangential strain (and radial stress) at the borehole wallthat the ice couples to the water. An analytic solution for steady-state power-law creepin a thick-walled cylinder is given by Finnie and Heller [1959, pp. 182-187]) 3dee = A (3 \ (N+1)/2 Pb ^2^(rmazdt^4)^(r,,,azirb)2/N — 1 N^rwhere the flow law of ice isin uniaxia1 compression. For the case of N = 3 considered here (4.1) becomesdee^A 3 (^1^rmaz NI 2dt = -6Ps (r„„„Irb)211s1 - 1^r )wherede^.1= Ace.dtThis is the desired solution, except for the coefficient A defined in (4.3). Comparison of(4.3) to the steady-state viscous part of the Shyam Sunder and Wu rheology allows A tobe expressed in terms of V, and yield the desired result. Recalling (2.65) and (2.72)A,Cri (c)dt 03 ( 1 N N-ver^cr.^2 V) eq^dt )) 3 () 2(4.1)(4.2)(4.3)dt(4.4)Chapter 4. Model TestingWhen stress is uniaxial (in x1 direction)crij 65(TOO000000(4.5)and lajj is given by (2.45)(2/3)cr 0 0— 0 — (1/3)a 00 0 —(1/3)crWith (2.51) for the definition of effective stress2^26req = 6rHence for uniaxial compressiondes^3 (1)3 2,^(de)dt v a all Tit 0de'^3 ( 1 \3 22 1de)dt a i7Vilt)0= (1\3 3(d6\Comparing (4.8) to (4.3) reveals that A = (del dt)o /V3. Hence (4.2) becomes) 3 (rmaz)21deo^pl  (de (dt^6V3 dt )0 k (r,./Tb)2/N — 1 )^r )and at the borehole wall(4.6)(4.7)(4.8)(4.9)374, ( de) ( ^1^(7-max 2)6V3^) — 1)^rb ) •deodt (4.10)For infinite outer radiusdeodt/d\6V3 dt ) 0rb(4.11)Chapter 4. Model Testing^ 66E 20.E180z 16CC1—(f) 140- 12_J• 0L.L.101.0^2.0ANALYTIC SOLUTIONBFI13 SOLUTION3.0^4.0log(rmax/rb)000Fig. 4.1: Effect of maximum radius on borehole-wall tangential strain in Glen-Law viscousice, after 1 yr creep test with Pb = 10 kPa, M = 9, and R 7.-213. The thick walledcylinder solution differs from the infinite solution by 107% at r,.../rb = 10, 15% at 100,3.1% at 1000, 1.0% at 5000, and 0.65% at 104. The bf 13 solution differs from the thickwalled cylinder solution by —1.6% at rmaz /rb = 10, —6.9% at 100, —10% at 1000, and—11% at 5000.Chapter 4. Model Testing^ 67Model simulations and analytic solutions were computed for a range of rmaxlrb values(Fig. 4.1). The thick walled cylinder solution differs significantly from the infinite-iceone; the difference for rmazIrb = 1000 is 3.1%. The model results differ significantlyfrom the thick walled cylinder solution due to error introduced by finite differencing.This error increases as rmazirb is increased, because the node spacing is increased. Nomodel output was obtained for rmazIrb > 104 as time step errors occurred. It appearsthat extending the grid to larger radius increases the stiffness of the problem. Sincethe thick walled cylinder solution overestimates the true strain, and the model solutionunderestimates the thick walled cylinder solution, there is an optimum number of nodesfor which the model solution closely approximates the true one (at least in terms of theborehole-wall strain). This point is between 100 and 1000, for R = r-213. I thereforechose rmazIrb = 1000 for the balance of the testing.The choice of radial coordinate transformation has a large effect (Table 4.3, andFig. 4.2). Transformations with negative powers of two or more could not be run with-out causing step size errors. Significant spatial oscillations, of wavelength 2 AR, are mostnoticable in the stress solution for R r and R In r. Similar oscillations exist in theother solutions, but are of smaller amplitude. Similar numerical oscillations can appear insolutions of the one-dimensional advection equation and the one-dimensional wave equa-tion, when centred spatial finite-difference operators are used [Mesinger and Arakawa,1976, pp. 26-30, 43-44]. These parasitic waves can in some cases be eliminated by usingun-centred finite-difference approximations. This may be true of the current problem aswell. Nevertheless, this model generates adequate solutions for transformations on theorder of R = r-2/3. It appears that, at least for rmaz/rb = 1000, that R = 7-31' is theoptimum transformation. The error is a reasonable -0.96% for 9 nodes.The number of nodes also affects the accuracy of the solution (Fig. 4.3). As expectedthe error drops as the number of nodes is increased. No step size errors were encountereda^ ANALYTIC SOLUTION* BFI13 SOLUTION, R = rBFI13 SOLUTION, R = ln(r)111111111- 00a.642(1) 0(Y —2—4_J< —6<Ct—8—100 —2—4LU—6_Jcr—8—100^200^400^600^800NON—DIMENSIONAL RADIUSANALYTIC SOLUTION* BFI13 SOLUTION, R = r-2/3BFI13 SOLUTION, R =7.'175BFI13 SOLUTION, R =Chapter 4. Model Testing^ 684^8^12^16^20NON—DIMENSIONAL RADIUSFig. 4.2: Glen-Law viscous stress distributions for various radial coordinate transforma-tions. M = 9 and r,„./rb = 1000. The outermost nodes are not plotted. Analyticsolutions are those for infinite ice.Chapter 4. Model Testing^ 69Fig. 4.3: Effect of number of nodes on borehole-wall tangential strain in Glen-Law ice.R = 7.-2/3 and rmazirb = 1000. The relative errors in the bfil3 solutions are 15% for 5nodes, 7.3% for 9, and 2.1% for 17.Chapter 4. Model Testing^ 70Transformation (R .) fe,ri. Relative error (%)r 5.477253 x 10' —99.99ln r 4.949487 x 10' —47r-2/3 8.710883 x-10-6 —7.3r 07 8.938541 x 10-6 —4.8r 075 9.307882 x 10-6 —0.96r 9.709685 x 10' 3.3r-1.0 1.161782 x 10 24Analytic solution (infinite ice) 9.398034 x 10' —Table 4.3: Dependence of borehole-wall tangential strain in Glen-Law ice on radial coor-dinate transformation.with R 7.213. However with R 7.-314 it was found that increasing M> 13 causedstep size errors. Thus increasing the number of nodes increases the problem stiffness inagreement with the observations of Schiesser [1991, p. 22].After consideration of the Glen-Law ice modelling, I chose as optimum the parametervalues listed in Table 4.4. The value of M can be chosen as large as possible, subject tothe stiffness constraint discussed earlier. These parameters lead to a Glen-Law solutionaccurate to within 1%. It is apparent from the results presented in this subsection thatthe errors caused by finite outer radius and the finite-difference approximation are fargreater than those introduced by the numerical integration, for transformations in therange of r-3/4 < R < r-213.2 (A + p)(r,27,az - rb2)(rm2ax - r)r •Pbrt2,7*= pbri2,r,n2 ax (4.12)Chapter 4. Model Testing^ 71Parameter bf i 13 Name ValueMaximum integration order mxord 5Error control method specifier i error 4Relative error tolerance eps 1 x 10-2Problem zero for maq, B, R, and a ewt (1 , 3-6 ,10,11) 1 x 10-6Problem zero for It ewt (2) 1 x 10'Problem zero for ev ewt (7-9) 1 x 10'Ratio of rmaz to rb r_max/r_b 1000Number of nodes M 9+Table 4.4: Optimum parameters for Glen-Law ice.4.2.2 Linear Elastic IceThe elastic behaviour of the model was tested by running it with power-law viscoelasticice (i.e., transient part of rheology disabled), and observing the solution at the end ofthe pressure ramp-up. The parameters of Table 4.1 , a pressure ramp-up to 10 kPa in 1 s,and a background pressure po = 0 were used. The solution after 1 s was taken to be theelastic solution; of a total borehole-wall tangential strain of 10-5, only 10-14 is viscous.The ddriv error control parameters of Table 4.2 were used, except that eps =This was done because the error tolerance could be tightened without causing step sizeerrors. No step size errors ever occurred in these tests, indicating that the viscoelasticproblem is much less stiff than the purely viscous one.First I assess the effect of the maximum radius on the solution. Jaeger and Cook[1969, p. 128] give the solution for a linearly elastic thick-walled cylinder undergoingboth internal and external pressure. For internal pressure only it reduces toChapter 4. Model Testing^ 72At the borehole wall2Pbrl ^Pbrbrmaz u,.(rb) =2 (A + p) (rm2 22 — 71) 2p, (rm2 ax —Pbrb ^I  r^+r„.42.z).2(r2 —^+ ft)^ttUsing the definition of tangential strain (2.188) this becomesPb^r2 2)2 (rLax — r)2 (A + it)Pb ^( 14. (rmaz/rb)2)2 ((rmaz/rb)2 — 1) k(A A)For infinite ice (4.14) reduces toee ,rb(4.13)(4.14)667,rb = 211• (4.15)Model simulations and analytic solutions were computed for a range of rmax/rb values(Fig. 4.4). The thick walled cylinder solution differs from the infinite-ice solution, butnot as much as in the Glen-Law case. Hence limiting r,,,./rb does not degrade the elasticsolution as severely as it does the Glen-Law one. Thus the r„,„z/rb = 1000 chosen duringthe Glen-Law testing is also acceptable as far as the elastic solution is concerned.The choice of radial coordinate transformation has a large effect on the solution(Table 4.5, Figure 4.5), as was the case for Glen-Law ice. The elastic solution is quiteaccurate for a wide range of transformations. The Glen-Law optimum transformation(R = 7-'14) results in a modest —0.01% error in the elastic solution, which is quiteacceptable.The effect of the number of nodes on solution accuracy was tested (Fig. 4.6). Asexpected, the error decreases as the number of nodes is increased. The error for M = 9 is—0.03%, again better than the Glen-Law solution with the same parameters. No stiffnesslimit on M was found; thus there is no reason not to use the maximum allowable numberof nodes (M = 20), except if computation time is a concern.0 ANALYTIC SOLUTIONBFI13 SOLUTION1.0^2.0^3.0^4.0^00log(rmairb)Chapter 4. Model Testing^ 731.5351.5301.525_1 1.5201.1.1z 1.515F-Fig. 4.4: Effect of maximum radius on borehole-wall tangential strain in linear elasticice. M = 9 and R = 7.-2/3. The thick walled cylinder solution differs from the infinitesolution by 1.4% at rmaz I rb = 10, 0.01% at 100, 0.0001% at 1000, and 0.000001% at10000. The bf 113 solution differs from the thick walled cylinder solution by —0.79% atrmax/rb = 10, —0.002% at 100, 0.007% at 1000, and 0.007% at 104.8—2Chapter 4. Model Testing^ 74-_-.-_--e+ +0 o+o a0---IF +0IIIIIIIII0.00 R = r-4- R = ln(r)0^200^400 600 800NON—DIMENSIONAL RADIUS4^8^12^16^20NON—DIMENSIONAL RADIUS1.2^1.6^2.0^2.4^2.8NON—DIMENSIONAL RADIUSFig. 4.5: Linear elastic stress distributions for various radial coordinate transformations.M = 9 and rmax/rb = 1000. The solid lines are the infinite-ice analytic solutions. Theoutermost node is omitted for clarity.Chapter 4. Model Testing^ 75--,E 1.5150E 1.5148010 1.5146,....1.5144z•TC 1.51421---v)o 1.5140rz2___J 1.5136<i=z 1.5134LJJ(2) 1.5132<H-ANALYTIC SOLUTION 0BFI13 SOLUTIONS^L......004^6^8^10^12^14^16^18 20NUMBER OF NODESFig. 4.6: Effect of number of nodes on borehole-wall tangential strain in linear elastic ice.R = r-2/3 and rmaxlrb = 1000. The relative errors in the bfil3 solutions are: —0.12%for M = 5, —0.03% for 9, and —0.007% for 17.Chapter 4. Model Testing^ 76Transformation (R =---) E8,rb Relative error (%)r —4.893174 x 10-7 —132ln r 1.188106 x 10-5 —21r-2/3 1.514491 x 10-6 —0.03r-3/4 1.514727 x 10-6 —0.01r-1 1.514925 x 10-6 0.0002T-2 1.514925 x 10-5 0.0002r-3 1.514923 x 10-5 0.00007Analytic solution (infinite ice) 1.514922 x 10-5 —Table 4.5: Dependence of borehole-wall tangential strain in linear elastic ice on radialcoordinate transformation.It appears, then, that the viscoelastic problem is better behaved than the Glen-Lawone. It is less stiff. The elastic component of the solution is less affected by the choiceof rmax I rb and suffers less finite-differencing error, compared to the Glen-Law viscoussolution. I find the optimum parameters determined for the Glen-Law problem to beacceptable for the elastic problem, except that eps = 1 x 10-5 be used, and that up toM = 20 nodes be used.4.3 Basal Groundwater Flow ModelTo test the accuracy of the basal part of the model, a step pressure change was simulatedusing bfil3; the number of nodes was varied and the results compared to the analyticsolution. The analytic solution can be obtained, by analogy, from the solution to theproblem of heat conduction in a whole space with imposed step temperature change atChapter 4. Model Testing^ 77a spherical inner boundary given by Carslaw and Jaeger [1959, p. 247]. The solution toOh^(1 a (20h= D^r -(w.)), r re^ (4.16)h(r,O) =h(co,t) =h(re,t) =0, r > rc0, t > 0hforce = Pforcel PRIM t > 0(4.17)(4.18)(4.19)isreit 4erh(r, t) = ^' " erfc (2r r^e)-Vr3t(4.20)where D = K I pwg(cr 70).To correspond with the analytic solution, bii.I.3 was run with po = 0 and values ofrc, K, and D representative of Trapridge Glacier used (K and D values from T. Murray,unpublished data, 1992). The rheological parameters of Table 4.2 and ddriv3 parametersof Table 4.3 were used, except that the relative error parameter was set to eps = 1 x 10-5.It was found that the accuracy of the solution increases with increasing number of nodes,as expected (Fig. 4.7). At ra = 0.129 (equivalent node 3) the relative errors in head are:Ma = 5 nodes, 4.9%; Ma = 9, 3.4%; and Ma = 17, 0.9%. Even with just 10 nodes thetime-dependent solution is well reproduced (Fig. 4.8). Note that the model solution atfirst diverges and later converges to the analytic one. No stability problems are evidentfor M < 20. It appears that 9 nodes provide a sufficiently accurate solution for mostpurposes, although more nodes would give greater accuracy, especially at early time.4.4 Coupled ModelThe coupled ice-water-bed model was tested quantitatively as well as qualitatively.Quantitative testing was restricted to the case of sudden pressurization of a sealed bore-hole, and subsequent pressure decay, with Glen-Law ice. This is the only coupled problem1.0.8ANALYTIC SOLUTION5 NODES9 NODES^BFI13 SOLUTIONS17 NODES.2.0Chapter 4. Model Testing^ 781^2^3^4^5^6^7^8^9EQUIVALENT NODEFig. 4.7: A comparison of model predictions with analytic solution for bed response to aborehole pressure step change for various numbers of nodes at t = 12h. The equivalentnode is that which has the same rlrb as for the 9-node case.1.0.8E2 .4=.2. 0Chapter 4. Model Testing^ 791^2^3^4^5^6^7^8^9NODEFig. 4.8: Comparison of model predictions and analytic solution for bed response toborehole pressure step change at various times. An essentially steady state has beenreached after 100yr).Chapter 4. Model Testing^ 80for which an analytic solution was found. Next the effect of adding the elastic and tran-sient viscous parts of the rheology was checked by simulating sealed-borehole suddenpressurization/decays. The results can be explained qualitatively. Lastly the coupledmodel was checked qualitatively by comparing sudden pressurization/decays for variouscombinations of ice and bed properties.4.4.1 Sealed Borehole with Glen-Law IceThe coupled ice—water model can be checked by simulating the sudden pressurization ofa sealed borehole and its subsequent decay for Glen-Law ice, and comparing results tothe analytic solution. The analytic solution is derived from the strain-rate solution ofSection 4.2 and borehole mass conservation equation of Chapter 3 as follows.Consider a cylindrical borehole with no cavity at the bottom. The borehole is short,so that water pressure (and hence density) is uniform throughout. For a sealed boreholewith no freezing, the mass of water is constant and(4.21)dt =0.At any given time the mass of water ismw(t) = p(t) dV(t)wheredmw(4.22)dV (t) = 7r7.(t) dz^ (4.23)andPw(t) = Po exP [0 (Pb(t) —^ (4.24)Substituting (4.23) and (4.24) into (4.22),rnw (t) Po exp [f3 (pb(t) — po)J Irq,(t) dzirpor(t)L exp fi (pb(t) —^. (4.25)Chapter 4. Model Testing^ 81Differentiating (4.25) with respect to t= 1rPo(t)L 0(43 exp [ig (pb(t) — RI)] ddPit'-1-2rb(t)71drib exp v3 (pb(t) — N)]}rPo(t)Lrb(t)exp [s (pb(t) — po)] (orb(t)ic 1 b 2drb)dt^dt ) •Conservation of mass demands that this be equal to zero, thusdpb^drb0 = igrb(t)— 4- 2dt^dt *dmwdt (4.26)(4.27)(4.28)Now recall the analytic solution for the strain rate of Glen-Law ice around a cylindricalhole in infinite ice (4.11)deo ^p(t) 'de)dt^6V3 dt )rbUsing the definition of tangential strain (2.188),drb^r5(t)p(t) (de)dt^6V3^kdt^) 0*(4.29)(4.30)Substituting (4.30) into (4.28) and rearranging,01@rbtIdpb clrb(t)pg(t) (dc) dt^6V3^dt ) 0dPb p(t) ( de)0dt^3V3 kdt ) 0dpb^p(t) (de"dt^3i3V3 dt )^ (4.31)01^ ( de)— p1,3(t)dpb =3f3V3^dt (4.32)0Integrating (4.32) from 0 to t (assuming instantaneous pressurization at t^0) withPb(0) — pforce, and simplifying, givesf Pb(t) _p_3 dp^1  (dc) ftPforce^3/3V3 dt ) 0^dr (4.33)0Chapter 4.^Model TestingP-2pb(t) t^(d\82(4.34)2 3)3V3 kdt ) 0Pforcep2(t)^Pior2 „ t^(de)(4.35)2^2 30V3 kdt ) 02t^( de\TOO^Pf—or2 ce (4.36)33y3 kdt ) 02t^( de\131,2(0= PiLe (4.37)3,QV3 kdtpb(t)^(601.2 ce 2t^(d€\ \- 1/2(4.38)3/3V3 ( dt)0)The program bfil3 was run with the optimum ice model parameters, rc = 0.001 m,and L = 0.01 m, to correspond to the assumptions of the analytic solution. The resultsvery closely approximate the analytic solution (Fig. 4.9). For 9 nodes the error in pressureafter 1 yr is 0.4%, and for 13 nodes -0.04%.4.4.2 Effect of Viscoelastic and Transient Viscoelastic IceSo far only the steady-state viscous and pure elastic behaviour of the ice model have beenchecked. The behaviour of viscoelastic and transient viscoelastic ice must be assessedqualitatively, as no analytic solutions exist. Three sealed borehole sudden pressuriza-tion/decay simulations were run with identical parameters: one with steady-state viscousice, one with power-law viscoelastic ice, and one with transient viscoelastic ice. Attemptswere also made to run the same test with transient viscous ice (elasticity turned off), butstep size errors occurred, even when the error tolerance was relaxed (eps = 1 x 10-i,ewt(1,3-6,10,11) = lx 10-2, ewt (2) = lx 10-5, and ewt (7-9) = lx 10-20). I do notconsider a solution obtained using a one-step error tolerance any greater than this (10%!)to be useful. It therefore appears that the transient viscous problem is even stiffer thanthe steady-state viscous one.Chapter 4. Model Testing^ 8310a_ 9w 8t27cc)1.1.1a. 6LiJo5LU40cri30^50^100^150^200^250^300^350TIME (DAYS)Fig. 4.9: Sudden pressurization and relaxation of a sealed borehole surrounded byGlen-Law ice.GLEN LAW VISCOUS— POWER LAW VISCOELASTIC----- TRANSIENT VISCO ELASTIC------------------Chapter 4. Model Testing^ 840^50^100^150^200^250^300^350TIME (DAYS)Fig. 4.10: Effect of ice rheology on pressure decay in a sealed borehole.It was found that the power-law viscous solution decays faster than the Glen-Lawviscous one (Fig. 4.10). This appears counter-intuitive at first. One might expect that asviscous creep causes pressure to drop, the ice would elastically relax, reducing boreholeradius; hence pressure would remain high for longer than would happen for Glen-Lawice. This is in fact what happens for the flow solution with linear viscoelastic ice arounda borehole and for power-law viscoelastic ice in uniaxial compression. The importantdifference between those two cases and the thesis problem is that they involve uniformstress distributions; the ratio cr,./cro/az is the same everywhere, as is ce/ev. Hence in thesecases the total strain solutions can be obtained by linear superposition of the elastic andviscous solutions. But superposition fails for the deformation of power-law viscoelasticice around a borehole. The stress solutions for linear elastic ice and power-law viscous iceChapter 4. Model Testing^ 85are different; hence total strain is not equal to the sum of the strains found from separateelastic and viscous solutions. For a given loading the stress and strain state depends onthe loading history. The initial state of stress is purely elastic, evolving later toward theviscous stress distribution [Murat et al., 19861. Hence intuition based on superpositionof elastic and viscous solutions fails, and cannot be used to predict the relative pressuredecay rates in the borehole. Thus the observed result is not ruled out, although I cannotexplain it conceptually.The accelerated pressure decay with transient viscoelastic ice compared to that withpower-law viscoelastic ice is consistent with the fact that transient creep is faster thansteady-state creep.4.4.3 Ice—Water—Bed ModelThis test was done to verify that the fully coupled model behaves as expected, in aqualitative sense. Sudden pressurization and decay simulations were run for the followingfour cases: deformable ice and impermeable bed, deformable ice and permeable bed (bothwith transient viscoelastic ice theology), Glen-Law deformable ice and permeable bed,and rigid ice and permeable bed. As expected, pressure decays much more quickly whenthe bed is permeable than when it is impermeable (Fig. 4.11). The situation is not soclear-cut for the case of deformable ice. Glen-Law ice results in a very slightly fasterdecay than rigid ice (by an amount imperceptible in Figure 4.11), which is expected.With transient viscoelastic ice the pressure decays first more quickly then more slowly,compared to that with rigid ice. Initially rapid transient creep allows the ice to relaxquickly, relieving the pressure on the water. Later, as transient creep slows down, iceelasticity and recovery of transient creep dominate, and slow pressure decay.(a) TRANSIENT VISCOELASTIC ICE, IMPERMEABLE BED'(b) TRANSIENT VISCOELASTIC ICE, PERMEABLE BED(c) GLEN LAW ICE, PERMEABLE BEDd) RIGID ICE, PERMEABLE BED\_\-^a\i x\-^\ -. bN --,c,d ,---,„_Chapter 4. Model Testing^ 86100a_8uJ(/)(1)6uJa_LU 40LU02.0^.2^.4^.6^.8^1 .0TIME (DAYS)Fig. 4.11: Pressure decay for various ice/bed combinations. The permeable bed hasK = 1 x 10-12 ms-1 and D = 1 x 10-10 m2 s-1Chapter 4. Model Testing^ 87Parameter bfii3 Name ValueMaximum integration order mxord 5Error control method specifier ierror 4Relative error tolerance eps 1 x 10-2Problem zero for maq, B, R, and a ewt (1,3-6,10,11) lx 10-6Problem zero for h. ewt (2) 1 x 10-6Problem zero for ev ewt (7-9) 1 x 10-40Ratio of 7.,„„^to rb r_max/r_b 1000Number of nodes (ice) M 9+Number of nodes (bed) M_a 20Table 4.6: Optimum parameters for coupled problem.4.5 ConclusionsIt is concluded that the model adequately simulates the problem. The elastic solutioncan be reproduced to within 0.01%. The Glen-Law viscous solution can be reproducedto within 1%. The bed solution can be reproduced to within 5%. The behaviour of thecoupled system is qualitatively correct. The optimum control parameters are summarizedin Table 4.6.Chapter 5Model Results5.1 IntroductionIn this chapter the bfil4 computer model is used to infer Trapridge Glacier ice and bedproperties. Of primary interest are the basal hydraulic properties of hydraulic conduc-tivity and bed compressibility. These are determined by using the model to interpretthe observations obtained from three boreholes drilled in the summer of 1992, as well asthree boreholes from earlier years. In 1992 two boreholes were drilled 5 m apart on thesame day. The first, 921124, was an unconnected hole. The second, 921126, was a blindhole. Both holes were instrumented with a pressure sensor and a thermistor string. Anal-ysis of the thermistor data provides information on the borehole-length function L(t).The pressure sensors provide pressure time series which can be fit by iterative forwardmodelling using the model developed in Chapters 2 and 3. First the ice properties areinferred from forward modelling of the 921126 pressure record, then the bed propertiesare inferred from forward modelling of the 921124 pressure record using the ice propertiesinferred in the first step. The sensitivity of the resulting basal parameters to changesin some of the more uncertain input parameters (such as ice parameters and boreholefreezing forcing) are tested by repeating the fitting with altered parameters. In additionthree other unconnected holes from other years (891150, 901154, and 911146) are modelledusing the 1992 freezing forcing and ice parameters determined from 921126, and a secondblind hole (921145) is modelled.88Chapter 5. Model Results^ 89Finally, the model is used to answer the question of to what extent the pressurerecords of unconnected holes reflect the "true" basal conditions, and to what extent theyare altered by the response of the borehole itself. This is done by examining the responseof the model to sudden pressurization, with the ice and bed properties obtained for 921E24and 921126.5.2 Freezing Forcing L(t) Obtained from Thermistor RecordsBoreholes through Trapridge Glacier are drilled with a hot-water drill. In all years except1992 the holes were drilled to their final diameter in one pass with the drill in roughlyone hour of drilling. In 1992, however, the thermal efficiency of the drill was impaired,so that the hole after one pass was too small, and a second pass was required. The firstpass took about one hour, and the second about 10 minutes.In 1992 two thermistor strings were built and installed: 92T01 in unconnected hole921124, and 92T02 in blind hole 921126. The thermistors were calibrated in 1990 byP. Langevin and zero readings checked in 1992 in an ice bath before cabling. All ther-mistors were stable, although the calculated temperatures were between —0.03 °C and—0.25°C rather than the expected 0 °C. For this study the absolute temperature isnot critical, so this is acceptable. The thermistors were calibrated and strings assem-bled as described in Jarvis [1973], except that Thermometrics P10 DA 402L thermistorswere used. Measurements were taken using a Fluke 8060A digital multimeter, using thetwo-polarity method of Clarke and Blake [1991] to average out spurious electro-chemicalpotentials. Readings were taken once a day for nine days following installation. Leadresistance was calculated and subtracted from the readings before the calibrations wereapplied.Trapridge Glacier is coldest near the surface^—5 °C at the surface if the winterChapter 5. Model Results^ 90cold wave is neglected), and below freezing throughout except at the bed where it is atthe pressure melting point [Clarke and Blake, 1991]. I therefore expect that a boreholefreezes shut first at the highest point with water in it, then downward from there. Thetime of the first freeze-over can be found from the borehole water-pressure record (It isthe time when the pressure begins to rise, usually about 12h after drilling). Hole lengththen decreases with time, asymptotically approaching zero, as the freezing time at thebed where ice is at its melting point is essentially infinite (neglecting melting due togeothermal or frictional heating).The temperature versus time after cessation of drilling was plotted for the thermistorsof the two strings (Figs. 5.1 and 5.2). For each thermistor a smooth curve was extrap-olated back until it intersected the initial temperature (measured soon after installationwhen all thermistors were in water); this was taken to be the time when the thermistorfroze in (Table 5.1). If freezing is assumed to be from the top down only, as in the model,then thermistor position versus freeze-in time gives the borehole length as a function oftime required by the model.The records from the two holes were then stacked together to reduce noise. I assumethat because the holes are nearby the ice surrounding both holes has the same tempera-ture stratification, and furthermore that because the drilling method was approximatelythe same, the two boreholes should freeze at the same rate. A point was added at timezero, corresponding to the initial water level in 921124. This was done under the assump-tion that freezing would happen first at the water surface. Also a point was added atlong time (L = 0.5 m at t = 2 yr). Pressure sensors located 0.5m above the bed appearto remain in water after one year, but be isolated by freezing after two years (Fig. 5.3).Therefore the lengths and times from both thermistor strings, and two extra pointsat zero and two years (weighted 10 times and 20 times respectively to reflect greaterconfidence in these points, and force the solution away from local minima) were fitted■■•-x- -x-^-x^8.11^1^1^1^1^1^1^11^1^1^11^I^1—\ N\ \\ \ \\ ■\\^ \\ \\^\ \^\ \\ \ \ \\^\ \\^\\\ \ \\\^ \\38.1THERMISTORDEPTH (m)^-x^48.1-x- -x\28.1-x^m-x^23.1■-x- --x --*-x^18.1HOLE FIRSTFREEZES SHUT13.1••■•■• ..1.0 NEN100 101ELAPSED TIME SINCE HOLE DRILLED (DAYS)Fig. 5.1: 921124 freeze-in temperature record (92T01). Solid lines are smooth extrap-olations used to find time of freeze-in. Thermistors at 58.1m and 68.1m displayed notemperature drop and are omitted for clarity.01—67THERMISTORDEPTH (m)45.9-x--K 35.9ik\25.9X- )4- 15.9--- \)4• +(- 10.9-9(I1-1 -41-4r4i 5.9$k\-• —\\\\\\\\ \\\ N\ \\ \tin OMB)4-HOLE FIRSTFREEZES SHUT100^ 101ELAPSED TIME SINCE HOLE DRILLED (DAYS)Fig. 5.2: 92H26 freeze-in temperature record (92T02). Solid lines are smooth extrap-olations used to find time of freeze-in. Thermistors at 55.9m and 65.9m displayed notemperature drop and are omitted for clarity.Chapter 5. Model Results^ 93ThermistorStringDepth(m)Borehole Length(m)Time AfterDrilling (days)Time After Startof Pressure Rise(days)92T0192T028.0613.0618.0623.0628.0638.0648.0658.0668.065.9210.9215.9225.9235.9245.9255.9265.9261.456.451.446.441.431.421.411.41.463.658.653.643.633.623.613.63.60.580.720.951.> 9.0> 9.00.510.671.> 8.9>> 8.4>> 8.4> 8.4Table 5.1: Freeze-in times for 92T01 and 92T02. The borehole length is glacier depth at921124 (69.48m) minus the thermistor depth below the glacier surface.- 19901991 ^- 19928060402080604020Chapter 5. Model Results^ 940^100^200^300DECIMAL DAYFig. 5.3: Multi-year pressure record for 90P02. Note that after one year there is consid-erable activity, indicating the sensor remains in water, whereas after two years only slowchanges occur, indicating that the sensor is frozen in and observes only slowly varyingice stress.92T01(92H24)X 92T02 (92H26)* INITIAL WATER LEVEL IN 92H24Chapter 5. Model Results^ 956050040LLI_J0 30UJ0CO 20100^1^3^4^5^6^7•^ 8^9^10TIME (DAYS)Fig. 5.4: 921124 length versus time: observed and fitted.^with the exponential series L(t) =^a exp(bit) (form required by the model) byminimizing the least-squares misfit using the down-hill simplex method of Press et al.[1987, pp. 289-293]. The bis were constrained to be negative, to prevent the length fromincreasing with time. The result is L(t) of the bottoming hole 921124, and the invertedfunction t(L) is the freezing time as a function of height above the bed in that borehole(Fig. 5.4). Since the thermal conditions in other boreholes are assumed to be the sameas those for hole 921124, it is also the freezing time as a function of height above thebed in other nearby boreholes. Hence for these nearby boreholes the length as a functionof time can be obtained from L(t) for 921126. To adapt the estimate to the blind hole921126 two changes must be made: L must be made smaller by the distance between thebottom of the hole and the bed (by adding a negative ao), and L must be reduced soChapter 5. Model Results^ 96that Lt=0 is consistent with the initial water level in the borehole. Since freezing timeversus depth is assumed to be the same, this is done by shifting the time zero (i.e., startlater on the same curve). So the new as are set to the old ai exp(bi At), where At is thetime at which L = Lt.0 in the new hole. This was done for 921126 and the other holesconsidered in this chapter.5.3 Acceptable Parameter RangesBefore the observed pressure records are fit by adjusting input parameters, it is necessaryto place reasonable limits on these parameters so that, for example, an ice temperatureof +5 °C may not be used to obtain a good fit. Some parameter ranges were obtainedfrom the literature (for example ice and bed properties), whereas others were obtainedfrom Trapridge Glacier conditions (for example temperature and geometry).5.3.1 GeometryThe freezing forcing L(t) was obtained from the 1992 thermistor data as outlined inSection 5.2. The borehole radius is that of Trapridge Glacier boreholes. These holes arenominally 5.0 cm diameter at the time of drilling, but are usually larger near the surface.Hence 0.02 m < rb < 0.03m, with a nominal value of 0.025m appears reasonable. Butit is quite unlikely that the borehole still has this radius when it begins to pressurize.Because freezing to the wall occurs more quickly in the colder ice near the glacier surfacethe borehole will be smaller at the top. Suppose the borehole always has a parabolicvertical cross-section (Fig. 5.5). The equation of this parabola isz = L (r/rb)2.^ (5.1)Chapter 5. Model Results^ 97rb->.Fig. 5.5: Parabolic borehole.Chapter 5. Model Results^ 98The volume of water in this borehole is found by evaluating the volume integral 0 (L z) r dr de.Substituting (5.1) for z1022.^rde f.') L (1 —r de [L frb r dr1021r de Lr21: —102w de^—21r 1deLr?, I: de—Lr22 bComparing this to the volume V = irrl of a cylindrical borehole of length L and radiusrb, it is evident that the parabolic borehole has half the volume of the cylindrical one.Therefore to approximate the parabolic borehole by a cylindrical one of the same volume,as required for the model, it is necessary to reduce the radius by 1/N5. Hence the boreholeradius, for modelling purposes, will be taken to be 0.018m, with minimum of 0.014m.The maximum will remain 0.03m, to allow for uncertainty in the above assumption.The radius of the subglacial cavity excavated by the drill at the bottom of the bore-hole is poorly constrained. Scratches in the paint of a dummy ploughmeter insertedand recovered in the summer of 1992 indicate that the area of disturbed and weakenedsediment (or no sediment) is as much as 20 cm deep (U.H. Fischer, unpublished data,1992). Although a large volume of material is disturbed by the drill, it is unlikely that allthis material is completely removed from the bed in an unconnected borehole in whichthe water (and sediment) has nowhere to go. It is likely stirred up and loosened, thenV (r/rb)2) r dr— (L/re2,) r r3drirbl±r44rLrd(5.2)(5 .3)Chapter 5. Model Results^ 99allowed to settle back into place. Hence the cavity radius might be considerably smallerthan the ploughmeter result suggests. The borehole radius will be taken to be the lowerbound. Hence 0.025 < rc < 0.20m is assumed for bottoming holes. In blind holes thereis no basal cavity, so a cavity radius of zero should be used. This would cause numericaldifficulties in the model though, so a small radius of 1 x 10-5m will instead be used forblind holes.5.3.2 Ice PropertiesThe minimum elastic parameters found in the literature were those computed from theseismic velocities of Paterson [1981, p. 78]; L = 3.3005 x 109 Pa, and A = 6.3608 x 109 Pa,corresponding to a Young's Modulus E = 8.77 GPa (where E = (3A 2p.) / (A -I- A)).Shyam Sunder and Wu [1989a, 1989b, 19901 use E = 9.5 GPa, and Gold [1977] gives arange of 9-11 GPa. I will therefore use the range 8.77 < E < 11 GPa. Assuming thatthe change in E is shared equally by and A, this gives the ranges 3.3005 x 109 < A <4.1 x 109Pa, and 6.3608 x 109 < A < 8.0 x 109Pa.The activation energies QL and QH will be assigned the values 'of Paterson [1981,p. 34], as was done in Chapter 4; QL = 67000J mori and QH = 126000J morl. Iconsider it unnecessary to allow these to vary, since they affect only V, which is alsoaffected by the variables Vo and T.The viscous parameter Vo is widely discussed in the literature. Shyam Sunder andWu [1990] give a range of values 5390 Pa sh/3 < Vo < 7860 Pa s1/3. These correspondto 73 x 105Pa s1/3 < V(-1°C) < 106 x 105Pa s1/3 (where V(T) is given by (2.24)).Values of V(-1 °C) were computed corresponding to the viscous parameters of Pa-terson [1981, p. 39] (102 x 106 Pa s1/3), Hooke [1981] (101 x 105Pas1/3), Budd andJacka [1989] (99 x 105 Pa s1/3), and the high-temperature experiments of Morgan [1991](131 x 105 Pa s1/3). Of these Morgan's is the highest. Hence I assume 73 x 105 <Chapter 5. Model Results^ 100A ft BoFresh water columnar (S-2)d = 3 mm [Sinha, 1978]Fresh water isotropic polycrystallined = 1 7 mm [Jacka, 1984]0.330.01420.4540.020.08650.286Table 5.2: Transient parameters from Shyam Sunder and Wu [1990]; d is the averagegrain size of the ice tested.V(-1 °C) < 131 x 106 Pa s1/3. Inverting (2.24) and rounding the results gives the rangefor V0 as 5400 <V0 <9700 Pa s1/3.The transient parameters A, H, and B0 are obtained from Shyam Sunder and Wu[1990]. They give values for A, H, and B0 for two kinds of ice derived from test datareported by others in laboratory tests (Table 5.2). The transient behaviour of ice hasbeen shown to be grain-size dependent [e.g., Sinha, 1979]. Both the tests quoted weredone on relatively fine grained ice. Glacier ice can have much larger grains, anywherefrom one to several millimetres in polar ice to hundreds of millimetres in temperate ice[Paterson, 1981, pp. 233-239]. Hence the test values must be extrapolated to larger grainsizes. If the grain-size dependence of Sinha [1979] is accepted, the transient parametersare related to the grain size by [Shyam Sunder and Wu, 1989a]A^dl .41 (5.4)dw+i)IN (5.5)Bo = gadliN (5.6)Chapter 5. Model Results^ 101A' (mm) ill (mm-113) .13L (mm-113)From data of Sinha [1978] 120 0.010 0.24From data of Jacka [1984] 9 0.105 0.06Table 5.3: Grain-size-independent transient parameters.A il Bod (mm) Sinha Jacka Sinha Jacka Sinha Jacka1 0.1 0.008 0.1 0.01 0.06 0.2100 10 0.8 50 5 0.3 1.0Table 5.4: Transient parameters computed from the data of Sinha [1978] and Jacka[1984].Inverting these and applying them to the parameters of Table 5.2 gives the grain-size-independent parameters (Table 5.3). Now assuming a minimum grain size of 1 mm and amaximum grain size of 100 mm, and applying (5.4)-(5.6), gives the grain-size-dependenttransient parameters (Table 5.4). Taking the extreme values from Table 5.4 (and round-ing) gives the ranges of the transient parameters as 0.01 < A < 10, 0.01 <H < 50, and0.05 < Bo < 1.In the absence of more detailed information, ice temperature will be assumed con-stant in time and space. The representative temperature will be some sort of averagewith depth in the borehole and distance from the borehole. Ice cannot be warmer thanthe pressure melting point, approximately 0°C. The coldest reasonable effective averageChapter 5. Model Results^ 102ice temperature must be somewhat warmer than the average background ice tempera-ture, because ice surrounding the borehole is warmed by drilling, and deformation willconcentrate in the warmer ice. Since the coldest temperature observed is —6.5°C(92T02, at 5.9m depth, Fig. 5.2), I will adopt a lowest reasonable average temperatureof —2°C. Hence —2 < T < 0°C.5.3.3 Bed PropertiesThe range of hydraulic conductivity and matrix compressibility for glacial till depositsare [Freeze and Cherry, 1979, pp. 29, 56]) 1012 < K < 10-8m s-1 and 10-8 < a <10-6 Pa.-1. Hence 10-8 <a + n13 < 10-6 Pa-1 (because ,8 10-1° Pa-1). Consolidationtests of Trapridge Glacier forefield till yield values of K = 2.2 x 10-8 m s-1 and a += 6.4 x 10_6 Pa' (T. Murray, unpublished data, 1992). Hence I will use the ranges1.0 x 10-12 < K < 1.0 x 10-6ms-1 and 1.0 x 10-8 < a -I- <1.0 x 10-5Pa-1. Theseand other parameter ranges are summarized in Table Ice Parameters from Blind Hole FitsBecause blind holes do not reach the bed, the pressure signal observed in them must bea function of ice properties, geometry, and freezing rate only. Hence it should be possibleto discover the ice properties by forward modelling and fitting the observed pressuresignals, assuming that the geometry and freezing forcings are chosen correctly.The first record to be modelled is that of 92P02, the pressure sensor in 921126 whichalso contained the thermistor string 92T02. The freeze-in record of this sensor shows avery rapid pressure rise to high pressure, followed by a sudden drop, and further repeatedsteep rises and sudden drops until, after several days, the pressure eventually stabilizesand begins to decline (Fig. 5.6). The rapid rise results from the compression applied byChapter 5. Model Results^ 103aN= 180E•"'' 160wc:D0 140owrxa_ 120r:wF_ 100<w 80___I0=w 60cK0co 202^204^206DECIMAL DAY, 1992208^210Fig. 5.6: Freeze-in pressure record for 92P02 (in borehole 921126).Chapter 5. Model Results^ 104Parameter Minimum Maximum Unitsrb 0.014 0.03 mrc 0.025 0.20 mA 3.3005 x 109 4.1 x 109 Pa'A 6.3608 x 10 8.0 x 109 Pa-1QL 67000 J moriQH 126000 J moriVo 5400 9700 Pa s1/3A 0.01 10ii. 0.01 50 —Bo 0.05 1 —T —5 0 °CTable 5.5: Summary of input parameter ranges.the expansion of freezing water. The sudden drops presumably reflect cracking of theice, allowing water to exit the borehole, thus lowering the pressure. It is likely that thiscracking occurs primarily at the top of the borehole where the sharp curvature of thesurface and the presence of the instrument cables causes the greatest stress concentration.To determine ice properties only the initial part of the record, from first pressure riseto first fracture, will be modelled. This is necessary because the SS&W rheology has nofracture mechanics embedded in it. The computer model cannot be started at the endof the drilling, as would be ideal, since it is not flexible enough to deal with a constantpressure followed by freezing-induced compression.Closer examination of the early part of the 92P02 record reveals that the initialpressurization is not one smooth curve (Fig. 5.7). I am unsure what causes this inflectionin the curve. All freeze-in pressure records examined for this thesis displayed this feature,Chapter 5. Model Results^ 105cL 1600Lu1400(I)LLI1200a_1-LI 1000I-L-1^80006000CO.0^.5^1.0^1.5^2.0^2.5^3.0^3.5TIME (HOURS)Fig. 5.7: Initial part of 921126 freeze-in from 92P02).1809,416014001200 100uJ0 8060Chapter 5. Model Results^ 106.0^.5^1.0^1.5^2.0^2.5^3.0^3.5TIME (HOURS)Fig. 5.8: Freeze-in pressure records considered in this thesis: (a) borehole 921E45 (pressuresensor 92P08), (b) 921126 (92P02), (c) 891150 (89PO4), (d) 901154 (90P10), (e) 921E24(91P12), and (f) 911146 (91P05). Dashed lines (a) and (b) were blind holes. Solid lines(c), (d), (e), and f) were believed to be drilled to the bed, although (c) has the characterof a blind hole. Broken lines indicate interpolated missing data.with the sole exception of 89PO4 (Fig. 5.8). One possible explanation involves waterentry into instrument cables. Pressurized water has been observed to flow up insidecables (between the outer jacket and the insulated conductors) and into data loggerenclosure boxes. Thus initially water could flow into the cables and force its way upwards.Eventually the water in the cable above the freezing front would freeze, preventing furtherwater entry. (The volue of water displaced by the assumed freezing before the inflectionpoint is 0.120 m3, somewhat less than the volume of air in the two cables of approximately0.400 m3.) This would cause the inflection point on the pressure records. After thatChapter 5. Model Results^ 107the system would behave as a water-filled borehole, as modelled in this thesis. (It isnot possible that the compressibility of cable materials, PVC and copper, plays a partsince these materials are less compressible than water [Van Krevelen and Hoflyzer, 1972,pp. 150-152].) Hence I will attempt to model only the curve after the inflection point. Todo so, I will extrapolate this curve back until it reaches the initial pressure (Fig. 5.7) andtime shift the borehole L(t) (as was done in Section 5.2) so that the simulation beginsat that point. This is yet another compromise, necessitated by the fact that the modeldoes not include all the physics involved.The model ice parameters were varied until the best visual fit with the pressure recordwas obtained. It was found that the pressure rise was scaled up by decreasing T, orincreasing 'Vo, it, or A. The rate of pressure rise was increased by increasing the transientparameters A, n, and Bo. Changing rb had no significant effect. The best fit, whichis not particularly good (Fig. 5.9), was obtained with T = —2 °C (coldest allowable);Vo = 9700 Pas", ji = 4.1 x 109 Pa.-1, and A = 8.0 x 109 Pa-1 (highest allowable); andA = 0.02, fi = 0.02, and B0 = 0.05 (rather close to the lowest allowable). The poorfit suggests that the model is not adequately describing the system. One likely cause isthe poorly constrained, approximate, top-down freezing forcing. Another is the use ofa one-layer model to simulate deformation in ice with depth-varying temperature, andhence rheological properties. A third effect is that new, unstrained ice is continuouslybeing added to the borehole wall by freezing, which could be expected to depress theapparent values of the transient parameters.Another possible cause of the poor fit is that the ice strains get large enough that theinfinitesimal-strain assumption fails. However, at the time of first fracture, the model-predicted borehole-wall tangential strain is only 2.5 x 10 (with best-fit ice parameters),which is still quite small. The strain elsewhere in the ice is lower. Hence the straindefinition is unlikely to be the source of error.Chapter 5. Model Results^ 108a_16001400(I)a_ 1200100080006000.0^.5^1.0^1.5^2.0^2.5^3.0^3.5TIME (HOURS)Fig. 5.9: Best fit to 921=126 freeze-in pressure record. Solid line is observed pressurerecord, dashed line is bfil4 model simulation.Chapter 5. Model Results^ 109uJ1800(f)16001400w 1200F-uJ100008000co 600.0^.5^1.0^1.5^2.0^2.5^3.0^3.5TIME (HOURS)Fig. 5.10: Best fit to 921145 freeze-in pressure record. Model ice parameters areT = —2°C, Vo = 9700Pas1/3, A = 4.1 x 10913a-1, A = 8.0 x 109 Pa-1, A = 0.01,= 0.015, and B0 = 0.05.Fitting the freeze-in pressure record of 921145, which contained pressure sensor 92P02and a geophone but no thermistor string, using the L(t) derived from 921124 and 921126(suitably time shifted) gave similar results (see Figure 5.10), and yielded even softertransient parameters.The 921126 record was also fit using a slower decaying L(t). This was done to inves-tigate the effect of off-centre thermistors on the results. If the thermistors are off-centrethey will freeze in sooner than if they were on the centre-line of the borehole, as I haveassumed. Hence the actual borehole length would be longer than assumed. To modelthis possibility I halved all bis in the L(t) expression and repeated the fit of 921126. TheChapter 5. Model Results^ 1101600LUD 1400LU1200a_1000800_J01.1.1 6000.0^.5^1.0^1.5^2.0^2.5^3.0^3.5TIME (HOURS)Fig. 5.11: Best fit to 921126 freeze-in pressure record, with slower-decaying L(t).ice properties for the best fit are somewhat different (T = —2°C, the coldest allowable;Vo = 9700Pa s1/3, = 4.1 x 10 Pa-1, and = 8.0 x 10Pa-1 (highest allowable); andA = 0.1, k = 0.2, and Bo = 0.07), although the fit is still poor (Fig. 5.11).5.5 Bed Parameters from Unconnected Hole FitsThis investigation concentrates primarily on 921124, one of the two 1992 boreholes fromwhich the length function L(t) was derived. It was also the only unconnected hole in 1992for which a pressure record of the freeze-in was obtained. This hole was unconnected atthe time of drilling, as indicated by the high and constant pressure in the hole betweendrilling and the freezing-induced pressure rise (Fig. 5.12). This hole did later connectChapter 5. Model Results^ 1110c,i 952E 9085LiCCD 80ocn1-1-1 750_70(2LiI— 65<60L.Li_J0 55iLu500CO 202^204^206DECIMAL DAY, 1992208^210Fig. 5.12: Freeze -in pressure record of 921124 (from 91P12).Chapter 5. Model Results^ 112to the subglacial drainage, as indicated by the onset of large-amplitude diurnal signalhalfway through day 201. However I believe that the hole was unconnected during thefreeze-in period, until approximately 0900h on day 201.The method used to obtain bed-parameter estimates was similar to that used, inthe previous section, to obtain ice-parameter estimates from blind holes. The freezingforcing was time shifted to account for the initial water level difference (none in the caseof 921124, because it is the hole for which L(t) was initially computed), and the delaybetween the first pressure rise and the start of the extrapolated post-inflection freeze-incurve.Since the basal-cavity radius re is poorly known, fits were done by adjusting K anda + ni3 for several values of re spanning its range. It was found that the height of thepressure peak was decreased by increasing K, a + 743, or re. The shape of the pressurepeak was affected differently by changing the different parameters. Increasing K causeda sharper peak with a more rapid subsequent pressure decline, whereas increasing a-En/3or re had the opposite effect of causing a rounder peak with a slower decline. Hence fora given re the records were fit by trading off K and a + rt,8 to obtain the correct heightand shape of curve.First 921124 was fit using the "standard" L(t) and ice properties determined fromthe fit of 921126 using the same L(t). Fits were obtained for re of 0.025m, 0.10m, and0.13m (Fig. 5.13 and Table 5.6). A good fit could not be obtained for re > 0.13mwithout decreasing a ni3 below the minimum allowable 1.0 x 10-8Pa-1. This suggeststhat the basal cavity is not larger than 0.13m in radius, consistent with the reasoning ofSection 5.3.1 that the material disturbed by the drill settles back into place.6509000850w 800750cr)I-1-1 700a_9000a__y 850800750cnI-1-1 7006509000850w 800CCcn 750(/)7000_650aChapter 5. Model Results^ 1130^1^2^3^4^5^6TIME (HOURS)Fig. 5.13: Best fits to 921124 freeze-in pressure record, with various assumed values of rc:(a) 0.025m, (b) 0.10m, and (c) 0.13m.Chapter 5. Model Results^ 114rc (m) K (ms-1) a + nfl (Pa-1)0.025 7.0 x 10-8 1.6 x 10-80.10 1.7 x 10-8 3.0 x 10-80.13 1.35 x 10-8 1.0 x 10-8Table 5.6: Best-fit bed parameters from fit to 921124 freeze-in pressure record, usingstandard L(t).5.5.1 Effect of Ice PropertiesAn important question is how much influence the ice properties have on these results. It isnecessary to know which parameters most strongly affect the results in order to assess thereliability of results obtained using uncertain input parameters. To do this the behaviourof the model was investigated with the softest allowable ice, the stiffest allowable ice, andrigid ice. The softest ice is that with all ice properties minimum, except temperaturewhich was the maximum. The stiffest ice is that with all ice properties maximum, excepttemperature, which was —5 °C. Rigid ice was simulated by running bfil4 with thetransient part of the rheology disabled, I/0 = 9700 x 1010 Pas", ji = 4.1 x 1030 Pa, andA = 8.0 x 1088 Pa (higher values resulted in step size errors). The best-fit bed propertiesfor r, = 0.025 m from Table 5.6 were used. The ice does indeed influence the simulationresults (Fig. 5.14), although not as strongly as K does. The pressure rise changes bymany orders of magnitude when K is varied over its range.The 921124 freeze-in record was then fit, for rc = 0.025m, with softest ice and rigidice (see Figure 5.15 and Table 5.7). The best fit obtained with the softest ice is not verygood. A good fit could not be achieved without reducing a + 70 below 1.0 x 10-813a-1,its minimum reasonable value. A proper fit could be obtained with a slightly higher KChapter 5. Model Results^ 115'E..; 950LLI 900v) 850Uia_ 800<•- 7500 700UJ0 650CD0^1^2^3^4^5^6TIME (HOURS)Fig. 5.14: Freeze-in pressure record of 921124, modelled using softest, stiffest, and rigidice.K (ms-1) a +743 (Pa-1)Softest Ice 7.2 x 10-9 1.0 x 10-6921126 Best Fit Ice 7.0 x 10-9 1.6 x 10-6Rigid Ice 6.5 x 10-9 4.0 x 10-8Table 5.7: Best-fit bed parameters from fit to 921124 freeze-in pressure record, withvarious ice properties. Standard L(t), rc = 0.025m.60 5a_900LU850Lii800F_ 750w 7000w 6500 2^3^4TIME (HOURS)Chapter 5. Model Results^ 116900Liiv) 850(I)1.1Jcca_ 800F-- 750w 7000iii 6500CD a0^1^.2^3^5^6Fig. 5.15: Best fit to 921124 freeze -in pressure record with (a) softest ice and (b) rigidice, for rc = 0.025m.Chapter 5. Model Results^ 117and a slightly lower a + nO. Nevertheless the softest ice values in Table 5.7 are probablywithin 20% of the best-fit values. It is evident that the hydraulic conductivity estimate isquite insensitive to the assumed ice properties, but the compressibility estimate is moresensitive. Also, it appears that assumed ice properties have less effect on the hydraulicconductivity estimate than does the assumed basal-cavity radius (Tables 5.6, 5.7). Thisis explained by the fact that for bed properties in this range the volume of water expelledinto the bed is much greater than the increase in volume of the borehole due to icedeformation.5.5.2 Effect of Freezing ForcingPerhaps the largest uncertainty surrounds the freezing forcing. To assess the effect of thefreezing forcing on the inferred bed properties, 921126 was separately fit with a weakerfreezing forcing than the standard L(t) and a stronger one.The weaker forcing was the (appropriately time shifted) one used in Section 5.4, exceptthat all bis in the L(t) expression were doubled (hence the borehole length decreases moreslowly). The ice properties derived from the fit of 921124 with this forcing were used,as was an rc of 0.025m. The best fit is fairly poor (Fig. 5.16). This suggests that thisfreezing forcing is not accurate.The stronger forcing was achieved by increasing the borehole radius rb to 0.025 m,from 0.018 m. This effectively increases the forcing since a greater volume of water mustbe expelled into the bed. 921126 was fit with 7-, of 0.025 m, 0.10 m, and 0.20 m (Fig. 5.17).The best-fit bed properties for both weaker and stronger freezing forcings, summarizedin Table 5.8, show little change in K, but large change in a + ni3.This investigation only scratches the surface of the freezing forcing issue. Although itappears that the hydraulic conductivity estimate is not overly sensitive to perturbationsto the freezing forcing, the fact that freezing is so crudely modelled suggests that thisChapter 5. Model Results^ 118-c:"ia. 900_Y..-..."LL.1CtD 850cn(I)Lif: 800a_ccwF— 750<w_1 7000=Li(: 6500m0^1^2^3^4^5TIME (HOURS)Fig. 5.16: Best fit to 921124 freeze-in pressure record, with weaker freezing forcing.rc = 0.025m.K (m s-1) a + nf3 (Pa-1)Weak Forcing 4.5 x 10 3.5 x 10-7Standard Forcing 7.0 x 10' 1.6 x 10-6Stronger Forcing 1.3 x 10-8 4.0 x 10-6Table 5.8: Best-fit bed parameters from fit to 921324 freeze-in pressure record, with weak,strong, and standard freezing forcings. re = 0.025m.; ; 1 1 ;-- ---------- ___--_1bChapter 5. Model Results^ 119.-. 9000CL__Nc 850..._...,w 800CCop 750U)700CCa_,—, 9000a_-- 850.......-I-Li 800OCD(f) 750U)1-LI 700CCa_, 9000a___y 850Li 800r:D 750(i)U)t_Li 700a_6506506500^1^2^3^4^5^6TIME (HOURS)Fig. 5.17: Best fits to 92H24 freeze-in pressure record with rb = 0.025m and variousvalues of rc: (a) 0.025m, (b) 0.10m, and (c) 0.20m.Chapter 5. Model Results^ 120remains an area of large uncertainty.5.5.3 Other BoreholesTo check the results inferred from the freezing of the single 1992 unconnected borehole,the freeze-in pressure records of three other boreholes from previous years were also exam-ined: 891150 (pressure sensor 89PO4), 901154 (90P10), and 911146 (91P05). These holesdid not connect when drilled. 901124 and 911146 both have freeze-in pressure curves withcharacters similar to the unconnected 921124, while 891150 looks suspiciously like a blindhole (Fig. 5.8). Nevertheless all three were modelled as bottoming unconnected holes,with the (suitably time-shifted) freezing forcing determined from the 1992 thermistor ob-servations and ice properties determined from the fitting of 921126. It is certain that thisfreezing forcing is at least somewhat in error, owing to the differing drill characteristicsand method of these years compared to 1992.The fits are not especially good (Figs. 5.18, 5.19, and 5.20). Hydraulic conductivitiesdetermined from fits to 901154 and 911146 agree quite well with those determined from921124, while that from 891150 is quite a bit lower (Table 5.9). The compressibilityestimates span the entire allowable range. The low hydraulic conductivity from 891150is no surprise, since the character of its pressure record suggests that it is a blind holeanyway. Therefore these results do not contradict the results obtained from 921E24.5.5.4 Summary of Bed PropertiesThe foregoing forward modelling has produced reasonably consistent estimates of hy-draulic conductivity and wildly varying estimates of bed compressibility. The hydraulicconductivity falls in the range 1.35 x 10' m s-1 (for 7-, = 0.13m) to 1.3 x 10-8m s-1 (forTb = 0.025m, r, = 0.025 m), inferred from the 1992 unconnected borehole 921124. Thisresult is supported by the results from two boreholes (901154 and 91H46) drilled in otherChapter 5. Model Results^ 121'-c:-a..1600wc:D0 14000wa.1200wi—<1000w0= 800w(2C0cn.0^.5^1.0^1.5^2.0^2.5^3.0^3.5TIME (HOURS)Fig. 5.18: Best fit to 891150 freeze-in pressure record with rc = 0.025m.K (m s-1) a + n f3 (Pa-1)891150 4.0 x 10-19 1.0 x 10'901154 6.0 x 10 2.0 x 10-7911146 7.0 x 10' 1.0 x 10'921124 7.0 x 10-9 1.6 x 10'Table 5.9: Best-fit bed parameters from fits to various borehole freeze-in pressure records.Standard L(t), rc = 0.025m.o__y 105011-1 1000950LUa_ 9008508007500Li-I 7000m 6503.0 3.5.0^.5^1.0^1.5^2.0^2.5TIME (HOURS)Chapter 5. Model Results^ 122Fig. 5.19: Best fit to 901154 freeze-in pressure record with re = 0.025 m.Chapter 5. Model Results^ 1230-▪ 900LL.1ct 850(/)tn• 800jrX 7507000650CC0Ca0^1^2^3^4^5^6TIME (HOURS)Fig. 5.20: Best fit to 911146 freeze-in pressure record with rc = 0.025m.Chapter 5. Model Results^ 1242 7570LUCC65(/)uJCL 60< 55500UifX 450220^240^260^280^300^320^340^360DECIMAL DAY, 1990Fig. 5.21: Pressure record of 90P07. Note the numerous sudden pressure changes andsubsequent decays at various time scales.years (using a somewhat different drilling method). The very low estimate from 891150is ignored as I conclude that it was a blind hole. This range is in good agreement withthe values obtained by other investigators (upper bound of 3.6 x 10'm s-1 [Smart andClarke, submitted], 1.27 x 10-8m s-1 [Clarke, 1987], and 2.2 x 108ms-1 (T. Murray,unpublished data, 1992). The estimates obtained for bed compressibility are so wildlyvarying that I attach no significance to them.5.6 Response to Sudden Pressure ForcingSudden pressure rises (or drops) are frequently seen on subglacial pressure records (e.g.,Fig. 5.21). Are these events intrinsic to the basal hydrological system, or could they beChapter 5. Model Results^ 125influenced by the response of the borehole itself? To answer this question I examine theresponse of an unconnected borehole to a sudden pressurization. If the system respondsslowly, and the excess pressure decays with a time constant similar to those observedin pressure records, then it is possible that the borehole is responsible for the observeddecays. If however the excess pressure decays away very quickly, then the observed slowdecays must be due to some other processes unrelated to the borehole.Simulations were run with the ice properties found from consideration of 921124 freeze-in (see Section 5.4). To place an upper bound on the length of the decay, the longestreasonable borehole (921324 after only one day, L = 45.3 m), and the smallest pressureforcing of interest (pfarce = 10 Pa, or about 1 m1120) were chosen.When bed properties found from 921124 were used (rc = 0.025 m, K = 7.0 x10' m s' ,and cx ni3 = 1.6 x 10' Pa-1), the decay curve has approximately the right shape, butthe decay is much too rapid (Fig. 5.22). The pressure drops to 1/e of its initial value injust 56 s. This response is so rapid that such an event would likely go unnoticed on apressure record sampled every 2 min, as taken at Trapridge Glacier during the summerfield season. It almost certainly would be missed by the 20 min overwinter sampling.Hence it appears that the borehole response is essentially instantaneous, and cannot beresponsible for the decay events observed in Trapridge Glacier pressure records.The simulation was repeated with a hydraulic conductivity of zero, to simulate theresponse of a blind hole (Fig. 5.23). The decay is considerably slower, with a 1/e timeof 2.6 days. The decay would be even slower for stiffer ice.Chapter 5. Model Results^ 1260^20^40^60^80^100^120TIME (SECONDS)Fig. 5.22: Sudden pressurization at t = 1 day of hole 921124 (simulation).Chapter 5. Model Results^ 127.0^.^1.0^2.0^3.0^4.0TIME (DAYS)Fig. 5.23: Sudden pressurization at t = 1 day of hole 92E124 assuming impermeable bed(simulation).Chapter 6SummaryIn this thesis I have presented a model for the mechanical response of water-filled bore-holes and applied it to interpret borehole water-pressure records from Trapridge Glacier.The model consists of single-layer plane-strain ice around a water-filled borehole, con-nected to a uniform half-space permeable bed through a hemispherical cavity, which canbe forced by an approximate top-down freezing forcing, or by sudden pressurization. Themodel was formulated using the non-linear transient viscoelastic rheology of Shyam Sun-der and Wu. This is the first' time that this theology has been used to model boreholeresponse or, to my knowledge, any other aspect of glacier physics.The resulting system of non-linear partial differential equations and algebraic equa-tions was solved using the method of lines, with integration by the ddriv3 integratorof Kahaner and Sutherland. The model was tested by comparing model results to theanalytic solutions for Glen-Law ice, linear elastic ice, and the response of a sphericallysymmetric aquifer to a step head change on its inside boundary. Satisfactory solutionswere found, although the stiffness of the problem did limit accuracy, especially for (in-compressible) Glen-Law ice. Spatial oscillations of the ice solution were observed forsome choices of spatial coordinate transformation, possibly as a result of using a centred-difference approximation of the spatial derivatives. The model gave qualitatively correctresults for various tests of the coupled system.The model was then used to infer the hydraulic properties of the bed under TrapridgeGlacier, by forward modelling of the freeze-in pressure records of six boreholes. The128Chapter 6. Summary^ 129model freezing forcing was inferred from the freeze-in temperature records from ther-mistor strings in two boreholes. Ice properties were then inferred from the modelling ofthe freeze-in pressure records from two blind holes. These results were not consideredreliable, as the model results did not fit the pressure records well. This is probably dueto the simple freezing forcing used in the model.The freeze-in pressure records of four unconnected boreholes were then modelled,using the ice properties from the blind-hole modelling, to obtain estimates of the prop-erties of the low-permeability layer. Good fits were achieved to the 1992 pressure record(the year for which the freezing forcing was estimated) whereas the fits to records fromother years were poorer. The hydraulic conductivity was found to be in the range1.3 x 10-9-1.3 x 10-8 ms', which agrees well with previous estimates (upper boundof 3.6 x 10-7 m s-1 [Smart and Clarke, submitted], 1.27 x 10-8 m s-i [Clarke, 1987], and2.2 x 10-8 m s' (T. Murray, unpublished data, 1992). The hydraulic conductivity wasfound to be most strongly influenced by the assumed freezing forcing and size of the cav-ity in the bed at the bottom of the borehole. It was found to be rather weakly dependenton the assumed ice properties. The bed compressibility was found to vary so much thatno useful estimate of it was obtained.In addition, the response of blind and unconnected boreholes to sudden pressurizationwas investigated. It was found that the pressure decays very quickly in unconnected holes,much more quickly than any decay observed in pressure records. In fact the decay is sorapid that it is unlikely to be resolved by the sampling rates used for Trapridge Glacierpressure sensors. It can therefore be concluded that the response of the borehole systemis essentially instantaneous compared to the signals being measured; hence the observedsignals reflect the true basal conditions, rather than the influence of the borehole (exceptduring the initial freeze-in period when the freezing forcing dominates).Some recommendations can be made, based on this work. The approach used in thisChapter 6. Summary^ 130thesis to obtain estimates of the hydraulic conductivity of the basal sediment is sound.The simple one-layer ice model is adequate. The added complexity of a multi-layer icemodel is not justified, since the ice deformation was found to be much less important thanwater flow into the bed from unconnected holes. Improvements could, however, be made.The approximate top-down freezing forcing is a major weak point of the model. Thisshould be replaced by a proper thermal freeze-in model, such as that of Humphrey andEchelmeyer, to accurately model the important freezing forcing. Then perhaps usefulestimates of ice properties could be made, as well as a better hydraulic conductivityestimate. The hydraulic conductivity estimate could also be improved if the size of thebasal cavity at the bottom of the borehole could be better constrained. The use of aformal inversion scheme instead of the trial and error forward modelling approach usedhere would improve the reliability of the results.The finite-difference approximation used for spatial derivatives in the ice model shouldbe re-examined. It is possible that replacement of centred-differences by an asymmetricdifferencing scheme might eliminate the observed spatial oscillations of the stress solution,and hence improve accuracy.The ice deformation part of the model could be used for other purposes, such asinterpreting pressuremeter results, which have previously been interpreted only in termsof the Andrade relation [Kjartanson et al., 1988; Murat et al., 1986; Shields et al., 1989].This would allow the Shyam Sunder and Wu theological parameters to be obtained fromin-situ pressuremeter tests.Finally, the pressure signals observed at Trapridge Glacier indeed reflect true basalconditions, and the influence of the borehole on the signal need not be considered, exceptduring the initial freeze-in period.ReferencesAbramowitz, M., and I.A. Stegun (Eds.), Handbook of Mathematical Functions, 3rded., 1046 pp., Dover, New York, 1972.Ashby, M.F., and P. Duval, The creep of polycrystalline ice, Cold Reg. Sci. Technol.,11, 285-300, 1985.Bear, J., and A. Verruijt, Modeling Groundwater Flow and Pollution, 414 pp., D. Rei-del, Dordrecht, Holland, 1987.Blake, E.W., The deforming bed beneath a surge-type glacier: Measurement of me-chanical and electrical properties, Ph.D. Thesis, Univ. of B.C., Vancouver, 1992.Blake, E., G.K.C. Clarke, and M.C. Gerin, Tools for examining subglacial bed defor-mation, J. Glaciol., 38, 388-396, 1992.Budd, W.F., and T.H. Jacka, A review of ice theology for ice sheet modelling, ColdReg. Sci. Technol., 16, 107-144, 1989.Clarke, G.K.C., Subglacial till: a physical framework for its properties and processes,J. Geophys. Res., 92, 9,023-9,036, 1987.Clarke, G.K.C., and E.W. Blake, Geometric and thermal evolution of a surge-typeglacier in its quiescent state: Trapridge Glacier, Yukon Territory, Canada, 1969—1989, J. Glaciol., 37,158-169, 1991.Clarke, G.K.C., S.G. Collins, and D.E. Thompson, Flow, thermal structure, and sub-glacial conditions of a surge-type glacier, Can. J. Earth Sci., 21, 232-240, 1984.Carslaw, H.S. and J.C. Jaeger, Conduction of Heat in Solids, 3rd ed., 510 pp., Claren-don, Oxford, 1959.Duval, P., Anelastic behaviour of polycrystalline ice, J. Glaciol., 21, 621-628, 1978.Engelhardt, H., N. Humphrey, B. Kamb, and M. Fahnestock, Physical conditions atthe base of a fast moving Antarctic ice stream, Science, 248, 57-59, 1990.Finnie, I., and W.R. Heller, Creep of Engineering Materials, 341 pp., McGraw-Hill,New York, 1959.Freeze, R.A., and J.A. Cherry, Groundwater, 604 pp., Prentice-Hall, Englewood Cliffs,NJ, 1979.131References^ 132Fung, Y.C., A First Course in Continuum Mechanics, 340 pp., Prentice-Hall, Engle-wood Cliffs, NJ, 1969.Gold, L.W., Engineering Properties of Fresh-Water Ice, J. Glaciol., 19, 197-212, 1977.Humphrey, N., Estimating ice temperature from short records in thermally disturbedboreholes, J. Glaciol., 37, 414-419, 1991.Humphrey, N., and K. Echelmeyer, Hot-water drilling and bore-hole closure in coldice, J. Glaciol., 36, 287-298, 1987.Jaeger, J.C., and N.G.W. Cook, Fundamentals of Rock Mechanics, 515 pp., Chapmanand Hail, London, 1969.Jarvis, G.T., and G.K.C. Clarke, Thermal effects of crevassing on Steele Glacier, YukonTerritory, Canada, J. Glaciol., 13, 243-254, 1974.Kahaner, D., C. Moler, and S. Nash, Numerical Methods and Software, 495 pp.,Prentice-Hall, Englewood Cliffs, NJ, 1989.Kamb, B., Glacier surge mechanism based on linked cavity configuration of the basalwater conduit system, J. Geophys. Res., 92, 9,083-9100, 1987.Kamb, B., C.F. Raymond, W.D. Harrison, H. Engelhardt, K.A. Echelmeyer, N.Humphrey, M.M. Brugman, and T. Pfeffer, Glacier surge mechanism: 1982-1983surge of Variegated Glacier, Alaska, Science, 227, 469-479, 1985.Kjartanson, B.H., D.H. Shields, L. Domaschuk, and C.S. Man, The creep of ice mea-sured with the pressuremeter, Can. Geotech. J., 25, 250-261, 1988.LeGac, H., and P. Duval, Constitutive relations for the non-elastic deformation of poly-crystalline ice, in IA TAM Symposium: Physics and Mechanics of Ice, Copenhagen,pp. 51-59, Springer-Verlag, Berlin, 1980.Ude, D.R., ed., Handbook of Chemistry and Physics, 71st ed., 2324 pp., CRC Press,Cleveland, 1991.Mase, G.E., Theory and Problems of Continuum Mechanics (Schaum's Outline Series),221 pp., McGraw-Hill, New York, 1970.Mesinger, F., and A. Arakawa, Numerical methods used in atmospheric models, Vol. 1,64 pp., GARP Publications Series No. 17, 1976.Morgan, V.I., High-temperature creep tests, Cold Reg. Sci. Technol., 19, 295-300,1991.Morland, L.W., and U. Spring, Viscoelastic fluid relation for the deformation of ice,References^ 133Cold Reg. Sci. Technol., 4, 255-268, 1981.Murat, J.R., P. Huneault, and B. Ladanyi, Effects of stress redistribution on creepparameters determined by a borehole dialatometer test, in 5th. International Sym-posium of Offshore Mechanics and Arctic Engineering, Tokyo, Vol 4, pp. 58-64,ASME, 1986.Nye, J.F., The flow law of ice from measurements in glacier tunnels, laboratory exper-iments and the Jungfraufirn borehole experiment, Proc. Roy. Soc. London, Ser. A,219, 477-489, 1953.Paterson, W.S.B., The Physics of Glaciers, 2nd ed., 380 pp., Pergamon, Oxford, 1981.Prager, W., Introduction to Mechanics of Continua, 230 pp., Dover, New York, 1961.Press, W.H., B.P. Flannery, S.A. Teukolsky, and S.T. Vetterling, Numerical Recipes,230 pp., Cambridge University Press, Cambridge, 1987.Raymond, C.F., How do glaciers surge? A review, J. Geophys. Res., 92, 9,121-9,134,1987.Schiesser, W.E., The Numerical Method of Lines, 326 pp., Academic Press, San Diego,1991.Shields, D.H., L. Domaschuk, F. Azizi, and B.H. Kjartanson, Primary Creep Parame-ters for ice measured in-situ, Cold Reg. Sci. Technol., 16, 281-290, 1989.Shyam Sunder, S., and M.S. Wu, A differential flow model for polycrystalline ice, ColdReg. Sci. Technol., 16, 45-62, 1989.Shyam Sunder, S., and M.S. Wu, A multiaxial differential model of flow in orthotropicpolycrystalline ice, Cold Reg. Sd. Technol., 16, 223-235, 1989.Shyam Sunder, S., and M.S. Wu, On the constitutive modeling of transient creep inpolycrystailine ice, Cold Reg. Sci. Technol., 18, 267-294, 1990.Sinha, N.K., Short-term rheology of polycrystalline ice, J. Glaciol., 21, 457-473, 1978.Sinha, N.K., Grain boundary sliding in polycrystalline materials, Philos. Mag. A, 40,825-842, 1979.Smart, C.C., and C.K.C. Clarke, Formation, propagation and release of a subglacialwater cavity beneath a surge-type glacier, J. Glaciol. (submitted).Spring, U., and L.W. Morland, Viscoelastic solid relations for the deformation of ice,Cold Reg. Sci. Technol., 5, 221—P, 234.1982References^ 134Spiegel, M.R., Mathematical Handbook (Schaurn's Outline Series), 271 pp., McGraw-Hill, New York, 1968.Stone, D.B., Characterization of the basal hydraulic system of a surge-type glacier:Trapridge Glacier, 1989-92, Ph.D. Thesis, Univ. of B.C., Vancouver, 1993.Stone, D.B., and G.K.C. Clarke, Estimation of subglacial hydraulic properties frominduced changes in basal water pressure: a theoretical framework for boreholeresponse tests, J. Glaciol., in press, 1992.Streeter, V.L., and E.B. Wylie, Fluid Mechanics, 7th ed., 562 pp., McGraw-Hill, NewYork, 1979.Van Krevelen, D.W., and P.J. Hoflyzer, Properties of Polymers: Correlation withChemical Structure, 427 pp., Elsevier, Amsterdam, 1972.Appendix ASummary of Model EquationsThis appendix presents the complete system of equations for the coupled ice—bed—borehole problem, in a form that can conveniently be programmed. Water density, pre-viously termed po, is written in this appendix as pw. The flow law exponent N 3, andthe expressions containing N have been simplified accordingly.Dependent Variableshi (for j = 2, Ma — 1)MaqBi,^ e, ^and crz,i (for i = 1,M)(for i = 1, M —There are a total of Ma + 9M — 2 dependent variables.Initial Conditionsh; = polPwg (for j = 1, Ma — 1)maq = 0= Bo,^= Re,i =-- Rz,i = 0,=^Cuz,i = 07 crr,i = crz,i = —PoDifferential Equationsam^aq 27r-pw exp (pb(0 — po)1Katr,(A.1)1 (for i = 1, M)135Appendix A. Summary of Model Equations^ 136For j = 2, Ma — 1For j = 1,MUIi^K^a2 hat — pwg (a + ni3)14 a R2 (A.2)(A.3)(A.4)(A.5)(A.6)(A.7)(A.8)(A.9)(A.10)(A.11)(A.12)(A.13)a Bi = fr E (  1  )3 d 2a -at^' k,Billi^eq'sORr,i . 2 Aza4.,iat^3^ata R,,i^2^aet i. 3 AE  at ,ata Rz,i^2^aftz -^ =at^3 AE  at ''• aei^act ,i.^a 6,..,iat^at + ataq,i^a et, i^a esi; i^ _^• +^'at^at^ata eLaetz i aezi^ ,  ^•' + '^at^at^atAlgebraic Constraint EquationsFor i = 1, M0^ezi •For i = 2, M — 1ace0 = Wi er,i •RFor i = M0 = ar,m +Po.For i = 1, when no pressure forcing is applied,0 = mw — mo maq m.f.For i = 1, when sudden pressurization is applied and t < tramp,0 = crr,1 pramp + Po. (A.14)Appendix A. Summary of Model Equations^ 137For i = 1, after application of sudden pressurization (t > tramp),0 = MYSTAuxiliary RelationsMramp Maq^Maq,ramp +(7r,1m1 — raf,rarnp• (A.15)(A.16)11 1PwgPo= (A.17)PwgL = E aiebit (A.18)Er^= (A.19)69,i^=^ee ^+ e (A.20)cz,i^=^Cez,i + eu (A.21)1 A(Crr,i2/./.°." az,i) (A.22)2/L^+1 A67),i 2 p.u9 's + z,i) (A.23)2p. (3A + 210 (crr,iEe1 0.^.2p.^z'sAcro,i az,i) (A.24)  2/./ (3A + 2/t) (ar,iGrr,i1 ,Cfr^kCfrCre,i CrZ (A.25)I. 3 •Cre ti^=crz,i1 ,^(re^— — (cr,^, 3^.1 ,— — k^09,i3^'cr.,)Crz,i)(A.26)(A.27)Creq,ifar,i2^ +Icrei (A.28)=^(—V11 30-2'1^2^,^eq'i(A.29)Appendix A. Summary of Model Equations^ 138(A.30)(A.31)(A.32)(A.33)(A.34)(A.35)(A.36)(A.37)(A.38)(A.39)(A.40)(A.41)(A.42)(A.43)(A.44)(A.45)AAs,j Ur,i494,iat^= A,,i'cre,iae.^  =atdt 0de)-d—t ode)(Tli I 0=^Crr,i^Rr,i_CrO,icrc!,i=^cre,i —=^— Rz,i= crd^1 ( „rd^j_— —3 Vr,i^'JO,:^'z,s1 j d dur 9 oi Cr° ,^Vr  i^Cr 0,i + Cr z,i—3^,d a 1 (^ d— Cri• i^4i + Cr z,iaz,i = az , 3^,3 (^ )3credAt 2 BiVi^d 2^, d 2^z,i, d 2)Crr,i^Cre,i^Cr/ d (CIEAticr —at• r's dt) 0afte,i^d^deAt^(—)at °A dt 0aez,,^deAt^—at z'i dt) 0a rxr arO,iCr^vv s OR1^for R = in r transformation=X Ri for R = rx- transformation3eq,iAppendix A. Summary of Model Equations^ 139Vo exp (Qi/3Rga„Ti)^ for T < 263.12K=^ (A.46)Vo exp [(QL — QH)/NRgasT263] exP [QH/NRga,T] for T> 263.12 Kraw2.ffrbexP (-0r,1 — Po)] [1 — exp (-0Pwg.0]2+yrr:Pw exP [13 (—aro. — po)]2^ 2 3irrb o mo =^' [1 exp (-13pwgLo)] -dirrcpwOg(A.47)(A.48)mf = 7r71,0 (L — Lo)^ (A.49)mf,ramp = 771,0 (-tramp — Lo) (A.50)Maq,ramp = rnaq(t =- tramp)^ (A.51)rb = rbro (1 + 694)^ (A.52)^(de) = 1s (A.53)dt 0f;+1 — fi-i for i = 2, M —1^(A.54)2AR—fi+2 4,fi+1 — 3fi ^for i = 1^ (A.55)2AR3fi -^+ fi-2 for i = M (A.56)2AR—fi+2^— 3f.; ^for j = 1^ (A.57)2ARafa+i — 2f; for j = 2, Ma — 1^(A.58)AR!where f is the dependent variable (e.g., cr,., fe, 10.dfdRdfdRdfdRdfdRad2 fdR!


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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


Related Items