Thermo-Mechanical Regime of the Greenland Ice Sheet and Erosion Potential of the Crystalline Bedrock

: Past glaciation is known to have caused a substantial morphological change to high latitude regions of the northern hemisphere. In the assessment of the long-term performance of deep geological repositories for radioactive wastes, future glaciation is a critical factor to take into consideration. This study develops a thermal-mechanical model to investigate ice sheet thermal evolution and the impact on bedrock erosion. The model is based on comprehensive ﬁeld data resulting from international collaborative research on the Greenland Analogue Project. The ice sheet model considers surface energy balance and basal heat ﬂux, as well as the temperature-dependent ﬂow of ice that follows Glen’s law. The ice-bedrock interface is treated with a mechanical contact model, which solves the relative velocity and predicts the abrasional erosion and meltwater ﬂow erosion. The numerical model is calibrated with measured temperature proﬁles and surface velocities at different locations across the glacier cross-section. The erosion rate is substantially larger near the glacier edge, where channel ﬂow erosion becomes predominant. The abrasional erosion rate is averaged at 0.006 mm/a, and peaks at regions near the ridge divide. The mean meltwater ﬂow erosion rate in the study area is estimated to be about 0.12 mm/a for the melted base region.


Introduction
In Canada and many other countries, deep geological disposal is being considered for the long-term management of radioactive waste. The objective is to contain and isolate the waste in a deep geological repository (DGR) at depths of hundreds of meters in a suitable geological formation for periods of up to one million years. Protection of the near-surface environment is achieved by a multi-barrier system, which should be sufficiently robust to resist extreme geological events that can occur within that long time frame.
In the past million years, the Northern hemisphere was subjected to many glaciationdeglaciation cycles. In Canada, there were nine such cycles, each one lasting approximately 120,000 years [1]. At the glacial maximum, the ice sheet thickness reached up to 5 km, imposing loads of 50 MPa on the earth's crust and depressing its surface by hundreds of meters. The advance and retreat of the ice sheet erodes the earth's surface and is accompanied by permafrost development and thawing. At the glacial retreat, the meltwater creates new lakes in depressions and recharges the groundwater. Glacial cycles, therefore, profoundly altered the surface geomorphology and subsurface environment, perturbing their thermal (T), hydrologic (H), mechanical (M), and chemical (C) regimes (e.g., the THMC processes).
Glacial cycles that occurred in the past will cyclically occur again in the future, with the next cycle predicted to occur within 50,000-100,000 years [1]. A DGR in Canada, or in similar Northern environments, would be subjected to the tremendous THMC perturbations that could potentially affect the structural integrity and performance of the multi-barrier system, reactivate fault zones near the DGR and erode the bedrock surface. Modelling The GAP field investigations consisted of studies to understand interactions between ice sheet hydrology and subglacial hydrology [3,4]. The ice sheet surface water production was measured. Mechanisms of water flow from the ice surface to its basal interface were also investigated by remote sensing, and spatial ice displacement measurement. The project also studied basal conditions to understand ice sheet hydrology and the characteristics of the basal boundary conditions. The ice sheet was drilled at several locations to measure water pressures at the ice-bedrock interface. This information provided important input for the conceptualization of hydraulic pressure during glacial conditions for groundwater models applicable in Fennoscandia and Canada.
Glacier surface meteorological conditions were measured [3,4]. Seasonal fluctuation of air temperature, wind speed, and humidity were monitored by on-site meteorological stations. Downward and upward irradiations, both shortwave and longwave, were measured throughout several years at different elevations. The surface albedos at different sites were also characterized against melting/precipitation conditions [3,4]. These data provide sufficient information to determine the upper boundary conditions for the glacier model.
The major findings from the GAP study [3,4] are summarized in Figure 1b that shows the main THMC processes, and Table 1 shows the characteristics of the ice sheet and underlying bedrock. Transient meltwater processes on the ice sheet surface are seasonal phenomena limited to 3-4 months in summer, lowering the ice thickness by up to 3-4 m, and completely dominating the meltwater production. At least the top 70-80 m of the ice sheet is firn, a permeable mixture of snow and ice. Under the firn is solid ice it is practically impermeable. However, the ice could be fractured under uneven deformation, which can cause abrupt draining of surface melt, and establish surface-to-bed pathways. Borehole drilling data [3,4] showed that the entire outer flank of the study area has a melted bed, with liquid water present, while further inland frozen basal conditions prevail. The location of the boundary between melted and interior frozen conditions is highly sensitive to the heat flux from surface melting through surface-to-bed pathways, but is relatively insensitive to the heat generated by sliding at the bedrock interface. It is implied that much of the bed extending 250 km inward of the margin is covered by water. The water head measured below the ice sheet is generally equivalent to the height of the overburden mass of the ice sheet, e.g., the subglacial water pressure equals 92% of the annual average ice thickness [3,4].
A fundamental characteristic of the subglacial environment is temperature. Once the bed of the glacier reaches the pressure melting point, it can produce meltwater, and the ice begins to slide on its bed, producing frictional heat that can further contribute to ice melting. The thermal regime is dependent on stresses within the ice, atmospheric temperature, geothermal heat flux, and the flow of ice. Therefore, in zones of the wet warm base, the ice sheet is supported by buoyancy pressures, and even with physical contact, transmits limited effective pressure to the bedrock. In addition, the rock at depths has low permeability, limiting the infiltration of the meltwater, thus allowing mounding of water within the ice sheet.
The stable isotope profiles indicate penetration of glacial meltwaters at this site to the subglacial bedrock, which is likely facilitated by the local structural geology and fracture distribution and the presence of high hydraulic gradients. Direct sampling of bed materials and ice borehole photography show that the bed consists of bedrock covered by a relatively thin veneer of sand and gravel-sized sediment [3,4]. The sediment cover appears to be patchy, with some areas having minimal sediment coverage and others with a sediment cover reaching tens of cm in thickness. Fine sediment was noticeably lacking, implying that such grain sizes were washed away by subglacial water flow. The ice itself often carries debris in a layer extending to, at most, approximately 10 m above the bed. The size of the debris may extend up to very large boulders plucked from the bed, based on observations of glacial erratics (large boulders deposited by ice) at the margin.
In the glacier margin, several boreholes were drilled to explore the sedimentation characteristics. Permafrost subsoil was observed in exposed zones outside the glacier margin. At depths below approximately 300 m in DH-GAP04, the rock type changes from mafic to felsic/intermediate gneisses [4], with an accompanied change in abundance of gypsum in fractures, suggesting a dominant mechanism of rock-fluid interaction. The hardness of the gneiss is 7, elasticity is 40~60 GPa, and its compressive strength is reported to be 158 MPa [5].
The above information was used to construct a generic thermal-mechanical coupled model for the GAP study area.

Mathematical Modeling of the Thermal-Mechanical (TM) Regime of the Greenland Ice Sheet
The thermal and mechanical processes of ice sheet are closely related. To appropriately estimate the glacial erosion rate on the bedrock, we need to study the interaction at the interface between the ice sheet body and the bedrock. The slip of the ice sheet governs the bedrock erosion rate, and is also governed by creep that is strongly temperature-dependent. Therefore, the thermal regime of the ice sheet has to be considered in the model.

Conceptual Model
A conceptual thermal-mechanical (TM) model [3] was developed based on the above observations, as illustrated in Figure 2. The ice sheet is conceptualized as a viscoelastic body that deforms and slides at its base, under the force of gravity. The ice sheet deformation and basal slip depend on the temperature distribution within the ice body influencing the viscoelastic properties of the ice, and its state. The thermal regime in the ice sheet depends on the heat exchange with the atmosphere and solar irradiation and the geothermal flux from the geosphere. The heat exchange with the atmosphere depends on the air temperature, which is warmer near the ice edge and becomes colder at higher altitudes inland. The net solar irradiation is uneven on the surface of ice sheet, as the albedo of surface depends strongly on its composition, i.e., ice allows more solar irradiation than snow. Warmer temperature and more ice coverage near the edge promote ice melting at shallow depths, and the meltwater infiltrates through cracks to the base. Therefore, wet base conditions were found up to approximately 250 km inland, resulting in higher basal flow velocity. In the area of active drainage, that extends more than 100 km inland, erosion, due to meltwater flow combined with abrasion from basal slip, results in a higher erosion rate than the drier base.

Heat Equations
The energy balance equation with the consideration of phase change and assumptions of local thermal equilibrium among the ice and liquid (meltwater) phases and heat conduction as the only mechanism of heat transfer is given as where U is specific internal energy, λe is the bulk thermal conductivity of the ice-liquid Based on this conceptual model, a mathematical model was developed to simulate the thermo-mechanical regime in the ice sheet. The main model outputs are the temperature and stress distributions in the ice sheet and its displacement velocity. Those outputs will be compared to the measured field data. The abrasion rate of the geosphere underneath the ice sheet is heavily dependent on the temperature and stress distribution at the interface. A flow model was also developed to calculate the discharge rate in the area of active drainage.
From this discharge rate the meltwater erosion rate was estimated. The total erosion rate is due to a combination of abrasion and meltwater flow. Both types of erosion rates were calculated using empirical equations. The mathematical model and erosion rate equations are described next.

Heat Equations
The energy balance equation with the consideration of phase change and assumptions of local thermal equilibrium among the ice and liquid (meltwater) phases and heat conduction as the only mechanism of heat transfer is given as where U is specific internal energy, λ e is the bulk thermal conductivity of the ice-liquid mixture, Q T is the thermal source term. The specific internal energy is given as δ TT c is the Kronecker delta; L f is the specific latent heat of phase transition of ice at 333 kJ/kg.
where T c is the melting point of ice at 273 K. When density varies during phase transition, the parameter α m is where S i , S l , ρ i , ρ l are respectively the volume ratio and density of ice and water. The thermal source term Q T is composed of frictional heating (q f ), and the geothermal flux (q g ).
Q T = q f + q g (5) Frictional energy dissipation will occur at the sliding ice-rock interface and contribute to a heat flux of similar magnitude as the geothermal one (q f = 0.03-0.3 J m −2 s −1 [6]). The frictional heating is defined as a product of frictional force and relative velocity: where µ is the interfacial dynamic frictional coefficient, σ is the normal stress, v is the relative velocity [7]. If a frictional contact interface is assumed, the relative velocity can be replaced with basal sliding velocity u b . The bulk thermal conductivity is given by a volumetric average approximation 3.

Mechanical Equations
The equation of equilibrium is applied for stress evaluation, where C is the tensor of material stiffness, u is the displacement vector, g is the acceleration vector, due to gravity. The total strain is given as and the stress-strain relationship, assuming that ice is a viscoelastic material, is given as where ε is the total strain, σ is the stress, γ is the viscous strain, and D is the elastic stiffness tensor. The viscosity of the ice sheet have been extensively studied [8]. Note that temperature is a factor that strongly affects the rheology of ice. Nonlinear behaviors are, thus, expected to play an important role in the constitutive models. Plasticity may not be significant, since viscous flow is believed to predominate the mechanical behavior. The nonlinear Glen's law was adopted to model the rheological behavior of the ice sheet.
Glen's law is applied for pressure-dependent creep modeling of ice sheet where A is model constant for viscosity, and n ≈ 3. A typical value of glacier ice viscosity (e.g., under 10 kPa shear stress) is 1 bar year, equivalent to a viscosity of some 10 12 Pa s, about 10 15 times that of water [9]. Note that A is a function of temperature. We adopted a power law as below, which is comparable to that of [10] A = A 0 10 0.1(T−273) (12) where A 0 is the rate constant at 0 • C, and is a parameter that is scale-dependent. Another form of temperature-dependence takes the form of Arrhenius equation where R is the ideal gas constant, T is absolute temperature, and Q is the creep activation energy. Q is reported to be dependent on temperature, in the range of 60-80 kJ/mol below 263 K, and raised to between 120-167 kJ/mol above 263 K [11]. Our modeling shows that Q = 100 kJ/mol gives the most consistent results with respect to the measured horizontal ice velocity. However, a nonlinear equation of the activation energy would better reflect the physics of the rheology of ice. Ice flow velocity varies with depth, with contributions from both the basal slip rate, u b , and the internal deformation rate u i [12] The basal slip ratio is defined as u b /u i . This parameter is hard to measure directly, and thus, was estimated through fitting of ice sheet morphology and surface motion data with the appropriate assumption of ice viscosity [13]. Variability in the estimates of the slip ratio is high from different sources, ranging from 10-30 to less than 5 within 100 km of glacier margin at the same site as the GAP study area [14]. In distances beyond 100 km inward the glacier, the basal slip ratio decreased substantially to around unity.

Erosion Laws
The dominant processes of bedrock erosion under an ice sheet are: Abrasion, plucking and quarrying, ripping, and meltwater erosion. For the above processes to take place beneath an ice sheet, at least two conditions must be met: (i) The base must be warm and wet, to allow meltwater to move rapidly and to allow basal slip; (ii) the base must be free of sediment to allow ice/rock contact. Glacial ripping is a newly recognized erosion process that involves the jacking, disruption, and transport of large, angular boulders [15,16]. The process of ripping requires groundwater overpressure, and is discovered beneath the retreating margin of the Fennoscandian Ice Sheet that operates to depths of several meters [16]. Ripping is similar to plucking or quarrying, but can relocate a large collection of joint blocks locally. Therefore, it may be a prominent driving force to modify the near-surface bedrock. This newly recognized process remains to be developed in a wellcalibrated erosion law, and is outside the scope of the present study.
The total erosion rate (E r ) would then be a combination of the above three processes where E a is the abrasion rate, E q is the quarrying rate, and E m the meltwater erosion rate. As a result of glacial erosion, sediment is transported by sliding ice and meltwater discharge, and deposits to form eskers near the edge [28]. Sediment discharge increases with warming climate, influences melting/collapse of glacier drainage channel, impacts downstream aquatic habitats and river infrastructure [28][29][30][31]. Sediment production and transport are important processes to be considered in future studies. However, the present modeling focuses on the three processes of abrasion, quarrying, and meltwater erosion.

Abrasion Rate
Subglacial abrasion operates regarding shearing at a certain velocity through the contact forces between the bed and clasts embedded in the ice. The abrasion rate is also affected by availability and concentration of entrained sediment. This is expressed in a well-known model by Hallet [17] as where K a is a model constant, S f is entrained sediment concentration, S c is critical slope reported at 1.5, u b is the basal slip rate, and F c is the contact force between the bed and the entrained clasts.
More recently, rock-on-rock wear under known normal loads was measured by Lee and Rutter [30] for a range of silicic rocks of porosities varying between 7% and 28%. The wear rate is related to the contact normal stress and the porosity (as shown in Figure 3). An empirical wear law was given as [30] where dς/dx is a dimensionless wear rate (meter wear (ς) per meter shear displacement (x)), σ is contact normal stress (MPa), φ is porosity, and B, n, and m are empirical constants. Multiplied by basal slip rate the above equation gives another form of erosion rate in m/s, i.e., Assuming that the entrained clasts mainly contribute to frictional abrasion erosion, then where B = 10 −8.11 , n = 8.33, m = 4.5, σ is in Mpa, c is volumetric concentration of entrained clast. These parameters are best fit of the experimental results on various types of rocks, as shown in Figure 3. Assuming that the entrained clasts mainly contribute to frictional abrasion erosion, then = (19) where B = 10 −8.11 , n = 8.33, m = 4.5, is in Mpa, c is volumetric concentration of entrained clast. These parameters are best fit of the experimental results on various types of rocks, as shown in Figure 3. Compared to the empirical Hallet abrasion law, which is a function of slip rate, this empirical wear law provides important perspectives on measurable properties of the bedrock beneath a sliding glacier. In addition, both the in situ normal stress and the horizontal velocity are available in either field measurements or numerical simulations. Uncertainty corresponds to the determination of the normal stress for a temperate melting bed. As pointed out by Lee and Rutter [30], normal stress concentration is expected to develop at clast/bed contacts from drag against ice flowing toward the bed.
These experimental data are consistent with wear rates inferred from silt outflow in subglacial streams, provided a tenfold stress concentration can occur relative to the ice overburden pressure. In this study, we use the Lee and Rutter model [29] to estimate the abrasion rate. The use of the empirical relationship between long-term erosion rates and basal power together with modeled basal energy dissipation for the ice sheet provides a sound basis for estimating erosion rates in the GAP study area.

Quarrying
Quarrying involves fracturing and crushing bedrock into blocks. The quarrying of bedrock favors low effective stress, high slope, and high sliding velocity. Hence, a ratelimiting relation for quarrying was given by Reference [27] where Kq is a model parameter about 10 times of Ka, Γ(β) is a function of the bed slope  in the direction of sliding, and diminishes when slope tilts upward. Quarrying needs a high slope and high sliding velocity. The GAP glacier is slowmoving. Subglacial sliding velocity, although not reported in the GAP study, is not ex-  Compared to the empirical Hallet abrasion law, which is a function of slip rate, this empirical wear law provides important perspectives on measurable properties of the bedrock beneath a sliding glacier. In addition, both the in situ normal stress and the horizontal velocity are available in either field measurements or numerical simulations. Uncertainty corresponds to the determination of the normal stress for a temperate melting bed. As pointed out by Lee and Rutter [30], normal stress concentration is expected to develop at clast/bed contacts from drag against ice flowing toward the bed.
These experimental data are consistent with wear rates inferred from silt outflow in subglacial streams, provided a tenfold stress concentration can occur relative to the ice overburden pressure. In this study, we use the Lee and Rutter model [29] to estimate the abrasion rate. The use of the empirical relationship between long-term erosion rates and basal power together with modeled basal energy dissipation for the ice sheet provides a sound basis for estimating erosion rates in the GAP study area.

Quarrying
Quarrying involves fracturing and crushing bedrock into blocks. The quarrying of bedrock favors low effective stress, high slope, and high sliding velocity. Hence, a rate-limiting relation for quarrying was given by Reference [27] where K q is a model parameter about 10 times of K a , Γ(β) is a function of the bed slope in the direction of sliding, and diminishes when slope tilts upward. Quarrying needs a high slope and high sliding velocity. The GAP glacier is slowmoving. Subglacial sliding velocity, although not reported in the GAP study, is not expected to be significant compared to the surface ice motion. In addition, the borehole drilling sample shows only thicknesses of centimeters to decimeters of fine sand sediments in subglacial interface with bedrock. The absence of blocks of sediments, together with other evidence, suggests that quarrying might not be a significant source of erosion for this study. By eliminating quarrying from the consideration, our erosion model can be focused on frictional abrasion and flow erosion. This further simplifies the problem and reduces uncertainties in model parameters.

Basal Water Flow and Erosion
According to the GAP study ( [4], , the surface runoff ends at the elevation of 1600~1700 m above the sea level. Borehole pressure observations confirmed that meltwater flows in a basal drainage system at the ice/bed interface at least 30 km from the ice margin. The basal drainage system is found to have a large capacity to transmit high volumes of water. The basal drainage system is also found to be well connected to the surface meltwater through a highly permeable glacier body, as suggested by the fast response of diurnal water pressure variations. Borehole sampling showed that the bedrock is largely exposed or covered with a thin layer of sand sediment.
The fast-moving basal channel flow may mobilize the sand sediment by saltation when the velocity surpasses a critical value. The drifted particle trajectory follows the fluid flow under viscous skin friction and flow pressure, and eventually drops back to the bed under gravity, at a distance of hop length from the original position on the bed. The settling particle impacts the bed at a certain velocity, and the dissipated momentum causes erosion [31]. As a critical mechanism governing sediment transport, saltation has been well investigated in the literature. This study adopts some of the classic concepts to characterize the flow regime and erosion behavior of the glacier's base drainage system.
Beaud et al. [32] and Lamb et al. [33] developed a theoretical model to predict this type of erosion. The saltation-abrasion model takes the following form: where I I is the particle impact rate against the bed, V i the volume of bedrock removed upon impact, and F e is the fraction of the bed that was sediment-free and exposed. The above equation can be further written as where L b is particle hop length, k v is dimensionless erodibility (k v = 10 6 ), E is Young's modulus (E = 50 Gpa), σ T is tensile strength (7 Mpa), q b is volumetric sediment flux per unit channel width (m/s), u i is particle impact velocity, and F e is the bed exposure fraction (0.5). The particle suspension concentration C b (0.001 kg/m 3 ) and average flow velocity u can be used to calculate the volumetric sediment flux q b = C b u/ρ s (~4 × 10 −7 u). The particle hop length L b can be calculated from the formula given in the Appendix A, which covers both the turbulent and viscous flow regime. With all these parameters accounted, the erosion equation becomes This equation is based on fundamental material properties that are easily accessible. The only variable to be determined is the particle impact velocity. Therefore, a computational fluid dynamic (CFD) model was developed separately in this study to investigate the behavior of basal flow and its effect on bedrock erosion. The numerical model assumes a 2D thin layer of the channel of height H, and length L (L = 100 km approximately the limit of melting distance from ice margin). Appropriate boundary conditions were assigned to simulate the pressure head distribution along the drainage flowline.
Cowton et al. [34] conducted a detailed investigation of the subglacial drainage and erosion of the Greenland ice sheet at the same location as the GAP study. The discharge and catchment area are used to calculate the average channel flow velocity, which is in turn used in the 2D CFD model to determine the channel height. The catchment is about 600 km 2 and extends 80 km along the flow line. Therefore, a cross-section of 7.5 km is a reasonable simplification. The field measured discharge peaks at 300 m 3 /s in the summer season. The channel is found to be approximately 1 cm high, which results in a flow rate of 5 m/s by the CFD model vs. 4 m/s as observed. The channel is comparatively shallow with respect to the sediment grain size of 1 mm. Therefore, the classic model of Engelund and Fredsoe [35] for particle velocity that correlates particle size and bed shear velocity is inapplicable [31], because the predicted saltation height according to the formula (estimated jump height of particles δ D ≈ 41.5) is too high for the channel in the dimension of H/D = 10. The derivation and calculation are given in the Appendix A.
Due to the turbulent nature of the subglacial channel flow, the fine particles can be regarded as moving in suspension. The experimental observation of Staben and Davis [36] may be applicable to extrapolate the particle velocity distribution. Fine particles flowing in a channel are found to be excluded from the wall region and focused near the center of the channel, resulting in an average velocity that is similar to the average fluid velocity [36]. Therefore, a simplified treatment would be to assume that the particle velocity is equivalent to the average fluid flow velocity.
Based on these rationales, and by relating the above formula to the three months of summer peak flow rates of 5 m/s, the annual subglacial erosion rate, due to meltwater flow can be calculated as 0.85 mm/a. Note that this value is an upper bound estimate, and is not representative of the whole catchment as the erosion law is a cubic power function of basal flow velocity, which increases from negligible in the inner part of the glacier to the peak value at the estuary outlet boundary. This value is comparable with the reported value of 4.9 mm/a by Cowton et al. [34], who based their calculation on the observed sediment flux from the catchment. This finding provides confidence in our methodology for estimating the subglacial erosion associated with meltwater flow, by first calibrating the CFD model of the channel flow with field data, and then estimating the erosion rate with the mechanistic erosion model. It also stresses the efficacy of subglacial drainage of glacier meltwaters on the rapid erosion around the margins of the ice sheet.

Numerical Model of Ice Sheet Thermo-Mechanical Regime
The Finite Element computational platform COMSOL (Ver 5.3, COMSOL Inc., Burlington, MA, USA) was used in this study to simulate the thermo-mechanical regime in the ice sheet. COMSOL has versatile functions that enable easy adaptation and solution of mathematical equations with high nonlinearity in parameters and strong coupling of multiple physical processes.

Glacier Model
The closest ice core at the ice divide is at the Dye-3 site and provides an estimate of the vertical temperature profile at the ice divide. The width of the ice sheet at the same latitude of 67N is about 700 km. The ice divide at this site is about 500 km from the GAP site. Therefore, the 2D FE model has a dimension of 500 km × 3 km, representing the study area, as shown in Figure 4. The long horizontal span extends almost to the divide, a permafrost region where basal slip is insignificant. The glacier surface elevation is the best-fit approximate to the data reported in Reference [4].
The model is composed of two parts, i.e., the ice body and the bedrock. The upper boundary is subjected to seasonal air temperature, net surface irradiation, and albedo of the air-ice interface. Adiabatic and zero normal displacements are assumed at the right boundary, and the bedrock bottom is fixed. The ice-bedrock interface is simulated with a contact model that accounts for compression and friction mechanisms. The frictional coefficient of the interface adopts a nonlinear expression to simulate ice/rock abrasion, as further elaborated later on. Heat fluxes, due to friction, geothermal gradient, and meltwater, are imposed at the ice-bedrock interface.
A Glen's law is applied to the ice body, while the bedrock is assumed to be linear elastic. The relevant model input parameters are given in Table 1.

CFD model for Subglacial Channel Flow
The regional-scale groundwater flow in the GAP study area has been investigated in detail [37]. The pressure head distribution is found to follow the trend of glacier elevation, while the flux is centralized around subglacial drainage networks [37]. Here we developed a simplified 2D CFD model for the subglacial channel flow, by solving the Navier-Stokes equations of the following form, where ρ is the density of the flowing fluid, u is the velocity vector, p is pressure, τ is the viscous stress tensor, and g is the gravity. A rectangular channel in the dimension of 80 km × 1 cm is considered for the 2D numerical simulation. The channel height is varied in the model to search for the ideal value so that the exit discharge rate agrees with the field measurement.
The channel features impermeable bottom, and right-side wall boundaries. The left side boundary is assigned with a zero pressure head to serve as a discharge exit. The upper boundary is imposed with water pressure equivalent to the overburden weight of the ice sheet.

Air Temperature
Monitored data at the edge of the ice sheet can be fitted with a sine function, where t is time in month. The above temperature variation ( Figure 5) can be used as a boundary condition at the edge of the ice sheet. For upper elevations, ambient air temperature decreases with an increase of altitude at a rate of 0.007 K/m. As the variation range of average air temperature in Greenland is about 2 °C in the last 2000 years, it is thus justifiable to use the decades-long weather station data as a boundary condition for the simulation of the past 1000 years of evolution.

CFD Model for Subglacial Channel Flow
The regional-scale groundwater flow in the GAP study area has been investigated in detail [37]. The pressure head distribution is found to follow the trend of glacier elevation, while the flux is centralized around subglacial drainage networks [37]. Here we developed a simplified 2D CFD model for the subglacial channel flow, by solving the Navier-Stokes equations of the following form, where ρ is the density of the flowing fluid, u is the velocity vector, p is pressure, τ is the viscous stress tensor, and g is the gravity. A rectangular channel in the dimension of 80 km × 1 cm is considered for the 2D numerical simulation. The channel height is varied in the model to search for the ideal value so that the exit discharge rate agrees with the field measurement.
The channel features impermeable bottom, and right-side wall boundaries. The left side boundary is assigned with a zero pressure head to serve as a discharge exit. The upper boundary is imposed with water pressure equivalent to the overburden weight of the ice sheet.

Air Temperature
Monitored data at the edge of the ice sheet can be fitted with a sine function, where t is time in month. The above temperature variation ( Figure 5) can be used as a boundary condition at the edge of the ice sheet. For upper elevations, ambient air temperature decreases with an increase of altitude at a rate of 0.007 K/m. As the variation range of average air temperature in Greenland is about 2 • C in the last 2000 years, it is thus justifiable to use the decades-long weather station data as a boundary condition for the simulation of the past 1000 years of evolution.

Solar Radiation
The annual variation of wave radiation of the Greenland glacier has been investigated in detail [3,4]. Solar radiation in shortwave form (Figure 6a) can be modeled with the following equation, where t is time of year in month. Meanwhile, the longwave radiation consists of mainly two parts, i.e., the incoming and outgoing radiation. The net value is always negative, suggesting a net loss of energy emitted from the ice body. Although ice surface temperature fluctuates with the season in a range of 30~40 °C, the net longwave radiation (Figure 6b) has a narrow distribution with an annual mean of The altitude (H) dependence has been reported in Reference [4], and is averaged using a linear function as

Solar Radiation
The annual variation of wave radiation of the Greenland glacier has been investigated in detail [3,4]. Solar radiation in shortwave form (Figure 6a Sensible heat flux over ice is a function of wind speed and air temperature, and is given as Meanwhile, the longwave radiation consists of mainly two parts, i.e., the incoming and outgoing radiation. The net value is always negative, suggesting a net loss of energy emitted from the ice body. Although ice surface temperature fluctuates with the season in a range of 30~40 • C, the net longwave radiation (Figure 6b) has a narrow distribution with an annual mean of The altitude (H) dependence has been reported in Reference [4], and is averaged using a linear function as

Supraglacial Sensible Heat Flux
Sensible heat flux over ice is a function of wind speed and air temperature, and is given as where h is the heat transfer coefficient, T ext and T are external air average temperature at about 2 m above ground and surface temperature, respectively. The effect of wind speed is reflected by the parameter h, which is variable according to fluid dynamics. Oerlemans and Grisogono [38] studied the parameter h from analysis of katabatic flow over several melting glaciers, including West Greenland, where the GAP is located, and found that h increases quadratically with (T ext − T), and is independent of the surface slope. From their dataset, the sensible heat flux can be regressed into a polynomial equation of temperature difference, and the formula is This is consistent with another report of q 0 ≈ 50 W/m 2 in summer month (July) of Greenland west where peak ∆T ≈ 10 • C [39]. With this observation, we adopted a heat transfer coefficient value h = 5 W/m 2 K.

Radiation Albedo
One key factor governing the thermal balance of ice sheets is the radiation heat. Solar irradiance is the primary source of heat inflow, in terms of both shortwave and longwave radiation, while the heat loss at night and cooling circumstance take place regarding longwave emission. For dielectric materials like snow and polycrystalline ice, the transmittance of the electromagnetic wave is largely dependent on the wavelength and density of impurity (i.e., air bubble and dust particle). The wavelength dependence of the various forms of solid water ubiquitously shows a valley shaped curve [40]. For snow and white ice on the glacier surface, the absorption coefficients (i.e., 1-albedo) are found to be in the range of 0.09-0.4 for visible light, which is orders of magnitude higher than bubble-free ice usually present at depth. Longwave radiation has a higher absorption coefficient and has less penetration depth than shortwave radiation in the visible range of the spectrum.
Using the averaged absorption coefficient, 8 m of radiation penetration is estimated, which is close to the observed depth where ice temperature varies seasonally. This observation suggests that the thermal profile of the shallow ice is more prone to solar radiation than the direct heat exchange on the ice-air interface. Our modeling results confirmed this hypothesis. Therefore, the radiation budget is considered in the ice modeling in terms of a depth-dependent heat source distribution to reflect the temperature fluctuation.
Net radiation is the major source of ablation energy, and turbulent fluxes (a.k.a. sensible heat) are smaller energy sources by about three times. The radiation albedo defines the proportion of sunlight illumination reflected the downward input, and varies with elevation, and is dependent on ice sheet surface ablation state. For higher elevations above the equilibrium line, which is mainly snow-covered all season, the albedo can be as high as 0.9. In contrast, lower elevations where snow accumulated in winter easily melts in summer can have a minimum albedo of 0.4. The GAP study compiled relevant data that were measured on-site at different elevations. The GAP study also measured the net radiation in both longwave and shortwave forms, at several different elevations along the ice sheet flowline. This study took a simplified approach by averaging the albedo coefficient of each radiation component into a single value and correlating it to the elevation. Albedo is set at 0.4 for y < 750 m, and 0.9 at y > 1750 m, with a gradual transition from 0.4 to 0.9 in elevations between 750 and 1750 m. This parameter is used in our numerical model to adjust the radiation-related heat source on top of the ice sheet.

Basal Heat Flux
Geothermal heat flux in the GAP region is reported to be 30 mW/m 2 . Our modeling showed that this magnitude of geothermal heat flux is insufficient to heat up the glacier base to the observed value. Creep dissipation provides another heat source for the ice sheet bed. Due to the dissipation of deformation energy, a substantial amount of heat is generated, as indicated by Funk et al. [41], where A(T) is the temperature-dependent rate factor, n is the flow-law exponent, and . ε I I the second invariant of the strain-rate tensor.
According to our simulation, the internal heat generation at equilibrium state is in the order of 10 −4~1 0 −3 W/m 3 . This is equivalent in quantity to the geothermal heat flux, as identified in the GAP report by Harper et al. [4].
Various sources of evidence suggested that the exchange of meltwater from surface glacier to the base through tunnels and fissures in the body of the ice sheet is believed to be the major contributor to the melting base. A simplified procedure was implemented here, where heat (induced by supraglacial discharge) was input in proportion to the surface net radiation. That proportion was determined at approximately 5%, by a calibration procedure to fit the measured temperature profile. This procedure may overlook some detailed features of the heat transport phenomenon in the deep region, however, it ensures the melting condition for the base, which is a critical boundary condition for the mechanical model.
In summary, three sets of basal heat flux were applied in the model, consisting of geothermal heat flux at 30 mW/m 2 , creep dissipation heat at 30 mW/m 2 , and supraglacial discharge induced heat flux at approximately 5% of the surface net radiation.

Initial Condition for Temperature Profile
For the ice divide boundary, the temperature measurements from the Dye-3 station were applied. The Dye-3 station is located at the ice sheet crest southern GAP site, in a region with mountainous bedrock topography. Deep core drilling was performed to a depth of 2037 m, and the temperature was recorded at that whole depth. The basal temperature was −13.22 • C, and was well below the melting point at the bed of the ice sheet. The minimum temperature, around −20.5 • C, was located at around 200 m depth. The temperature gradient for the lower 800 m is lower than the gradient seen at Camp Century (0.012 K/m vs. 0.018 K/m), and the small reduction in gradient near the bottom may reflect the variability of the geothermal flux. An exponential function can best fit the Dye-3 temperature profile, and is applied in the model.
For regions near the ice sheet edge, the shallow portion has gradually increasing temperature, while the middle portion is at the minimum temperature. A polynomial function would best fit the temperature distribution at the edge.
Above the equilibrium line altitude at about 1550 m ASL [4], the surface of the ice sheet is covered with loose firn that has an exponentially decreasing porosity with depth. For regions that are covered with melting ice below the equilibrium line altitude (ELA), a firn layer of about 2 m thick exists. Therefore, our model assumes that the initial porosity distribution follows the equation given in Table 1 for domains above the ELA, while a constant initial porosity is assumed for the remaining domains. The thermal conductivities of firns at various densities are given as an empirical equation [42] λ = 0.138 − 1.01ρ + 3.233ρ 2 (33) Table 2 shows the boundary and initial conditions for the heat flow simulation.
, where x is distance from th edge of the ice sheet in km, and ReDepth(x,z) is relative depth between 0 and 1.

Basal Friction Model
The ice sheet bed is in contact with the bedrock, the interaction of which is essentially governed by the temperature distribution along with the interface. Under conditions of permanent frost, the glacier bed can either be bonded to the bedrock or undergo shear deformation. In regions with a melting base, the ice sheet separates from the bed, with the existence of drainage channels created by flowing meltwater. Therefore, the interaction between glaciers and bedrock will be affected by the variation of water head. The motion of the ice sheet is driven by gravity, but is resisted by the bottom friction effect. To model the subglacial boundary condition, a contact model is necessary to take friction into account. Here we develop a contact model with a temperature-dependent frictional coefficient in the range of 0.1 < µ < 0.5 at freezing temperatures down to −20 • C. A systematic iceon-rock friction experiment was conducted under varying temperatures, as reported by Reference [26]. Figure 7 summarizes the observed steady-state friction coefficient as a function of ice temperatures in the range of −20~0 • C. An exponential function was found to best fit the experimental data with high correlations, 15 20 0.693 (34) by Reference [26]. Figure 7 summarizes the observed steady-state friction coefficient as a function of ice temperatures in the range of −20~0 °C. An exponential function was found to best fit the experimental data with high correlations, = 0.5exp −0.5 −ln − − 273. 15 20 . (34) Figure 7. Temperature-dependence of steady-state friction coefficient of ice-rock (Barre granite) interface ( [25], normal stress at 100 kPa, R 2 = 0.98).

Parameter Determination
The Glen's law was fine-tuned in the model to reproduce the level of the observed ice velocity. The commonly accepted value of n = 3, as inferred from laboratory creep tests on ice [43,44]. A previous study [45] discussed the effect of polycrystalline preferential orientation on the creep of ice, and confirmed the variation in the value of power index n, e.g., which can be as high as 4 [46]. A recent study further demonstrated that Glen's law n = 4 instead of n = 3 based on modeling of the northern region of Greenland Ice sheet with frozen base [43], suggesting the ice is softer than usually believed. The grain size was found to influence the Glen's stress exponent n [43]. The mechanism of the rheology of ice involves dislocation creep (n ≈ 4), basal slip limited creep (n ≈ 2.4), and grain boundary sliding (n ≈ 1.8) under different stress levels and grains sizes, with n ≈ 3 over the range of strain rates found in most natural systems [11,43]. Caution needs to be taken in dealing with model input obtained from lab tests for large scale field models. Our calibration of the model with field data further emphasizes this.

Parameter Determination
The Glen's law was fine-tuned in the model to reproduce the level of the observed ice velocity. The commonly accepted value of n = 3, as inferred from laboratory creep tests on ice [43,44]. A previous study [45] discussed the effect of polycrystalline preferential orientation on the creep of ice, and confirmed the variation in the value of power index n, e.g., which can be as high as 4 [46]. A recent study further demonstrated that Glen's law n = 4 instead of n = 3 based on modeling of the northern region of Greenland Ice sheet with frozen base [43], suggesting the ice is softer than usually believed. The grain size was found to influence the Glen's stress exponent n [43]. The mechanism of the rheology of ice involves dislocation creep (n ≈ 4), basal slip limited creep (n ≈ 2.4), and grain boundary sliding (n ≈ 1.8) under different stress levels and grains sizes, with n ≈ 3 over the range of strain rates found in most natural systems [11,43]. Caution needs to be taken in dealing with model input obtained from lab tests for large scale field models. Our calibration of the model with field data further emphasizes this.
Studies of direct measurement of ice-rock abrasion indicated limited information. Shamsutdinova et al. [25] conducted laboratory studies on concrete-ice abrasion, and showed a low abrasion of concrete with the maximum abrasion coefficient of 10 −7 . This is consistent with the results of [47]. Laboratory ice-rock abrasion tests usually could not consider the effect of entrained fine sediments, which is known to enhance the erosion rate proportionally.

Computational Fluid Dynamics (CFD) Model for Subglacial Channel Flow
A 2D thin layer of the rectangular domain was developed for the simulation of channel flow. The fluid properties are assigned with those of water at zero degree. The length of the channel is assumed to be equal to the length of the flowline of the catchment drainage system of the GAP site (approximately 80 km) [33]. The Navier-Stokes equations for laminar incompressible flow were used to describe flow in the channel, and were numerically solved with the COMSOL code, with the following boundary conditions: The upper boundary of the channel was assigned with a pressure head equivalent to the overburden glacier loading; the channel bottom and the right boundary are treated as impermeable; and the left boundary is a discharge outlet with zero pressure head.
The channel height H was varied in a parametric study to calibrate the calculated outlet discharge with the measured discharge at the GAP site. The best match was obtained with a channel height of 1 cm. The distribution of flow velocity and shear rate can then be determined-providing input to calculate the subglacial erosion rate, due to meltwater.

Mathematical Modeling Results and Discussion
The glacier TM model is a 2D study with a representative profile taken along the glacier flowline. The subglacial channel flow is a 2D model as well. The meshing size is constrained to a maximum of 3 m in Z and 300 m in X directions. The solution time step is adjusted to be less than a day to reflect the daily variation of the upper thermal boundary. Convergence is automatically achieved by the COMSOL solver without manual set up. Figure 8 shows the vertical cross-sections with monitoring results available for calibration. The horizontal distances from the ice edge are 20, 40, 60, and 500 km, respectively. The cross-section at 60 km is from the Jakobshavn Glacier, 2 • north of the GAP site. The GAP study does not cover this distance from the ice sheet edge. Therefore, the temperature profile of the Jakobshavn Glacier is used as a qualitative calibration at this distance.   Figures 9-11 shows the comparison of modeling results with field measurements. The model is found to reproduce the observed trend of temperature distribution with depth. Basal melting is well identified from the modeling results, with about 10~20% thickness of the ice sheet being melted. Ice is in its phase transition stage as the temperature reaches the melting point. The simulated melting zone reaches about 150 km inward, and up to the 1100 m elevation, which is close to the inference from Harper et al. [4]. Figure 12A,B shows the estimated temperature field in the GAP report [4]. Figure 13 shows the simulated temperature field in the present study. The shape and extent of temperature contours as predicted by our model resembles the ones in the GAP report. This behavior is also confirmed by ice core temperature profiles. The base is generally warmer than the body of ice, and the basal temperature is consistent with the reported distribution of GrIS [48]. The melted base ( Figure 14) is predicted to be smaller in our model than in Reference [4]. Figure 8. Cross-section of the ice sheet model, with vertical boreholes highlighted in red that are compared with available monitoring data. The horizontal axis is directed towards the ridge divide, while the vertical axis indicates elevation. The three red vertical lines are respectively located at 20, 40, and 60 km from the edge.                 Figures 15 and 16 show the comparison of modeled surface ice movement with field observations. The observations were obtained by GPS monitoring. Since we do not consider in the model the ablation/accumulation that affects the elevation change, only the horizontal velocity is considered for calibration. The modeling results are consistent with the GAP site measurements. Note that the EGIG profile shown in Figure 15a is located at Jakobshavn Glacier and spans northeastward. This is not an exact comparison, but rather a trending characterization. The horizontal velocity peaks at about 80-100 km inland and decreases gradually to the level of 10 m/s at the ice divide. Bartholomew et al. [49] reported an average ice velocity of ~100 m/a extending to 80 km from the ice margin. The GAP study reveals a similar range of supraglacial horizontal velocity of 100-400 m/a. It seems,  Figures 15 and 16 show the comparison of modeled surface ice movement with field observations. The observations were obtained by GPS monitoring. Since we do not consider in the model the ablation/accumulation that affects the elevation change, only the horizontal velocity is considered for calibration. The modeling results are consistent with the GAP site measurements. Note that the EGIG profile shown in Figure 15a is located at Jakobshavn Glacier and spans northeastward. This is not an exact comparison, but rather a trending characterization. The horizontal velocity peaks at about 80-100 km inland and decreases gradually to the level of 10 m/s at the ice divide. Bartholomew et al. [49] reported an average ice velocity of~100 m/a extending to 80 km from the ice margin. The GAP study reveals a similar range of supraglacial horizontal velocity of 100-400 m/a. It seems, therefore, that the model proposed in this study can adequately reproduce the thermal and deformational behaviors of the ice sheet at the GAP site.

Surface Ice Movement in Horizontal and Vertical Cross-Sections
Maier et al. [46] measured shear strain rates in boreholes at a field site 30 km inland from the Western Margin of the GrIS (67.2N 49.6W). During the 2015/2016 winter, the mean surface velocity was 114.8 m/a. Ice deformation is reported to contribute 4.6 m/a, and basal sliding is responsible for 110.2 m/a, which is termed as plug flow that the basal motion dominates the surface motion. Our modeling shows the surface velocity is half of the above measurement, and basal sliding below 10~15% of ice thickness contributes most of the motion. This behavior is consistent in general with the observation of [46]. The discrepancy is likely attributed to (1) regional-scale velocity variability [47], and (2) the application of the ice-bedrock friction model that uses a larger than real friction coefficient.
The modeling results suggest a gradual transition of ice velocity from the glacier edge inland (Figure 16). The glacier edge zone where surface melting starts to penetrate to the ice bed moves the fastest, at the level of 100 m/a, which is consistent with the reported annual average value of the GPS measurements. The seasonal variation of temperature enhances the mobility of the ice sheet in summer, which is not considered in our model. Another important observation is the vertical variation of ice mobility, as shown in Figure 17. The lower 10% of the ice sheet contributes to~60% of the horizontal surface velocity, which is in good agreement with other studies [10,40].  Figures 15 and 16 show the comparison of modeled surface ice movement with field observations. The observations were obtained by GPS monitoring. Since we do not consider in the model the ablation/accumulation that affects the elevation change, only the horizontal velocity is considered for calibration. The modeling results are consistent with the GAP site measurements. Note that the EGIG profile shown in Figure 15a is located at Jakobshavn Glacier and spans northeastward. This is not an exact comparison, but rather a trending characterization. The horizontal velocity peaks at about 80-100 km inland and decreases gradually to the level of 10 m/s at the ice divide. Bartholomew et al. [49] reported an average ice velocity of ~100 m/a extending to 80 km from the ice margin. The GAP study reveals a similar range of supraglacial horizontal velocity of 100-400 m/a. It seems, therefore, that the model proposed in this study can adequately reproduce the thermal and deformational behaviors of the ice sheet at the GAP site.

Surface ice Movement in Horizontal and Vertical Cross-Sections
Maier et al. [46] measured shear strain rates in boreholes at a field site 30 km inland from the Western Margin of the GrIS (67.2N 49.6W). During the 2015/2016 winter, the mean surface velocity was 114.8 m/a. Ice deformation is reported to contribute 4.6 m/a, and basal sliding is responsible for 110.2 m/a, which is termed as plug flow that the basal motion dominates the surface motion. Our modeling shows the surface velocity is half of the above measurement, and basal sliding below 10~15% of ice thickness contributes most of the motion. This behavior is consistent in general with the observation of [46]. The discrepancy is likely attributed to (1) regional-scale velocity variability [47], and (2) the application of the ice-bedrock friction model that uses a larger than real friction coefficient.   The modeling results suggest a gradual transition of ice velocity from the glacier edge inland ( Figure 16). The glacier edge zone where surface melting starts to penetrate to the ice bed moves the fastest, at the level of 100 m/a, which is consistent with the reported annual average value of the GPS measurements. The seasonal variation of temperature enhances the mobility of the ice sheet in summer, which is not considered in our model. Another important observation is the vertical variation of ice mobility, as shown in Figure 17. The lower 10% of the ice sheet contributes to ~60% of the horizontal surface velocity, which is in good agreement with other studies [10,40].  The modeling results suggest a gradual transition of ice velocity from the glacier edge inland (Figure 16). The glacier edge zone where surface melting starts to penetrate to the ice bed moves the fastest, at the level of 100 m/a, which is consistent with the reported annual average value of the GPS measurements. The seasonal variation of temperature enhances the mobility of the ice sheet in summer, which is not considered in our model. Another important observation is the vertical variation of ice mobility, as shown in Figure 17. The lower 10% of the ice sheet contributes to ~60% of the horizontal surface velocity, which is in good agreement with other studies [10,40].

Stress Distribution
The stress contour map as shown in Figure 18 indicates that the mean pressure is dependent on the ice sheet thickness, and the shear stress is mostly negligible (<0.1 MPa) in the ice body but is concentrated (0.2-0.25 Mpa) on the ice bed near the most mobile region where the peak of horizontal velocity takes place. Figure 19 shows the calculated stress components on the ice sheet base. Friction stress grows quickly from zero at the edge to a peak of 125 kPa within tens of km from the edge, and gradually decreases towards the ice divide, with slope changes near the boundary of the melting line. This indicates that the base provides the necessary frictional force that helps stabilize the ice sheet. The supraglacial ice sheet exerts tensional drag force that causes a downward slip, which is reflected by the shear stress that peaks at the location coincident with the peak horizontal velocity. The shear stress is within the reported range of driving stress of the ice sheet

Stress Distribution
The stress contour map as shown in Figure 18 indicates that the mean pressure is dependent on the ice sheet thickness, and the shear stress is mostly negligible (<0.1 MPa) in the ice body but is concentrated (0.2-0.25 Mpa) on the ice bed near the most mobile region where the peak of horizontal velocity takes place. Figure 19 shows the calculated stress components on the ice sheet base. Friction stress grows quickly from zero at the edge to a peak of 125 kPa within tens of km from the edge, and gradually decreases towards the ice divide, with slope changes near the boundary of the melting line. This indicates that the base provides the necessary frictional force that helps stabilize the ice sheet. The supraglacial ice sheet exerts tensional drag force that causes a downward slip, which is reflected by the shear stress that peaks at the location coincident with the peak horizontal velocity. The shear stress is within the reported range of driving stress of the ice sheet for the western coast of Greenland [46]. The stress contour map as shown in Figure 18 indicates that the mean pressure is dependent on the ice sheet thickness, and the shear stress is mostly negligible (<0.1 MPa) in the ice body but is concentrated (0.2-0.25 Mpa) on the ice bed near the most mobile region where the peak of horizontal velocity takes place. Figure 19 shows the calculated stress components on the ice sheet base. Friction stress grows quickly from zero at the edge to a peak of 125 kPa within tens of km from the edge, and gradually decreases towards the ice divide, with slope changes near the boundary of the melting line. This indicates that the base provides the necessary frictional force that helps stabilize the ice sheet. The supraglacial ice sheet exerts tensional drag force that causes a downward slip, which is reflected by the shear stress that peaks at the location coincident with the peak horizontal velocity. The shear stress is within the reported range of driving stress of the ice sheet for the western coast of Greenland [46].

Erosion Rate Estimation
From the ice sheet thermal-mechanical (TM) model we developed and calibrated against various field data, the erosion rate can be calculated from the empirical wear law of Lee and Rutter [29]. One parameter of the law is the volumetric concentration of entrained clast (c). The other required parameter is the rock porosity. We assumed c = 1 as a conservative scenario. The GAP study compiled a collection of borehole rock parameters, including the porosity distribution with depth [4]. The average density of the gneiss bedrock is 3208 kg/m 3 . In general, the porosities of both felsic and mafic rocks range from 0.2 to 0.6%, and average around 0.45%. All boreholes have samples that plot in the range 0.7-0.9%, regardless of depth. Another distinct group of felsic rocks, occurring between 500 and 650 m borehole length in DH-GAP04, have matrix porosities ranging from 0.9% to 1.32%.
The TM model results include the basal slip rate and contact normal stress, which are used to calculate the erosion rate from the Lee and Rutter model [29]. The TM model re-

Erosion Rate Estimation
From the ice sheet thermal-mechanical (TM) model we developed and calibrated against various field data, the erosion rate can be calculated from the empirical wear law of Lee and Rutter [29]. One parameter of the law is the volumetric concentration of entrained clast (c). The other required parameter is the rock porosity. We assumed c = 1 as a conservative scenario. The GAP study compiled a collection of borehole rock parameters, including the porosity distribution with depth [4]. The average density of the gneiss bedrock is 3208 kg/m 3 . In general, the porosities of both felsic and mafic rocks range from 0.2 to 0.6%, and average around 0.45%. All boreholes have samples that plot in the range 0.7-0.9%, regardless of depth. Another distinct group of felsic rocks, occurring between 500 and 650 m borehole length in DH-GAP04, have matrix porosities ranging from 0.9% to 1.32%.
The TM model results include the basal slip rate and contact normal stress, which are used to calculate the erosion rate from the Lee and Rutter model [29]. The TM model results also include velocity and stress distributions in the ice sheet. The spatial distributions of the erosion rates are shown in Figure 20. within 20 km from the ice margin. It is orders of magnitude higher than the predicted abrasion erosion. This observation is consistent with the findings of Cowton et al. [33] and Smith et al. [49], which emphasized the important role of drainage flow on the base erosion of a temperate glacier.  The erosion on the base is found to peak near the ice divide, which is mainly due to the extremely high overburden stress. The flow velocity of ice, even with a friction contact boundary, decreases in magnitude with depth, and reaches a minimum at the bedrock interface. A rapid shift in ice-flow mode occurs at the 10% lower part of the ice sheet, with greatly enhanced shear deformation observed in various studies [10,40]. Our results suggest the same inference, i.e., the majority of the shear strain occurs in the lower 10% of the ice sheet, as shown in Figure 18.

Implications and Limitations
The average erosion rate along the ice sheet flowline can be calculated by dividing the integral of E r by distance as where L is the length of the bed along X coordinate. The average erosion rate for the bed is 0.0057 mm/a. This is in good agreement with the reported value of 0.01 mm/a by Reference [4]. Note that the known enhancement effect of entrained fragments in ice in temperate glaciers that exerts extra normal contact stress against the rock bed, is not considered in the current model. As the basal layer of ice melts against the bedrock under geothermal and frictional heating, it imposes a frictional drag force upon the entrained clast in addition to the overburden stress under self-weight of the overlain ice sheet.
The subglacial erosion resulted from meltwater flow is calculated to be 0.12 mm/a on average, with a peak of 0.9 mm/a at the ice margin. This is the average value over the 80 km flowline of the melted base, with the majority of the erosion occurring in the region within 20 km from the ice margin. It is orders of magnitude higher than the predicted abrasion erosion. This observation is consistent with the findings of Cowton et al. [33] and Smith et al. [49], which emphasized the important role of drainage flow on the base erosion of a temperate glacier.

Implications and Limitations
The developed model addresses the thermal and mechanical processes in a moving glacier. The temperature profiles and horizontal velocity distributions were found to agree with the measurements from the GAP study and literature reports. The application of a frictional contact model between glaciers and bedrock successfully simulated the kink-type rapid contrast in shear rate near the glacier bed. The model correctly simulated a faster sliding rate at higher ice temperature, particularly for regions with a temperature close to the melting point. The simulated bedrock erosion rate was within reported ranges of variation.
The model is based on classical theories of heat transport and viscoelasticity. All model parameters were justified with appropriate range of values reported from the GAP publications or general literature. For instance, the ice-rock abrasion constant and friction coefficient were taken from separate experimental data previously reported. A trial and error procedure was adopted to calibrate the parameters, within the reported range of values, to best fit the observed TM behavior of the ice sheet and bedrock abrasion. It would be worthwhile to perform a systematic uncertainty analysis in the future. One additional limitation of the current model lies in the uncertainties related to the heterogeneous behavior of the subglacial drainage channel. The model is in 2D along the glacier flowline; therefore, the observed regional-scale velocity variability cannot be reproduced [47]. Future efforts will be devoted to developing 3D models with consideration of drainage networks that govern the sediment transport and heat flux regimes.
It is noteworthy that there are some recent development on perspectives about Glen's law parameters for ice sheet that shed light on the nonlinear nature, as well as the dependency of the parameters on temperature, depth, and grain size. This remains to be explored numerically in-depth in the framework of the current TM coupled theories.

Conclusions
In this study, efforts were made to develop a well-calibrated heat transfer model to study the temperature field of the ice sheet in the GAP site. The current model can reproduce the ice core temperature profile consistently. The obtained temperature field, especially on the ice sheet base, provides a solid foundation for the success of a subsequent mechanical model of the ice sheet/bedrock system, where the viscous properties of the ice are temperature-dependent.
The modelling suggests a warm base extending to about 150 km from the edge. The basal heat flux from both geothermal and frictional heating plays a negligible role in the overall temperature distribution. The main contribution to a warm condition is from meltwater percolating from the glacier surface. An extra heat flux of 0.5 W/m 2 on the ice sheet base is sufficient to cause base melting.
The ice velocity is governed by the rheology of the ice sheet, which is driven by self-weight, thermal field, and topographic slope. A thermal-mechanical coupled model was developed, and successfully calibrated to simulate the effect of temperature on ice sheet mobility based on the supraglacial and basal elevations in the GAP study area.
With an inland distance of 500 km and glacier height of 3.0 km, the simulated horizontal ice velocity peaked at the level of 100 m/a at 80~100 km to the edge, in agreement with the GPS measurements. The vertical distribution of horizontal ice velocity suggests that about 60% of deformation takes place in the lower 10% of ice thickness, in agreement with various field observations. The ice-bedrock interface is characterized by an increasing friction force, and decreasing shear stress in the ice body inland towards the ice divide. The peak shear stress in the ice body coincides with the location of the maximum horizontal velocity.
The abrasional erosion rate of the GAP study area is estimated at the level of 0.006 mm/a. This value is in agreement with the finding of a well-cited review paper by Hallet et al. [14], which indicated the same level of erosion for polar glaciers and thin temperate plateau glaciers on crystalline bedrock.
However, further analysis of the meltwater erosion emphasizes its significant contribution to the subglacial erosion at locations within 20 km from the melted ice margin. A model was developed to model the subglacial channel flow, and was calibrated with field data of the catchment discharge. The predicted maximum channel erosion is 0.85 mm/a, and corresponds well to a recent field investigation of the site. It is noteworthy that this channel erosion depends strongly on the flow velocity that concentrates at the estuary outlet. It drops quickly to less than 0.1 mm/a in regions further than 20 km from the edge. The average value is 0.12 mm/a over the 80 km flowline of the studied subglacial drainage system.

Data Availability Statement:
The data presented in this study are available on request from the Canadian Nuclear Safety Commission.