A cross-validated three-dimensional model of an englacial and subglacial drainage system in a High-Arctic glacier

Abstract

Recent speleological surveys of meltwater drainage systems in cold and polythermal glaciers have documented dynamic englacial and in some cases subglacial conduits formed by the ‘cut-and-closure’ mechanism. Investigations of the spatial distribution of such conduits often require a combination of different methods. Here, we studied the englacial drainage system in the cold glacier Longyearbreen, Svalbard by combining speleological exploration of a 478 m long meltwater conduit with a high-resolution ground penetrating radar (GPR) survey with two different centre-frequencies (25 and 100 MHz). The results yielded a 3-D documentation of the present englacial drainage system. The study shows that the overall form of englacial conduits can be detected from velocity−depth converted GPR data, and that the 3-D model can facilitate a method to pinpoint the reflections in a radargram corresponding with the englacial drainage system, although fine detail cannot be resolved. Visible reflections approximately parallel to the mapped englacial water drainage system likely result from sediment incorporated in the ice or from abandoned parts of the englacial drainage system.

1. Introduction

Historically, due to the inaccessibility of the interior of glaciers, fewer studies have been carried out on englacial and subglacial meltwater systems than on supraglacial meltwater systems. However, recently there has been a shift towards a greater focus on englacial hydrology, especially in the context of high magnitude drainage events (Greenwood and others, 2016). Investigations of englacial and subglacial drainage systems have been conducted on small and moderately-sized warm and polythermal glaciers, which has helped to improve the understanding of the mechanisms of surface meltwater drainage into these glaciers (Bingham and others, 2005; Irvine-Fynn and others, 2011) as well as large ice sheets (Zwally and others, 2002; Cowton and others, 2013). Direct methods, including speleological investigations and down-borehole imaging, have allowed a first-hand view of englacial and subglacial drainage systems and documented complex networks of fractures and conduits both in warm and cold ice (e.g. Vatne, 2001; Stuart and others, 2003; Fountain and others, 2005a, 2005b; Gulley and Benn, 2007; Benn and others, 2009; Gulley and others, 2009a; Naegeli and others, 2014).

The increasing amount of evidence from direct exploration does not support the broadly accepted conceptual ‘Shrevian’ model (Shreve, 1972), which has been used as a fundamental concept for predicting water movement in temperate glaciers. The model states that the spatial distribution of meltwater conduits in glaciers is strongly related to theoretical hydraulic gradients and evolves due to the primary permeability of warm ice. However, the primary permeability of ice is extremely low (Lliboutry, 1996; Jordan and Stark, 2001) and, consequently, englacial drainage systems should not be expected to evolve from water-filled vein networks (Gulley and others, 2009b). To address the issue of englacial conduit formation, new theories of routing supraglacial meltwater to the englacial and subglacial system have been proposed, one of them being the ‘cut-and-closure’ mechanism. This was first described by Vatne (2001) and later elaborated by Gulley and others (2009a), who observed that in crevasse-free areas of cold and polythermal glaciers (Longyearbreen, Svalbard and Khumbu Glacier, Himalaya) where sufficient meltwater generates an incision rate greater than the rate of surface ablation, englacial and occasionally subglacial channels can evolve from supraglacial streams that gradually cut through the glacier surface.

The cut-and-closure process results in meandering vadose canyons with prominent sutures marking the lines of closure along the channel roof. The channels typically have long, low-gradient sections separated by steep drop-offs (nickpoints) independent of theoretical hydraulic gradients (in contrast to Shreve, 1972). To better understand the spatial distribution of englacial and subglacial drainage networks within polythermal and cold glaciers, indirect methods such as ground penetrating radar (GPR) surveys have been introduced, involving the transmission of electromagnetic waves in radio frequency (referred to as a radar wave), from an antenna into the subjacent surface (e.g. Stuart and others, 2003; Bælum and Benn, 2011; Zhen and others, 2013; Temminghoff and others, 2018). In these studies, a single velocity, similar to the theoretical values for clean dry ice of 0.168 ± 0.003 m ns⁻¹ (Petterson and others, 2003), has been used to depth-convert the GPR data. Although a great variety of radar wave velocities is known in glacier ice and, especially, in drainage systems containing air- and water-filled conduits (e.g. Moorman and others, 2000; Lee and others, 2010), the method of single-velocity-to-depth-conversion is broadly used.

Stuart and others (2003) and Temminghoff and others (2018) have shown that a combination of GPR survey and direct observation through speleological exploration can validate the interpretations of GPR signals. Stuart and others (2003) describe a way to measure the depth to an englacial channel and the water table within it, by comparison of synthetical modelled travel times and waveforms within air, ice and water and measurement of the relative travel times within radargrams. Furthermore, when a radar wave travels through a medium with homogeneous properties such as cold ice, the amplitude will be reduced due to spherical spreading and attenuation, and the polarity of that wavelet will stay the same. However, when the radar wave encounters a boundary between two media, polarity and phase will change. For example, the reflection at an interface between ice and air will be of normal polarity and that between ice and water will be of reversed polarity. With this method, they determined the dimensions of englacial conduits and the depth of water within conduits in the central part of Austre Brøggerbreen, Svalbard.

A more visual approach was followed by Temminghoff and others (2018), who compared GPR surveys and direct speleological exploration of a conduit system in Scott Turnerbreen, Svalbard, to estimate the location and shape of englacial conduits from composite hyperbolae focused around the same reflector in GPR images.

Extending these approaches, we aim here to validate the GPR technique by documenting the englacial drainage conduits of Longyearbreen that occur in a complex setting where the drainage system is highly sinuous and located close to the glacier margin where debris is abundant. This is done by combining a speleological exploration and a detailed high-resolution radar survey (Hansen, 2016). This approach, involving two independent study methods, is expected to yield the most extensive cross-validated 3-D model of the drainage system in a cold glacier in order to (1) verify the consistency of these two methods, and (2) provide new insights into identifying englacial conduits from single-velocity-to-depth-converted GPR data.

2. Study area

Fieldwork was carried out in the early spring of 2015 on Longyearbreen in central Spitsbergen, Svalbard (78.17°N; 15.46°E). Longyearbreen is located in the valley Longyeardalen. The lithology in and above the valley is sandstone, and siltstone with shale coal seams at the bottom of the valley. The sediment is defined as easily erodible soft and fine-grained sedimentary rock (Etzelmüller and others, 2000). Longyearbreen is a small, mainly uncrevassed, valley glacier with an area of ~2.5 km² (measured in 2010 by Langford and others, 2014). The altitude range is from ~250 to ~1000 m a.s.l. and the glacier flows northward from a single-cirque accumulation basin at the head of the valley (Fig. 1). The glacier is located at the lee side of a mountain slope facing northwest and the terminus is at a wide ice-cored moraine. In 2004, the equilibrium line altitude was estimated at 615 m a.s.l. (Bringedal, 2004; Yde and others, 2008). Studies from 1977 to 1992 showed a negative specific mass balance of −0.55 m a⁻¹ (Hagen and Liestøl, 1990; Hagen and others, 2003). Furthermore, Nuth and others (2010) calculated that the glacier is thinning at a rate of 0.55 m a⁻¹. This is in agreement with observations made by the authors in summer 2014, when net ablation occurred over the entire glacier.

Fig. 1. Location of the study area on Longyearbreen in Svalbard. Entrance to the ice cave is marked with a red circle. Topographic base from the Norwegian Polar Institute (available at: http://toposvalbard.npolar.no/); contour lines in m a.s.l.

The cold thermal regime of Longyearbreen is well documented by several studies (e.g. Etzelmüller and others, 2000; Riger-Kusk, 2006). Using GPR data, Sevestre and others (2015) documented a predominance of transparent cold ice in the glacier. Furthermore, englacial studies of the drainage system on the lower western margin showed a basal temperature of −4 °C measured in 2001, 2002 and 2003 (Humlum and others, 2005). During the ablation season, numerous meltwater channels develop on the surface of Longyearbreen. The largest channels form in the hollows between the main body of the glacier and the ice-cored lateral moraines (Etzelmüller and others, 2000; Gulley and others, 2009a). The glacier cave system studied here is situated close to the steep western mountain slope, and the surrounding frost-sensitive sedimentary rock is slowly incorporated into the glacier ice through rock falls. The glacier ice in this area is therefore expected to contain a certain amount of debris.

3. Methods

3.1 Speleological survey

The glacial cave system of Longyearbreen was surveyed several times in March 2015 and once at the beginning of April 2015. The entrance to the cave is located at 520 m a.s.l. in the upper part of the glacier (Fig. 1). Standard speleological equipment and methods modified for glacier cave surveys (Dasher, 1997; Gulley and Benn, 2007) were used to descend into and map the cave system. The measuring equipment consisted of a Suunto Tandem liquid-filled precision compass and clinometer, and a handheld Leica DISTO D210 laser distance metre.

The survey started at the cave entrance (GPS coordinates 78.17347°N, 15.46945°E). While penetrating the cave, 60 survey stations were set up (numbered A1–A43, B0–B4, C1–C9; Fig. 2) at which measurements were made. Successive stations were placed at distances as large as possible while still in sight from the previous station enabling a direct bearing, plunge and distance reading along a straight line in the down-cave direction. Between the stations, 2-D line-plot geometrical observations were performed to obtain continuous longitudinal planar views (Fig. 3). At selected survey stations, cross-sectional shapes of the cave were determined with an additional reading of bearing, plunge and distance above and below the station. From geometrical calculations, each reading was plotted and then connected on the basis of a sketch of the morphology. The survey produced ten cross-sectional views and covered a total of 478 m of the cave path.

Fig. 2. Longitudinal profile of the ice cave system with every second speleological survey station labelled along the cave path (B4–B0, A1–A38 and C1–C9). Distances from the ice cave to the glacier surface and glacier bed are measured from the individual reflections in radargrams covering the area above the ice cave.

Fig. 3. Plan view of the ice cave system showing its most important characteristics. Ten cross-sections of the in situ mapped cave morphology within the upper cave part are plotted along the cave path. Four of these cross-sections indicated by red lines are aligned with GPR lines that cross the cave path approximately perpendicularly (GRP-1 to GRP-4).

3.2 Ground penetrating radar

In glaciological studies, transmitted radar waves will penetrate the substratum (air, snow, ice and sediment) and will be reflected back by different dielectric properties, which will be received by the antenna in the form of a hyperbola (Plewes and Hubbard, 2001). Because of the transparency of cold ice to radar signals and the strong velocity contrast between ice (0.168 m ns⁻¹) and air (0.300 m ns⁻¹) (Plewes and Hubbard, 2001), GPR surveys have shown great potential for determining the presence of englacial conduits within cold ice (Stuart and others, 2003).

In this study, the GPR system used was a Malå radar system composed of a Malå XV monitor, a ProEx Professional 235 Explorer control unit and a fixed unshielded ‘rough terrain’ antenna with elements oriented in-line. Two different centre frequencies were used, 25 and 100 MHz. The GPR system was towed behind a snowmobile on the glacier surface at speeds below 10 km h⁻¹ and equipped with a high-precision GPS Real Time Kinematic 30/30 mode Trimble SPS855.

The GPR data were collected by driving the GPR system over the glacier surface above the englacial drainage system. The survey covered an area ~400 m long and 250 m wide (Fig. 4). The area was crisscrossed with a spacing of 5–50 m between each line on 11 and 17 April 2015. A grid of 27 lines was driven over twice, first with the centre frequency of 25 MHz and then 100 MHz (Fig. 4). The best horizontal resolution of the ice cave, located close to the glacier surface, is achieved with the 100 MHz centre frequency, and therefore these lines are used to examine the reflection of the ice cave (cf. Reynolds, 1997). With the 25 MHz frequency, small objects will not be resolved as well as with the 100 MHz frequency. However, attenuation is less for the lower-frequency radar waves, hence the reflection from the glacier bed will appear stronger (Daniels and Utsi, 2013). Six of the 27 100 MHz lines, GPR-1 to GPR-6 (Fig. 4), have been selected to compare the speleological survey and the radargrams in detail. Lines with both frequencies are used to investigate and construct a 3-D model of the englacial cave system.

Fig. 4. A close-up view of the part of Longyearbreen investigated showing all GPR lines and the mapped cave path. 100 MHz lines are plotted in red and 25 MHz lines are plotted in dashed green. Radargrams shown and discussed in the paper are labelled GPR-1 to GPR-6. Topographic base from the Norwegian Polar Institute (available at: http://toposvalbard.npolar.no/); contour lines in m a.s.l.

The GPR data have been processed using the software ReflexW version 7.5 (Sandmeier Scientific Software), the main processing steps consisting of an algorithm to suppress low-frequency energy emitted by the field near the transmitter, and an exponential gain function to increase the strength of the signal with depth. Furthermore, a Kirchhoff migration (Hagedoorn, 1954; Kaufman and others, 2002) was performed to rearrange the reflections and focus reflections into their correct spatial locations. In order to get a complete migration and rearrangement of the reflection, all the diffractions need to be orientated orthogonally to the survey direction, which can vary slightly when towing the antennas over the glacier surface. The signals reflected from a small area within the ice cave will not only come from the top and base but also from the sides, so to facilitate the best comparison between the two methods, all radargrams were migrated. Finally, a topographic correction and conversion of the y-axis from nanoseconds to metres was made by assuming a typical constant velocity propagation of radar waves for Svalbard cold-ice glaciers of 0.167 m ns⁻¹ (Björnsson and others, 1996). The effects of each of the main processing steps are displayed in Figure 5.

Fig. 5. The main processing steps on radargram GPR-4. (a) Raw data. (b) Resetting time zero. (c) Dewow and exponential gain function. (d) Dynamic correction and interpolation of traces. (e) Kirchhoff migration. (f) Topographic correction.

3.3 Combined speleological and geophysical approach

Using Petrel software (Schlumberger, USA), the mapped cave path and all processed GPR data, which are geo-referenced in spatial coordinates, were used to construct a 3-D model. Each radargram was adjusted to the depth scale of the y-axis and compared with the calculated maximum depth from the two-way travel time. If the depth diverged by more than ±2 m from the calculated maximum depth, the sample interval was adjusted (within the minimum and maximum range of 0.12–0.2 s). This method was the most efficient way of uploading GPR data to Petrel, which is constructed for seismic data. All radargrams are presented with four times vertical exaggeration to provide a clearer image of the internal structures of Longyearbreen, which is very long in relation to its thickness.

After uploading the radargrams to Petrel software, the mapped ice cave entrance, which is located highest in the investigated area on the glacier surface, was chosen as the point 0 to simplify interpretations of the relative depths in the radargrams. The top reflector of radargram GPR-4, which intersects the cave entrance and is parallel to the glacier's longitudinal direction (Fig. 4), was used as a reference and aligned to be 0 at the entrance to the cave. The additional aligned radargrams and the radargrams transverse to the glacier's longitudinal direction (e.g. southeast to northwest) were then fitted to the top reflector of the radargram GPR-4 (Fig. 6). In this way, each radargram was topographically corrected to the glacier surface downglacier, resulting in a relative height for each profile. The topographic correction for longitudinal radargrams utilizes interpolated geometries of the glacier surface and the glacier bed (Figs 6a, b). These surfaces were generated by locating and picking the glacier surface reflector and the bed reflector in each radargram and then performing a convergent interpolation, resulting in a surface resolution of 50 × 50 m.

Fig. 6. A first-order interpretation of the processed GPR signals through radargrams GPR-4 and GPR-5 (centre frequency of 100 MHz). (a) Radargram GPR-4 located approximately along the direction of the cave path. (b) Radargram GPR-5 located perpendicular to the direction of the cave path. (c) Location of both radargrams with respect to other radargrams nearby. The grey arrow shows the sets of parallel surfaces made by the ‘ringing effect’ of the radar system. The upper part of the ice cave is outlined by a blue line and the descent into the lower part through the main nick point is outlined by a green line.

3.3.1 Methodological uncertainties

The accuracy of the cave mapping and the subsequent spatial models created in Petrel depend on the accuracy of the GPS measurements, cave mapping method and the uploading and picking of reflectors in Petrel.

The accuracy of the GPS system should reach a theoretical accuracy of 30 cm both horizontally and vertically. However, at high latitudes, GPS systems are less precise due to a lower amount of visible satellites and these satellites are at a low angle relative to the horizon. Therefore, the precision of the GPS data was estimated from a stationary GPS logger provided by the Kjell Henriksen Observatory in Longyearbyen. There, the GPS deviations in latitude and longitude were measured during April 2015, and showed a typical deviation of up to 4 m. The errors of the GPS positions are estimated to be between 30 and 400 cm in all directions.

Cave survey measurements were made with a precision of ±5 mm (distance) and ±0.5° (bearing and dip), the latter being equivalent to ±70 mm over the mean survey step of 8 m (Table 1). Because measurement errors are random and uncorrelated, this yields a maximum positional error (at the end of the survey) of ±54 cm. When uploading the GPR line into Petrel, an offset of ±2 m was accepted which may result in an error between the height of the profiles in Petrel and the reference point of the cave system mapped by the GPS point at the top of the ice cave entrance of up to 2 m in the vertical direction. The precision in picking the reflectors from the glacier surface and bed in Petrel depends on the subjective interpretation of the radargrams and on how sharply the bed reflectors appear in the radargram. In areas with multiple bed reflections, the first appearance of a continuous reflector was chosen as an indicator for the glacier bed, as proposed by Irvine-Fynn and others (2006).

Table 1. Spatial parameters of the ice cave (upper and lower parts) determined by in situ speleological mapping

*From survey station C9, three measurements were made into the cave to determine the dimensions of the large oval cave room (Chamber 2).

Altogether, the methodological precision between the GPS, the speleological data and the uploading of each radargram into Petrel results in an error of up to 6.5 m vertically and horizontally.

3.3.2 GPR resolution

Whether a reflective boundary will appear as a single reflector or as a cluster of reflections in the radargram is dependent on the resolution of the radargram. The resolution is the ability to differentiate between two adjacent reflectors and is divided into vertical and horizontal resolution. The vertical resolution depends on the centre frequency and the velocity of the radar wave through the medium. The velocity depends on the radar wavelength in the media under consideration. In ice, it depends on ice temperature, orientation of ice crystals and the amount of water and impurities. In temperate ice, the velocity changes rapidly with water content and can vary from 0.130 to 0.170 m ns⁻¹ (Endres and others, 2009), whereas in cold ice, the velocity is relatively constant at 0.168 m ns⁻¹ (Petterson and others, 2003). Given that Longyearbreen is a cold glacier and with a limited yearly snow layer, the vertical resolution can be estimated using the velocity through cold ice. The vertical resolution will then be 0.56 m with a 100 MHz centre frequency and 2.24 m with a 25 MHz centre frequency.

The horizontal resolution depends on the distance between two reflections on a line perpendicular to the radar wave's direction of propagation, also called the target area. The target area needs to be equivalent to, or exceeding, the area of the first Fresnel zone (Fowler, 2004). The area of the first Fresnel zone is dependent on the propagation velocity in the subsurface (0.167 m ns⁻¹) and the vertical distance to the reflective surface, thus the greater the distance the poorer the spatial resolution. The ice cave is located between 0 and 30 m below the glacier surface, and by measuring the distance from the glacier surface to the top of the ice cave, the horizontal resolution can be estimated. At 5 and 15 m depth, the resolution is 4.5 and 7 m, respectively, with a centre frequency of 100 MHz, and 9.5 and 15 m, respectively, with a centre frequency of 25 MHz.

4. Results

4.1 Speleological observations

The mapped cave representing the englacial drainage system is depicted in a longitudinal profile with every second survey station labelled along the cave path (Fig. 2) and in the plan view (Fig. 3). In addition, Figure 3 shows cross-sections through the cave measured at ten localities.

The ice cave is divided into upper and lower parts, separated by a large vertical step (frozen waterfall) with a drop of 15 m and a liquid water pool below (Fig. 7a), named the main nickpoint (Fig. 3). The geometric parameters of the cave are listed in Table 1. The upper part consists of a single narrow and meandering passage with three minor nickpoints. Apart from the entrance, which is a steep vertical descent (B0–B4), the upper part (A1–A38) is more than 300 m long and has a slope of 3.8°–5.6°. The floor elevation in the upper part is from 0 to −25 m relative to the entrance. The lower part begins at the main nickpoint (Fig. 7a) where the deepest point of the whole cave is located, with an elevation of −40 m relative to the entrance (with ~15 m of ice above the cave ceiling, see Fig. 2). From here, the cave turns to the northwest revealing two large chambers (Chamber 2 is seen in Fig. 7b), connected by a narrow passage and a small chamber in the far western corner where the valley floor was encountered. The elevation increases to −30 m and remains constant through the remaining part of the cave.

Fig. 7. Photographs of features within Longyearbreen ice cave. (a) Main nickpoint with a pool of water at the bottom of Chamber 1 (survey station A40). (b) Chamber 2 at the northern end of the cave (survey stations C4 to C8). Note the big blocks of ice fallen from the roof and the walls covering most of the cave floor. (c) Cross-section through an englacial v-shaped meltwater channel at survey station A13. The channel is filled with imbricated granular sediment. The channel is ~2 m wide and up to ~0.5 m deep. Note the nearly vertical banded ice layers exposed in the cave wall. (d) Section of a debris-filled crevasse in the upper part of the cave between survey stations A13 and A14.

Along the cave path, several structures appear in the ice walls. At survey station A13, a v-shaped channel structure is exposed in the wall of the cave, infilled with imbricated sediment and boulders (Fig. 7c). Between stations A13 and A14, a large body of coarse-grained, angular debris and boulders are situated in the upper wall of the cave (Fig. 7d). A common structure seen in the entire cave is banded ice consisting of light grey to almost black stripes dipping nearly vertically. Icicles, hoar frost and refrozen meltwater together with imbricated rocks in the ice walls and the floor occur in several places (Fig. 3).

Four types of cave morphologies were observed and classified according to the system of Gulley and others (2009a). In the upper part, the most common morphology was plugged canyons (Figs 8a, d). In these canyons, the roofs are high and consist of snow or firn infill. The canyons are narrow and often tilted towards the centre of meander bends. The second most commonly observed morphology was sutured canyons (Figs 8b, c). These passages typically have lower roofs than the plugged canyons. The walls are brought into contact by ice creep closing the roofs and these canyons are not infilled by snow. However, in some cases, the top parts of canyons were inaccessible to direct observation, thus determining the nature of the top closure in these areas was not possible. We classify these roofs as ‘false roofs’ (see partly false roofs in Fig. 8a). The third type of the observed cave morphology was tubular passages, with an approximately elliptical upper cross-section (Fig. 8e). In the lower part, the cave terminates in the fourth morphological type, closed canyons (Fig. 3), which are vertical sutures that had been filled with frozen snow (cf. Gulley and others, 2009a). The closed canyons are located at both ends of the two large chambers (Chamber 1 and 2 in Fig. 3). In the northern corner of Chamber 2, a possible passage continuation is observed below an area of fallen ice blocks, but due to the small size of the passage, it was not possible to explore it any further.

Fig. 8. Sections of the Longyearbreen cave system. (a) Plugged canyon with a partly false snow roof between survey stations A2 and A3. Passage height is 7.5 m and width ~1 m. (b) Sutured canyon between survey stations A34 and A35. Passage height is 6 m and width between 0.5 and 2 m. (c) Plugged canyon between survey stations A6 and A7. Passage height is 4.5 m and width ~1 m. (d) Sutured canyon between survey stations A14 and A15. In the middle of the passage, there is debris frozen to the cave floor. Passage height is 7 m and the width between 0.5 and 3 m. (e) Tubular canyon between survey station A26 and A27. Passage height is 5 m and width between 0.5 and 3 m.

4.2 GPR survey

4.2.1 First-order interpretation

A first-order interpretation of the processed GPR signals is done through the radargrams GPR-4 and GPR-5 (centre frequency of 100 MHz) which are located approximately along and perpendicular to the direction of the cave path, respectively (Fig. 4). Individual features interpreted from the radargrams are illustrated in Figures 6a, b. The interpolated glacier surface and bed are indicated by the thin black dotted lines. Just below the glacier surface, several sets of parallel surfaces (marked by the grey arrow) are seen. The parallel surfaces represent the ‘ringing effect’ of the radar system and are artificial reflections caused by the antenna. From the speleological survey, we know that the upper part of the englacial cave system consists of a single narrow and meandering passage with a low slope angle, and that the lower part begins at the main nickpoint with a 15 m drop. From the radargram, we can see the meandering form of the cave system from the discontinuous groups of reflections, which appear and disappear several times in the top part of the radargram (marked by a thin blue line). The steep slope from the main nickpoint into Chamber 1 (Fig. 3) is clearly discernible in the radargram (marked by a thin green line) and the connection between the supraglacial and subglacial system is visible in the reflections reaching the reflector at the ice/bed interface (the thin black dotted line).

4.2.2 Three-dimensional model of the drainage system

By mapping the glacier surface and bed, we gain an insight into the location of the mapped ice cave relative to the glacier surface and the ice/bed interface. The model generated from the cave surveys and migrated radar data is shown from three different perspectives in Figure 9. In Figure 9a, the cave system is plotted together with radargram GPR-4. The cave is viewed from the southeast, and thus the lower cave part extending to the northwest that includes the two chambers is not visible. In Figure 9b, the cave system is viewed up-glacier towards the southwest against the radargram GPR-6. This perspective shows how the cave turns to the left (northwest) after entering the lower cave through the main nickpoint. The cave floor reaches its maximum depth at the base of the main nickpoint; thereafter, it rises and intersects the interpolated glacier bed in the area where the valley floor was encountered during the cave survey. In Figure 9c, the cave system is viewed from above, and from this angle, the cave path is seen to follow and cross the radargram GPR-2 several times.

Fig. 9. Three-dimensional visualizations of the ice cave system generated in Petrel software. (a) Side-view including radargram GPR-4. (b) Side-view including radargram GPR-6. (c) View from above showing the cave path in black crossing the radargram GPR-2 at eight points. In panels (a) and (b) the cave roof is plotted as a dark purple line and the cave floor as a black line; both lines are derived from in situ mapping. The transparent light-orange surface is the interpolated glacier surface and the blue surface is the interpolated glacier bed. Due to interpolation smoothing, both surfaces appear less rough than they are in reality. Numbers 33, 48 and 63 are metres below the cave entrance.

4.2.3 Observations from radargrams compared with the mapped cross-sections

From the 3-D model, it is possible to further test the correspondence between direct observation and the radar survey. Four close-up views of 100 MHz radargrams (GPR-1 to GPR-4; Figs 10a–d) crossing the cave path perpendicularly are compared with the area where the mapped cave intersects the radargram in the 3-D model (Fig. 10e). All cross-sections are located at a depth of 5–16 m (Table 2). At this depth, the vertical resolution is 0.56 m and the horizontal resolution (which decreases with depth) is 4.5 m at 5 m depth, 6.3 m at 11 m depth and 7.5 m at 16 m depth.

Fig. 10. Four close-up views of 100 MHz radargrams with four times vertical exaggeration and an overview map of the area of each close-up view from the GPR lines which crossing the cave path perpendicularly. (a) Radargram GPR-1 between stations A2 and A3. (b) Radargram GPR-2 between stations A6 and A7. (c) Radargram GPR-3 between stations A14 and A15. (d) Radargram GPR-4 between stations A34 and A35. (e) Overview map with lines GPR-1 to GPR-4 plotted onto the mapped cave path.

Table 2. Radargrams (100 MHz) used for comparison with cave morphology as determined by in situ speleological mapping

In radargrams GPR-1 and GPR-2 (Figs 10a, b), which cross the cave path between stations A2 and A3 (Fig. 8a), and between stations A6 and A7 (Fig. 8c), respectively, the outline of the in situ mapped cave morphology (black outline) is placed where the mapped cave path intersects the radargram in the 3-D model. In Figure 10a, the mapped position of the cave corresponds with a wide region of reflections. In Figure 10b, there is an area of reflections in a horizontal location similar to the in situ mapped cave morphology and position, but the lower part of the cave outline is located below the reflections in an area without any clear reflections. Both areas of reflections are much wider than the in situ mapped cave width. In Figure 10a, continuous reflections are seen below the mapped cave floor, whereas in Figure 10b, there is no reflection below the mapped cave outline. In both radargrams, small areas of reflections are seen below and next to the mapped cave outline.

In both radargrams, there are areas of reflections labelled regions of complex reflections not accessed in the cave survey in Figure 10 (outlined by yellow broken lines). The reflections are within relatively focused areas that continue up close to the glacier surface. In Figures 10a, b, the region of complex reflections is clearly discernible from the reflections located at the cave position.

Radargrams GPR-3 and GPR-4 (Figs 10c, d) cross the cave path between stations A14 and A15 (Fig. 8d), and between stations A34 and A35 (Fig. 8b), respectively. In Figure 10c, there is a small amount of scattering in the upper part of the observed cave position, whereas in the lower part, there is a distinct large area of reflections which appears to continue below the cave floor. In Figure 10d, the observed position of the cave corresponds to a distinct and width region of reflections. In both radargrams, as in Figures 10a, b, distinct and relatively focused areas of complex reflections that continue up close to the glacier surface are outlined by yellow broken lines.

5. Discussion

5.1 Speleological observations

At station A13, the sediment-filled v-shaped channel exposed in the cave wall records a former channel position, oriented approximately perpendicular to the current channel (Fig. 7c). Here we can observe that the current channel has developed independently of the sediment-filled channel. The large body of coarse-grained and angular debris between stations A13 and A14 (Fig. 7d) is interpreted as a debris-filled surface crevasse because of the debris angularity indicating its supraglacial origin (cf. Boulton and Eyles, 1979; Small, 1987; Benn and others, 2003; Benn and Evans, 2010). Here, the cave developed parallel to the crevasse, possibly because of the supraglacial debris that entered the glacier through the crevasse and weakened the ice. This may be similar to the observations from some debris-covered glaciers in Nepal, where large amounts of debris have weakened the ice and influenced the flow direction of englacial meltwater (Gulley and Benn, 2007). These englacial structures observed during the speleological survey indicates a dynamic englacial meltwater system that both cut across and alongside pre-existing englacial structures.

The englacial channels comprise straight, meandering and nickpoint (waterfalls and pools) patterns formed by the ‘cut-and-closure’ mechanism (Gulley and others, 2009a). These geometries are transient in time and space adapting to changes in glacier hydrology and dynamics, in particular water recharge at the ice surface. Therefore, the present geometry of the Longyearbreen drainage system is a snapshot in time, and its future evolution may not be a simple progression of the down-cutting process. Gulley and others (2009a) described how upstream migration of nickpoints can occur at times of increased melting rates and create long, straight passages, and Vatne and Irvin-Fynn (2015) found nickpoints to be the main mechanism of elevation change in the englacial systems of Austre Brøggerbreen – an entirely cold-based glacier twice the size of Longyearbreen. The process of upstream migration of nickpoints has been suggested as the early stage of moulin formation in cold uncrevassed glaciers (Myreng, 2015; Vatne and Irvine-Fynn, 2015). During the speleological survey, the valley floor was encountered at the very end of the cave path after descending from the upper meandering part of the cave from the main nickpoint (Fig. 7a) into the two large chambers in the lower cave (Fig. 3). Whether or not the main nickpoint in Longyearbreen would start to migrate upstream forming a direct and more effective transfer of meltwater from the glacier surface to the englacial and maybe the subglacial regimes is unknown. However, this study shows that there is a direct connection between the glacier surface and the glacier bed, and that there in the northern corner of Chamber 2 (Fig. 3), is a possible passage continuation below an area of fallen ice blocks. This is most likely where the meltwater flows out and may leave a signature in the geomorphological and sedimentary record (cf. Greenwood and others, 2016). If so, glacial landsystems consisting of subglacial channels (e.g. Evans, 2003), traditionally interpreted as evidence of warm-based ice, may require a new scrutiny and critical re-evaluation leaving room for cold-ice regimes to be considered.

The structures inside the ice cave, such as the crevasse-infill and the nearly vertical banded ice, interpreted as annual snow accumulation layers tilted from the originally horizontal position due to differential ice flow rates (Hambrey and others, 1999; Goodsell and others, 2005; Gulley and Benn, 2007), indicates a time characterized by drastically different mass balance and dynamic state than the present. It is clear that speleological exploration can contribute with valuable information, which can help to interpret the previous thermal regime of the glacier.

5.2 Three-dimensional model of the drainage system

In the constructed 3-D model (Fig. 9), the speleological survey and the GPR data are for the first time plotted together making it possible to view the system from different angles and to examine intersections between the cave survey and the individual radargrams. By migrating all radargrams before uploading them into Petrel, the reflections are moved into their theoretical correct spatial locations. Despite the complexity of the englacial system, characterized by high and narrow passages with multiple surfaces, the reflections within the radargrams show clear alignment with the in situ mapped location of the ice cave. In the radargrams orientated along the glacier length, it is possible to recognize the upper and lower parts of the englacial meltwater system, and in the transversely orientated radargrams, the locations of the valley side relative to the ice cave. The fact that the interpolated glacier bed corresponds to the point where the valley floor was encountered during the speleological survey provides further verification of the comparability of the two methods.

5.3 Comparison between the radargrams and the mapped cave cross-sections

The direct cave survey revealed four different cave morphologies: plugged canyons, sutured canyons, tubular passages and closed canyons. When comparing the mapped cave morphologies with radargrams that cross the cave path vertically in that exact area (Figs 10a–d), no difference in cave morphologies is discernible in the reflections of the radargrams. This is not surprising because a GPR survey is an indirect method where the radar waves are reflected from the surfaces inside the ice cave in all directions, which may result in a complex pattern of reflections and scatters. One example is the strong reflections observed below the mapped cave floor in Figure 10a. These reflections cannot be assigned to specific features within the ice cave since the speed of radar waves through air is 0.300 m ns⁻¹ (Plewes and Hubbard, 2001), which is approximately twice the speed in ice (0.167 m ns⁻¹), and the depth conversion by a single velocity of 0.167 m ns⁻¹ will underestimate the propagation velocity through the conduit by a factor of ~2, and therefore overestimate its thickness by a factor of ~2. With this in mind, this method can only show the depth to the top of the conduit and identify the underlying reflections as coming from the base of the conduit, but not the exact base depth of the conduit.

In many areas, the reflections are seen to be much wider than the in situ measured width of the cave (Figs 10a, d), and in other places, there is an absence of reflections in the lower (Fig. 10b) or upper (Fig. 10c) part of the cave outline. In Figure 10b, the offset between the in situ mapped cave roof and the reflections above is <6 m apart in the vertical direction and it may be explained by the estimated error of up to 6.5 m between the different methods.

In Figure 10c, there is an absence of reflections at the top of the in situ mapped cave (black outline) and an area of reflections is particularly strong at the cave base. The reflections could be identified as coming from the actual mapped cave, because it falls under the error of 6.5 m in the vertical direction. Another explanation could be that at this place on the cave path, a large amount of debris was frozen to the ground in the middle of the cave tract between stations A14 and A15 (Fig. 8d). The sediment below Longyearbreen consists of sand and siltstone and shale coal seams which radar waves would travel through with a velocity of 0.13–0.09 and 0.14 m ns⁻¹, respectively (GPR rental supplies, 2019). This would result in a strong reflection due to the velocity contrast between ice and the sediment entrained in the ice. The debris may contribute to the resulting large area of reflections and the depth to which this debris continues below the cave floor is unknown.

Figure 10d is the cross-section farthest away from the cave entrance and with the deepest position in the glacier among the four radargrams. From this radargram, it is estimated that the top of the conduit is located ~12 m below the glacier surface. At this depth, the horizontal resolution is ~6.6 m and the radar wave will therefore not be able to distinguish between two reflective surfaces with a distance <6.6 m. The reflections surrounding the cave outline are wider than 6.6 m. This may be due to the location of the GPR line within a meander bend (Fig. 10e), and the extensive reflections may mirror the meandering bend where the multiple reflections from the surfaces within the bend overlap.

In all four radargrams (Figs 10a–d), there are focused regions of complex reflections (yellow broken lines) that are seen in the area where the mapped cave path intersects the radargrams, but which do not correspond to the mapped cave locations. In Figures 10a, c, the mapped cave path is only intersected once along the close-up view of the radargram, and during the speleological survey, there were no signs of additional branches (Fig. 10e). The areas of complex reflections are located at a distance from the mapped cave outline greater than the estimated error of up to 6.5 m in the horizontal direction. Since the reflections are located parallel to the mapped cave position (20 m in Fig. 10a and 15 m in Fig. 10c), these reflections may indicate abandoned parts of the cave system. Gulley and others (2009a) observed that the deeper parts of the englacial conduits in Longyearbreen can become blocked by a combination of ice accumulation and creep. They theorized that blockage could cause abandonment of deeper parts of the system, and re-routing of the meltwater to higher levels. With re-routing of the meltwater, tabular morphology is formed as the water reaches phreatic stages which were evident during the speleological explorations (Fig. 8e). The regions of complex reflections outlined in yellow may therefore represent such abandoned parts of the englacial meltwater system.

In Figures 10b, d, the GPR line crosses the mapped ice cave through meandering bends to the southwest, but to the northeast of the in situ mapped cave, there is no mapped cave path and the regions of complex reflections outlined in yellow may therefore also represent abandoned parts of the englacial meltwater system.

From observations in the walls of the ice cave and its location adjacent to a steep mountain slope, it is expected that the glacier ice contains a relatively large amount of debris. To determine if a reflector is from a boundary between ice and air rather than ice and sediment analyses of a single trace, polarity can be used (e.g. Stuart and others, 2003; Bælum and Benn, 2011). In this approach, the area above the reflector needs to consist of relatively clean glacier ice to prevent the signal from being disturbed before reaching the reflection of interest, hence changing the interpretation (or identification) of the medium. Stuart and others (2003) used single trace analyses on a radargram in the middle of Austre Brøggerbreen, where the glacier mainly consists of clear ice. However, in our study area, it would be problematic to use this method due to the proximity of the steep mountain slope, and the likelihood of extensive amount of sediment within the glacier ice, such as between stations A13 and A14 (Fig. 7d). The amount of sediment can appear undetected in the radargram, but that can easily change the polarity of the wavelet and the phase of the source wavelet would therefore be unknown.

Areas of reflections, such as the reflections below the observed cave position at 18 m depth in radargram GPR-2 (Fig. 10b), seem too restricted to be part of the englacial meltwater system and may indicate sediment in the ice. However, by making a 3-D model, we avoid the use of single trace analyses for each reflection and characterize the reflectors as coming from the englacial meltwater system by the proximity to the mapped ice cave path.

The comparison of the cave cross-section and the 3-D model shows that the speleological mapping and the GPR survey are to a large extent compatible and thus the cave tract can often be identified from a radargram. The reflections where the cave is located will often be more extensive and more diffuse than the in situ measured morphology of the cave due to the complex spatial structure of the cave. The horizontal resolution will diminish with depth and the area of reflections will increase as more reflections cluster together instead of resolving the different parts of the cave path. It also shows that the method can be used to pinpoint reflections in a radargram corresponding with the speleological exploration and that without a speleological survey, differentiating between reflections from a part of the englacial system, incorporated sediment and regions of complex reflections, possible abandoned part of the cave system, may be very challenging.

To improve the quality and detail of the 3-D model, the next step would be to develop a full 3-D model of the reflections corresponding with the mapped cave path to view the 3-D propagation of the reflections (e.g. Murray and Booth, 2010). A way to compare the mapped cave base with the reflections within the radargram would be to calculate the travel time of radar waves in air (instead of cold ice) in the area where the ice cave is located from the non-depth converted data. This would allow the actual depth of the cave base to be estimated and not to overestimating it. However, because it is not always possible to directly map the englacial meltwater system, the method of single velocity–depth conversion may be the only way to depth convert the data.

Our observations confirm the existence of complex meltwater drainage systems in cold glaciers and that the Shrevian model (Shreve, 1972) cannot account for the drainage system of the entire glacier in all its complexity (cf. Gulley and others, 2009b; Greenwood and others, 2016).

6. Conclusions

By combining direct speleological exploration with GPR survey, we document a dynamic englacial meltwater drainage system in Longyearbreen formed by the ‘cut-and-closure’ mechanism. The 478 m long cave consists of sections with a single narrow and meandering passages with a low slope angle (upper part Fig. 2) and steep waterfall-pool areas ending in large chambers (lower part Fig. 2). The cave extends from the glacier surface to its bed and reaches a maximum depth of 30 m, establishing a hydrologic connection between the ice surface and its bed. The hydrological connection may impact the geomorphological and sedimentary record by eroding or depositing sediment as the meltwater reaches the subglacial environment. Our data contribute to current research on cold glaciers and suggest that cold glacier hydrology is more complex than typically assumed and emphasize the need to abandon the Shrevian model of englacial conduit formation.

Combination of the two independent methods enabled us to generate for the first time a cross-validated 3-D model of the drainage system in a cold glacier. The 3-D model can be studied from every angle and the correspondence between the reflections within the migrated radargrams and the in situ mapped location of the ice cave can be evaluated. The model verifies that the two methods are largely consistent. It shows that it is possible to detect englacial conduits within Longyearbreen using a single-velocity–depth conversion despite the complexity of the ice cave morphology. In addition, the model yields a possibility to pinpoint which reflections in a radargram are to be used when evaluating a GPR survey results.

Several areas of reflections are visible in the radargrams which were not accessed during the cave survey. Some of these cannot be explained by the sinuosity of the cave path, and may result from sediment incorporated in the ice or from abandoned parts of the englacial drainage system. The latter would show that the meltwater system in Longyearbreen has undergone several cycles of incision and abandonment (Gulley and others, 2009b).