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 SEDIMENT DETERMINED FROM THE MECHANICAL RESPONSE OF WATER-FILLED BOREHOLES By Brian Stewart Waddington B.A.Sc. (Mechanical Engineering) University of British Columbia  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF GEOPHYSICS AND ASTRONOMY We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA June 1993 © Brian Stewart Waddington, 1993  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  (Signature)  Department of  Geopky5iG5 and Asir°notttY  The University of British Columbia Vancouver, Canada  Date JUne  DE-6 (2/88)  ae., 19F3  Abstract  The hydraulic properties of basal sediments must be known if the surge mechanism of soft-bedded glaciers is to be understood. Freezing of water-filled boreholes drives water into the subglacial bed and the associated pressure effects give information about subglacial hydraulic properties. A numerical model describing the mechanical response of an unconnected borehole and the bed beneath it to this freezing forcing was developed, using a non-linear transient viscoelastic ice rheology and an approximate model of top-down freezing. The resulting system of equations was solved using the method of lines. Results agree well with analytic solutions, if parameters are correctly chosen. Forward modelling of pressure records from three 1992 boreholes and three from other years indicates that the till underlying Trapridge Glacier has a hydraulic conductivity of 1.3 x 10-9-1.35 x 10' m s-1, in agreement with previous indirect estimates. The model was also used to investigate the response of a borehole to sudden pressure changes. The response is very fast compared to pressure sensor sampling rates: thus the "true" basal signal is essentially unaffected by the presence of the borehole, except during the initial freeze-in. The decay events seen in Trapridge Glacier pressure records are caused by the response of the basal drainage system, not the kprehole.  Table of Contents  Abstract List of Tables^  vi  List of Figures^  viii  List of Symbols^  xii  Acknowledgements^  xvii  1  Introduction  1  2  Ice Deformation Model Development  7  2.1  Introduction ^  7  2.2  Governing Equations ^  8  2.2.1^Rheology of Shyam Sunder and Wu ^  9  2.2.2^Balance Equations  2.3  ^  15  2.2.3^Definition of Strain and Strain Rate ^  15  2.2.4^Compatibility of Strains ^  15  Application of Limiting Assumptions ^  16  2.3.1^Isotropic Ice ^  16  2.3.2^Very Slow Flow ^  23  2.3.3^Radial Symmetry ^  23  2.3.4^Small Vertical Gradients (Except Hydrostatic Pressure) ^  26  ili  2.3.5^Infinitesimal Strains  ^  2.4^Initial and Boundary Conditions ^  3  4  31 34  2.4.1^Initial Conditions ^  34  2.4.2^Boundary Conditions ^  35  Bed Equations, Forcings, and Model Synthesis  37  3.1^Basal Groundwater Flow ^  38  3.2^Borehole Equations ^  42  3.3^Forcings ^  46  3.3.1^Freezing ^  46  3.3.2^Sudden Pressurization  ^  47  3.3.3^Pressure Step Change  ^  48  3.4^Summary of Coupled Initial Value Problems ^  48  3.5^Method of Lines Solution ^  49  3.5.1^Spatial Coordinate Transformation ^  50  3.5.2^Finite Differencing ^  55  3.6^Computer Programming of Solution ^  58  Model Testing  60  4.1  Introduction ^  60  4.2  Ice Deformation Model ^  61  4.2.1^Glen-Law Viscous Ice ^  61  4.2.2^Linear Elastic Ice ^  71  4.3  Basal Groundwater Flow Model ^  76  4.4  Coupled Model ^  77  4.4.1^Sealed Borehole with Glen-Law Ice ^  80  4.4.2^Effect of Viscoelastic and Transient Viscoelastic Ice ^  82  iv  4.4.3 Ice—Water—Bed Model ^ 4.5 Conclusions ^ 5 Model Results^  5.1 Introduction ^ 5.2 Freezing Forcing  85 87 88  88  L(t) Obtained from Thermistor Records ^ 89  5.3 Acceptable Parameter Ranges ^ 5.3.1 Geometry ^ 5.3.2 Ice Properties 5.3.3 Bed Properties ^  96 96 ^99 102  5.4 Ice Parameters from Blind Hole Fits ^  102  5.5 Bed Parameters from Unconnected Hole Fits ^  110  5.5.1 Effect of Ice Properties ^  114  5.5.2 Effect of Freezing Forcing ^  117  5.5.3 Other Boreholes ^  120  5.5.4 Summary of Bed Properties ^  120  5.6 Response to Sudden Pressure Forcing ^  124  6 Summary^  128  References^  131  Appendices^  135  A Summary of Model Equations ^  135  List of Tables  2.1  Notation and physical definitions of Shyam Sunder and Wu rheology. . . .  10  4.1  Standard rheological parameters used for model testing ^  62  4.2  Optimum ddriv3 error control parameters ^  63  4.3  Dependence of Glen-Law strain on coordinate transformation .  70  4.4  Optimum parameters for Glen-Law ice ^  71  4.5  Dependence of linear elastic strain on coordinate transformation ^  76  4.6  Optimum parameters for coupled problem ^  87  5.1  Freeze-in times for 92T01 and 92T02 ^  93  5.2  Transient parameters from Shyam Sunder and Wu [1990] ^  100  5.3  Grain-size-independent transient parameters ^  101  5.4  Computed transient parameters ^  101  5.5  Summary of input parameter ranges ^  104  5.6  Bed parameters from fit of 921124 freeze-in ^  114  5.7  Bed parameters from fit of 921124 freeze-in, with various ice properties  115  5.8  Bed parameters from fit of 921124 freeze-in, with various freezing forcings  118  5.9  Bed parameters from fits of various freeze-ins ^  121  vi  List of Figures  1.1 Types of Trapridge Glacier boreholes and freeze-in pressure records . . . 1.2 Response to freezing forcing  ^4  3.1 Bed geometry and coordinate system ^  39  3.2 Ice finite-difference grid  ^55  3.3 Bed finite-difference grid  ^56  4.1 Effect of maximum radius on strain in Glen-Law ice ^  66  4.2 Glen-Law stress distributions for various coordinate transformations . ^68 4.3 Effect of number of nodes on strain in Glen-Law ice ^  69  4.4 Effect of maximum radius on strain in linear elastic ice ^ 73 4.5 Elastic stress distributions for various coordinate transformations .^74 4.6 Effect of number of nodes on strain in elastic ice ^  75  4.7 Aquifer response for various numbers of nodes ^  78  4.8 Aquifer response at various times ^  79  4.9 Sudden pressurization of a sealed borehole surrounded by Glen-Law ice^83 4.10 Effect of ice rheology on pressure decay in a sealed borehole ^ 84 4.11 Pressure decay for various ice/bed combinations ^  86  5.1 92H24 freeze-in temperature record (92T01) ^  91  5.2 921126 freeze-in temperature record (92T02) ^  92  5.3 Multi-year pressure record for 90P02 ^  94  5.4 921124 length versus time: observed and fitted 5.5 Parabolic borehole ^ vii  ^95 97  5.6 Freeze-in pressure record for 92P02 (in borehole 921126) ^ 103 5.7 Initial part of 921126 freeze-in (from 92P02) ^  105  5.8 Freeze-in pressure records considered in this thesis ^ 106 5.9 Best fit to 921126 freeze-in pressure record ^  108  5.10 Best fit to 921145 freeze-in pressure record ^  109  5.11 Best fit to 921126 freeze-in pressure record, with slower-decaying L(t)^110 5.12 Freeze-in pressure record of 921124 (from 91P12) ^ 111 5.13 Best fits to 921124 freeze-in pressure record, with various r, ^ 113 5.14 Comparison of freeze-in with softest, stiffest, and rigid ice ^ 115 5.15 Best fit to 921124 freeze-in pressure record with softest and rigid ice . .^116 5.16 Best fit to 921124 freeze-in pressure record, with weaker freezing forcing . 118 5.17 Best fits to 921124 freeze-in, with rb = 0.025m and various rc ^ 119 5.18 Best fit to 891150 freeze-in pressure record ^  121  5.19 Best fit to 901154 freeze-in pressure record ^  122  5.20 Best fit to 911146 freeze-in pressure record ^  123  5.21 Pressure record of 90P07 ^  124  5.22 Sudden pressurization at t = 1 day of hole 921124 ^ 126 5.23 Sudden pressurization at t = 1 day of hole 921124 impermeable bed. . . ^ 127  List of Symbols  a()^a5^coefficients of L(t) expression (m) A^Glen Law flow law parameter (Pa's1/3)  A^kinematic hardening parameter grain-size independent kinematic hardening parameter. b1 • • • b5^exponents of L(t) expression (s-1) scalar drag stress associated with isotropic hardening B0^initial scalar drag stress grain-size independent initial scalar drag stress c1^c6^anisotropic parameters  c*^normallizing coefficient, c* = c2c3 c3c1 c1c2 ice grain size (mm) hydraulic diffusivity (m2s1) Dijki^orthotropic  elastic stiffness tensor (Pa)  Young's modulus (Pa) 60^unit vector function (of which finite difference is taken)  f.^body force vector (N) gravitational acceleration, g = 9.806 m s'  Gipet^fourth-order stress transformation tensor hydraulic head (m) ho^initial hydraulic head (m) Hijkl^fourth-order  strain-rate transformation tensor 1X  temperature-independent isotropic hardening parameter  H'^grain-size independent temperature-independent isotropic hardening parameter hydraulic 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 nodes Ma^number of bed nodes porosity of basal sediment flow-law exponent hydrostatic 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 ice Ra^transformed radial coordinate in bed Rgaa^gas  constant, Rgas = 8.314 J mo1-1K-1  elastic 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 the deviatoric stress tensor) (Pa) time (s) tramp^ramp-up time for pressure forcing (s) temperature (K) T263^constant  temperature,  T263 =  263.12 K  ui^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 temperatures  xi  T> —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 derivatives xi^spatial coordinates (m)  X^power of radial coordinate transformation vertical spatial coordinate, z = x3 (m) compressibility of basal sediment matrix (Pa-1) compressibility of water (Pa-1)  7^engineering shear strain (s-1) Kronecker delta total 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 stress cylindrical 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)  xii  lor^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 derivative  a/Ot^spatial time derivative  Acknowledgements  My 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 those involved 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 subglacial pressure data, and for their insightful suggestions. I also thank my family and friends for their support during the inevitable ups and downs of this project. This work has been supported by grants from the Natural Sciences and Engineering Research Council and the University of British Columbia.  xi v  Chapter 1  Introduction  The question of how glacier surges work is a major problem in glaciology. It is now well established that subglacial hydrology plays a pivotal role — surges result from accelerated basal sliding induced by sustained high water pressure [Raymond, 1987]. How this high water pressure comes about and how it is sustained are still unclear. The surge of Variegated Glacier has been modelled in terms of a linked cavity drainage system between the ice and an impermeable rock bed [Kamb, 1987], but not all surge-type glaciers overlie a 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 Western Antarctica [Engelhardt et al., 1990]. Previous investigations have established that Trapridge Glacier, a small sub-polar surge-type glacier currently in its quiescent phase, is underlain by permeable, deforming sediment [Clarke et al., 1984; Blake, 1992; Blake eta., 1992]. Stone [1993] suggests that Trapridge Glacier rests on a layer of relatively low-permeability till, perhaps one metre thick, which overlies till and alluvial deposits of much higher permeability. Smart and  Clarke [submitted] propose that in the winter, and at other times of low surface-water input, water generated at the bed by frictional and geothermal heat percolates through the low-permeability layer, into the higher-permeability layers beneath (from which it easily escapes by groundwater flow). There also exists during the summer a high-permeability (although spatially discontinuous) drainage system immediately under ice, above the low  1  Chapter 1. Introduction^  2  permeability-layer [Stone, 1993]. This drainage system has been investigated using borehole 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. Typical freeze-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 displays diurnal oscillations as shown in Fig. 1.1). Unconnected holes experience a pressure rise starting 6-12 h after drilling, when the hole freezes over. Blind holes experience a much more rapid pressure rise with repeated sudden drops, believed to be caused by hydraulic fracturing of the ice.) Knowledge of the hydraulic properties of the low-permeability layer is important to understanding the surge mechanism of Trapridge Glacier because this layer controls water escape 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 is low [Kamb et al., 1985]. This is precisely the time of year when the low-permeability layer controls the Trapridge drainage.) The only existing in-situ estimate of the hydraulic conductivity of this layer is an upper bound of 3.6 x 10-7 m s-1 obtained by Smart and  Clarke [submitted] from consideration of the minimum observable water level drop in an unconnected borehole. Clarke [1987] has estimated the hydraulic conductivity to be 1.27 x 10-8 m s-1, based on the granulometry of samples taken from the bed, and assumptions about porosity. Tavi Murray (unpublished data, 1992) obtained a value of 2.2 x 10-8 m s-1 from consolidation tests of forefield till, which is presumed to be similar to the basal material. None of these estimates is completely satisfactory and this fact motivates my attempt to establish a more direct, in-situ measurement of subglacial hydraulic properties. In this thesis I examine the freeze-in of unconnected boreholes and use this analysis  Chapter 1. Introduction^  3  a  100  50 0  1  < 200 0  cc  D 100 — I cn cn 0^1 cc 300  200  ^(77// c  100 0^ 1 TIME (DAYS)  Fig. 1.1: Types of Trapridge Glacier boreholes and freeze-in pressure records. Boreholes are 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) unconnected 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 hole does not reach the bed. Note the differing time and pressure scales.  Chapter 1. Introduction^  4  a  Fig. 1.2: Response to freezing forcing. These cartoons illustrate the two ways that excess water displaced by freezing can be accommodated: (a) by ice deformation, and (b) by water flow into the bed. In blind holes only ice deformation operates. In connected holes water flow into the bed dominates. It is unclear which mechanism dominates in unconnected holes. to estimate the hydraulic conductivity of this layer. When unconnected boreholes freeze shut, water in the borehole is compressed and its pressure rises. The magnitude of the rise depends on the hydraulic properties of the bed, as well as the rate of freezing of the borehole, 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 describes the mechanical aspects of a freezing borehole, including the ice deformation and groundwater flow of water into the bed, and an approximate freezing forcing obtained from temperature measurements made as boreholes freeze. The inclusion of a proper thermal model, such as that of Jarvis and Clarke [1974] or Humphrey and Echelmeyer [1990], instead of the approximate freezing would improve the model, but was precluded by limited time. I make many other simplifications to render the problem more tractable, including using only a single layer of deforming ice, neglecting non-hydrostatic background stress,  Chapter I. Introduction^  5  and 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 the use of a transient viscoelastic rheology [Shyam Sunder and Wu, [198913] to describe the ice; the inherently transient nature of the problem precludes the use of a simpler GlenLaw description. Compressible water is used in the borehole so that blind holes can be modelled. The resulting system of non-linear partial differential and algebraic equations is 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 estimate ice and bed properties, by iterative forward modelling until observed freeze-in pressure records are satisfactorily fit. The sensitivity of the results to uncertainties in inputs such as freezing rate, ice properties, and basal-cavity size (excavated by the hot water drill at the base of the borehole) was investigated. Hydraulic conductivity proved to be most sensitive to the freezing rate and the basal-cavity size, and essentially insensitive to the ice properties. I estimate the hydraulic conductivity of the low-permeability layer to be 1.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 for this parameter were so highly variable. Boreholes in Trapridge Glacier freeze from the top down, since the coldest ice is near the surface [Clarke and Blake, 1991]. Pressure sensors located 0.5 m above the bed continue to indicate water pressure for 1-2 yr after insertion, before they become isolated by freezing. Hence above every pressure sensor there is a water column which varies in length over time. The question arises of whether the response of this borehole could influence the pressure observed in unconnected holes, essentially filtering the "true" basal pressure signal. This question must be resolved before the pressure signals in unconnected holes (and those holes that, during winter, become unconnected) can be confidently interpreted.  Chapter 1. Introduction^  6  To settle this issue I use the model to investigate the effect of the borehole response  on the pressure observed in unconnected boreholes after the initial freeze-in period. The modelling results indicate that the borehole has a negligible effect on the observed pressure (except during the initial period of rapid freezing).  Chapter 2  Ice Deformation Model Development  2.1 Introduction To predict the pressure response of an unconnected borehole as it freezes shut after hotwater drilling, the mechanical response of the ice—bed—borehole system must be modelled together with the freezing forcing itself. The model must describe the deformation of ice in response to elevated pressure, as well as the flow of water from the bottom of the borehole into the bed. In this chapter the ice deformation model is developed. In Chapter 3 the bed and borehole equations are developed and the numerical solution to the coupled problem is described. The problem in this thesis is clearly transient, so an appropriate transient ice rheology must be used. This rheology must describe the response to time-varying pressure at the borehole wall, including pressure which first rises, then falls. Therefore the strain recovery behaviour of ice [Duval, 1978] must be properly described. The simple Andrade-Law exponential form [e.g. Paterson, 1981, p. 26] adequately describes the behaviour of ice under a suddenly applied steady load (i.e., in a creep test), but does not model strain recovery. The viscoelastic fluid relation of Morland and Spring [1981] and the viscoelastic solid relation of Spring and Morland [1982] model only the response to constant load or constant strain rate. A number of formulations seek to model full response, including strain recovery [Sinha, 1978; LeGac and Duval, 1980; Ashby and Duval, 1985; Shyam Sunder and Wu, 1989a; Shyam Sunder and Wu, 1989b]. Of these the physically based  7  Chapter 2. Ice Deformation Model Development^  8  phenomenological model of Shyam Sunder and Wu [1989a, 1989b] appears to best predict test results [Shyam Sunder and Wu, 1990]. For the present work I accept the model of  Shyam Sunder and Wu [1989b]. It must be noted that none of the foregoing transient rheologies describes ice deformation at large strains (tertiary creep), where recrystallization is important. Instead they describe the behaviour from zero strain up to the time when the minimum strain rate is attained (at between 0.5% and 2% [Budd and Jacka, 1989]). Thus the model that I develop is valid only for small strains. The freezing-borehole problem is inherently three-dimensional, but I simplify it to make it tractable. I assume the background ice stress to be purely hydrostatic. I consider the borehole to be a long cylindrical hole, so that end effects can be ignored and ice deformation assumed radial. Furthermore variations in ice properties (e.g., those caused by temperature change) and stress state with depth (the hydrostatic pressure in the borehole changes more quickly with depth than the background hydrostatic pressure in the ice due to the different densities of ice and water) are ignored, so that the deformation can be considered constant with depth. This allows the ice deformation problem to be reduced 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. The development is lengthy and complex, an inescapable consequence of the complexity of the physical model. Readers who have little appetite for mathematical details are encouraged to proceed directly to Chapter 4.  2.2 Governing Equations The governing equations fall into three basic categories: balance equations (mass and momentum), equations defining strain, and the theological equations. These are the  ^ ^  Chapter 2. Ice Deformation Model Development ^  9  starting point for the derivation of the ice deformation model.  2.2.1 Rheology of Shyam Sunder and Wu The Shyam Sunder and Wu (SS86W) rheology can be described as a transient, nonlinear, viscoelastic rheology with internal state variables. Shyam Sunder and Wu [1989b] present their 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 polycrystalline ice [Shyam Sunder and Wu, 1989b, p. 224], only primary and secondary creep are 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]. The exclusion of recrystallization effects, and hence the inability to describe tertiary creep, make this rheology unsuitable for modelling large-strain deformations. Because I am primarily interested in small strains, this limitation is acceptable. No existing rheologies can model the full range of primary, secondary, and tertiary creep. The SSSEW rheology is summarized on p. 231 [Shyam Sunder and Wu, 1989b], and restated, in my notation (which differs slightly from theirs), below. Notation and physical definitions are given in Table 2.1. Strain Rate Decomposition clei;^de7i^de:i dt^dt^dt^dt  (2.1)  Equations of State (in rate form) dui;^deli dt^Mid dt dafij =- dui;^dRi; dt^dl^dt  (2.2) (2.3)  Chapter 2. Ice Deformation Model Development  ^  SS&W  thesis  parameter description  A  A  al ... a6  c1...c6  B  B  D  Dijkl  E  E  G H  Giiki  fourth-order stress transformation tensor  Hijki  fourth order strain-rate transformation tensor  k  .i--/  temperature-independent "isotropic hardening parameter" (my terminology)  R S*  Ri; S.i'i  elastic back stress tensor  Sa*  S'ciii  pseudo-deviatoric reduced stress tensor (for isotropic materials this reduces to the deviatoric stress tensor)  SR*  Sliij  back stress measure tensor  0  77  E  Ei;  total strain tensor  "kinematic hardening parameter" (my terminology) anisotropic parameters scalar drag stress associated with isotropic hardening orthotropic elastic stiffness tensor Young's modulus  applied loading tensor, pseudo-deviatoric stress tensor  constant which scales the equivalent stress  <17  elastic strain tensor  ct  ,ii  transient viscous strain tensor  c,  c:i  steady-state viscous strain tensor  A  A,  "steady-state viscous softness parameter" (my terminology)  Ad  At  "transient viscous softness parameter" (my terminology)  cr  a-ii  stress tensor  ad  ,ij  stress difference tensor, reduced stress tensor  Ee Et  Table 2.1: Notation and physical definitions of Shyam Sunder and Wu rheology.  10  Chapter 2. Ice Deformation Model Development ^  11  where dR•;^deki ^ — AE 'kJ^• dt^3^3 dt  (2.4)  Evolution Equations deYildt  (del dt)0  de:j I dt  AtSli  (de I dt)0 dB dt  HE ( de) crd,eq kdt )  The 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 strain be ^v  ^t  = fij  gives ^dey-^de5 def. ^1 ^1^ + ^dt ^dt^dt and ==E  75  "  Equations (2.11) and (2.12) replace (2.1). The preceding equations and those appearing throughout this thesis are in tensor notation. Unfortunately many of the same variables are defined in Shyam Sunder and Wu [198913] in "engineering" vector notation. In this convention the symmetric stress and strain tensors are represented as six-element vectors jr. = E^  [eyy,  Cryy) Crzz, Crxy, Cryz, Crzxl  (2.13)  7zxi T  (2.14)  Eyy, CZZ 7xy, 7yz  ^  Chapter 2. Ice Deformation Model Development ^  12  where -y = 2E are the engineering shear strains, equal to twice the conventional shear strains. 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-rank transformation 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^--3C1C2C3 c*2 C*2^  C*2  3 (c1 c3) 4^--3c1c2c3  H  c*2^c*2 3 (ci^c3) c*2  0  (2.16)  1 (2c4) 1  SYM.  (24) 1 (2c6)  [Shyam Sunder and Wu, 1989b, p. 229], where c* c2c3 + c3c1 + c1c2. The term (de/dt)0 may 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 [Shyam  Sunder and Wu, 1989b, p 228] N —1 = (1 )N eq  v^  2^3  Creq^ (id TG  (2.17)  (2.18)  Chapter 2.  Ice Deformation Model Development  13  = 1^ (ci +c3 )  ^  (2.19)  — -31 C1^— -3 C1  I kci^c3)^--3c1  51 (C1 +  G=  0 C3)  (2.20) 2c4  SYM.  2c5 2c6  where  77  = Vo exp (Q/NRg„,T)  (2.21)  =^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 state  of stress is uniaxial (i.e., a = (0, cryy, 0, 0, 0, 0)). This will result in creq being the usual equivalent uniaxial stress. Equation (2.21) is the commonly used Arrhenius Equation for the temperature dependence of the viscous stress factor V (which is loosely equivalent to the flow law parameter A in the Glen flow law), and is generally accepted to be valid for temperatures below —10°C [e.g., Shyam Sunder and Wu, 19894 At higher temperature the viscous stress factor changes more rapidly than (2.21) predicts, due to the presence of 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- and low-temperature expressions must, however, give the same value at V at —10° C, which would 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 using  ^ ^  Chapter 2. Ice Deformation Model Development ^  14  the low-temperature expression. Hence V(-10 ° C) V0,11 exP (QH/NR9as71263) — Vo exP (QL/NRga5T283) where T263 = 263.12 K. Thus VO,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 At  = 3 ( 1 )N,TN-1 d,eq  (2.25)  3  crd2,eq = —CrdGe/d^ 71  (2.26)  :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 stress  evolution equation (2.7) is defined as (c/c/dt)d,eq (de/dt)0 [Shyam Sunder and Wu, 1989b, p. 231].  (Crcl,eq) N  BV  (2.29)  Chapter 2. Ice Deformation Model Development^  15  2.2.2 Balance Equations The relevant balance equations here are those of mass and momentum. Energy balance is not important to the mechanical model. In Cartesian coordinates we have mass balance [e.g., Mase, 1970, p. 126] op a --(97 5Ti(pvi). 0  (2.30)  and momentum balance [e.g., Mase, 1970, p. 128]  Ph+  &xi;^dv; xi  = P dt  (2.31)  2.2.3 Definition of Strain and Strain Rate The rheology of Shyam Sunder and Wu [1989b] was derived under the assumption of small strains [Shyam Sunder and Wu, 1989b, p. 225]. Hence it is acceptable to use the infinitesimal 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 axi  (2.32)  The strain rate is (for any state of strain) [e.g., Mase, 1970, p. 113] deij^1 (Ovi  di^2Oxi+ Oxi  (2.33)  2.2.4 Compatibility of Strains Arbitrarily chosen strains will not in general yield a single-valued, continuous displacement field (because the definition of strain provides six equations to determine the three components of displacement). A real displacement field is guaranteed to exist if the components of strain satisfy the compatibility equations of strain. These are [e.g., Mase,  •▪  ^  Chapter 2. Ice Deformation Model Development ^  16  1970, p. 92J: u612 , 402 f 22 ^2 n 8X1^(9X1^UXiCIX2 ^02 623 02E22 02633 1-^2 04 04^10M20X3 .92633^02611^az 2 8 aX? OX3OX1 02611  34  az eii a^8623 8631^0E12^ .  az, ( oxi ax2 ax3^8x28x3 8623 _ 0631+0612^ a (^ 82622 _ . 8x1 ax2 8z3^axoxi ^ax2 49^8623 , ,8631^.8612 -t- r,  ax3  (  ax,  ox2 ax3  .  ^  a8.2:83:2 •  (2.34) (2.35) (2.36) (2.37) (2.38) (2.39)  Similar expressions for the compatibility of strain rates also exist, but are not required because compatibility of strain rates is assured by the compatibility of strains.  2.3 Application of Limiting Assumptions The following assumptions will be applied in this section to simplify the problem equations: isotropic ice, very slow flow, radial symmetry, small vertical gradients, and infinitesimal strains.  2.3.1 Isotropic Ice It will be assumed that the ice surrounding the borehole is initially composed of randomly oriented crystals, and hence is macroscopically isotropic. While this is unlikely to be true in glacier ice, the assumption of isotropy leads to enormous simplification of the rheology and no information exists that would justify a more complicated assumption. For an isotropic material, (2.8) describing the relationship between stress and elastic  Chapter 2. Ice Deformation Model Development^  17  strain 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 function of stress [e.g., Muse, 1970, p. 143] 1 e7 = 2/1cYij  ^A 2p(3A -I- 2)  akkk.  (2.41)  For an isotropic material the coefficients c1...c6 in (2.16), (2.20), and (2.22) are all equal to unity [Shyam Sunder and Wu, 1989b, p. 228]. Substituting these into (2.22) yields  n = 2. Furthermore the transformation matrices G and H  become  2/3 —1/3 —1/3 2/3 —1/3^0 2/3  G=  (2.42)  2 SYM.  2 2  and 2/3 —1/3 —1/3 2/3 —1/3^0 2/3  11=  (2.43)  SYM. 1/2 These transformation matrices can now be used to evaluate the pseudo-deviatoric stress vectors .C"`,  §:ri,  and fill".  Chapter 2. Ice Deformation Model Development^  18  Substituting for G in (2.19) gives 2/3 —1/3 —1/3  Crxx  2/3 —1/3^0  0" VY  2/3  azz  =  2  Crzy  SYM.^2  Cfyz  2  CrZX  2^1^1 axx 3 0Jl - azz  + 0. ___  -  3 ^  -  1^1^2 5rCrxx - ayy  3 yy 3 zz  5crzz  2o-xy  (2.44)  2tryz 2o-zz 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,  1  Crij = C•ij -01k8ij  •3  and expands to the individual terms 1, cru. =^— —^+ 0'22 + ass) 3 1 ^ G-22 = 1722^( a22 033) 3 1, C733 = 0.33 - - Vrii (r22 + 0.33 ) 3 V12 = (712 U23 = 0-23  Cr31 = U31-  (2.45)  Chapter 2. Ice Deformation Model Development^  19  These can be translated into engineering vector notation, using the definition (2.15) 1 (ri„  40-•^■-• yy^t'zz)  —1 (cr.. -I- 2o Cr  =  -  ^crzz)  (crz. — crw 2crzz)  (2.46)  2crzi, 2o.yz 2cr z.  Note that (2.46) is identical to (2.44), the expression just derived for ,§*• Therefore it is apparent that for isotropic materials kr = CI, and that the matrix G is equivalent to a deviatoric operator (because G = 1-6) . Thus  = = = or, in tensor notation, Crij  (2.47)  d Crij A"ci^ S  (2.48)  * SR  (2.49)  Now the definition of the equivalent stress can be simplified. Substituting for 77 and GO in (2.18) using 77^2, (2.19), and (2.47), eq —  3 —T CT —cr 2  Chapter 2. Ice Deformation Model Development ^  20  (2Crxx — Cryy Crzz ) —1 2o-yy 3(a^ xx „.2 " eq  —  3r —2 La.., airy, zz, Crzy, (ryz, o]  =er _ :x —  — Cize  —1 3(a xx— cryy 2U)  1^1^1  1  —2o-xxoyy — —2o-zzcr.. io-z„ayy ti — —2 Crvy trzy  1^1 — — crzzo-zz — — cryyc•zz  2  3  crz2z  2azy 2cryz  2crzz  1 (0.2. _^4_ 0.2 ) 4_ 1 _ _ (0.2 _ 20.yy0.zz 0.z2z) 2^ 2 YY  1^, +— (T" — 2crzzaz. cr:x)^3o-yz 3azx 2^zz  —1 ( Cr2 — 2 CrxxCryy cr2^—1 (cr2 — 2cryycrzz o)  2 k "^ 2 1111 1 +— (0-2 — 20-zza= cr:z) + 3 [crzy Cryz crzz] 2^ zz  1 — [(azx — ayy) 2 + (Cfyy — Cree) 2 + ((Tye — crxx) 2 } + 3 [o-zy  2  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] creq 2 = 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 to 1  _ =— 6  I  —^ 0-22)^ -r v722  — u33) 33)  ,^  k0.33  ,^  kai2-1-  0.23 +  0'30-  Hence  1 2^  Cre  q  ——  2  [kali —  ,^ , f^N^ 2 .3 I Cr22)^k(022 — Cr33) --r- “733 —^ } -I- k0-12+ Cr23  a31)  (2.52)  21  Chapter 2. Ice Deformation Model Development ^  Note that (2.52) is identical to (2.50) just derived. Hence (2.51) may be used for the definition of  creq•  Therefore, from (2.5), (2.17), (2.47), and (2.51), the steady-state viscous  component of the strain rate is given by deyi dt jLJjJ  (del dt)0  3 (  \N N-11  3^)N  eq ^(N —1)12  (2.53)  —2 V  This is identical in form to the commonly used Glen Flow Law for "steady-state" viscous deformation of ice [see e.g., Paterson, 1981, p. 31]. Similarly, it can be shown that 2^3 1 c Cr d,eq^crip'ij•  (2.54)  The evolution equation (2.4) for R must be examined. The vector equivalent of hod de! I dt is Hdetl dt. Thus in vector form (2.4) becomes, with  77 =  2,  d;t 2 = — AEH dt^3^dt  (2.55)  Expanding H det/ dt using (2.43) for the isotropic definition of H gives 2/3 —1/3 —1/3 2/3 —1/3^0  H  t  dt  2/3  =  de/dt den t /dt detzz/dt  2 detxy/dt SYM.  2 detyz/dt 2 dezz / dt  Chapter 2. Ice Deformation Model Development^  22  t Idt — dezzldt) 23- (2 dexxl dt — dcyy t /dt — detzz/dt) (—dexzldt + 2 deyy d-;t dt  —  (—det../dt — /dt + 2 detzz/dt)  (2.56)  (Idyl dt deytzldt de/dt  Because de!a/dt is proportional to 'crej, for which the first invariant vanishes, the first t Idt dezzldt = 0 and invariant of de:i/dt also vanishes. Hence dext.Idt deyy dd.! dt deyy t dt  dt  =  dezzldt  Zet  2 dexyl dt  dt  (2.57)  2 d4zIdt 2 detz.ldt or, in tensor notation,  de de^de • Hi •'kl^= 8i1c87^= dt^dt^dt •  (2.5 8)  Therefore (2.4) becomes  dRii 2 AEd4i dt^3^dt •  (2.59)  In summary, the rheology for isotropic ice is described by the equations ==^dh EZi^  (2.60)  de' V deBj 2i^,J^ (2.61) dt^dt^dt^ 1^A akkSia^ (2.62) e • =^Cli3^ 2/.z (3A + 2f.t) — v —^  (2.63)  ^  Chapter 2. Ice Deformation Model Development^  dRi; 2 AE detii dt^3^dt  dt de!i  =^d  dt dB  (2.64)  ( —(16^ dt)0  (2.65)  (dE)^  (2.66)  -1 1t o  E 1 crd,eq)N cra,eq BV )  dt  23  (1-3-117 )N a N-1 d,eq^  (2.67) (2.68)  where ( de)  = 1s-1^ 3,  / =^ creq^ a i- Jai.; 2  ,....2^3/ di d ' d,eq = — Cri• Cri .  2 33  A, =^( 1 )N crit-i 2 Vl^q N 3(1 N-1 At = crd ,eq •^ 2 ./3 .17)^  (2.69) (2.70) (2.71) (2.72) (2.73)  2.3.2 Very Slow Flow This assumption simplifies the momentum balance equation by eliminating the inertial term, thus reducing it to a force balance equation. Slow flow implies that accelerations are very small, and p dvildt < pfj,acrii/Oxi. Thus momentum balance (2.31) becomes acri • Pf, + ^3 = 0-^ uf Xi  (2.74)  2.3.3 Radial Symmetry Here I convert from Cartesian to cylindrical coordinates, and apply the assumption of radial symmetry.  Chapter 2. Ice Deformation Model Development^  24  First consider the definition of strain (2.32). The gradient of a vector field in terms of cylindrical coordinates is [Prager, 1961, P. 36] Oui oxi  -  (r) (r) °Ur^(r) (9)^  (r) (z) aUz ei ei^- F ei ej^+ ei ei or  CU?^)^  (0) (r) 1 (0) 0)1 (au 0 + ei e^ `•^— — us + e• e • — — + u r ao^r 00^r  (e) (z)1 auz ei ej — r 00  +e(z)e(r) OUr ecoe(e)Oue e(.)e(z)auz  (2.75)  3 az^3 az^3 az  With radial symmetry alae 0 and ue = 0 so that the gradient of the displacement field becomes au  Oxi  eloe(r)2t-`1. + e(r)e(:9) (0) + e(r)e(z)auz 3 Or^3^2  3 Or  er (0) + er^. e(8) + e• e 3 (0) r^3  ^+er  +elz)er -b--+ elz)er (0) + e1 e 5(z) au z .  az  (2.76)  Substituting this definition into the definition of strain (2.32) yields  Our  err =^  (2.77)  fee =^  (2.78)  Or  err  auz  Oz^  (2.79)  ere = 0^  (2.80)  ee. = 0^ 1 (auz au?) erz =^ + az 2 ar^  (2.81) (2.82)  Similarly the strain-rate components become de,.,.^Dv,.  at^Or  de00vr ^ =^ dt  (2.83) (2.84)  ^  Chapter 2. Ice Deformation Model Development^  de„^avz dt^az dere =0 dt dee. = 0 dt derz^1 (avz avr) dt^2 ar + az )  25  (2.85) (2.86) (2.87) (2.88)  Consider the mass balance equation (2.30). Expanding it by carrying out the spatial differentiation gives  op^Ova^Op +19:9X; + vi azi = 0.  (2.89)  Noting that Op/at vi 8p/Oxj = dpldt, the material time derivative of density [ .g.,  Prager, 1961, P. 74], we have avj dp P axj^dt^°.^  (2.90)  The divergence of a vector field expressed in cylindrical coordinates is [Prager, 1961, p. 36]  avj av,  ( ay,^avz ax,^Or + -7: a° +vrj^az • With radial symmetry vo = 0 and  1  alae  (2.91)  0 so that mass balance (2.90) becomes  av, vr avz\ dp P^+ 7 +^+ Tit  °•  (2.92)  Next consider the force balance equation (2.74). The divergence of a tensor field expressed in cylindrical coordinates is [Prager, 1961, p. 37]  ^3  + 1 acrer + (Tr,^coo + aux,Or^r ae^az  00r,. acrii^ ^ ,- e(r) ^  ax,^3  +e  (^ Ocrre +1 ace° + aro + as,^au + ze ^ + 59  ^  ar^r  ae^az  (aurz +1 (threz 4_ crTz) + OCTzz) +e^ Or 09^az • 3^ar^r  (2.93)  Chapter 2. Ice Deformation Model Development ^  With radial symmetry  26  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 the  stress field reduces to  ao-i; ax,  Carr 1 ,^ = e.(r) ^ + ((Tr^cree)+ aaz r) 3^ ar^r^az +e(z) Car.. + crr. acrzz Or^r^az ) •  (2.94)  Thus the force balance condition yields two non-trivial equations  Pfr  acrrr^1 ^ +(crrr aeo) acrzr) = 0 ^(2.95)  ar^r^az  oar..^CrTZ^aCr2.2  +^+^= pfz+ ^ Or^r^az  0.^(2.96)  Gravity is the only body force acting on the ice. Therefore the force balance equations (2.95) and (2.96) become  (aarr + 1 (arr ow) Or^r  pg +  (9;zr)  (aarz arz acrzz) 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 independent of hydrostatic pressure. For linear elasticity, hydrostatic pressure does not affect the solution, except to cause a small change in ice density, which is second order, and  hence insignificant. We can therefore neglect the effects of vertical gradients in crx,. and CTZT  •  Applying this simplification the mass balance relation (2.92), becomes (2.99)  Chapter 2. Ice Deformation Model Development^  27  and the force balance relations, (2.97) and (2.98), become aurr , 1 f -r- kurr — Or^r  Cr t 9 0) = 0  aa. = 0. pg + ^ az  ( 2.1 0 0) (2.101)  Since crzz = —p 'crzz, and a'azziaz = 0, (2.101) becomes  Pg  OP oz  = 0,  ( 2. 1 02 )  the standard equation for hydrostatic pressure. Since the hydrostatic pressure has an insignificant effect on the flow solution, this equation has no influence, and I shall henceforth 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 Our  Err^  Or  (2.103) (2.104)  coo E rr  =0  (2.105)  Er  =  (2.106)  fez = 0  (2.107)  - 0  (2.108)  EZT  and the strain-rate expressions (2.83)—(2.88) become derr  (2.109)  dt (1E09  yr  dt 0  (2.111)  0  (2.112)  0  (2.113)  = 0.  (2.114)  dt clEro dt dEez dl (kyr  dt  (2.110)  =  ^ ^  28  Chapter 2. Ice Deformation Model Development^  Hence, 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 of the 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 the  diagonal terms remain. The notation for the diagonal terms of the tensors can now be simplified, so that o„ becomes or, and so on. In the simplified notation the ice rheology becomes E,. = E +^  (2.115)  €8 = fee + Eto'^  (2.116)  0 =  cez + evz^  ^de:).^de:.^de: =^ ^dt ^dt + dt dee'^dee de: = ^dt ^dt + dt del^det^des =^z + z^ ^dt ^dt^dt 1^A  21i  E,e. = — ar  1  2p, (3A + 2/t) A  2p (3A + 2/4) A  =^ ez^ — Cr z ^ 2p^2p, (3A + 2) (a'. + a° ± d CI r = ar —  (2.118) (2.119) (2.120)  (ar ± ao + az )^(2.121)  Ee^—a, =^ 211,^^ (ur ± CO ± Cr z) 1  (2.117)  (2.122) (2.123)  gz)  Rr  (2.124)  d ^go^as = -Re  (2.125)  d ^Rz uz -^ a=z  (2.126)  29  Chapter 2.^Ice Deformation Model Development  dt  dR9 dt dRz dt  2  de'  3 2  dt dct  dt  3 2  dEt  3  dt de) 0  de,'9 dt  ^ ^(1) Asicrs^0  desz^  dt der  (2.128) (2.129)  (2.130) (2.131)  de) 0^ Azicrz (7i  (2.132)  At' ad cA ()^  (2.133)  dt^dt  dee = Azi4 ((A)^ dt^dt detz dt  (2.127)  At =^'acz‘  (1)  dB^1 y cr d dt^BV^eq  (2.134)  (2.135)  (2. 13 6 )  where  1s'  (2.137)  3 ( N N_i 2\ V) "  cr„ =  (To Crz  -23- (iur2^'1792^icrz2) 1 ar — — ar ao ± az ) 3 1 — — (ar^+ az)  3  1 o-z — — (q,.^aer + az) 3  (2.138) (2.139) (2.140) (2.141) (2.142)  ^^  Chapter 2. Ice Deformation Model Development^  3( 1^ )N d N-1 = At^creg 2 BV  ad  eg  i d  3 ( IF^ .^  2 + isre + ,012)  2\  1 (ard + 4 ..^ ad — — 3  ar^r  + azd)  /ad^ = erect _ 1_ (0.7.d + crcj + azd) e 3  / d^ad =^— 3_1 (ard + 4 + 0.zd) cr.^.^  30  (2.143) (2.144)  (2.145) (2.146)  (2.147)  It 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 equation (2.99), we have  (de.^deo dez)^dp = 0. P^ dt^dt^dt )  (2.148)  The term in brackets, dekkldt, is the cubical dilation rate, that is the rate of change of volume of an element. Because the mass of an element (m = pV) is constant  ^dm^dV dt = ° = P dt^dt and  dp^p dV dt^V dt.  (2. 14 9)  Substituting (2.149) into (2.148) gives  ^P  ( de,.^deo^dez^p dV^0. dt^dt^dt ) V dt  Simplifying and rearranging,  ( de,.^dco^dez)^1 dV dt + di + dt ) V dt which is just the definition of the cubical dilation rate. Hence it is apparent that the mass 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^  31  The compatibility of strain equation can now be developed. One possibility would be to 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 strain definitions in cylindrical coordinates. From the definition of tangential strain (2.104) it follows that V.,. =Te o.^  (2.150)  Substituting this into the definition of radial strain (2.103), and differentiating, gives , er^— krce)  ar  r  569  Or  Or 569  r^— 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 Strains The assumption of infinitesimal strains has already been invoked a number of times in the preceding derivation. Here I use it once more to set the advective part of spatial time derivatives to zero, and allow the preceding equations to apply at spatially fixed nodes. As a result of this assumption all the material time derivatives may be replaced by spatial time derivatives. Hence the ice rheology becomes (2.152) (2.153) (2.154)  ^  ^r  Chapter 2. Ice Deformation Model Development ^  aft,:^ae^SE. at^at^at (34 = afte .94 at at + & af: _ aft 06: at — at + &  32  (2.155)  =^r + r  (2.156) (2.157)  1 A = —cr, ^ fry + cre + crz)  (2.158)  crz)  (2.159)  (cr, cre + az)  (2.160)  c^  2p^2p (3A + 2p) A e^1 =^ ee — (ar^+  ^2p^2p (3A + 2p)  ^1 ^A =  2tt az 2(3A + 2) "r  =  (2.161)  crea^era —  (2.162)  or^crz — Rz  (2.163)  OR,.— =  ^at^ToT  (2.164)  ^at^3^at 20 4 aRz^—AE a 3^at ^at^  (2.165)  ^Me^2^ete —AEa  (2.166)  (2.167) (2.168) (2.169)  (2.170) (2.171) (2.172)  Chapter 2. Ice Deformation Model Development^  aB — = E()N d Cr^. N1 BV " at^  33  (2.173)  where (de 0 =  1 s-1  (2.174) (2.175) (2.176)  / c•r  =  (2.177)  /ao  =  (2.178)  =  (2.179)  =  (2.180)  At  = \/_3 (42 + /act + ,(712) 2 4) a.r = '4 3k  (2.181)  0'0 = '4  (2.183)  d o-eq  4 — 1 (4 + 4 +  4— _1 (4 + 4 + 4) 3\ 1^1 (crd 4_ 4 + 4) •„,s1- ' z = crd.^ _ _13^ T .  (2.182)  (2.184)  The the other relevant equations are summarized below: Compatibility of Strains  ace - c,. = 0^ co+ r^ Dr  (2.185)  Momentum Balance  par 1 +- o-e) 0 Dr r^  ( 2. 1 86)  Definition of Strain  Er  (2.187)  Chapter 2. Ice Deformation Model Development^  34  Ur es = —  ( 2. 1 88 )  ez = 0.^  (2.189)  r  The above equations are in no particular order. They can be grouped as follows to form the field (differential) equations of an initial value problem, with algebraic constraint equations. 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 other equations are used to calculate the intermediate variables needed for the differential and constraint equations.  2.4 Initial and Boundary Conditions Here I will state the initial conditions and boundary conditions, which together with the field equations just developed, fully specify the initial value problem in the ice. These follow directly from the assumptions of the model, as outlined in section 2.1.  2.4.1 Initial Conditions The 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 scalar drag stress, which depends on the grain size of the ice [Shyam Sunder and Wu, 1989a, p. 56].  Chapter 2. Ice Deformation Model Development^  35  2.4.2 Boundary Conditions The condition on the inside boundary is fairly obvious. The radial compressive stress must be equal to the water pressure in the borehole, and the displacement at the borehole wall must equal the change in borehole radius. Thus cr,. rb , t = —Pb  ^  u,. (rb,t) = Arb.^  (2.194)  (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), yielding E9  (re,, t) =  Arb rb  (2.196)  At the outside boundary there are numerous possibilities. Although the outside boundary in the problem as formulated so far is at infinite radius, in the numerical problem this boundary will have to be placed at some large, but finite radius. The choice of 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 compressive stress only (equal to hydrostatic pressure). The problem with zero strain is that if the boundary is moved to finite radius there will be no flow unless the ice is elastically compressible. It would be useful if the model could be run with elasticity turned off, to facilitate comparison with analytic solutions. The second and third options are similar. The specified-radial-stress option does have the advantage that some analytic solutions 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; Jaeger and Cook, 1969, p. 128]. Therefore  the outside boundary.  I will apply the specified-radial-stress condition to  Chapter 2. Ice Deformation Model Development ^  36  Thus, to summarize the ice boundary conditions:  (rb, t)^—Pb  (2.197)  Arb rb  (2.198)  cr,. (oo , t) = —Po.  (2.199)  Es (rb, t) =  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 deformation model developed in this chapter, to form the complete coupled model.  Chapter 3  Bed Equations, Forcings, and Model Synthesis  The flow of water in a thin sheet aquifer when flow is induced by borehole response tests has been modelled by Stone and Clarke [1992] and Stone [1993] using a methodof-lines solution on a radially symmetric grid. Smart and Clarke [submitted] consider the same problem that I do here; namely, to determine the hydraulic conductivity of the low-permeability layer beneath Trapridge Glacier by considering the behaviour of water in an unconnected borehole. They find an upper bound to the hydraulic conductivity by considering the bed to be a uniform half-space bounded by impermeable ice, except at the bottom of the borehole where a hemispherical water-filled cavity exists. They consider the steady-state analytic solution for Darcian groundwater flow. In this model I also approximate the bed as a homogeneous half-space, with hydraulic properties equal to those of the low-permeability till layer. This is a reasonable approximation as long as drainage channels are not too close to the hole, and the radius of the hole-bottom cavity is small relative to the thickness of the low-permeability layer. I further assume that the flow is described by the standard diffusion equation for transient groundwater 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 that the only motion of the soil matrix is the small vertical movement caused by elastic deformation of the matrix as effective pressure is changed. The ice and bed are coupled by conservation of mass of compressible water in the borehole. The borehole freezing thermal problem has been considered thoroughly in connection  37  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  38  with the interpretation of thermistor measurements and hot water drilling in cold ice [Jarvis and Clarke, 1974; Humphrey and Echelmeyer, 1990; Humphrey, 1991]. Models such as the finite-element model of Humphrey and Echelmeyer [1990] describe the freezing of 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 the freezing forcing. Unfortunately time constraints prevented me from implementing such a 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 as a function of time is described as a series of decaying exponentials, having coefficients that can be determined from temperature measurements made as a borehole freezes. I also allow for arbitrary sudden pressurization of the borehole, to permit the modelling of sudden pressurization and decay events. In this chapter I develop the bed model, borehole model, and description of freezing and other forcings, then synthesize these with the ice deformation model to form the complete coupled model, and formulate the problem for numerical solution using the method of lines.  3.1 Basal Groundwater Flow Consider a homogeneous, isotropic, permeable half-space bounded on top by impermeable ice, and by water in a hemispherical cavity at the bottom of a borehole. It is postulated that such a cavity is formed during the hot water drilling process. Let the coordinate system 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 [Blake et al., 1992], it will be assumed that this deformation is slow enough that the standard  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  Fig. 3.1: Bed geometry and coordinate system.  39  ^  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  40  equations for groundwater flow in a porous medium still apply. Taking the till deformation into account would introduce much more complexity than is justified in such a rough model. The standard water mass balance equation in a saturated porous medium [e.g., Bear  and Verruzjt, 1987, p. 63] gives  ^, Oh^ ^-v^Pwg (a +^7113)-07  (3.1)  where is specific discharge (or volume flux) of water, pw is density of water, a is soil matrix 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 general form of Darcy's Law (which applies to low-Reynolds-Number flow through a stationary matrix) [e.g., Bear and Verruijt, 1987, p. 38] is written  Oh = —Kij(zi) Q^  (3.2)  where Kij is the second rank tensor of hydraulic conductivities. I will further assume that the basal material is hydraulically homogeneous and isotropic, in which case Darcy's Law simplifies to  qi  = —K  Oh Ox;  (3.3)  or, 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-known governing equation for transient flow in a compressible porous medium  Oh ICC1211 = pwg (a + ni3)^  (3.5)  ^ ^  ^r2  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  41  This equation can be expressed in spherical coordinates and simplified to take advantage of symmetry. First express Viz in spherical coordinates [e.g., Spiegel, 1968, p. 126] v2h 1 a  Or  (7. 2  k  With spherical symmetry  a 1. 8ah\ + 1 02h (3.6) Or^r2^e ao^ao^r 2 2 9 802.^  ah\ +  moo  1  0 and moo 0, hence (3.6) simplifies to give  10 al ) .^ V2h =^(r2-±. r2 Or^Or Thus (3.5) becomes  ,^(2 1 a r- —) ah —^Oh r2 ar^ar^Pwg(a+n13)Wi  (3.7)  (3.8)  which can be rearranged to give Oh ^K ^1 a ( 2ah\ at - pw g (a + nfl) r2 Or kr Or ) •  (3.9)  The initial conditions and boundary conditions for the basal groundwater flow problem follow directly from the assumptions that the bed is initially in a no-flow (uniform hydraulic head) state, and that the background head in the bed changes much more slowly than the head in the borehole. The head at the cavity wall must also be equal to the borehole head. Hence h (r, 0) = po I pwg, r, <  < 00  h (oo,t) = polpwg h (re, t) = pb(t)Ipwg, t > 0 where re is cavity radius (the interface with the borehole), polpwg = ho is background head 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 a uniform permeable half-space, initially at constant head, to an arbitrary pressure forcing  Chapter 3. Bed Equations, Forcings, and Model Synthesis^  42  at a hemispherical interface is given by  Oh ^K r < 00 etpwg(a+ =(2ah nt3) r2 -6; r^  (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 Equations The borehole is filled with compressible water, which acts to couple the ice and bed. The coupling is provided by the water mass balance equation for the borehole. The mass of water in the borehole plus the mass of water lost by flow into the bed or freezing (see section 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 which has flowed into the bed between time 0 and t, mf is the mass of water which has frozen between time 0 and t (see section 3.3), and mo is the mass of water in the borehole at time 0: The mass of water in the borehole at any instant is derived as follows. The compressibility of water must be taken into account; hence the density of water may in general vary in space as well as time. Thus  mw = f w pw dV v where Vw is the sum of the volume of the borehole VB and the volume of the basal cavity  Vc. The depth of the cavity is small, so that the density of water in it can reasonably  Chapter 3. Bed Equations, Forcings, and Model Synthesis^  43  be considered constant. Density in the borehole cannot, however, be considered constant since the pressure change with depth is of similar magnitude to the pressure change with time. First, consider water mass of the cavity  m(t) = pc(t)Vc = pc(t)17rr!  ^  (3.18)  where pc is cavity water density, given by the assumed equation of state for water  Pc(t) = Po exP {/(3 (Pb (t) 1)0)]  (3.19)  2 m(t) = 1rrpo exp (Pb (t) — po)]  (3.20)  Thus  Next, consider the water mass of the borehole itself  mb(t) = iv pw (z ,t) dV^  (3.21)  where  dV = rr(Odz Pw(z,t) = po exP (P(z, t) — Po)1.  ^ (3.22)  Here p(z,t) is  p(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 pressures considered are never especially large, pw can be considered a constant Po in (3.23). Hence  p(z,0 = pb(t)— pogz.^  (3.24)  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  44  Substituting (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 simplifying mb(t) =  rL(t)  jo1  Po exP (Pb(t) — Po)] exP [—PPogz]rrl(t)  dz  1(t)  I^exP [—PPogz] dz Po exP b6 (Pb(t) — po)] irrl(t) Jo 1 (t)PoexP [P(Pb(t) — 790)] [_ opog exP (—Opogz)]  rr(t) g  exp [p(pb(t)— p0)][1  Note that the height of the water-filled borehole  —  L(t)  exp (—PpogL(t))1 .  L is  (3.26)  time dependent. This allows for  the approximate freezing forcing developed in Section 3.3. Borehole radius  rb(t) can  be  found from the tangential strain in the ice at the borehole wall. The ice strain boundary condition at the borehole wall (2.198) leads to fe(rb,t) _  Arb(t) rb(t) rb(t) — rb,o(t) ?b(t)^•  For infinitesimal displacements rb(t)— rb,o(t)^rb(t)— rb,o(t) rb(t)^rb,0 and rb(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 in  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  mw =  45  R-7-1(t)  exp V) (pb(t) — AO] [1 — exp (—,3 pog LW)] 13 g 2 ,^r„ ^„ + —Irr-po exP E (Pb (t) — P0).1 • 3 '  (3. 28)  Evaluating this at time t = 0 gives mo, hence MO =  „ 2 3 IrT/'° exp [,8 (po — RI)] [1 — exP (-0^ PogLo)1+-3- 1'7'W° exP [fl (Po — Po)] 13g 2  =- rrb,o,0 r1 _ exp FfipogL0)} + i2irr:poe0 )3g^1 2 2 --= 7rrb'0 [1 — exp (-13pogLo)] + '71- 7.3cP0i fig  (3 .29)  Now let us consider the flow of water from the borehole into the bed. First, Darcy's Law must be converted to spherical coordinates and simplifying assumptions applied. In spherical coordinates the gradient can be written [e.g., Spiegel, 1968, p. 125]  ^a Or^r ao e r sine ao and for spherical symmetry 0/89 = 0 and a 1 ao = 0, so that v _  e(r) (9 + e(9)-1 a + „(0) 1  V' =  e(r) a  (3.30)  (3.31)  Or  Substituting (3.31) into Darcy's Law (3.4) yields the radial specific discharge qr =  ,Oh -Or .  (3.32)  —A  Equation (3.32) gives the volume flow rate per unit of area (the flux) across a hemispherical surface centred on the origin. The total volume flow rate through such a surface is simply the flux multiplied by 277-,2, the surface area of a hemisphere. Therefore the volume flow rate into the bed (the flow rate at the cavity surface) is Qw (r,) = —271-7.1C :-711  (3.33) r,  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  46  and the mass flow rate is  Q, =  —27rperif  ^  ( 3.34)  rc  Substituting (3.19) for pc gives Qm  :K — = —21-po exp (pb(t) — po)lr^  (3.35)  re  The total mass of water lost by the borehole to the bed is the mass flow rate integrated over time. For compatibility with the rest of the problem, this can be•expressed as a differential equation thrtaq^Oh  = —2rpo exP [fi (Pb(t) — P0)]r2K — ar at^  (3.36)  subject to maq(0) = 0. 3.3 Forcings I shall consider three distinct types of forcing: volume forcing caused by' freezing to the top of the borehole, sudden pressurization, and step pressure change. 3.3.1 Freezing Freezing 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 arbitrary rate, with no freezing to the borehole walls. Under this simplifying assumption freezing forcing enters the model by specifying L(t), the height of the water-filled borehole as a function of time. The most important place L(t) appears is in the term M1 in the borehole mass-conservation equation (3.17). Let the height of the water-filled borehole as a function of time be given by some arbitrary function having coefficients that can be assigned to approximate the desired  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  47  freezing rate. As a functional form, I have chosen a series of decaying exponentials which displays the appropriate behaviour of initially fast decay followed by slower decay  L(t) = a() aieb't a2eNt a3eb3t a4eNt a5e where  5  t^  (3.37)  ai and bi can be arbitrarily chosen. The mass of water Mf which has frozen to ice  is then the volume of ice produced multiplied by the density of ice. The volume of newly frozen ice V1 is  dL  Mt) = irrl(t)— dt. dt Because rb changes much more slowly than  (3.38)  L, rb(t) may be replaced by rb,o, a constant  V1(t) = irr [L(0) — L(t)1^ (3.39) thus mf = rrt,o [L(0) —  L(t)]^  (3.40)  3.3.2 Sudden Pressurization The sudden-pressurization forcing is included to allow modelling of sudden pressurization and decay events observed in Trapridge Glacier pressure records. In principle a sudden pressure change could be modelled by appropriate choice of initial conditions. The borehole pressure would be something other than  Po, and stresses and strains in the ice would  be those corresponding to the analytic elastic solution. The difficulty with this approach is that, in general, the analytic elastic stresses and strains will be inconsistent with the system 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 numerical integration schemes fail when started with inconsistent initial conditions. An alternate approach which avoids this difficulty is to start with normal equilibrium initial  48  Chapter 3. Bed Equations, Forcings, and Model Synthesis^  conditions, 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 is Po + pforce) and tramp is ramp-up time (half the cosine period).  3.3.3 Pressure Step Change It is desirable to be able to simulate a creep test, in which the borehole pressure is suddenly changed, then remains constant. This allows comparison with analytic solutions (even though the infinitesimal-strain assumption of this model will cause inaccurate results at large strains). Creep tests are simulated by using (3.41) to ramp up the pressure, and then holding borehole pressure constant pb(t) = Po pforce for t > tromp.  ^  (3.42)  3.4 Summary of Coupled Initial Value Problems I now have all necessary equations to describe the coupled initial value problems: ice deformation, groundwater flow in the bed, and coupling through the water mass balance in the borehole. The initial value problem in the ice is specified by (2.152)-(2.189), with initial conditions (2.190)-(2.193) and boundary conditions (2.197)-(2.199). The initial value problem in the bed is given by (3.13), with initial condition (3.14) and boundary conditions (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^  49  3.5 Method of Lines Solution In the preceding sections the equations of the coupled ice and bed initial value problems have been derived. The ice and borehole equations are non-linear, and there are both partial differential and algebraic equations to be simultaneously solved. There is no hope of an analytical solution, so numerical methods must be turned to. Few numerical methods can handle such a problem. The finite element method has been used to solve problems involving the ice rheology of Shya.m and Wu, at least in one dimension [Shyam  Sunder 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 algebraic constraints (if an appropriate numerical integration routine is selected), and coupled initial value problems also present no special difficulty. I have chosen to use the method of lines, primarily as a result of my familiarity with the method and its comparative simplicity. Briefly, the method of lines involves discretizing the region into a grid, replacing continuous functions of the spatial variables with discrete functions at the grid nodes, replacing spatial derivatives in these equations with finite-difference approximations, and then solving the resulting system of ordinary differential equations, and possibly algebraic equations 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 regions replaced by systems of ordinary differential equations. The equations that apply at the interior and boundary nodes will be stated. Finally the resulting system of equations for the coupled problem will be presented in a form suitable for programming.  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  50  3.5.1 Spatial Coordinate Transformation  The method of lines requires that the ice and bed regions be discretized, so that finitedifference 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 small radius, close to the borehole. It is readily apparent that at infinite distance from the hole, all these variables will become spatially constant. Intuitively one would expect a more accurate approximation to the true result if the nodes are more numerous where variables change rapidly, namely at small radius. Some possible ways to achieve a closer node spacing at small radius include: (1) increasing the total number of nodes without altering their distribution, (2) applying a spatial coordinate transformation and spacing nodes 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 increases the 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 differing  time scales, requiring the integrator to use very small time steps to achieve reasonable accuracy.) Hence the problem may become too stiff to be solved. The adaptive grid method is very good for problems in which the solution shape is not known a priori, or where the solution shape changes with time as in the present problem. Unfortunately it has the disadvantage of being complex to program. As  a result of the drawbacks of the  other two alternatives given above, I chose to transform the spatial coordinate r, in the ice and bed. One approach would be to transform the spatial variable such that the field variables  Chapter 3. Bed Equations, Forcings, and Model Synthesis^  51  would change by comparable amounts between each pair of adjacent nodes. In this case the resulting solutions for the field variables would approximate straight lines in the transformed space. Although I do not in general know what the ice and bed solution shapes 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 proportional to lir (see [e.g., Carslaw and Jaeger, 1959, p. 247] for the solution to the equivalent heat conduction problem). 2. For steady-state or transient flow of linear (elastic, viscous, or viscoelastic) ice surrounding 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, stress is 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 be allowed: R = ?X where x is a finite nonzero real number, and R = In r. The logarithmic transformation is included because the steady-state solution to the diffusion equation with cylindrical symmetry is logarithmic [e.g., Carslaw and Jaeger, 1959, p. 189], and one 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, or  €9 +^— Er = 0 and force balance (2.186) r  our  or  ^  —  179  = .  Chapter 3. Bed Equations, Forcings, and Model Synthesis^  52  The initial and boundary conditions are dealt with by direct substitution. The only term in the compatibility of strains and force balance equations that needs to be transformed is ra/ar. The power-law transformed expression is derived as follows:  R = rx so that r = R-i.^  (3.43)  Hence  a^aR  T— = ar^r Or  aR a (r^ x) aR a  7*-5-;  rxrx-i  xrxa OR =xR a  a OR  aR  (3.44)  The equivalent logarithmically transformed expression is derived as follows: R = ln r so that r = eR.  (3.45)  Hence  a^aR a r Or OR  =  a^a r-ar (ln r) aR r aR a 0R  (3.46)  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  53  A convenient way to express the two transformation forms is to let  a , '  a  (3.47)  where W is a weighting function given by  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 the initial and boundary conditions. Hence compatibility of strains (2.185) becomes fe — — = 0 WOR  (3.50)  and force balance (2.186) becomes acr  aR  + cr — cre = O.  (3.51)  The R = 1/r transformation of the bed problem is similar In the bed equations r appears only in initial condition (3.14), boundary conditions (3.15) and (3.16), and the field equation (3.13)  Oh _^10 ( ,ah) at pwg(cx+ ni3) r2 Or e. Or  The initial and boundary conditions are dealt with by direct substitution. In the field equation the only term to transform is  1  a  ( 20h)  r2 -87r^Or First note that  •  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  54  and  OhOh OR Or = OR Or _2 0h TOR'  Therefore  1 a ( oh) -Or)  1 0  (3.52)  r 21 _2ahm  7 ;ir -  8 ( Oh\ 1 aR a ph) r2 8r OR 1:1) 1 _2) 02h -- —r  r2 OR2 _4 02h r 0112 4 02h R 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. The field equation (3.13) for head becomes Oh _^-02h ^ at — pwg (a + nfl)R4 OR2.  (3.54)  The borehole equations also require transformation. Equation (3.35) for the flow rate of water into the bed must be changed as it includes a spatial derivative term. Applying the R = 1/r transformation and using (3.52) to convert (3.35) yields  Q. =  —271-po exP (pb(t)  = —27rpo exp [fi (pb(t)  =  —  —  ah) aR ah  po)] r:K^ Po)ir:K  ah  2rpo exP (Pb(t) — po)) K oR  aR  rc  rc  (3.55)  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  R =Rb  55  R Rma),.  ^BOREHOLE rARH^ WALL^  N  1-1^2^M-1^M Fig. 3.2: Ice finite-difference grid.  3.5.2 Finite Differencing The next step in the method-of-lines formulation is to establish grids in the ice and bed, on which the partial differential equations in It and t can be replaced by systems of ordinary differential equations in t only. The differentials with respect to R will be replaced by finite-difference approximations. First consider the ice. The innermost node is at the It value corresponding to the borehole radius rb. I will let the outermost node be at Rmaz, the R value corresponding to some TmaI, to be be chosen later. This approach preserves flexibility, and retains the ability to model the response of thick-walled cylinders. The remaining nodes will be spaced 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, there are M nodes, identified by subscript i, with i = 1 at R Rb and i = M at R Rmaz. At all 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 normalstress condition (2.197) applies, unless a pressure forcing is applied in which case or  —Pb^  (3.56)  Chapter 3. Bed Equations, Forcings, and Model Synthesis^  Ra= Rc CAVITY WALL  rARal  56  Ra=0  N^•^I 2^Ma-1 Ma  Fig. 3.3: Bed finite-difference grid. where Pb is given by (3.41) for sudden pressurization, or by (2.5) for constant borehole  pressure. 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 node is Ra = 0, corresponding to 7.. = co (Fig. 3.3). Henceforth the subscript, a will be used to identify quantities pertaining to the basal grid, as opposed to the ice grid. As shown in Figure 3.3 there are M. nodes, identified by subscript j, with j = 1 at R. = 1/re and  j = 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 ice and bed equations listed above, as well as the water mass flow expression (3.16) which appears in the borehole equation. Examination of these equations reveals that finitedifference approximations are required for the first derivative at all ice nodes (innermost, interior, and outermost) and at the innermost basal node. Finite-difference approximations of the second derivative are required at bed interior nodes. There are numerous possible choices of finite-difference operators; I use second-order (3 point) operators for  Chapter 3. Bed Equations, Forcings, and Model Synthesis^  57  both first and second derivatives. Although accuracy could be improved by taking higherorder approximations, the stability of the resulting method of lines solution would be reduced [Schiesser, p. 110]. A general expression for the finite-difference approximations of first and higher derivatives, along with tabulated coefficients, is given by Abramowitz and Stegun [1972, p. 914]  die f (x)  k!  dxk  m!hk  E Aif('i)  (3.57)  where k is the order of differentiation, m is the order of the finite-difference approximation, h is the step size (or node spacing), j is the node at which the derivative is being approximated, and Ai are the (tabulated) coefficients. From these, I obtain the following expressions for second-order finite differences: 1. First derivative, centred difference (applies at interior nodes) df  -  dx  2Ax  (3.58)  2. First derivative, forward difference (applies at inner boundary node) df dx  -fi+2 + 4fi+i 2 Ax  - 3fi  (3.59)  3. First derivative, backward difference (applies at outer boundary node) df  3f= --- 4A-1 -I- A-2  dx  2 Ax  (3.60)  4. Second derivative, central difference (applies at interior nodes). (Second derivatives are never taken at boundary nodes.) d2 f dx2  f=1-1. (Ax)2  (3.61)  This completes the method-of-lines formulation of the coupled ice-bed-borehole problem. The resulting system of equations is summarized in Appendix A.  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  58  3.6 Computer Programming of Solution  The final step in the problem solution is the solution of the system of equations of Appendix A, using a numerical integration routine. I have chosen ddriv3 of Kahaner and Sutherland (public domain software, available from Los Alamos National Laboratory),  an implicit integrator which uses the Backward Differentiation Formula method, an extension of the Euler method [Schiesser, 1991, p. 65]. ddriv3 is similar to the stiffly stable integrator sdriv2 described in Kahaner et al. [1989, pp. 295-297, 341-346], except that it is double precision and capable of handling systems of differential and algebraic equations. The solution has been implemented in the program bfi.I.4 (Borehole Flow Infinitesimal, 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 themselves when the background pressure Po 0 0. The tests of Chapter 4 are valid because all the tests 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 followed  by free relaxation, sudden pressurization followed by freezing forcing, or freezing forcing 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, and  transient viscous parameters. 4. Ice temperature at each node.  Chapter 3. Bed Equations, Forcings, and Model Synthesis ^  59  5. Basal hydraulic properties including hydraulic conductivity and compressibility + 743.  6. Water compressibility. 7. Grid parameters including number of ice and bed nodes, and spatial transformation in 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 ice stresses, basal hydraulic heads, and debug output.  11. ddriv3 control parameters Parameters 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 transient 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 be achieved by suitable choice of viscous parameters (very stiff). The ice-borehole system can be modelled independently of the basal system by setting hydraulic conductivity and compressibility a ± nig to zero. The bed-borehole system (i.e., ice assumed rigid) can be modelled by suitable choice of ice parameters (very stiff). The outputs of the program are: 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 all ice nodes at all time steps ( str file), and the head at all aquifer nodes at all time steps ( .hed file).  Chapter 4  Model Testing  4.1 Introduction Any computer model must be tested to ensure that its results are correct. This is done by comparing model output to available known solutions (usually analytic), and by qualitative assessment when no know solutions exist. The performance of a methodof-lines solution is limited by both the accuracy and stability of the solution. In general accuracy is increased by decreasing the grid spacing (increasing the number of nodes), increasing the finite-difference approximation order, increasing the integration-method order, and decreasing the temporal step size [Schiesser, 1991, pp. 19-31]. In this problem the finite-difference approximation order has already been fixed, so that only the grid spacing, 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 the estimated error at each time step within a user-specified tolerance. Hence in this case the error tolerance replaces the temporal step size as the user-specified parameter. Stiff problems may also require the integrator to take very small time steps to meet the error tolerance. Stability can also be a problem, causing the model solution to diverge from the true solution, or causing*temporal numerical oscillation [Schiesser, 1991, pp. 122-129]. The model testing reported on here has two purposes: to determine the control parameter values which give optimum stability and accuracy, and to determine if the optimum solutions are acceptably close to known solutions. These questions are first addressed  60  Chapter 4. Model Testing^  61  for the ice part of the model, then for the bed part, and last, in a qualitative way, for the coupled model. The ice and bed models are tested and results compared to analytic solutions. It is shown that although numerical oscillation of the ice solution can be a problem and stiffness limits accuracy for certain choices of ice properties, acceptable solutions can in general be obtained with the correct choice of input parameters.  4.2 Ice Deformation Model The ice model was tested by running creep tests and systematically varying input parameters. These parameters include: ddriv3 error control parameters, spatial coordinate transformation, ratio of maximum radius to borehole radius rmazIrb, and number of nodes M. The effect on the Glen-Law viscous and linear elastic solutions was investigated. The  values of ice parameters used for all tests in this chapter are listed in Table 4.1. 4.2.1 Glen Law Viscous Ice -  These tests were run by modelling one-year creep tests, with transient and elastic parts of the rheology disabled. First I ran a number of tests to determine the optimum ddriv3 error control parameters. The ddriv3 integrator can be run using several different onestep error control strategies: absolute error control, relative error control, or relative error control with user-specified definitions of "problem zero". The last case is identical to relative error control, except that any dependent variable having a value less than its problem zero is considered equal to zero when its relative error is computed. This has 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 accurately computing numbers which are essentially zero. The integrator allows the problem zeros to be set separately for each variable. I have taken advantage of this capability because the  Chapter 4. Model Testing ^  Parameter  Value  Vo  6590 Pa s113  A  0.0142  H  0.02  0.286  Source  Shyam Sunder and Wu [1989a] estimated from data of Sinha [1978] Shyam _Sunder and Wu [1989a] fitted from data ofJacka [1984] Shyam Sunder and Wu [1989a] fitted from data of Jacka [1984] Shyam Sunder and Wu [1989a] fitted from data of Jacka [1984] Shyam Sunder and Wu [1989a]  Rga,  8.314 J mo1-1K-1  QL  67000J mori  Paterson [1981, p. 34]  QH  Paterson [1981, P. 34]  it  139000J mol-1 3.3005 x 109Pa  A  6.3608 x 109Pa  ice temperature  0 °C  water temperature  0 °C  Pw  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]  calculated from seismic velocities of Paterson [1981, p. 78] calculated from seismic velocities of Paterson [1981, p. 78]  Table 4.1: Standard rheological parameters used for model testing.  62  63  Chapter 4. Model Testing^  Parameter  bf i 13 Name  Maximum integration order Error control method specifier Relative error tolerance Problem zero for m,,,q, B, R, and a Problem zero for it Problem zero for et'  mxord ierror eps _ ewt (1 ,3 6 ,10 ,11) ewt (2) ewt (7 9) -  -  Value  5 4 1 x 10 1 x 10 1 x 10 1 x 10  2 6 10 40  -  -  -  -  Table 4.2: Optimum ddriv3 error control parameters. dependent variables have widely differing characteristic magnitudes (for instance stresses have typical magnitude 100-104, whereas viscous strains have typical magnitude i010-s.) It was found that the error tolerance (eps) and problem zeros (ewts) were limited by the step size required to keep the error within the tolerance. ddriv3 terminates with an error message whenever the step size has to be reduced more than 50 times in the attempt to take one time step, or if the time step size falls to less than the roundoff error in time. All step size errors observed during the model testing occurred during the pressure ramp-up. It was found that the solution changed by less than one part in 107 when eps was changed 1 x 10' to 1 x 10. The strictest parameters that could be run under most conditions are summarized in Table 4.2. It was found that the Glen-Law ice deformation problem is very stiff; many combinations of parameters resulted in problems so stiff that a step size error resulted, even with the loose error tolerance of eps = 1 x In contrast, when elasticity was included all trial runs with eps = 1 x 10-5 succeeded. The ratio of maximum node radius to borehole radius (r,„,,z/rb) has a large effect on how well the model solution approximates the true solution for infinite ice. This is  ^2  Chapter 4. Model Testing^  64  because bfil3 really solves the problem of flow in a thick-walled cylinder, which has a solution different to that of infinite ice, especially for small (rmazIrb). Hence even if the numerical solution were perfect it could not be expected to do any better than the analytic solution to the thick walled cylinder problem. In this section I focus on tangential strain at the borehole wall (comb) as a measure of solution accuracy. It is through tangential strain (and radial stress) at the borehole wall that the ice couples to the water. An analytic solution for steady-state power-law creep in a thick-walled cylinder is given by Finnie and Heller [1959, pp. 182-187] dee =  A (3  )3  \ (N+1)/2  Pb ^2^(rmaz  dt^4)^(r,,,azirb)2/N  —1  N^r  (4.1)  )2  where the flow law of ice is  in uniaxia1 compression. For the case of N = 3 considered here (4.1) becomes  dee^A  ) 3 ( 3 (^1^rmaz NI 2  (4.2)  dt = -6Ps (r„„„Irb)211s1 - 1^r ) where de^ = Ace. .1  (4. 3)  dt  This 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 to be expressed in terms of V, and yield the desired result. Recalling (2.65) and (2.72)  dt  A,Cri (c) dt 0  N-v 3 ( 1 N er^cr. V)  eq^dt )  (4.4)  Chapter 4. Model Testing  65  When stress is uniaxial (in x1 direction) (TOO  crij  (4.5)  000 000  and  lajj  is given by (2.45)  —  (2/3)cr  0  0  0  — (1/3)a  0  0  0  —(1/3)cr  (4.6)  With (2.51) for the definition of effective stress 2^2  (4.7)  6req = 6r Hence for uniaxial compression  des^3 (1)3 2,^(de dt^v a all Tit )0 de'^3 ( 1 \3 22 1de) dt^a i7Vi lt)0 = (1\3 3(d6\  (4.8)  Comparing (4.8) to (4.3) reveals that A = (del dt)o /V3. Hence (4.2) becomes ) (rmaz)2 deo^pl (de ( 1 dt^6V3 dt )0 k (r,./Tb)2/N — 1 )^r ) 3  (4.9)  and at the borehole wall  deo dt  3  74, ( de) ( ^1^(7-max ) 2 6V3^)^— 1)^rb ) •  (4.10)  For infinite outer radius  deo dt  rb  /d\ 6V3 dt ) 0  (4.11)  Chapter 4. Model Testing^  E .E 0  66  20 18 ANALYTIC SOLUTION  z 16  BFI13 SOLUTION  CC  1— (f) 14 0 -  12  _J •  L.L.1 0  0  0 1.0^2.0  3.0  ^  4.0  00  log(rmax/rb)  Fig. 4.1: Effect of maximum radius on borehole-wall tangential strain in Glen-Law viscous ice, after 1 yr creep test with Pb = 10 kPa, M = 9, and R 7.-213. The thick walled cylinder 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 thick walled cylinder solution by —1.6% at rmaz /rb = 10, —6.9% at 100, —10% at 1000, and —11% at 5000.  Chapter 4. Model Testing^  67  Model 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-ice one; the difference for rmazIrb = 1000 is 3.1%. The model results differ significantly from 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. No model output was obtained for rmazIrb > 104 as time step errors occurred. It appears that extending the grid to larger radius increases the stiffness of the problem. Since the thick walled cylinder solution overestimates the true strain, and the model solution underestimates the thick walled cylinder solution, there is an optimum number of nodes for which the model solution closely approximates the true one (at least in terms of the borehole-wall strain). This point is between 100 and 1000, for R = r-213. I therefore chose rmazIrb = 1000 for the balance of the testing. The choice of radial coordinate transformation has a large effect (Table 4.3, and Fig. 4.2). Transformations with negative powers of two or more could not be run without causing step size errors. Significant spatial oscillations, of wavelength 2 AR, are most noticable in the stress solution for R r and R In r. Similar oscillations exist in the other solutions, but are of smaller amplitude. Similar numerical oscillations can appear in solutions of the one-dimensional advection equation and the one-dimensional wave equation, 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 using un-centred finite-difference approximations. This may be true of the current problem as well. Nevertheless, this model generates adequate solutions for transformations on the order of R = r-2/3. It appears that, at least for rmaz/rb = 1000, that R = 7-31' is the optimum 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 expected the error drops as the number of nodes is increased. No step size errors were encountered  68  Chapter 4. Model Testing^  0  a.  6  a  4 2  (1)  0  (Y  —2  _J < < Ct  —4 —6 —8 —10  -  0  ^ ANALYTIC SOLUTION * BFI13 SOLUTION, R = r BFI13 SOLUTION, R = ln(r) 111111111  0^200^400^600^800  NON—DIMENSIONAL RADIUS 0  —2  —4 LU  —6  ANALYTIC SOLUTION * BFI13 SOLUTION, R = r 2/3 BFI13 SOLUTION, R =7.'175 BFI13 SOLUTION, R = -  _J  —8  cr —10 4^8^12^16  ^  20  NON—DIMENSIONAL RADIUS Fig. 4.2: Glen-Law viscous stress distributions for various radial coordinate transformations. M = 9 and r,„./rb = 1000. The outermost nodes are not plotted. Analytic solutions are those for infinite ice.  Chapter 4. Model Testing^  69  Fig. 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 5 nodes, 7.3% for 9, and 2.1% for 17.  Chapter 4. Model Testing^  Transformation (R .)  r ln r r-2/3 r 07 r 075 r r-1.0  Analytic solution (infinite ice)  70  fe,ri.  Relative error (%)  5.477253 x 10' 4.949487 x 10' 8.710883 x-10-6 8.938541 x 10-6 9.307882 x 10-6 9.709685 x 10' 1.161782 x 10  —99.99 —47 —7.3 —4.8 —0.96 3.3 24  9.398034 x 10'  —  Table 4.3: Dependence of borehole-wall tangential strain in Glen-Law ice on radial coordinate transformation. with R 7.213. However with R 7.-314 it was found that increasing M> 13 caused step size errors. Thus increasing the number of nodes increases the problem stiffness in agreement with the observations of Schiesser [1991, p. 22]. After consideration of the Glen-Law ice modelling, I chose as optimum the parameter values listed in Table 4.4. The value of M can be chosen as large as possible, subject to the stiffness constraint discussed earlier. These parameters lead to a Glen-Law solution accurate to within 1%. It is apparent from the results presented in this subsection that the errors caused by finite outer radius and the finite-difference approximation are far greater than those introduced by the numerical integration, for transformations in the range of r-3/4 < R < r-213.  Chapter 4. Model Testing^  71  Parameter  bf i 13 Name  Maximum integration order Error control method specifier Relative error tolerance Problem zero for maq, B, R, and a Problem zero for It Problem zero for ev Ratio of rmaz to rb Number of nodes  mxord i error eps ewt (1 3-6 ,10,11) ewt (2) ewt (7-9) r_max/r_b M  Value  ,  5 4 1 x 10-2 1 x 10-6 1 x 10' 1 x 10' 1000 9+  Table 4.4: Optimum parameters for Glen-Law ice. 4.2.2 Linear Elastic Ice The elastic behaviour of the model was tested by running it with power-law viscoelastic ice (i.e., transient part of rheology disabled), and observing the solution at the end of the 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 the elastic 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 size errors. No step size errors ever occurred in these tests, indicating that the viscoelastic problem 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 undergoing both internal and external pressure. For internal pressure only it reduces to  =  Pbrt2,7*  pbri2,r,n 2 ax  2 (A + p)(r,27,az - rb2)(rm2ax  -  r)r •  (4.12)  Chapter 4. Model Testing^  72  At the borehole wall u,.(rb) =  2 Pbrl ^Pbrbrmaz  2 (A + p) (rm 2 22 — 71) 2p, (rm 2 ax — Pbrb  ^I  2.z). r^+r„.4  2(r2 —^+ ft)^tt  (4.13)  Using the definition of tangential strain (2.188) this becomes ee  ,rb  Pb^r2  2)  2 (rLax — r)2 (A + it) Pb  ^( 14.  (rmaz/rb)2)  2 ((rmaz/rb)2 — 1) k(A A)  (4.14)  For infinite ice (4.14) reduces to 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, but not as much as in the Glen-Law case. Hence limiting r,,,./rb does not degrade the elastic solution as severely as it does the Glen-Law one. Thus the r„,„z/rb = 1000 chosen during the 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 quite accurate 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 quite acceptable. The effect of the number of nodes on solution accuracy was tested (Fig. 4.6). As expected, 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 stiffness limit on M was found; thus there is no reason not to use the maximum allowable number of nodes (M = 20), except if computation time is a concern.  Chapter 4. Model Testing^  73  1.535  1.530  0 ANALYTIC SOLUTION BFI13 SOLUTION  1.525  _1 1.520 1.1.1 z 1.515 F-  1.0  ^  2.0  ^  3.0  ^  4.0  ^  00  log(rmairb)  Fig. 4.4: Effect of maximum radius on borehole-wall tangential strain in linear elastic ice. M = 9 and R = 7.-2/3. The thick walled cylinder solution differs from the infinite solution by 1.4% at rmaz I rb = 10, 0.01% at 100, 0.0001% at 1000, and 0.000001% at 10000. The bf 113 solution differs from the thick walled cylinder solution by —0.79% at rmax/rb = 10, —0.002% at 100, 0.007% at 1000, and 0.007% at 104.  Chapter 4. Model Testing^ 8  _-  0  o  --IF  _  a  o  -.  ++  74  +  +  00  --e  0  0. 0 R -4- R  =r = ln(r)  IIIIIIIII  0^200^400 600 800  NON—DIMENSIONAL RADIUS —2  4^8^12^16^20  NON—DIMENSIONAL RADIUS  1.2^1.6^2.0^2.4^2.8  NON—DIMENSIONAL RADIUS Fig. 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. The outermost node is omitted for clarity.  Chapter 4. Model Testing^  --, E 1.5150  ANALYTIC SOLUTION 0  E 1.5148 0  1 0 1.5146 , .... 1.5144 z •TC 1.5142  75  BFI13 SOLUTIONS  1---  v) o 1.5140 rz 2  ___J 1.5136 < i= z 1.5134 LJJ (2) 1.5132 < H-  ^L......  4^6^8^10^12^14^16^18 20  00  NUMBER OF NODES  Fig. 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^  Transformation (R =---)  76  E8,rb  Relative error (%)  —4.893174 x 10-7  —132  r-3  1.188106 x 10-5 1.514491 x 10-6 1.514727 x 10-6 1.514925 x 10-6 1.514925 x 10-5 1.514923 x 10-5  —21 —0.03 —0.01 0.0002 0.0002 0.00007  Analytic solution (infinite ice)  1.514922 x 10-5  —  r ln r r-2/3 r-3/4 r-1 T-2  Table 4.5: Dependence of borehole-wall tangential strain in linear elastic ice on radial coordinate transformation. It appears, then, that the viscoelastic problem is better behaved than the Glen-Law one. It is less stiff. The elastic component of the solution is less affected by the choice of  rmax I rb  and suffers less finite-differencing error, compared to the Glen-Law viscous  solution. I find the optimum parameters determined for the Glen-Law problem to be acceptable for the elastic problem, except that eps = 1 x 10-5 be used, and that up to  M = 20 nodes be used. 4.3 Basal Groundwater Flow Model To test the accuracy of the basal part of the model, a step pressure change was simulated using bfil3; the number of nodes was varied and the results compared to the analytic solution. The analytic solution can be obtained, by analogy, from the solution to the problem of heat conduction in a whole space with imposed step temperature change at  Chapter 4. Model Testing^  77  a spherical inner boundary given by Carslaw and Jaeger [1959, p. 247]. The solution to  Oh^(1 a (20h , r re^ (4.16) = D^r -(w.)) h(r,O)  =  0, r > rc  (4.17)  h(co,t)  =  0, t > 0  (4.18)  h(re,t)  =  hforce = Pforcel PRIM  t>0  (4.19)  is  reit 4er  re) h(r, t) = ^' " erfc (2r^ -Vr 3t  (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 of rc, 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 parameters of 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 the time-dependent solution is well reproduced (Fig. 4.8). Note that the model solution at first diverges and later converges to the analytic one. No stability problems are evident for M < 20. It appears that 9 nodes provide a sufficiently accurate solution for most purposes, although more nodes would give greater accuracy, especially at early time.  4.4 Coupled Model The 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 borehole, and subsequent pressure decay, with Glen-Law ice. This is the only coupled problem  78  Chapter 4. Model Testing ^  1.0  .8 ANALYTIC SOLUTION 5 NODES 9 NODES^BFI13 SOLUTIONS 17 NODES  .2  .0 1  ^  2^3^4^5^6^7  ^ ^ 9 8  EQUIVALENT NODE  Fig. 4.7: A comparison of model predictions with analytic solution for bed response to a borehole pressure step change for various numbers of nodes at t = 12h. The equivalent node is that which has the same rlrb as for the 9-node case.  79  Chapter 4. Model Testing^  1.0 .8  E  2 .4 = .2 .0 1  ^  2^3^4^5^6 NODE  ^ ^ ^ 9 8 7  Fig. 4.8: Comparison of model predictions and analytic solution for bed response to borehole pressure step change at various times. An essentially steady state has been reached after 100yr).  Chapter 4. Model Testing^  80  for which an analytic solution was found. Next the effect of adding the elastic and transient viscous parts of the rheology was checked by simulating sealed-borehole sudden pressurization/decays. The results can be explained qualitatively. Lastly the coupled model was checked qualitatively by comparing sudden pressurization/decays for various combinations of ice and bed properties.  4.4.1 Sealed Borehole with Glen-Law Ice The coupled ice—water model can be checked by simulating the sudden pressurization of a sealed borehole and its subsequent decay for Glen-Law ice, and comparing results to the analytic solution. The analytic solution is derived from the strain-rate solution of Section 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 borehole with no freezing, the mass of water is constant and  dmw dt =0.  (4.21)  mw(t) = p(t) dV(t)  (4.22)  At any given time the mass of water is  where  dV (t) = 7r7.(t) dz^  (4.23)  = Po exP [0 (Pb(t) —^  (4.24)  and Pw(t)  Substituting (4.23) and (4.24) into (4.22), rnw (t)  Po exp [f3 (pb(t) — po)J Irq,(t) dz  irpor(t)L exp fi (pb(t) —^.  (4.25)  Chapter 4. Model Testing^  81  Differentiating (4.25) with respect to t dmw dt  = 1rPo(t)L 0(43 exp [ig (pb(t) — RI)] ddPit' -1-2rb(t)7 drib 1 exp  (4.26)  v3 (pb(t) — N)]}  rPo(t)Lrb(t)exp  1 b 2drb) [s (pb(t) — po)] (orb(t)icdt^dt )•  (4.27)  Conservation of mass demands that this be equal to zero, thus dpb^drb 0 = igrb(t)— 4- 2 dt^dt *  (4.28)  Now recall the analytic solution for the strain rate of Glen-Law ice around a cylindrical hole in infinite ice (4.11) ^p(t) 'de) dt^6V3 dt )  deo  (4.29)  rb  Using the definition of tangential strain (2.188), drb^r5(t)p(t) (de) dt^6V3^kdt^) 0*  (4.30)  Substituting (4.30) into (4.28) and rearranging, 0 0  clrb(t)pg(t) (dc) 1@rbtIdpb dt^6V3^dt ) 0  p(t) ( de) dPb dt^3V3 kdt ) 0 dpb^p(t) (de" dt^3i3V3 dt )^ 0 — p1,3(t)dpb =  ^ 1( de) 3f3V3^ 0  dt^  (4.31) (4.32)  Integrating (4.32) from 0 to t (assuming instantaneous pressurization at t^0) with Pb(0) — pforce, and simplifying, gives f Pb(t) _p_3 dp^1 (dc) ft Pforce^  0 3/3 V3 dt ) 0^  dr^ (4.33)  Chapter 4.^Model Testing  82  P  -2  2  pb(t)  t^(d\ 3)3V3 kdt ) 0  Pforce  p2(t)^Pior2 „  (4.34)  2^2  t^(de) 30V3 kdt ) 0  (4.35)  TOO^Pf—or2 ce  2t^( de\ 33y3 kdt ) 0  (4.36)  2t^( de\  131,2(0= PiLe 2 ce pb(t)^(601.  (4.37)  3,QV3 kdt 2t^(d€\ \3/3V3 ( dt)0)  1/2  (4.38)  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 results very closely approximate the analytic solution (Fig. 4.9). For 9 nodes the error in pressure after 1 yr is 0.4%, and for 13 nodes -0.04%.  4.4.2 Effect of Viscoelastic and Transient Viscoelastic Ice So far only the steady-state viscous and pure elastic behaviour of the ice model have been checked. The behaviour of viscoelastic and transient viscoelastic ice must be assessed qualitatively, as no analytic solutions exist. Three sealed borehole sudden pressurization/decay simulations were run with identical parameters: one with steady-state viscous ice, one with power-law viscoelastic ice, and one with transient viscoelastic ice. Attempts were also made to run the same test with transient viscous ice (elasticity turned off), but step 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 not consider 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 than the steady-state viscous one.  Chapter 4. Model Testing^  83  10 a_  9  w 8 t2  cc)  7  1.1.1  a. 6 LiJ  o5 LU  0 cri  4 3 0^50^100^150^200^250  ^  300  ^  350  TIME (DAYS)  Fig. 4.9: Sudden pressurization and relaxation of a sealed borehole surrounded by Glen-Law ice.  Chapter 4. Model Testing^  84  GLEN LAW VISCOUS — POWER LAW VISCOELASTIC ----- TRANSIENT VISCO ELASTIC  -----------------  0^50^100^150^200^250  ^  300  ^  350  TIME (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-Law viscous one (Fig. 4.10). This appears counter-intuitive at first. One might expect that as viscous creep causes pressure to drop, the ice would elastically relax, reducing borehole radius; hence pressure would remain high for longer than would happen for Glen-Law ice. This is in fact what happens for the flow solution with linear viscoelastic ice around a borehole and for power-law viscoelastic ice in uniaxial compression. The important difference between those two cases and the thesis problem is that they involve uniform stress distributions; the ratio cr,./cro/az is the same everywhere, as is ce/ev. Hence in these cases the total strain solutions can be obtained by linear superposition of the elastic and viscous solutions. But superposition fails for the deformation of power-law viscoelastic ice around a borehole. The stress solutions for linear elastic ice and power-law viscous ice  Chapter 4. Model Testing^  85  are different; hence total strain is not equal to the sum of the strains found from separate elastic and viscous solutions. For a given loading the stress and strain state depends on the loading history. The initial state of stress is purely elastic, evolving later toward the viscous stress distribution [Murat et al., 19861. Hence intuition based on superposition of elastic and viscous solutions fails, and cannot be used to predict the relative pressure decay rates in the borehole. Thus the observed result is not ruled out, although I cannot explain it conceptually. The accelerated pressure decay with transient viscoelastic ice compared to that with power-law viscoelastic ice is consistent with the fact that transient creep is faster than steady-state creep.  4.4.3 Ice—Water—Bed Model This test was done to verify that the fully coupled model behaves as expected, in a qualitative sense. Sudden pressurization and decay simulations were run for the following four cases: deformable ice and impermeable bed, deformable ice and permeable bed (both with transient viscoelastic ice theology), Glen-Law deformable ice and permeable bed, and rigid ice and permeable bed. As expected, pressure decays much more quickly when the bed is permeable than when it is impermeable (Fig. 4.11). The situation is not so clear-cut for the case of deformable ice. Glen-Law ice results in a very slightly faster decay 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 relax quickly, relieving the pressure on the water. Later, as transient creep slows down, ice elasticity and recovery of transient creep dominate, and slow pressure decay.  Chapter 4. Model Testing^  (a) TRANSIENT VISCOELASTIC ICE, IMPERMEABLE BED' (b) TRANSIENT VISCOELASTIC ICE, PERMEABLE BED (c) GLEN LAW ICE, PERMEABLE BED d) RIGID ICE, PERMEABLE BED  10  0 a_ 8 uJ  (/) (1)6 uJ  a_ LU  0  4  LU  02  86  \ _  \-^ \i x \ -^\ -. b N --,  _  c,d  a  ---, „  .0^.2^.4^.6  , ^.8^1 .0  TIME (DAYS)  Fig. 4.11: Pressure decay for various ice/bed combinations. The permeable bed has K = 1 x 10-12 ms-1 and D = 1 x 10-10 m2 s-1  Chapter 4. Model Testing^  87  Parameter  bfii3 Name  Maximum integration order Error control method specifier Relative error tolerance Problem zero for maq, B, R, and a Problem zero for h. Problem zero for ev Ratio of 7.,„„^to rb Number of nodes (ice) Number of nodes (bed)  mxord ierror eps ewt (1,3 6,10,11) ewt (2) ewt (7 9) r_max/r_b -  -  Value 5 4 1 x 10 lx 10 1 x 10 1 x 10  M M_a  2 6 6 40  -  -  -  -  1000 9+ 20  Table 4.6: Optimum parameters for coupled problem. 4.5 Conclusions It is concluded that the model adequately simulates the problem. The elastic solution can be reproduced to within 0.01%. The Glen-Law viscous solution can be reproduced to within 1%. The bed solution can be reproduced to within 5%. The behaviour of the coupled system is qualitatively correct. The optimum control parameters are summarized in Table 4.6.  Chapter 5  Model Results  5.1 Introduction  In this chapter the bfil4 computer model is used to infer Trapridge Glacier ice and bed properties. Of primary interest are the basal hydraulic properties of hydraulic conductivity and bed compressibility. These are determined by using the model to interpret the observations obtained from three boreholes drilled in the summer of 1992, as well as three boreholes from earlier years. In 1992 two boreholes were drilled 5 m apart on the same day. The first, 921124, was an unconnected hole. The second, 921126, was a blind hole. Both holes were instrumented with a pressure sensor and a thermistor string. Analysis 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 forward modelling using the model developed in Chapters 2 and 3. First the ice properties are inferred from forward modelling of the 921126 pressure record, then the bed properties are inferred from forward modelling of the 921124 pressure record using the ice properties inferred in the first step. The sensitivity of the resulting basal parameters to changes in some of the more uncertain input parameters (such as ice parameters and borehole freezing forcing) are tested by repeating the fitting with altered parameters. In addition  three other unconnected holes from other years (891150, 901154, and 911146) are modelled using the 1992 freezing forcing and ice parameters determined from 921126, and a second blind hole (921145) is modelled.  88  Chapter 5. Model Results ^  89  Finally, the model is used to answer the question of to what extent the pressure records of unconnected holes reflect the "true" basal conditions, and to what extent they are altered by the response of the borehole itself. This is done by examining the response of the model to sudden pressurization, with the ice and bed properties obtained for 921E24 and 921126. 5.2 Freezing Forcing L(t) Obtained from Thermistor Records  Boreholes through Trapridge Glacier are drilled with a hot-water drill. In all years except 1992 the holes were drilled to their final diameter in one pass with the drill in roughly one 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 first pass took about one hour, and the second about 10 minutes. In 1992 two thermistor strings were built and installed: 92T01 in unconnected hole 921124, and 92T02 in blind hole 921126. The thermistors were calibrated in 1990 by P. Langevin and zero readings checked in 1992 in an ice bath before cabling. All thermistors 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 is not critical, so this is acceptable. The thermistors were calibrated and strings assembled as described in Jarvis [1973], except that Thermometrics P10 DA 402L thermistors  were used. Measurements were taken using a Fluke 8060A digital multimeter, using the two-polarity method of Clarke and Blake [1991] to average out spurious electro-chemical potentials. Readings were taken once a day for nine days following installation. Lead  resistance was calculated and subtracted from the readings before the calibrations were applied. Trapridge Glacier is coldest near the surface^—5 °C at the surface if the winter  Chapter 5. Model Results^  90  cold wave is neglected), and below freezing throughout except at the bed where it is at the pressure melting point [Clarke and Blake, 1991]. I therefore expect that a borehole freezes shut first at the highest point with water in it, then downward from there. The time of the first freeze-over can be found from the borehole water-pressure record (It is the time when the pressure begins to rise, usually about 12h after drilling). Hole length then decreases with time, asymptotically approaching zero, as the freezing time at the bed where ice is at its melting point is essentially infinite (neglecting melting due to geothermal or frictional heating). The temperature versus time after cessation of drilling was plotted for the thermistors of the two strings (Figs. 5.1 and 5.2). For each thermistor a smooth curve was extrapolated back until it intersected the initial temperature (measured soon after installation when all thermistors were in water); this was taken to be the time when the thermistor froze 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 of time required by the model. The records from the two holes were then stacked together to reduce noise. I assume that because the holes are nearby the ice surrounding both holes has the same temperature stratification, and furthermore that because the drilling method was approximately the same, the two boreholes should freeze at the same rate. A point was added at time zero, corresponding to the initial water level in 921124. This was done under the assumption that freezing would happen first at the water surface. Also a point was added at long time (L = 0.5 m at t = 2 yr). Pressure sensors located 0.5m above the bed appear to 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 points at zero and two years (weighted 10 times and 20 times respectively to reflect greater confidence in these points, and force the solution away from local minima) were fitted  ^  ••■•■• ..1.0 NEN  \  THERMISTOR DEPTH (m)  —  N \ \  -x^48.1  \ \ \  \  ■\  \  \^ \ \^ \ \ \^ \^\ \ \^\^\ \ \^\^\  -x- -x  38.1  28.1  \^\\  \^\ \ \ \\ \^  -x^m-x^23.1  ■  -x- --x --*-x^18.1  13.1  HOLE FIRST FREEZES SHUT  ■■•  -x- -x-^-x^8.1 1^I^1  1^1^1^1  1^1^1^1^1^1^1^1  ^  100 ELAPSED TIME SINCE HOLE DRILLED (DAYS)  101  Fig. 5.1: 921124 freeze-in temperature record (92T01). Solid lines are smooth extrapolations used to find time of freeze-in. Thermistors at 58.1m and 68.1m displayed no temperature drop and are omitted for clarity.  0 tin  $k\-• — 1  THERMISTOR DEPTH (m)  OMB  \\  45.9  \  -x--K  \\ \\  35.9  \ \  \  \ N  \  \  \  25.9  \  ik  \ X- )4-  )4-  —6  7  ---  \ )4• +(-  15.9  10.9  HOLE FIRST FREEZES SHUT  I1-1 -9( -41-4r4i  100^ ELAPSED TIME SINCE HOLE DRILLED (DAYS)  5.9  101  Fig. 5.2: 92H26 freeze-in temperature record (92T02). Solid lines are smooth extrapolations used to find time of freeze-in. Thermistors at 55.9m and 65.9m displayed no temperature drop and are omitted for clarity.  Chapter 5. Model Results^  93  Thermistor String  Depth (m)  Borehole Length (m)  Time After Drilling (days)  Time After Start of Pressure Rise(days)  92T01  8.06 13.06 18.06 23.06 28.06 38.06 48.06 58.06 68.06  61.4 56.4 51.4 46.4 41.4 31.4 21.4 11.4 1.4  0.58 0.72 0.95 1.3 2.1 4.2 6.0 > 9.0 > 9.0  0.02 0.16 0.39 0.7 1.5 3.6 5.4 > 8.4 > 8.4  92T02  5.92 10.92 15.92 25.92 35.92 45.92 55.92 65.92  63.6 58.6 53.6 43.6 33.6 23.6 13.6 3.6  0.51 0.67 1.0 1.8 3.7 6.4 > 8.9 > 8.9  0.01 0.17 0.5 1.3 3.2 5.9 > 8.4 > 8.4  Table 5.1: Freeze-in times for 92T01 and 92T02. The borehole length is glacier depth at 921124 (69.48m) minus the thermistor depth below the glacier surface.  Chapter 5. Model Results^  94  80 60 40 20 - 1990  80 60 40 20  1991  ^  - 1992 ^ 300 0^100^200 DECIMAL DAY Fig. 5.3: Multi-year pressure record for 90P02. Note that after one year there is considerable activity, indicating the sensor remains in water, whereas after two years only slow changes occur, indicating that the sensor is frozen in and observes only slowly varying ice stress.  ^  Chapter 5. Model Results^  60  95  92T01(92H24) X 92T02 (92H26) * INITIAL WATER LEVEL IN 92H24  50  0 40 LLI _J  0 30 UJ  0 CO 20  10 0  ^^ 1 3^4^5^6^7 ^  •  8  ^ ^ 9 10  TIME (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) by minimizing 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 from increasing with time. The result is L(t) of the bottoming hole 921124, and the inverted function 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 same  as those for hole 921124, it is also the freezing time as a function of height above the bed in other nearby boreholes. Hence for these nearby boreholes the length as a function of time can be obtained from L(t) for 921126. To adapt the estimate to the blind hole 921126 two changes must be made: L must be made smaller by the distance between the bottom of the hole and the bed (by adding a negative ao), and L must be reduced so  Chapter 5. Model Results^  96  that Lt=0 is consistent with the initial water level in the borehole. Since freezing time versus depth is assumed to be the same, this is done by shifting the time zero (i.e., start later on the same curve). So the new as are set to the old ai exp(bi At), where At is the time at which L = Lt.0 in the new hole. This was done for 921126 and the other holes considered in this chapter.  5.3 Acceptable Parameter Ranges Before the observed pressure records are fit by adjusting input parameters, it is necessary to place reasonable limits on these parameters so that, for example, an ice temperature of +5 °C may not be used to obtain a good fit. Some parameter ranges were obtained from the literature (for example ice and bed properties), whereas others were obtained from Trapridge Glacier conditions (for example temperature and geometry).  5.3.1 Geometry The freezing forcing L(t) was obtained from the 1992 thermistor data as outlined in Section 5.2. The borehole radius is that of Trapridge Glacier boreholes. These holes are nominally 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. But it 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 surface the borehole will be smaller at the top. Suppose the borehole always has a parabolic vertical cross-section (Fig. 5.5). The equation of this parabola is  z = L (r/rb)2.^  (5.1)  Chapter 5. Model Results^  rb->.  Fig. 5.5: Parabolic borehole.  97  Chapter 5. Model Results^  98  The volume of water in this borehole is found by evaluating the volume integral  0  (L z) r dr de.  (5.2)  Substituting (5.1) for z  V  r  de [L frb r dr — (L/re2,)  1021r de Lr21: 102w  r  de f.')r L (1 — (r/rb)2) r dr  1022.^  de^—  21r 1  —  ±r4 4r  r3dri  rbl  Lrd  de  Lr?, I: de —Lr2 2 b  (5 .3)  Comparing this to the volume V = irrl of a cylindrical borehole of length L and radius  rb, 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 borehole radius, 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 borehole is poorly constrained. Scratches in the paint of a dummy ploughmeter inserted and recovered in the summer of 1992 indicate that the area of disturbed and weakened sediment (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 all this material is completely removed from the bed in an unconnected borehole in which the water (and sediment) has nowhere to go. It is likely stirred up and loosened, then  Chapter 5. Model Results^  99  allowed to settle back into place. Hence the cavity radius might be considerably smaller than the ploughmeter result suggests. The borehole radius will be taken to be the lower bound. Hence 0.025 < rc < 0.20m is assumed for bottoming holes. In blind holes there is no basal cavity, so a cavity radius of zero should be used. This would cause numerical difficulties in the model though, so a small radius of 1 x 10-5m will instead be used for blind holes. 5.3.2 Ice Properties The minimum elastic parameters found in the literature were those computed from the seismic 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 a range of 9-11 GPa. I will therefore use the range 8.77 < E < 11 GPa. Assuming that the 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. I  consider it unnecessary to allow these to vary, since they affect only V, which is also affected by the variables Vo and T. The viscous parameter Vo is widely discussed in the literature. Shyam Sunder and Wu [1990] give a range of values 5390 Pa sh/3 < Vo < 7860 Pa s1/3. These correspond to 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 Paterson [1981, p. 39] (102 x 106 Pa s1/3), Hooke [1981] (101 x 105Pas1/3), Budd and Jacka [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^  Fresh water columnar (S-2) d = 3 mm [Sinha, 1978]  Fresh water isotropic polycrystalline d = 1 7 mm [Jacka, 1984]  100  A  ft  Bo  0.33  0.454  0.0865  0.0142  0.02  0.286  Table 5.2: Transient parameters from Shyam Sunder and Wu [1990]; d is the average grain size of the ice tested. V(-1 °C) < 131 x 106 Pa s1/3. Inverting (2.24) and rounding the results gives the range  for 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 data reported by others in laboratory tests (Table 5.2). The transient behaviour of ice has been shown to be grain-size dependent [e.g., Sinha, 1979]. Both the tests quoted were done on relatively fine grained ice. Glacier ice can have much larger grains, anywhere from 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 grain sizes. If the grain-size dependence of Sinha [1979] is accepted, the transient parameters are related to the grain size by [Shyam Sunder and Wu, 1989a]  A^dl .41 dw+i)IN Bo = gadliN  (5.4) (5.5)  (5.6)  Chapter 5. Model Results^  101  A' (mm)  ill (mm-113)  .13L (mm-113)  120 9  0.010 0.105  0.24 0.06  From data of Sinha [1978] From data of Jacka [1984]  Table 5.3: Grain-size-independent transient parameters.  d (mm) 1 100  Sinha 0.1 10  A  il Jacka  Sinha  Jacka  0.008 0.8  0.1 50  0.01 5  Bo Sinha Jacka 0.06 0.3  0.2 1.0  Table 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-sizeindependent parameters (Table 5.3). Now assuming a minimum grain size of 1 mm and a maximum grain size of 100 mm, and applying (5.4)-(5.6), gives the grain-size-dependent transient parameters (Table 5.4). Taking the extreme values from Table 5.4 (and rounding) gives the ranges of the transient parameters as 0.01 < A < 10, 0.01 <H < 50, and 0.05 < Bo < 1. In the absence of more detailed information, ice temperature will be assumed constant in time and space. The representative temperature will be some sort of average with depth in the borehole and distance from the borehole. Ice cannot be warmer than the pressure melting point, approximately 0°C. The coldest reasonable effective average  Chapter 5. Model Results^  102  ice temperature must be somewhat warmer than the average background ice temperature, because ice surrounding the borehole is warmed by drilling, and deformation will concentrate 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 temperature of —2°C. Hence —2 < T < 0°C.  5.3.3 Bed Properties The range of hydraulic conductivity and matrix compressibility for glacial till deposits are [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). Consolidation tests 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 ranges 1.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. These and other parameter ranges are summarized in Table 5.5.  5.4 Ice Parameters from Blind Hole Fits Because blind holes do not reach the bed, the pressure signal observed in them must be a function of ice properties, geometry, and freezing rate only. Hence it should be possible to discover the ice properties by forward modelling and fitting the observed pressure signals, 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 which also contained the thermistor string 92T02. The freeze-in record of this sensor shows a very rapid pressure rise to high pressure, followed by a sudden drop, and further repeated steep rises and sudden drops until, after several days, the pressure eventually stabilizes and begins to decline (Fig. 5.6). The rapid rise results from the compression applied by  Chapter 5. Model Results  ^  103  a  N  = 180  E  •"'' 160 w c: D 0 140 o w rx 120 a_ r: w F_ 100  <  w 80 ___I 0 = w 60 cK 0 co  202^204^206  208  ^  DECIMAL DAY, 1992  Fig. 5.6: Freeze in pressure record for 92P02 (in borehole 921126). -  210  104  Chapter 5. Model Results^  Parameter  rb rc A A QL QH Vo  A ii. Bo T  Minimum  Maximum  0.014 0.03 0.025 0.20 3.3005 x 109 4.1 x 109 6.3608 x 10 8.0 x 109 67000 126000 5400 9700 0.01 10 0.01 50 1 0.05 —5 0  Units  m m Pa' Pa-1 J mori J mori Pa s1/3 —  —  °C  Table 5.5: Summary of input parameter ranges. the expansion of freezing water. The sudden drops presumably reflect cracking of the ice, allowing water to exit the borehole, thus lowering the pressure. It is likely that this cracking occurs primarily at the top of the borehole where the sharp curvature of the surface 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 rise to first fracture, will be modelled. This is necessary because the SS&W rheology has no fracture mechanics embedded in it. The computer model cannot be started at the end of the drilling, as would be ideal, since it is not flexible enough to deal with a constant pressure followed by freezing-induced compression. Closer examination of the early part of the 92P02 record reveals that the initial pressurization is not one smooth curve (Fig. 5.7). I am unsure what causes this inflection in the curve. All freeze-in pressure records examined for this thesis displayed this feature,  Chapter 5. Model Results ^  105  cL 1600 Lu 1400 (I) LLI  a_ 1-LI  1200 1000  I-L-1^800  0 0  600  CO  .0^.5^1.0^1.5^2.0^2.5  ^  3.0  TIME (HOURS)  Fig. 5.7: Initial part of 921126 freeze-in from 92P02).  ^  3.5  Chapter 5. Model Results  ^  106  180 9,4160 140 0 120 0 100 uJ  0  80  60  .0^.5^1.0^1.5^2.0^2.5  ^  3.0  ^  3.5  TIME (HOURS) Fig. 5.8: Freeze-in pressure records considered in this thesis: (a) borehole 921E45 (pressure sensor 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 character of a blind hole. Broken lines indicate interpolated missing data. with the sole exception of 89PO4 (Fig. 5.8). One possible explanation involves water entry into instrument cables. Pressurized water has been observed to flow up inside cables (between the outer jacket and the insulated conductors) and into data logger enclosure 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 further water entry. (The volue of water displaced by the assumed freezing before the inflection point is 0.120 m3, somewhat less than the volume of air in the two cables of approximately 0.400 m3.) This would cause the inflection point on the pressure records. After that  Chapter 5. Model Results^  107  the system would behave as a water-filled borehole, as modelled in this thesis. (It is not possible that the compressibility of cable materials, PVC and copper, plays a part since 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. To do so, I will extrapolate this curve back until it reaches the initial pressure (Fig. 5.7) and time shift the borehole L(t) (as was done in Section 5.2) so that the simulation begins at that point. This is yet another compromise, necessitated by the fact that the model does not include all the physics involved. The model ice parameters were varied until the best visual fit with the pressure record was obtained. It was found that the pressure rise was scaled up by decreasing T, or increasing 'Vo, it, or A. The rate of pressure rise was increased by increasing the transient parameters A, n, and Bo. Changing rb had no significant effect. The best fit, which is 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); and  A=  0.02, fi = 0.02, and B0 = 0.05 (rather close to the lowest allowable). The poor  fit suggests that the model is not adequately describing the system. One likely cause is the poorly constrained, approximate, top-down freezing forcing. Another is the use of a one-layer model to simulate deformation in ice with depth-varying temperature, and hence rheological properties. A third effect is that new, unstrained ice is continuously being added to the borehole wall by freezing, which could be expected to depress the apparent values of the transient parameters. Another possible cause of the poor fit is that the ice strains get large enough that the infinitesimal-strain assumption fails. However, at the time of first fracture, the modelpredicted 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 strain definition is unlikely to be the source of error.  Chapter 5. Model Results  ^  108  a_ 1600 1400 (I)  a_ 1200 1000  0 0  800 600  .0^.5^1.0^1.5^2.0^2.5  ^  3.0  ^  3.5  TIME (HOURS)  Fig. 5.9: Best fit to 921=126 freeze-in pressure record. Solid line is observed pressure record, dashed line is bfil4 model simulation.  Chapter 5. Model Results  ^  109  1800 uJ  1600 (f)  1400 w 1200 FuJ  0 0 co  1000 800 600  .0^.5^1.0^1.5^2.0^2.5  ^  3.0  ^  3.5  TIME (HOURS) Fig. 5.10: Best fit to 921145 freeze-in pressure record. Model ice parameters are  T = — 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 92P02 and 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 softer transient parameters. The 921126 record was also fit using a slower decaying L(t). This was done to investigate the effect of off-centre thermistors on the results. If the thermistors are off-centre they will freeze in sooner than if they were on the centre-line of the borehole, as I have assumed. Hence the actual borehole length would be longer than assumed. To model this possibility I halved all bis in the L(t) expression and repeated the fit of 921126. The  Chapter 5. Model Results  ^  110  1600 LU  D 1400 LU  a_ 1200 1000 _J 0 1.1.1  0  800 600 .0^.5^1.0^1.5^2.0^2.5  ^  3.0  ^  3.5  TIME (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); and  A = 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 Fits This investigation concentrates primarily on 921124, one of the two 1992 boreholes from which the length function L(t) was derived. It was also the only unconnected hole in 1992 for which a pressure record of the freeze-in was obtained. This hole was unconnected at the time of drilling, as indicated by the high and constant pressure in the hole between drilling and the freezing-induced pressure rise (Fig. 5.12). This hole did later connect  Chapter 5. Model Results ^  111  0c,i 95 2  E  90  Li 85 CC  D 80 o cn  1-1-1 75 0_ 70 (2 Li I— < 65  L.Li  _J  60  0 55  i Lu 0  CO  50 202^204^206  208  ^  DECIMAL DAY, 1992  Fig. 5.12: Freeze - in pressure record of 921124 (from 91P12).  210  Chapter 5. Model Results^  112  to the subglacial drainage, as indicated by the onset of large-amplitude diurnal signal halfway through day 201. However I believe that the hole was unconnected during the freeze-in period, until approximately 0900h on day 201. The method used to obtain bed-parameter estimates was similar to that used, in the previous section, to obtain ice-parameter estimates from blind holes. The freezing forcing was time shifted to account for the initial water level difference (none in the case of 921124, because it is the hole for which L(t) was initially computed), and the delay between the first pressure rise and the start of the extrapolated post-inflection freeze-in curve. Since the basal-cavity radius re is poorly known, fits were done by adjusting K and a + ni3 for several values of re spanning its range. It was found that the height of the pressure peak was decreased by increasing K, a + 743, or re. The shape of the pressure peak was affected differently by changing the different parameters. Increasing K caused a sharper peak with a more rapid subsequent pressure decline, whereas increasing a-En/3 or re had the opposite effect of causing a rounder peak with a slower decline. Hence for a given re the records were fit by trading off K and a + rt,8 to obtain the correct height and shape of curve. First 921124 was fit using the "standard" L(t) and ice properties determined from the fit of 921126 using the same L(t). Fits were obtained for re of 0.025m, 0.10m, and 0.13m (Fig. 5.13 and Table 5.6). A good fit could not be obtained for re > 0.13m without decreasing a ni3 below the minimum allowable 1.0 x 10-8Pa-1. This suggests that the basal cavity is not larger than 0.13m in radius, consistent with the reasoning of Section 5.3.1 that the material disturbed by the drill settles back into place.  113  Chapter 5. Model Results^  0  900 850  w 800  cr)  750  I-1-1 700  a_  0 a_ _y  a  650 900 850 800  cn  750  I-1-1 700 650 0  900 850  w 800 CC cn 750 (/) 700 0_ 650 0^1^2^3^4^5  ^  6  TIME (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^  114  rc (m)  K (ms-1)  a + nfl (Pa-1)  0.025 0.10 0.13  7.0 x 10-8 1.7 x 10-8 1.35 x 10-8  1.6 x 10-8 3.0 x 10-8 1.0 x 10-8  Table 5.6: Best-fit bed parameters from fit to 921124 freeze-in pressure record, using standard L(t). 5.5.1 Effect of Ice Properties An important question is how much influence the ice properties have on these results. It is necessary to know which parameters most strongly affect the results in order to assess the reliability of results obtained using uncertain input parameters. To do this the behaviour of the model was investigated with the softest allowable ice, the stiffest allowable ice, and rigid ice. The softest ice is that with all ice properties minimum, except temperature which was the maximum. The stiffest ice is that with all ice properties maximum, except temperature, which was —5 °C. Rigid ice was simulated by running bfil4 with the transient part of the rheology disabled, I/0 = 9700 x 1010 Pas", ji = 4.1 x 1030 Pa, and A = 8.0 x 1088 Pa (higher values resulted in step size errors). The best-fit bed properties for r, = 0.025 m from Table 5.6 were used. The ice does indeed influence the simulation results (Fig. 5.14), although not as strongly as K does. The pressure rise changes by many 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 rigid ice (see Figure 5.15 and Table 5.7). The best fit obtained with the softest ice is not very good. 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 K  •Chapter 5. Model Results  ^  115  'E..; 950  LLI  900  v) 850 Ui  a_ 800 < 750 0 700 UJ  0 650 CD  0  ^^ ^ ^ 2^3^4 1 5 6  TIME (HOURS) Fig. 5.14: Freeze-in pressure record of 921124, modelled using softest, stiffest, and rigid ice.  Softest Ice 921126 Best Fit Ice Rigid Ice  K (ms-1)  a +743 (Pa-1)  7.2 x 10-9 7.0 x 10-9 6.5 x 10-9  1.0 x 10-6 1.6 x 10-6 4.0 x 10-8  Table 5.7: Best-fit bed parameters from fit to 921124 freeze-in pressure record, with various ice properties. Standard L(t), rc = 0.025m.  116  Chapter 5. Model Results^  Lii  900  v) 850 (I) 1.1J  cc a_ 800 F-- 750  w  700  a  0 iii  650  0 CD  0  ^^ ^ ^ ^ 1 6 .2 5 3  0  2^3^4  a_ LU  900 850  Lii  800 F_ 750 w 700  0 w 650 0  5  6  TIME (HOURS) Fig. 5.15: Best fit to 921124 freeze in pressure record with (a) softest ice and (b) rigid ice, for rc = 0.025m. -  Chapter 5. Model Results^  117  and a slightly lower a + nO. Nevertheless the softest ice values in Table 5.7 are probably within 20% of the best-fit values. It is evident that the hydraulic conductivity estimate is quite insensitive to the assumed ice properties, but the compressibility estimate is more  sensitive. Also, it appears that assumed ice properties have less effect on the hydraulic conductivity estimate than does the assumed basal-cavity radius (Tables 5.6, 5.7). This is explained by the fact that for bed properties in this range the volume of water expelled into the bed is much greater than the increase in volume of the borehole due to ice deformation. 5.5.2 Effect of Freezing Forcing  Perhaps the largest uncertainty surrounds the freezing forcing. To assess the effect of the freezing forcing on the inferred bed properties, 921126 was separately fit with a weaker freezing forcing than the standard L(t) and a stronger one. The weaker forcing was the (appropriately time shifted) one used in Section 5.4, except that all bis in the L(t) expression were doubled (hence the borehole length decreases more slowly). 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 this freezing 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 must  be 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, summarized in 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 it appears that the hydraulic conductivity estimate is not overly sensitive to perturbations to the freezing forcing, the fact that freezing is so crudely modelled suggests that this  Chapter 5. Model Results  ^  118  -c:"i a. _Y 900 ..-..." LL.1  Ct D 850 cn (I) Li  f: a_ 800 cc w F 750 < —  w _1 700 0 = Li (: 650 0 m 0  ^ ^^ 5 2^3^4 1  TIME (HOURS) Fig. 5.16: Best fit to 921124 freeze-in pressure record, with weaker freezing forcing. rc = 0.025m.  Weak Forcing Standard Forcing Stronger Forcing  K (m s-1)  a + nf3 (Pa-1)  4.5 x 10 7.0 x 10' 1.3 x 10-8  3.5 x 10-7 1.6 x 10-6 4.0 x 10-6  Table 5.8: Best-fit bed parameters from fit to 921324 freeze-in pressure record, with weak, strong, and standard freezing forcings. re = 0.025m.  119  Chapter 5. Model Results^  .-. 900 0 CL __Nc 850 ..._..., w 800 CC op 750 U)  CC a_  700 650  -- ---------- _  _ _ I-Li 800 OC D (f) 750 U) 1-LI 700 _ CC  ,—, 900 0 a_ -- 850 .......-  a_  650  b 1  ;  ;  1  1  ;  , 900 0  a_  __y 850  Li 800 r: D 750 (i)  U)  t_Li 700 a_  650 0^1^2^3^4^5  ^  6  TIME (HOURS)  Fig. 5.17: Best fits to 92H24 freeze-in pressure record with rb = 0.025m and various  values of rc: (a) 0.025m, (b) 0.10m, and (c) 0.20m.  Chapter 5. Model Results^  remains an  120  area of large uncertainty.  5.5.3 Other Boreholes  To 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 examined: 891150 (pressure sensor 89PO4), 901154 (90P10), and 911146 (91P05). These holes did not connect when drilled. 901124 and 911146 both have freeze-in pressure curves with characters similar to the unconnected 921124, while 891150 looks suspiciously like a blind hole (Fig. 5.8). Nevertheless all three were modelled as bottoming unconnected holes, with the (suitably time-shifted) freezing forcing determined from the 1992 thermistor observations and ice properties determined from the fitting of 921126. It is certain that this freezing forcing is at least somewhat in error, owing to the differing drill characteristics and method of these years compared to 1992. The fits are not especially good (Figs. 5.18, 5.19, and 5.20). Hydraulic conductivities determined from fits to 901154 and 911146 agree quite well with those determined from 921124, while that from 891150 is quite a bit lower (Table 5.9). The compressibility estimates span the entire allowable range. The low hydraulic conductivity from 891150 is no surprise, since the character of its pressure record suggests that it is a blind hole anyway. Therefore these results do not contradict the results obtained from 921E24. 5.5.4 Summary of Bed Properties  The foregoing forward modelling has produced reasonably consistent estimates of hydraulic conductivity and wildly varying estimates of bed compressibility. The hydraulic conductivity falls in the range 1.35 x 10' m s-1 (for 7-, = 0.13m) to 1.3 x 10-8m s-1 (for Tb  = 0.025m, r, = 0.025 m), inferred from the 1992 unconnected borehole 921124. This  result is supported by the results from two boreholes (901154 and 91H46) drilled in other  Chapter 5. Model Results  '-c:a..  ^  121  1600  w c:  D  0 1400 0 w a. w i— <  1200  1000  w 0 = 800 w (2C 0 cn .0^.5^1.0^1.5^2.0^2.5  ^  3.0  ^  3.5  TIME (HOURS) Fig. 5.18: Best fit to 891150 freeze-in pressure record with rc = 0.025m.  891150 901154 911146 921124  K (m s-1)  a + n f3 (Pa-1)  4.0 x 10-19 6.0 x 10 7.0 x 10' 7.0 x 10-9  1.0 x 10' 2.0 x 10-7 1.0 x 10' 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.  Chapter 5. Model Results  ^  122  o_ _y 1050 11-1 1000 LU a_  950 900 850 800  0  750  Li-I  700  0 m  650  .0^.5^1.0^1.5^2.0^2.5  3.0  3.5  TIME (HOURS)  Fig. 5.19: Best fit to 901154 freeze in pressure record with re = 0.025 m. -  ▪ Chapter 5. Model Results  ^  0- 900 LL.1  ct 850 (/)  tn  • 800  rX j 750 700 0  650  CC 0 Ca  0  ^ ^ ^^ 5 6 2^3^4 1  TIME (HOURS)  Fig. 5.20: Best fit to 911146 freeze-in pressure record with rc = 0.025m.  123  Chapter 5. Model Results  ^  124  2 75  LU  70  CC  (/) uJ  65  CL 60 < 55  0  50  Ui  fX 45 0  220^240^260^280^300^320  ^  340  ^  360  DECIMAL DAY, 1990 Fig. 5.21: Pressure record of 90P07. Note the numerous sudden pressure changes and subsequent decays at various time scales. years (using a somewhat different drilling method). The very low estimate from 891150 is ignored as I conclude that it was a blind hole. This range is in good agreement with the values obtained by other investigators (upper bound of 3.6 x 10'm s-1 [Smart and Clarke, 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 wildly varying that I attach no significance to them. 5.6 Response to Sudden Pressure Forcing Sudden 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 be  Chapter 5. Model Results^  125  influenced by the response of the borehole itself? To answer this question I examine the response of an unconnected borehole to a sudden pressurization. If the system responds slowly, and the excess pressure decays with a time constant similar to those observed in pressure records, then it is possible that the borehole is responsible for the observed decays. If however the excess pressure decays away very quickly, then the observed slow decays must be due to some other processes unrelated to the borehole. Simulations were run with the ice properties found from consideration of 921124 freezein (see Section 5.4). To place an upper bound on the length of the decay, the longest reasonable borehole (921324 after only one day, L = 45.3 m), and the smallest pressure forcing 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, but the decay is much too rapid (Fig. 5.22). The pressure drops to 1/e of its initial value in just 56 s. This response is so rapid that such an event would likely go unnoticed on a pressure record sampled every 2 min, as taken at Trapridge Glacier during the summer field 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 be responsible for the decay events observed in Trapridge Glacier pressure records. The simulation was repeated with a hydraulic conductivity of zero, to simulate the response of a blind hole (Fig. 5.23). The decay is considerably slower, with a 1/e time of 2.6 days. The decay would be even slower for stiffer ice.  Chapter 5. Model Results  ^  0^20^40^60^80 TIME (SECONDS)  126  ^  100  ^  120  Fig. 5.22: Sudden pressurization at t = 1 day of hole 921124 (simulation).  127  Chapter 5. Model Results^  .0^.^1.0  ^  2.0  ^  3.0  ^  4.0  TIME (DAYS)  Fig. 5.23: Sudden pressurization at t = 1 day of hole 92E124 assuming impermeable bed (simulation).  Chapter 6 Summary  In this thesis I have presented a model for the mechanical response of water-filled boreholes 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, connected to a uniform half-space permeable bed through a hemispherical cavity, which can be forced by an approximate top-down freezing forcing, or by sudden pressurization. The model was formulated using the non-linear transient viscoelastic rheology of Shyam Sunder and Wu. This is the first' time that this theology has been used to model borehole response or, to my knowledge, any other aspect of glacier physics. The resulting system of non-linear partial differential equations and algebraic equations was solved using the method of lines, with integration by the ddriv3 integrator of Kahaner and Sutherland. The model was tested by comparing model results to the analytic solutions for Glen-Law ice, linear elastic ice, and the response of a spherically symmetric aquifer to a step head change on its inside boundary. Satisfactory solutions were found, although the stiffness of the problem did limit accuracy, especially for (incompressible) Glen-Law ice. Spatial oscillations of the ice solution were observed for some choices of spatial coordinate transformation, possibly as a result of using a centreddifference approximation of the spatial derivatives. The model gave qualitatively correct results for various tests of the coupled system. The model was then used to infer the hydraulic properties of the bed under Trapridge Glacier, by forward modelling of the freeze-in pressure records of six boreholes. The 128  Chapter 6. Summary^  129  model freezing forcing was inferred from the freeze-in temperature records from thermistor strings in two boreholes. Ice properties were then inferred from the modelling of the freeze-in pressure records from two blind holes. These results were not considered reliable, as the model results did not fit the pressure records well. This is probably due to 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 properties 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 from other years were poorer. The hydraulic conductivity was found to be in the range 1.3 x 10-9-1.3 x 10-8 ms', which agrees well with previous estimates (upper bound of 3.6 x 10-7 m s - 1 [Smart and Clarke, submitted], 1.27 x 10-8 m s-i [Clarke, 1987], and 2.2 x 10-8 m s' (T. Murray, unpublished data, 1992). The hydraulic conductivity was found to be most strongly influenced by the assumed freezing forcing and size of the cavity in the bed at the bottom of the borehole. It was found to be rather weakly dependent on the assumed ice properties. The bed compressibility was found to vary so much that no useful estimate of it was obtained. In addition, the response of blind and unconnected boreholes to sudden pressurization was 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 so rapid that it is unlikely to be resolved by the sampling rates used for Trapridge Glacier pressure sensors. It can therefore be concluded that the response of the borehole system is essentially instantaneous compared to the signals being measured; hence the observed signals reflect the true basal conditions, rather than the influence of the borehole (except during the initial freeze-in period when the freezing forcing dominates). Some recommendations can be made, based on this work. The approach used in this  Chapter 6. Summary^  130  thesis 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 ice model is not justified, since the ice deformation was found to be much less important than  water 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. This should be replaced by a proper thermal freeze-in model, such as that of Humphrey and Echelmeyer, to accurately model the important freezing forcing. Then perhaps useful estimates of ice properties could be made, as well as a better hydraulic conductivity estimate. The hydraulic conductivity estimate could also be improved if the size of the basal cavity at the bottom of the borehole could be better constrained. The use of a formal inversion scheme instead of the trial and error forward modelling approach used here would improve the reliability of the results. The finite-difference approximation used for spatial derivatives in the ice model should be re-examined. It is possible that replacement of centred-differences by an asymmetric differencing 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 as interpreting pressuremeter results, which have previously been interpreted only in terms of 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 from  in-situ pressuremeter tests. Finally, the pressure signals observed at Trapridge Glacier indeed reflect true basal conditions, and the influence of the borehole on the signal need not be considered, except during the initial freeze-in period.  References  Abramowitz, M., and I.A. Stegun (Eds.), Handbook of Mathematical Functions, 3rd ed., 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. Reidel, Dordrecht, Holland, 1987. Blake, E.W., The deforming bed beneath a surge-type glacier: Measurement of mechanical 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 deformation, J. Glaciol., 38, 388-396, 1992. Budd, W.F., and T.H. Jacka, A review of ice theology for ice sheet modelling, Cold Reg. 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-type glacier 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 subglacial 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., Clarendon, 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 at the 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. 131  References^  132  Fung, Y.C., A First Course in Continuum Mechanics, 340 pp., Prentice-Hall, Englewood 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 disturbed boreholes, J. Glaciol., 37, 414-419, 1991. Humphrey, N., and K. Echelmeyer, Hot-water drilling and bore-hole closure in cold ice, J. Glaciol., 36, 287-298, 1987. Jaeger, J.C., and N.G.W. Cook, Fundamentals of Rock Mechanics, 515 pp., Chapman and Hail, London, 1969. Jarvis, G.T., and G.K.C. Clarke, Thermal effects of crevassing on Steele Glacier, Yukon Territory, 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 basal water 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-1983 surge 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 measured with the pressuremeter, Can. Geotech. J., 25, 250-261, 1988. LeGac, H., and P. Duval, Constitutive relations for the non-elastic deformation of polycrystalline 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^  Cold Reg. Sci. Technol.,  133  4, 255-268, 1981.  Murat, J.R., P. Huneault, and B. Ladanyi, Effects of stress redistribution on creep parameters determined by a borehole dialatometer test, in 5th. International Symposium 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 experiments 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 Parameters 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, Cold Reg. Sci. Technol., 16, 45-62, 1989. Shyam Sunder, S., and M.S. Wu, A multiaxial differential model of flow in orthotropic polycrystalline ice, Cold Reg. Sd. Technol., 16, 223-235, 1989. Shyam Sunder, S., and M.S. Wu, On the constitutive modeling of transient creep in polycrystailine 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 subglacial water 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.1982  References^  134  Spiegel, M.R., Mathematical Handbook (Schaurn's Outline Series), 271 pp., McGrawHill, 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 from induced changes in basal water pressure: a theoretical framework for borehole response tests, J. Glaciol., in press, 1992. Streeter, V.L., and E.B. Wylie, Fluid Mechanics, 7th ed., 562 pp., McGraw-Hill, New York, 1979. Van Krevelen, D.W., and P.J. Hoflyzer, Properties of Polymers: Correlation with Chemical Structure, 427 pp., Elsevier, Amsterdam, 1972.  Appendix A  Summary of Model Equations  This appendix presents the complete system of equations for the coupled ice—bed— borehole problem, in a form that can conveniently be programmed. Water density, previously termed po, is written in this appendix as pw. The flow law exponent N 3, and the expressions containing N have been simplified accordingly. Dependent Variables hi (for j = 2, Ma — 1)  Maq Bi,^  e,  ^and crz,i (for i = 1,M)  (for i = 1, M — There are a total of Ma + 9M — 2 dependent variables. Initial Conditions h; = polPwg (for j = 1, Ma — 1) maq = 0  = Bo,^= Re,i =-- Rz,i = 0, =^Cuz,i  = 07 crr,i = crz,i = — Po  1  (for i =  1, M)  Differential Equations  am ^ aq 27r-pw exp (pb(0 — po)1K at^  135  (A.1) r,  Appendix A. Summary of Model Equations^  For j = 2, Ma — 1  136  UIi^K^a2 h at — pwg (a + ni3)14 a R2  (A.2)  a Bi = fr E ( 1 )3 d 2 a at^' k,Billi^eq's  (A.3)  For j = 1,M  ORr,i . 2 Aza4.,i  at^3^at 2^aet i a R,,i^ . AE , at 3 at a Rz,i^2^aftz AE ^ = at^3 at '' • aei^act ,i a 6,..,i .^ at^at + at aq,i^a et, i^a esi; i _^ ^ • +^' at^at^at i aezi a eLaetz ^ ^ , ' + ^' • at^at^at  (A.4) (A.5) (A.6)  (A.7) (A.8) (A.9)  Algebraic Constraint Equations  For i = 1, M 0^e z i  (A.10)  •  For i = 2, M — 1  0  =  Wi  ace  R  er,i •  (A.11)  For i = M 0 = ar,m +Po.  (A.12)  For i = 1, when no pressure forcing is applied, 0 = mw — mo maq m.f.  (A.13)  For i = 1, when sudden pressurization is applied and t < tramp, 0 = crr,1 pramp + Po.  (A.14)  ^  Appendix A. Summary of Model Equations  ^  137  For i = 1, after application of sudden pressurization (t > tramp), 0=  MYST  Mramp  Maq^Maq,ramp  + m1 — raf,rarnp•  (A.15)  Auxiliary Relations (7r,1  11 1  =  (A.19)  cz,i^=^Cez,i  67),i  2 p.u9 's 1 0.^  Ee  .  2p.^z's  A 2/L^  +e  (A.20)  + eu  (A.21)  +  ,  — — (cr, ^(re^ Cre = ,^3^. ti^  crz,i  Creq,i  + z,i)  cro,i  ,  3 •  kCfrCre,i Cfr^ I.^ 1  az,i)  (Crr,i  A 2p. (3A + 210 (crr,i A 2/./ (3A + 2/t) (ar,i  1  Grr,i  (A.18)  =  69,i^=^ee ^  1 2/./.°." 1  (A.17)  Pwg  E aiebit  L=  Er^  (A.16)  Pwg Po  1, — — k^09,i 3^' far,i2^Icrei  +  11 30-2 =^(— eq'i V '1^2^,^  az,i)  (A.22) (A.23)  (A.24)  CrZ  (A.25)  cr.,)  (A.26)  Crz,i)  (A.27)  (A.28) (A.29)  ^ ^ ^ ^  Appendix A. Summary of Model Equations ^  A  As,j  494,i  Ur,i  (A.30)  dt 0 de)  = A,,i'cre,i  at^ ae. ^ = at  138  (A.31)  -d—t o de)  (A.32)  (Tli I 0  =^Crr,i^Rr,i_  (A.33)  CrO,i  =^cre,i  (A.34)  crc!,i  =^— Rz,i  —  (A.35)  1 ( „rd^j_ — —3 Vr,i^'JO,:^'z,s  (A.36)  1 j d^d ur 9 oi ^Vrr Cr° , — 3^, i^Cr 0,i + Cr z,i  (A.37)  =  d  az,i  crd  a  =  az ,  3 eq,i  (^  d  + Cr z,i  3 (^ )3cred 2 BiVi  At ^d  1  Cri• i^4i 3^, —  , d 2) 2^, d 2^ Crr,i^Cre,i^Crz,i  —  / d (CIE Aticr  at•  r's  dt) 0  afte,i^d^de At^(—) at^°A dt 0  aez,,^de At^— z'i dt) 0 at^ a OR  rxr ar  s Cr^vv O,i  =  1^for R = in r transformation X Ri for R = rx- transformation  (A.38) (A.39) (A.40)  (A.41) (A.42) (A.43)  (A.44)  (A.45)  ^  139  Appendix A. Summary of Model Equations ^  for T < 263.12K Vo exp (Qi/3Rga„Ti)^ =^ Vo exp [(QL — QH)/NRgasT263] exP [QH/NRga,T] for T> 263.12 K .ffrb2  raw  (A.46)  exP (-0r,1 — Po)] [1 — exp (-0Pwg.0]  2 +yrr:Pw exP [13 (—aro. — po)]  (A.47)  2^ 2 3 irrb o mo =^' [1 exp (-13pwgLo)] -dirrcpw  (A.48)  Og  mf = 7r71,0  (L — Lo)^  (A.49)  mf,ramp = 771,0 (-tramp — Lo)^ (A.50) Maq,ramp = rnaq(t =- tramp) ^  rb = rbro (1 +  694)^  (de) = 1s  ^  dt 0  ^for  df dR df i dR df dR  f;+1 — fi-i 2AR  (A.52) (A.53)  for i = 2, M —1^(A.54)  —fi+2 4,fi+1 — 3fi1^ (A.55) = 2AR  3fi -^+ fi-2  for i = M^(A.56)  2AR  df —fi+2^— 3f.; ^for 2ARa dRa  d2 f dR!  (A.51)  fa+i — 2f; AR!  j = 1^ (A.57)  for j = 2, Ma  where f is the dependent variable (e.g., cr,., fe, 10.  —  1^(A.58)  


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