Mass-balance measurements

To calibrate a glacier melt model, we conducted mass-balance measurements over parts of the 2003 ablation season. In April 2003, we installed ten ablation stakes distributed over Gruvefonna and along Høganesbreen (Fig. 1) to capture the mass-balance distribution. The stakes are distributed over an altitudinal range of 280-880 m and the dominating aspect classes.

To assess the water equivalent of the snow cover at the beginning of the period over which the mass balance is to be determined, we measured thickness and density of the snowpack at each stake location. In addition, we probed the snow thickness manually in a regular grid of ∼250m spacing across Høganesbreen, and on the Gruvefonna plateau we used Ramac ground-penetrating radar with 500 MHz antennae. The radar data were later resampled to a 200 m resolution. The thickness measurements were spatially interpolated according to the description by Reference Liu, Jezek and LiLiu and others (1999) to construct a map of the initial snow cover. The snow thickness was converted to water equivalent using a mean density of 328 kg m–³. The standard deviation of the density measurements from this mean value is <12% and we estimate the total error associated with the map to be in the same order.

In August 2003, the winter snow cover had completely disappeared at all stakes and we measured the stake emergences in two instances on 9 and 10 August and 17 August.

Automatic weather station

Close to the stake at 480 m a.s.l., we operated an automatic weather station (AWS) to measure air temperature and surface lowering. We used an Aanderaa 2775 air-temperature sensor mounted in a radiation screen and installed it initially 2 m above the ice surface. An ultrasonic rangefinder ULS 3600 Aanderaa was installed to record a high-resolution time series of surface lowering. During the 2003 melt season the glacier surface had melted back such that the rangefinder was >3 m away from the glacier surface. The instrument was subsequently set to a lower position to reduce its vibrations with the wind and thus attain more stable readings. The station was operated until the end of the ablation season in November 2003 at intervals of 10 min.

To analyze the temperature variations with elevation, a second temperature sensor was installed at the top of Gruvefonna, at 880 m a.s.l. The sensor and instrument set-up was identical to that at the lower station, and the summit station was operated from 7 April to 10 August 2003 and yielded an average lapse rate of −0.0068Km–¹.

Discharge measurements in the mine

Water intrudes at many different places into the mine and flows in small creeks towards a pool of ∼ 1000m³ at the lowest part of the mine. Pumps conducting the water to the surface are triggered when the water level in this pool exceeds a critical value. The effective pump rate was measured at intervals of 10min using two ultrasonic Doppler flow meters (Endress prosonic flow 93 P) that were attached to the pipes through which the water was conducted to the surface. The producer quotes the accuracy of these measurements to be better than 5%.

Model description

We have applied a distributed temperature-index method developed by Reference HockHock (1999) that includes potential direct solar radiation. The latter accounts for shading of the glacier due to surrounding topography, and the effects of local slope and aspect of the glacier surface. The diurnal melt rate M (mm d^{- 1} ) is then calculated as

(1)

where MF is a melt factor (mmd–^1˚C–¹), a _snow/_ice is a radiation coefficient different for snow and ice surfaces (m²W–¹ mm d^-1 °C–¹) and T is air temperature (°C). The potential direct solar radiation at the glacier surface, I (Wm^–2) , is calculated as a function of solar geometry and topographic characteristics (Reference HockHock, 1999). In several different applications, the model formulation was able to simulate accurately meltwater production (e.g. Reference HockHock, 1999; Reference Flowers and ClarkeFlowers and Clarke, 2002; Reference Schuler, Fischer, Sterr, Hock and GudmundssonSchuler and others, 2002).

Model operation

Melt rates at the glacier surface were calculated using Equation (1) for each glacier gridcell of a digital elevation model (DEM). To this end, we used a 50 m resolution DEM of the Gruvefonna region as the computational domain. This DEM was constructed from 1990 aerial photographs. The outlines of the glacier-covered regions were derived from the same aerial photographs. The computation of potential clear-sky solar radiation requires further digital maps of slope and aspect angles, which were derived directly from the DEM. In addition, the effect of topographic shading at each gridcell was determined at hourly time-steps based on the effective horizon and the position of the sun. A further raster map contains the distribution of initial snow water equivalent and served as a starting situation for the calculations. Air temperature is required as input data for the model (Fig. 2b). Precipitation is also considered because the period over which the model is run includes times of snow accumulation, and also because it accounts for the contribution of rainfall to total water input. For this purpose, we used precipitation data from the meteorological station in Longyearbyen, approximately 37 km north of Gruvefonna (Fig. 2a).

Fig. 2.

Daily values of (a) precipitation, P, in Longyearbyen, (b) air temperature, T, at the AWS on Høganesbreen and (c) melt rate, M, averaged over the entire model domain.

Typically, model runs were performed for the period of the temperature record, from 8 April 2003 to 4 November 2003 when the recording of pump rates was stopped. The model was run with daily time-steps since precipitation data were available only at this interval. The air temperature recorded at the AWS was extrapolated and distributed to each gridcell using the linear lapse rate of –0.0068 Km^–1 that was derived from measurements. Precipitation was distributed spatially using a gradient with elevation for which we adopted a value of +15%(100m)^–1 in agreement with measurements of Reference Sand, Winther, Maréchal, Bruland and MelvoldSand and others (2003). A threshold temperature of 1˚C was used to decide whether the precipitation was snow or rain, and a correction factor was applied to account for undercatch of the precipitation gauge. The calculated melt rate is shown in Figure 2c.

The melt-factor and radiation coefficients for snow and ice are empirical coefficients, the values of which were adjusted iteratively to optimize the agreement between model results and mass-balance measurements. In doing so, attention was paid to a number of different criteria. Mainly, we aimed to minimize the deviation of modelled from observed values at individual stakes (Fig. 3a). An additional condition was to reproduce the variation of mass balance with elevation (Fig. 3b), and a third target was to replicate the high-resolution time series of ice ablation that was provided by the ultrasonic rangefinder after the winter snow cover disappeared (Fig. 3c). Calibrated parameter values are presented in Table 1. Compared to former applications of this model (Reference HockHock, 1999; Reference Schuler, Fischer, Sterr, Hock and GudmundssonSchuler and others, 2002), we find that the melt factor and the radiation factors are approximately twice as large, probably reflecting the persistent summer melting in the High Arctic.

Fig. 3.

Model performance. (a) Measured vs modeled mass balance at individual ablation stakes. (b) Mass balance vs elevation; squares represent values of the period 8 April–10 August 2003, and circles those of 8 April–17 August 2003. (c) Time series of cumulative melt at the weather station site.

Table 1.

Adjusted model parameters and their optimized values

Meltwater production

In general, the simulated melt corresponds fairly well to the measured melt (Fig. 3), and evaluating the model performance we find a discrepancy of ∼19%. This value is comparable to accuracy estimates in previous studies (e.g. Reference Schuler, Fischer, Sterr, Hock and GudmundssonSchuler and others, 2002). However, some deviations are recognized from Figure 3a, where the points off the 1:1 line at low melt rates indicate that the model fails to explain the observed variability on the plateau of Gruvefonna. We associate this low variability of our results with the weakly pronounced topography on the plateau and the nearly homogeneous distribution of solar radiation, two factors that control melt distribution in our model formulation. In contrast, the natural distribution of melt is strongly influenced by snow redistribution by wind, a process that is not accounted for by our simple model. The underestimation of melt at the stakes located between 400 and 600ma.s.l. and at the top of Gruvefonna (Fig. 3b) can be explained by the same mechanism. Erosion of snow due to higher wind speed at the summit and at the entrance into the narrow valley that confines Høganesbreen may lead to an earlier exposure of glacier ice and thus induce increased melting. The good match of measured and calculated time series of melt shown in Figure 3c demonstrates that apart from uncertainties due to the complex distribution of snow, the melt model performs well. The total meltwater production in the period 8 April–4 November 2003 varies from 1000 mm on Gruvefonna to ∼3000mm further down along Høganesbreen at lower altitudes.

The time series of meltwater production (Fig. 2c) in the entire model domain demonstrates that melt in early June occurred sporadically. Continuous melt production is marked by a sudden onset in late June. The main melt season lasted throughout July and large parts of August, and shows large melt variability, presumably reflecting local weather patterns. Melt production during that period reached peak values of 2.5 × 10⁶ m³ d–¹ and rarely declined below 1 × 1 0 ⁶m³d–¹ . After 20 August, melt production dropped suddenly and was almost terminated. A last pronounced melt event occurred at the beginning of September before melting at the glacier surface ceased.

Input–output analysis

To analyze the transfer of water from the glacier surface to the mine, we compare the total volumes of water production (including rain) and of pumped water out of the mine. However, the exact extent of the surface area that contributes to the mine discharge is a priori unknown; therefore we have estimated the contributing area considering two extreme situations. For a minimum contributing area, we assume that only the glacier surface located directly above the mining area contributes to mine discharge (scenario A). For a maximum contributing area, we assume that the entire surface located above the subglacial catchment drains to the mine (scenario B).

The subglacial catchment area was determined based on the distribution of hydraulic potential at the glacier bed as described above. The calibrated melt model was applied to each of these two surface areas.

The simulated, total melt production in the area considered in scenario A amounts to 1.69 × 106m³ (±20%), most of it originating from melt. The contribution of rain is <2%. In comparison, the total discharge volume pumped from the mine was about 2.83 × 10⁶ m³ (±5%) and exceeds the total water production on the glacier surface by 67%. This implies that the mine discharge receives contributions from a significantly larger surface area than in scenario A.

The subglacial catchment area of the current mine is approximately 1.7 times larger than the mining area itself, and the calculated runoff from the glacier surface directly above this area (scenario B) amounts to 2.85 × 1 0 ⁶m³ (±20%). Again, the bulk of the water stems from glacier melt, and the contribution of rain is 3.5%. The excellent agreement between calculated and measured water volumes indicates that the approach in scenario B is appropriate for estimating the catchment area of the mine, and we apply this procedure in our subsequent considerations.

Figure 5 displays the calculated meltwater flux from the surface directly above the subglacial catchment, and the measured mine discharge. The mine discharge shows a sudden onset similar to the melt production but with a delay of a few days. It appears that the variations of melt production during July and August are much larger than those of the mine discharge. After surface melt drops to low levels around 20 August, mine discharge continues at a high level for a few days before starting a gentle, exponential decline. This decline lasts approximately until the beginning of November when the logging of pump rates was stopped, while melt production had already ceased in September. This behaviour can be described using a linear reservoir approach (e.g. Reference Chow, Maidment and MaysChow and others, 1988). The method describes routing of water through a system, the outflow of which is proportional to the stored water volume. At the time nΔt, the response of a linear reservoir to a sequence of rectangular pulses is expressed by

(3)

Fig. 5.

The melt- and rainwater flux calculated for the catchment area of the actual mine (grey line), and the discharge in the mine simulated by the linear reservoir model (dotted line) and derived from pumping records (black solid line).

where I and Q denote the input and output discharge, k is a storage constant and Δt is the time-step of the calculation (e.g. Reference Chow, Maidment and MaysChow and others, 1988).

Applying this model, we have calibrated the storage constant and found that a value of 14 days yields minimum discrepancy between the modelled output and the observed mine discharge (Fig. 5). This simple model simulates the discharge in the mine surprisingly well (r ² = 0.95) considering that it describes in reality a complex system of water transfer at the glacier surface, through the glacier body and through the cracked bedrock to the mine. The discrepancies between model predictions and discharge measurements presumably originate from the simple structure of the linear reservoir model, which cannot account for the entire complexity of the considered system. Due to its overall good performance, we have selected this approach to produce scenarios of future water intrusions into a progressively enlarged mine.

Implications for future evolution

To predict future water intrusions into a mine that is progressively enlarged, we applied the coupled melt and transfer models to several stages of the planned extension. In detail, we have mapped the subglacial catchment area for each of the planned panels and applied the melt model to the corresponding glacier surface. The melt model was driven by the meteorological data available for 2003. Finally, the calculated melt- and rainwater fluxes served as input to the linear reservoir model to predict the rate of water intrusions into the mine for each considered step of enlargement.

The model results for eight steps of enlargement are presented in Figure 6. Since we have used the same meteorological data in all calculations, and the water intrusion is largely controlled by glacier melt, the shape of the hydrograph as simulated for the current mining area is preserved in the future scenarios. However, the predicted hydrographs appear vertically stretched relative to the reference curve, since the total discharge volume increases considerably with enlargement of the mine. Our predictions suggest that the total water volume will increase from 2.8 × 1 0 ⁶m³ (±20%) to about 9 . 8 × 1 0 ⁶m³ (±20%) (an increase of 350%) when the mining area is enlarged from the current to full size. The maximum diurnal discharge will increase from ∼54 × 10³ m³ d–¹ (±20%) in 2003 to roughly 190 × 10³ m³ d^–1 (±20%) as the maximum extent is reached (an increase of 350%).

Fig. 6.

The hydrograph of discharge in the mine as predicted for the current mine (2003) and eight steps of enlargement. The hydrograph is based on the meteorological input data from 2003.