Introduction

The Greenland ice sheet is experiencing mass loss at its margins (Reference ZwallyZwally and others, 2005; Reference Chen, Wilson and TapleyChen and others, 2006) due to the acceleration of outlet glaciers (Reference Joughin, Abdalati and FahnestockJoughin and others, 2004; Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006) and increased surface melting (Reference KrabillKrabill and others, 2000, Reference Krabill2004; Reference MoteMote, 2007). Rates of peripheral mass loss are much higher than can be explained from ablation alone, leading Reference KrabillKrabill and others (1999) to suggest that the imbalance is dynamic, possibly related to decreased basal friction as a result of meltwater drainage to the bed of the ice sheet (Reference AndreasenAndreasen, 1985; Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Joughin, Das, King, Smith, Howat and MoonJoughin and others, 2008). Observations by Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others (2002) suggest that seasonal acceleration occurs much further inland than previously expected, possibly originating at the equilibrium line near Swiss Camp (69.57° N, 49.33° W; Fig. 1). That work has stimulated new interest in englacial drainage systems, to understand how meltwater reaching the base of the ice sheet might affect the ice-sheet interior as the climate warms (e.g. Reference Parizek and AlleyParizek and Alley, 2004; Reference Alley, Dupont, Parizek and AnandakrishnanAlley and others, 2005; Reference Van der VeenVan der Veen, 2007; Reference DasDas and others, 2008). An alternative explanation for melt-induced acceleration is given by Reference Price, Payne, Catania and NeumannPrice and others (2008), who show that the velocity variations at Swiss Camp observed by Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others (2002) may also be explained by longitudinal coupling, due to melt-induced accelerations occurring closer to the ice-sheet margin.

Fig. 1. Landsat ETM (Enhanced Thematic Mapper) image from 1 August 2001, showing the locations of vertically stacked hyperbolae (VSH) patterns observed in radar profiles and described in the text (circles); weather stations maintained by Greenland Climate Network (GC-Net) (pink squares); approximate borehole locations drilled by Reference Thomsen, Olesen, Braithwaite and BøggildThomsen and others (1991) (stars); and radar profiles described in the text (green lines with labeled end-points). Circle diameter indicates the off-axis distance to the topmost hyperbola and has been doubled relative to the map scale for visibility. Surface-elevation contours are from Reference Bamber, Ekholm and KrabillBamber and others (2001). Inset shows location of study area.

Critical to the understanding of how melt-induced acceleration affects ice-sheet mass balance is the extent to which surface meltwater can quickly, and locally, access the base of the ice sheet through thick, cold ice. Englacial drainage is common on temperate glaciers and has recently been observed on polythermal glaciers despite the existence of a cold mantel of ice at the surface which might prevent supra-glacial water penetration (Reference Boon and SharpBoon and Sharp, 2003; Reference Bingham, Hubbard, Nienow and SharpBingham and others, 2008). In Greenland, surface meltwater penetration to the bed is thought to be more difficult because the ice is thicker and colder (∼1000 m; <0°C). At Swiss Camp, the ice is ∼1200 m thick and theory suggests that local connection of surface meltwater to the bed through water-filled crevasse propagation (or moulin formation) should be difficult; to offset refreezing and creep closure at depth, a crevasse must fill with meltwater at rates ∼1 m h⁻¹ (Reference Van der VeenVan der Veen, 2007).

Moulins and englacial conduits have been reported at and within ∼10 km of the ice-sheet margin in the Swiss Camp region (personal communication from L. Moreau, 2007; personal communication from A. Behar, 2007) but they have not been observed inland of the equilibrium line. Recent observations of Reference DasDas and others (2008) indicate that surface meltwater can penetrate thick ice (980 m) via lake drainage into an existing crevasse. In that case, an existing crevasse was enlarged to a moulin by a large meltwater supply, produced in the warm summer of 2007, at a location where ice temperature may be closer to the pressure-melting point (their study region is ∼20 km downstream of the equilibrium line).

An improved understanding of melt-induced acceleration requires characterization of the path water takes to the bed. This includes how the englacial drainage system may evolve over the summer melt season as well as an understanding of how quickly surface water can reach the bed. Ice-penetrating radar has contributed to the understanding of the hydrology (englacial and subglacial) of temperate valley glaciers (e.g. Reference Walford, Kennett and HolmlundWalford and others, 1986; Reference Arcone and YankielunArcone and Yankielun, 2000; Reference Murray, Stuart, Fry, Gamble and CrabtreeMurray and others, 2000; Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003; Reference Fountain, Jacobel, Schlichting and JanssonFountain and others, 2005) and polythermal glaciers (e.g. Reference Boon and SharpBoon and Sharp, 2003; Reference Bingham, Hubbard, Nienow and SharpBingham and others, 2008). We build on this work by examining ground-based, ice-penetrating radar data from near the equilibrium zone on the western margin of the Greenland ice sheet (near Swiss Camp; Fig. 1). Radar data are analyzed in order to characterize features of the englacial drainage system and to assess the evidence for connections between the surface and subglacial water systems at or near the equilibrium line.

Data Acquisition

We use a custom-built, low-frequency, short-pulse, ground-based radar system to image deep internal layers (∼100–1500 m) and the ice–bed interface. Our transmitter is designed by Kentech Instruments and produces a ±2 kV output pulse with a repetition frequency of 1 kHz. The center frequency of the radar system is determined by the length of the resistively loaded dipolar antennas. The data presented here were acquired using 1 MHz antennas with 40 m half-lengths, creating a pulse wavelength of ∼100 m. The large wavelength permits penetration through the entire ice column but prevents resolution of small-scale englacial features. Transmitter and receiver sleds were separated by ∼120 m. To improve the signal-to-noise ratio, our data consist of several hundred (∼800) stacked waveforms recorded over a horizontal spacing of ∼10 m. As a result, steeply dipping internal layers may be aliased.

Additional processing includes reduction of noise and instrumentation artifacts (through removal of the mean waveform from each radar profile), bandpass filtering and conversion of the two-way travel time to a depth coordinate. We assume that the wave speed in ice is 168 m μs^{− 1} and that higher wave speeds in the surface snow cover are negligible since the snow depth is only a few meters throughout our study area (snow cover varied between 0.1 and 5 m thick throughout the ablation region). A global positioning system (GPS) receiver is mounted on our radar system and allows latitude, longitude and elevation to be calculated at points defining each averaged radar waveform.

Over 250 km of radar data were acquired during our spring 2006 field campaign (Fig. 1). During this time period, no signs of significant surface melting were visible. In addition, we did not observe direct evidence of relict supraglacial streams beneath the snowpack. We did, however, observe several possible supraglacial lake basins (as indicated by relict shorelines and remnant icebergs remaining after lake drainage). Our radar profiles include three along-flow profiles that span the ice-sheet equilibrium zone, separating accumulation and ablation zones, and three across-flow profiles at different elevations within the accumulation, equilibrium and ablation zones. A simple two-dimensional migration routine was applied to the data but gave poor results, indicating that a large number of observed diffractors are from the radar-profile axis. For this reason, we present unmigrated data below.

Radar Observations: Stratigraphy

The bed is topographically rough throughout the region, and ice thickness decreases along each of our radar lines from ∼1300 m in the accumulation zone to ∼600 m in the ablation zone (e.g. Fig. 2). We crossed the equilibrium line approximately midway through each along-flow profile (elevation ∼1200 m; Fig. 1). Internal layers are continuous upstream of the equilibrium line and appear consistent with flow over rough basal topography (e.g. Fig. 2). Downstream of the equilibrium line, however, internal layers become deformed in places and are cross-cut by numerous hyperbolic diffractors.

Fig. 2. 1 MHz radar profile through Swiss Camp (at x = 0 km) from JAR1 (left) to CMPB (right) in Figure 1. Ice-flow direction is indicated by the arrow. The diffraction patterns observed here are described in the text.

Two dominant diffraction patterns are observed: (1) simple diffractors and (2) complicated diffractors that appear to be stacked vertically. Simple diffractors are comprised of only one or two hyperbolic reflections (e.g. Fig. 2 at 10 km and Fig. 3b at 31–35 km). We interpret these diffractors to be the result of surface or near-surface crevasses located either below or off-axis from the radar profile. We expect that most of the crevasses in the ablation region fill with available surface meltwater in the summer and, if they cannot drain englacially, the water pooled in these crevasses refreezes each winter. Direct evidence of this was seen in the field; large (1–2 m wide, up to 1 m tall) ‘humps’ of bubble-free ice were seen on the ice-sheet surface throughout the ablation region, and the locations of these features correlate well with simple-type diffraction patterns seen in the radar data.

The second type of diffraction pattern, which we refer to as vertically stacked hyperbolae (VSH), is more unusual, but is present throughout the ablation zone. This pattern consists of multiple diffraction hyperbolae, up to 2 km wide, which appear to be stacked vertically or near-vertically (several examples are shown in Fig. 3). A total of 33 VSH patterns are identified from our radar data (Fig. 1). Although it is certain that more of these features exist outside the areas covered by our radar profiles, we observe none in the accumulation region covered by radar profiles. While it remains possible that some VSH patterns could be observed upstream of the equilibrium line, they are much more frequent throughout the ablation region.

Fig. 3. 1 MHz radar profiles across several examples of VSH features (indicated by white arrows). Ice-flow direction is indicated at the bottom left of each profile. Irregular hyperbolic returns between the two VSH patterns in (d) are due to noise from instruments located at Swiss Camp. Additionally, one of the VSH patterns in (d) is due to a relict borehole drilled to 600 m in 1990.

A survey of the VSH patterns seen in our dataset allows us to compile a list of characteristics that are common among some or all of the features:

1. 1. all VSH patterns consist of a suite of bright, asymmetric hyperbolae (Fig. 3);

2. 2. some VSH patterns contain hyperbolae that diminish in amplitude with increasing ice thickness (e.g. Fig. 3a, c, d and g), while others have nearly constant hyperbolic amplitude with depth (e.g. Fig. 3b, e and f). In part, diminishing amplitude can result from attenuation of the radar signal in thicker ice; however, we observe VSH patterns in similar ice thickness with different attenuation characteristics (e.g. Fig. 3b compared with Fig. 3g);

3. 3. in some cases, diffraction hyperbolae contained in the VSH patterns are observed below the reflection of the ice–bed interface and often maintain large amplitudes with increasing ice thickness (e.g. Fig. 3b);

4. 4. some VSH patterns are coincident with downwarped internal layers (Fig. 3b and e). In both of these cases, the amplitude of internal layer downwarping increases with depth and the deepest layers appear to truncate at the bed.

While ice-penetrating radar has been used to examine the structure of englacial conduits (e.g. Reference Murray, Stuart, Fry, Gamble and CrabtreeMurray and others, 2000; Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003), diffraction patterns similar to the VSH patterns have rarely been observed, with two exceptions. Reference Raymond, Catania, Nereson and van der VeenRaymond and others (2006) observed a similar pattern due to scattering from a radio antenna erected on the ice-sheet surface in a remote camp (their fig. 2b). Reference Neumann, Conway, Waddington, Catania and MorseNeumann and others (2008) observed a similar pattern in radar profiles across the site of the Byrd ice core (their fig. 2a). While both these examples are related to man-made structures, they offer instructive information about the features causing the VSH pattern; in both these cases, diffraction is from an object that is much longer than it is wide with respect to the radar pulse length (∼100 m). In addition, these two observations suggest two possible orientations for the objects causing the VSH patterns: either horizontal along or normal to the ice-sheet surface. In the latter case, a VSH pattern would be produced when profiling in close proximity to the object, crossing it in any direction. In the former case, profiling along the feature would produce a different diffraction pattern than profiling across it.

As the radar system transmits energy in all directions, we can assume that many diffraction hyperbolae in the data are due to features located some distance off-axis from our radar profile. Further, all VSH patterns are located in the ablation region, so we assume that they are not due to buried structures (since they would eventually melt out to the surface). We measure the travel time to the topmost hyperbola and use the radio-wave velocity in ice (assuming there is no firn present) to determine the off-axis distance to each of the VSH patterns. We map increasing distance using circles with increasing diameters (Fig. 1); large circles indicate that our radar line crossed the feature far off-axis, and small circles indicate that we crossed very close to the feature. We detect VSH patterns up to ∼400 m off-axis from our profiles.

There are a few man-made structures in our study area that may appear in our radar profiles and thus require a non-glaciological interpretation. The approximate locations of several boreholes drilled near to, and downstream of, Swiss Camp in 1990 by Reference Thomsen, Olesen, Braithwaite and BøggildThomsen and others (1991) are shown in Figure 1 (yellow stars). These holes were drilled to depths of up to 600 m and were used to install thermistor strings to measure englacial temperature. Only one of these boreholes lies close enough to our radar-profile locations to appear in our radar data (Fig. 3d, at 4 km). In addition to this, several automated weather stations exist in the region, which might create diffraction patterns similar to the VSH patterns. These stations are part of the Greenland Climate Network (Reference Steffen and BoxSteffen and Box, 2001), and two of them (JAR1 and Swiss Camp) are in our study area (Fig. 1). While we tended to avoid profiling close to these stations, we believe that noise from these instruments can be seen in profiles near Swiss Camp (Fig. 3d, at 2.2–2.5 km). Other than these examples (which are not included in the analysis below), the VSH patterns observed in these profiles are assumed to be glaciological in origin. Our hypothesis is that these patterns represent drainage features, possibly moulins, which allow surface meltwater access to the subglacial drainage system.

It is difficult to determine the true orientation of the objects causing the VSH patterns because horizontally oriented surface features (e.g. supraglacial stream channels) and englacial conduits may appear similar to moulins in radar records (Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003). Despite this, we argue that most of the VSH patterns are due to vertical conduits; radar observations across and along englacial channels would show very different diffraction patterns (Reference Arcone and YankielunArcone and Yankielun, 2000; Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003) and several of the VSH patterns examined here were crossed by multiple profiles at nearly right angles, yet all profiles show the same VSH pattern.

Radar Simulation

In order to understand the origin of the VSH patterns and their variations (listed above), and to test the hypothesis that they are related to moulins, we use a publicly available, two-dimensional, electromagnetic waveform model, ‘GprMax2D’ (Reference GiannopoulosGiannopoulos, 2005). The model uses the finite-difference time-domain approach to the numerical solution of Maxwell’s equations and makes several assumptions including:

1. 1. all media are linear (dielectric properties are independent of the propagating electromagnetic field) and isotropic;

2. 2. the transmitting antenna is modeled as a line source;

3. 3. the constitutive parameters are frequency-independent;

4. 4. the boundaries of the model domain are regions where any waves impinging on them are absorbed, thus allowing for truncation of the computational domain.

Media in the model are described by their dielectric constant (relative permittivity) and electrical conductivity. Our model domain consists of three layers: a 4 km wide by 5 m thick layer of air overlying a 4 km wide by 1195 m thick layer of ice (with constant electrical properties throughout, i.e. no firn/snow cover), which overlies a 4 km wide by 200 m thick layer of bedrock (Fig. 4a). Typical values are taken for cold ice, fresh water, bedrock and air (Table 1; Reference HobbsHobbs, 1974; Reference GriffithsGriffiths, 1999). We use a 2 MHz line-source radar signal and separate the source and receiver points by a constant interval, equivalent to the separation between our actual radar system transmitter and receiver (120 m). The source and receiver points are advanced across the domain (on the surface of the ice) at an increment of 10 m. We record simulated waveforms at each increment for travel times up to 20 μs. We use a 2 MHz radar signal rather than a 1 MHz signal (as in the radar data) to better resolve near-surface features and to increase the resolution of modeled hyperbolae. The model domain is discretized by 2 m by 2 m gridcells except for some cases where a greater horizontal resolution, of 0.5 m by 2 m, was required.

Table 1. Typical electrical properties for materials used in the models (Reference HobbsHobbs, 1974; Reference GriffithsGriffiths, 1999; Reference Miners, Wolff, Moore, Jacobel and HempelMiners and others, 2002)

Fig. 4. (a) Geometry of layered model in the region of a 2 m wide, vertical conduit centered in the model domain with horizontal internal layers. (b–d) Simulated waveforms from (b) a water-filled conduit; (c) an air-filled conduit; and (d) an ice-filled conduit. (e) Model run with internal layers but no vertical conduit. (f–h) Difference between (e) and water-filled (f), air-filled (g) and ice-filled (h) cases.

Our goal is to replicate some of the distinct characteristics of the VSH patterns listed earlier, in order to understand their origin. We insert a 2 m wide by 1000 m long vertical rectangle into the center of the model domain to simulate a smooth-walled, vertical conduit that penetrates the entire ice thickness (i.e. a moulin). We acknowledge that this idealization does not reflect the complicated geometry of actual moulins, as observed by investigators in the field (e.g. Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003). Our approach is to begin simply and to increase moulin complexity only when necessary. Our initial model results using this geometry (not shown) produced only two diffraction hyperbolae: one at the top of the ice, where the conduit opens to the free surface, and one at the base of the ice, where the conduit connects to the bed. The lack of hyperbolae within the ice in this simulation is due to the absence of point reflectors within the ice sheet, which are created when the conduit truncates other objects in the ice column (such as internal layers). To simulate this, we introduce nine 2 m thick horizontal internal layers across the model domain at various depths (Fig. 4a). These layers have electrical properties approximating a mixture of volcanic ash and ice (Table 1; Reference Adams, Perger, Rose and KostinskiAdams and others, 1996; Reference Miners, Wolff, Moore, Jacobel and HempelMiners and others, 2002).

Simulations with different conduit shape

We also ran simulations to investigate how conduit shape affects the existence and amplitude of hyperbolae. First, we modeled a wide (30 m), ice-filled conduit (Fig. 5b) and compared the results to those from the narrow (2 m), ice-filled case described earlier (Fig. 4d). While we do not expect moulins to be this wide, this simulation aids in understanding the influence of changing conduit width on hyperbolic amplitude. We find that increasing conduit width also increases the amplitude of hyperbolic returns and we expect that this also holds true for moulins filled with different materials. Brighter hyperbolae are possible for wider conduits because their size approaches the horizontal resolution of the 2 MHz signal, which is 1/4λ or 37.5 m (Reference YilmazYilmaz, 1994). For this reason, we suggest that VSH patterns with very bright hyperbolae can also result from relatively wide conduits or from conduits containing large voids.

Fig. 5. (a) Model geometry and (b) results for a 30 m wide, ice-filled conduit centered at 2000 m. (c) Model geometry and (d) results for an irregularly shaped, water-filled conduit in ice with no internal layers.

Next, we investigate the effect of changing conduit shape on hyperbolic amplitude. The geometry used in the simulation consists of an irregularly shaped conduit (comprised of several rectangular elements) that narrows with depth (Fig. 5c). To investigate the effects of shape only, we removed the internal layers in this simulation. In the previous examples, hyperbolae are created from truncation of internal layers, but with the irregularly shaped conduit, hyperbolae also result at the edges of each rectangle (Fig. 5d). This indicates that some of the diffraction hyperbolae observed in the data may result from changes in conduit shape. Additionally, because of asymmetry in the conduit geometry used in this simulation, hyperbolae (especially those nearest the surface) also appear slightly asymmetrical compared to those created with a smooth, singular rectangular conduit (Fig. 4b). Thus, the complex diffraction hyperbolae seen in the observed radar profiles may result from complicated, asymmetrical conduit geometries. Observations from descents into moulins confirm that, while their geometry in the top half of the ice is generally straight, moulins become wider and connected to other englacial conduits deeper down (Reference HolmlundHolmlund, 1988). Further, moulins may contain plunge pools and failed englacial conduits, which may cause diffraction hyperbolae (Reference ReynaudReynaud, 1987; Reference HolmlundHolmlund, 1988; Reference SchroederSchroeder, 1998; Reference VatneVatne, 2001). This simulation also shows that when a conduit narrows with depth, hyperbolic amplitude decreases with depth and may explain why some VSH patterns show increased attenuation with depth.

Simulations with different penetration depth

Lastly, we investigate changes in conduit depth using a 2 m wide, water-filled conduit that penetrates only 300 m below the ice surface (Fig. 6c), and compare this to a fully penetrating conduit (modeled previously). In the unconnected case, the model domain contains several continuous (unbroken) internal layers below the unconnected conduit. Despite this, the simulation produced diffraction hyperbolae appearing at each layer (Fig. 6d). Differencing the results of the partially connected simulation from the fully connected simulation highlights the subtle differences in hyperbolae amplitude at depth (Fig. 6e); multiple reflections below the unconnected conduit have smaller amplitudes than those below the fully connected conduit. Thus, attenuation of hyperbolic amplitude with depth can also result from radar profiles over deeply penetrating crevasses rather than fully connected moulins. The borehole identified in Figure 3d represents an analogy to this model simulation because it only penetrates to a 600 m depth but exhibits attenuating multiples that continue below this depth.

Fig. 6. (a) Model geometry and (b) results for a 2 m wide, water-filled vertical conduit centered at 2000 m penetrating through the entire ice thickness. (c) Model geometry and (d) results for a 2 m wide, water-filled vertical conduit with penetration to 300 m depth. (e) shows the difference between (b) and (d).

Discussion

Based on our model results we re-examine the characteristics listed earlier to determine the likelihood that the VSH patterns are related to englacial conduits and moulins. Vertically stacked hyperbolic diffractors can result from simple, smooth-walled vertical conduits truncating internal layers (e.g. Fig. 4b) and also from conduits that have non-uniform shape with depth (e.g. Fig. 5d). In the latter case, hyperbolae are produced at each kink in the conduit. Brighter amplitude returns are more likely if the conduit infill consists of both water and air in places. Complex and asymmetrical hyperbolae probably result from conduits with complex shapes, including changes in conduit width and orientation with depth and connections to other englacial conduits. An alternative explanation for the complex hyperbolic returns is that the radar is detecting a set of individual, vertically connected cracks that may represent a moulin in its initial stage of development. Over time, water flowing into the cracks widens and connects them, smoothing the walls, and evolves into a single, dominant conduit.

Diminishing hyperbolic amplitude with depth may result for numerous reasons. The conduit may not penetrate the full ice thickness, in which case the VSH pattern results from a deeply penetrating crevasse (Fig. 6) and deep hyperbolic returns are simply multiple reflections. Alternatively, the conduit may change direction at depth (not modeled), turning off-axis from the radar profile such that reflected energy is diminished. Additionally, the conduit may be ice-filled at depth, due either to creep closure or to refreezing of un-drained water (Fig. 4), or the conduit may narrow with depth (Fig. 5d). These scenarios are possible for moulins.

Several of the observed VSH patterns have strong-amplitude, below-bed hyperbolic returns. Below-bed diffraction hyperbolae can result from off-axis structure indicating a change in orientation of the conduit as it nears the base of the ice sheet. Alternatively, they can result from multiple reflections. Radar multiples are more likely if the permittivity contrast between the conduit and the surrounding ice is large (e.g. water and ice) (Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003), validating our observations that the large-amplitude diffraction hyperbolae result from water pockets or a partial water fill. Further, multiples at large (below-bed) travel times may indicate primary reflections from objects with correspondingly large travel times (within the ice column). This suggests that below-bed multiples may result from conduits that penetrate deep into the ice and/or are connected to the bed.

Finally, VSH patterns coincident with locally dipping internal layers and layer truncation at the bed are a clear indication of localized basal melting (Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others, 2001; Reference Catania, Conway, Raymond and ScambosCatania and others, 2006; Reference Buchardt and Dahl-JensenBuchardt and Dahl-Jensen, 2007). These VSH patterns more obviously represent moulins because they indicate vertical orientation and a clear connection between the surface and the subglacial drainage system. In the case of Figure 3e, layer downwarping is limited to the region immediately surrounding the VSH pattern indicating a smaller amount of subglacial melting. In contrast, Figure 3b shows internal layers deformed as far away as ∼7.5 km from the VSH pattern, indicating widespread subglacial melt due to surface water introduction. While internal layer shapes can help to confirm the presence of a moulin that once connected to the subglacial water system, the interpretation of internal layers is not useful everywhere. In particular, the internal stratigraphy in the downstream region of our study area (e.g. along the profile from GPSD1 to GPSD4 in Fig. 1) is difficult to resolve due to a high degree of scattering from near-surface or off-axis crevasses.

Our goal is to replicate some of the distinct characteristics of the VSH patterns listed earlier, in order to understand their origin. We insert a 2 m wide by 1000 m long vertical rectangle into the center of the model domain to simulate a smooth-walled, vertical conduit that penetrates the entire ice thickness (i.e. a moulin). We acknowledge that this idealization does not reflect the complicated geometry of actual moulins, as observed by investigators in the field (e.g. Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003). Our approach is to begin simply and to increase moulin complexity only when necessary. Our initial model results using this geometry (not shown) produced only two diffraction hyperbolae: one at the top of the ice, where the conduit opens to the free surface, and one at the base of the ice, where the conduit connects to the bed. The lack of hyperbolae within the ice in this simulation is due to the absence of point reflectors within the ice sheet, which are created when the conduit truncates other objects in the ice column (such as internal layers). To simulate this, we introduce nine 2 m thick horizontal internal layers across the model domain at various depths (Fig. 4a). These layers have electrical properties approximating a mixture of volcanic ash and ice (Table 1; Reference Adams, Perger, Rose and KostinskiAdams and others, 1996; Reference Miners, Wolff, Moore, Jacobel and HempelMiners and others, 2002).

Our model results and radar observations do not allow us to conclude that all VSH patterns in the study region result from moulins. They do, however, allow us to argue that many of them are likely to be moulins, especially those coincident with locally dipping internal layers. We now explore possible linkages between moulins and parameters that might influence moulin formation, such as the local ice thickness and the proximity to supraglacial lakes.

Relationship to ice thickness

Ice-thickness contours from our radar data are plotted in Figure 7. In addition to showing a general thinning toward the margin of the ice sheet, these data show that the downstream portion of our study region is characterized by a broad bed-topographic high beginning at ∼800 m a.s.l., reaching ∼200–400 m in height, extending across our survey region and aligned approximately orthogonal to the ice-flow direction. This bed ridge is coincident with a dark band of ice on the surface visible in satellite imagery, which represents a crevassed region observed in the field and described by Reference Price, Payne, Catania and NeumannPrice and others (2008). Ice flowing over this ridge thins by ∼40–50% on the stoss side, leading to along-flow extension, and thickens on the lee side, leading to along-flow compression. Reference Nath and VaughanNath and Vaughan (2003) indicate that tensile stresses in excess of ∼100 kPa are required for initiation of surface crevassing. We use the two-dimensional flowline model of Reference Price, Waddington and ConwayPrice and others (2007, Reference Price, Payne, Catania and Neumann2008) to calculate the tensile stress along a portion of our central (along-flow) radar line (CMPB–JAR1). We find two regions of elevated tensile stress: a region of ∼50 kPa in the vicinity of Swiss Camp and a region of ∼100 kPa at a distance of 12–15 km downstream of Swiss Camp (Fig. 8). The latter region is coincident with the bedrock ridge described above. Radar data in this region (Fig. 2) also show numerous VSH patterns and other diffraction hyperbolae, which we interpret as resulting from heavy surface crevassing. These observations are consistent with the idea that the bed topography creates an environment favorable to crevasse formation and, assuming sufficient meltwater supply, moulin formation. Along the other radar lines we also observe an increased frequency of moulins in regions with ice thickness less than ∼800 m, coincident with ice thinning over this bedrock ridge. Indeed, the two moulins associated with dipping internal layers (b and e in Fig. 7) are located in the crevassed region associated with the bedrock ridge.

Fig. 7. Map showing the off-axis distance to moulins (size of circle), supraglacial lake area (red dashed lines), the location of previously drilled boreholes (yellow stars), radar profiles (green lines) and ice thickness (colored contours). Circle diameter has been doubled relative to map scale for visibility. Labels a–g correspond to detailed radar images in Figure 3.

Fig. 8. Along-flow tensile stress in the downstream region of the radar data shown in Figure 2. Basal topography has been smoothed for efficiency of model calculation. Contours of tensile stress are as indicated.

Meltwater supply increases downstream from the equilibrium line due to warmer surface temperatures and a larger upstream catchment area; this may partially explain why we observe a greater number of moulins downstream from the equilibrium line than upstream (Fig. 7). The exception to this is the GPSD3–GPSD1 line (Fig. 1), where a few moulins are visible within a few kilometers of the equilibrium line. This is probably due to the presence of thinner ice further upstream along this line than the other radar lines; thinner ice may be necessary to form moulins in a region where meltwater supply is not as plentiful. The recent observations of Reference DasDas and others (2008) were made in a region south of our field site ∼20 km downstream from the equilibrium line. Ice here is estimated to be ∼1000 m thick, however, and a much larger water supply (due to drainage of a supraglacial lake, and an increased distance downstream of the equilibrium line) and warmer ice may have permitted localized moulin formation despite the relatively thick ice here.

Conclusions

Our radar observations and modeling results allow us to identify and understand the characteristics of unique en-glacial diffraction patterns found in the ablation region of the Greenland ice sheet. We conclude that many of the observed patterns are due to the presence of moulins, which connect the supra- and subglacial drainage systems and provide a mechanism for seasonal increases in glacier velocity. The observed moulins probably have a complicated geometry, since the diffraction patterns are complex and asymmetrical. This is consistent with numerous reports on the shapes of moulins and englacial conduits through direct exploration (Reference ReynaudReynaud, 1987; Reference HolmlundHolmlund, 1988; Reference SchroederSchroeder, 1998; Reference VatneVatne, 2001; Reference Stuart, Murray, Gamble, Hayes and HodsonStuart and others, 2003). Since moulins may initially form from surface crevasses, adequate tensional stress is required for surface meltwater to enter the englacial environment via crevassing. Surface meltwater availability increases downstream from the equilibrium line, as a result of both warmer surface temperatures and a larger upstream catchment area, and the ice thins towards the ice-sheet margin. All of these conditions are likely to promote moulin formation. We propose that moulins are more likely to develop closer to the ice-sheet margin, where meltwater is more plentiful and where rough bed topography beneath relatively thinner ice is more likely to result in significant surface crevassing; bed topography has a larger influence on the longitudinal stress closer to the margin, where the ratio of bed topography to ice thickness is greater. Consequently, the threshold for crevasse formation is likely to be more easily reached in thinner ice.

Reference Price, Payne, Catania and NeumannPrice and others (2008) show that surface velocity variations measured at Swiss Camp (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002) and 15 km downstream, near JAR1 (Reference Rumrill, Neumann and CataniaRumrill and others, 2006), are adequately explained by increased sliding in the region 8–12 km downstream from Swiss Camp. This hypothesis is contrary to that in which sliding variations at Swiss Camp are the result of local surface meltwater penetration (via moulin formation) through 1200 m thick ice. Instead, Reference Price, Payne, Catania and NeumannPrice and others (2008) hypothesize that velocity variations at Swiss Camp result from longitudinal coupling within the ice, between slower-moving ice near the equilibrium line and faster-moving (faster-sliding) ice 8–12 km downstream from Swiss Camp. Although our radar lines under-sample the area, the data presented here show general support for this hypothesis; we observe nearly all moulins in the region of elevated tensile stress associated with a broad bedrock ridge (in ice that is thinner than ∼800 m). The larger tensile stress in this area allows crevassing, which provides an initiation point for moulin formation if meltwater supply is large enough.

While tensile stress and a sufficient meltwater supply may be realized at the bottom of a supraglacial lake (Reference Alley, Dupont, Parizek and AnandakrishnanAlley and others, 2005), we observe only a moderate correlation between supraglacial lakes and moulin location. According to Reference Joughin, Das, Howat and MoonJoughin and others (2006), lakes on the surface of the Greenland ice sheet re-form in the same place from year to year, with location determined by the local surface topography, which changes slowly through time. Although there is evidence that some lakes may drain locally to the base of the ice sheet (Reference DasDas and others, 2008), we infer that many drain over the ice surface or through the ice in such a way that no internal reflectors are generated; the former seems more likely than the latter.

We find little evidence for moulins at or above the equilibrium line of the ice sheet, despite the presence of lakes here. The importance of lake-drainage events on ice dynamics in the interior of the ice sheet depends on how the ice thickness and surface elevation evolve in the future. Given the transient nature of the observed response, lake-drainage events appear to represent only a moderately positive feedback on ice-thickness change. In order to significantly increase the spatial extent of moulins, and thus the coupling between surface melting and increased basal motion, we suggest that the ice thickness (and associated surface elevation) would need to change appreciably by a different means, such as surface melting or enhanced ice flow. It is possible that moulins could form in the equilibrium or accumulation zones of the ice sheet given enough localized meltwater (as appears to have been available in the 2007 melt season (Reference MoteMote, 2007)). Our data, however, suggest that moulins are much more common in the ablation zone, where there is more water on the surface, where the ice is thinner and where tensile stresses are higher.