Landforms produced beneath former ice sheets offer insights into inaccessible subglacial processes and present analogues for how current ice masses may evolve in a warming climate. Large subglacial channels cut by meltwater erosion (tunnel valleys [TVs]) have the potential to provide valuable empirical constraints for numerical ice-sheet models concerning realistic melt rates, water routing, and the interplay between basal hydrology and ice dynamics. However, the information gleaned from these features has thus far been limited by an inability to adequately resolve their internal structures. We use high-resolution three-dimensional (HR3-D) seismic data (6.25 m bin size, ∼4 m vertical resolution) to analyze the infill of buried TVs in the North Sea. The HR3-D seismic data represent a step-change in our ability to investigate the mechanisms and rates at which TVs are formed and filled. Over 40% of the TVs examined contain buried glacial landforms including eskers, crevasse-squeeze ridges, glacitectonic structures, and kettle holes. As most of these landforms had not previously been detected using conventional 3-D seismic reflection methods, the mechanisms that formed them are currently absent from models of TV genesis. The ability to observe such intricate internal structures opens the possibility of using TVs to reconstruct the hydrological regimes of former mid-latitude ice sheets as analogues for contemporary ones.

Throughout the Quaternary, the growth and retreat of ice sheets across high- and mid-latitude continental shelves drove major changes to topographic relief and global sea level (e.g., Batchelor et al., 2019). Subglacial water flow to the margins of these ice sheets excavated huge channels, kilometers wide and hundreds of meters deep, which are referred to as tunnel valleys (TVs) (ÓCofaigh, 1996; Huuse and Lykke-Andersen, 2000; Kehew et al., 2012). Multiple cross-cutting generations of TVs, potentially correlated with seven glacial cycles, are buried beneath the seafloor in the North Sea (Kristensen et al., 2007; Stewart and Lonergan, 2011). TVs provide a unique opportunity to investigate the subglacial plumbing system of the ice sheets that formerly covered northwestern Europe and to learn more about inaccessible basal processes that regulate ice-sheet flow and retreat. However, there is currently no consensus on the timescales and mechanisms through which glacial meltwater forms TVs or their impact on ice-sheet dynamics, and theories about their formation range between rapid erosion during catastrophic floods to gradual development through multiple incision events (Kehew et al., 2012). The infill of TVs is similarly complex; some models propose that TVs were cut and filled incrementally (e.g., Jørgensen and Sandersen, 2006), while others favor simultaneous cutting and infilling in a conveyor-like fashion tracking ice front recession (Praeg, 2003; Kristensen et al., 2008). Disagreement over TV infill reflects spatial heterogeneity in TV form and sedimentary composition in outcrops, a scarcity of borehole samples, and data resolution constraints that limit the extent to which their internal architecture can be analyzed in three dimensions (3-D; Huuse and Lykke-Andersen, 2000; Praeg, 2003). We use novel high-resolution 3-D (HR3-D) seismic data from the North Sea to examine the infill of TVs in unprecedented detail and discuss the implications for TV genesis.

We examined six HR3-D seismic data sets, covering ∼60 km2, from the central North Sea (Figs. 1A and 1B). The acquisition system comprised two 1200-m-long streamers towed 3 m beneath the sea surface with 96 hydrophone groups at 12.5 m spacing, a 6.25 m shot interval, and a 1 ms sample rate (Games, 2012). The seismic source consisted of two 160 in3 sleeve air gun clusters with a 20–250 Hz signal frequency. Data processing included swell noise attenuation, tide correction, multiple suppression, two passes of velocity analyses run at 250 × 250 m intervals, normal moveout correction, and bandpass filtering. The final processed data sets consist of time-migrated 3-D stacks with a 1 ms sample rate, a 6.25 × 6.25 m bin size, a vertical resolution of ∼4 m, and a detection limit along individual reflectors of ∼0.5 m. In comparison, the 3-D seismic data previously used to examine TVs in this region typically have bin sizes of 12.5–50 m and a vertical resolution of ∼8–16 m (Stewart et al., 2013). Depth conversions used a velocity of 1800 m s−1 (Stoker et al., 1985).

The six HR3-D seismic data sets image 19 cross-cutting incisions that are 300–3000 m wide, up to 300 m deep, and possess undulating thalwegs. We interpret these as TVs formed by subglacial meltwater based on their distinctive morphology (e.g., Stewart et al., 2013). The TVs contain between one to four discrete fill units, and most consist of two seismic facies. The infill is highly variable and lacks consistent patterning between each TV, although some common facies are present (Figs. 1C and 1D). Over 50% of the TVs contain a chaotic, largely homogenous fill package that often overtops the TV shoulders. We interpret this facies as subaqueous or subaerial outwash. In addition, the lowermost facies of ∼60% of TVs contains discontinuous sub-parallel reflections that are interpreted as clays and sands deposited in a subglacial or proglacial subaqueous setting. We estimate that twice the number of infill units can be resolved in our HR3-D seismic data compared to conventional 3-D seismic data; this improvement represents a step-change in our ability to observe fine-scale TV infill structures (Fig. 2).

Several TVs contain high-amplitude reflections between their upper and lower fill units. When mapped in 3-D, these reflections delineate sinuous ridges that are up to 2.8 km long, 30–150 m wide, and 5 m high on average (Fig. 3A). Based on their morphological similarity to features elsewhere (e.g., Storrar et al., 2014; see Fig. S1 in the Supplemental Material1), we interpret the ridges as eskers that were deposited in meltwater conduits present at the base of an ice sheet.

Small hollows containing well-layered fill are observed buried within the outwash unit that caps some TVs. The hollows are 100–300 m wide, 200–750 m long and have side slopes of 10–35° (Figs. 3B and 3F). We interpret the hollows as kettle holes that formed when stagnant ice blocks were stranded in shallow water (e.g., Ottesen et al., 2017). Their well-layered internal structure likely reflects infilling with fine-grained sediments once the ice blocks melted in place or floated off in a marine setting.

Subtle ridge patterns of two morphologies are buried within the TVs and appear as undulating reflections with high acoustic amplitude. The features form irregular networks of four- to six-sided polygons that are rhombohedral in planform morphology, 40–80 m in diameter, and slightly hollow in the center. Individual ridges are 20–250 m long, 20–30 m wide, and <4 m high from base to crest (Figs. 3C and 3G). Several processes have been documented to produce rhombohedral ridges, including fluid escape, diagenesis, glacier surging, and permafrost processes (e.g., Ottesen and Dowdeswell, 2006; Morley et al., 2017; Bellwald et al., 2019). We rule out fluid escape due to the absence of visible chimneys or pipes above the features and suggest that substantial diagenetic change is unlikely to occur at such shallow depths (<250 m below seabed). The fact that the ridges are confined to the TV rather than extending onto the surrounding banks, combined with their irregular rims, suggests these landforms do not result from permafrost processes. Rather, based on their morphological similarity to features on glaciated terrains elsewhere (Ottesen and Dowdeswell (2006); Fig. S2), we interpret these landforms as crevasse-squeeze ridges that formed through sediment injection upwards into basal fractures beneath grounded ice (Rea and Evans, 2011; Evans et al., 2016).

Groups of sub-parallel curvilinear ridges, commonly located in the upper third of TV infill, constitute the second type of ridge pattern (Fig. 3D). The symmetrical ridges are spaced 50–100 m apart, 100–300 m long, <3 m high, terminate sharply, and are oriented perpendicular to the long axis of the TVs (Fig. 3D). We interpret these as crevasse-squeeze ridges formed as grounded ice retreated through the TVs.

Distinctive features are also found at the base of the TVs in the form of networks of anabranching channels (80 m wide and ∼6 m deep on average) incised around streamlined bars 55–540 m long and 30–165 m wide (Fig. 3E). The thalweg of the TV displayed in Figure 3E undulates, which suggests a subglacial rather than subaerial fluvial origin. Several TVs also contain chaotic and displaced reflections along the sides and bases, which we interpret as evidence of slumping, faulting, and glacitectonic thrusting (Figs. 1D and 3H).

Over 40% of the TVs examined here contain buried landforms of glacial origin (Fig. 3). Most of these features are too small to interpret using conventional 3-D seismic data and would be difficult to detect using 2-D seismic data or boreholes alone (Fig. 2). Accordingly, glacial landforms, both erosional and depositional in nature, are likely far more common inside TVs than previously recognized.

The ice sheets formerly occupying the North Sea Basin were underlain by thick sequences of flat-lying unconsolidated sediments. In such settings, subglacial water is thought to be transported in networks of broad and shallow sedimentary channels (Walder and Fowler, 1994) that may resemble the braided channels observed at the base of some TVs (Fig. 3E). Braided channel systems may be eroded into the substrate surrounding a subglacial conduit when the latter drainage system is temporarily overwhelmed by pulses of meltwater supplied from the ice sheet surface (Fig. 4A; Lewington et al., 2020). This could be achieved by supraglacial lake drainage via hydrofracture (e.g., Das et al., 2008). Repeated transfer of surface meltwater to the bed, coupled with lateral enlargement by overriding ice, would gradually excavate a TV and induce glacitectonic deformation structures and mass movements along its sides (Fig. 3H); these processes have also been inferred from TVs exposed in sections of Late Ordovician glaciogenic rocks (Hirst et al., 2002; Le Heron et al., 2004).

Although eskers have been documented along the base of unfilled subaerial TVs in North America (e.g., Brennand and Shaw, 1994), these features are rarely reported inside filled TVs (van der Vegt et al., 2012). Traditionally, eskers and TVs are thought to seldom coexist, as eskers are associated mainly with hard bed substrates while TVs form in poorly consolidated sediments (Clark and Walder, 1994; Huuse and Lykke-Andersen, 2000). The eskers observed within the TVs here indicate that this association probably reflects the poor preservation potential of eskers on soft substrates (Storrar et al., 2019). The near-central stratigraphic position of the eskers within the TV infill is significant, as it demonstrates that grounded ice occupied the TV until the channel filled to approximately half of its accommodation space, which implies significantly greater longevity of ice occupation than was assumed in previous models of TV genesis (van der Vegt et al., 2012). It is possible that eskers in the TVs represent the final sedimentary cast of migrating meltwater conduits that were filled during the last stages of deglaciation under a thinning ice-sheet terminus (Storrar et al., 2014; Beaud et al., 2018). Lateral continuity and a high degree of esker preservation implies gradual ice retreat from these TVs and precludes reworking by ice readvances before the eskers were buried by outwash sediments. As esker continuity reflects the ice dynamics at the time they were formed (Storrar et al., 2014), further investigation into the degree of esker fragmentation within TVs may hint at the style and rate of past ice-sheet retreat (Livingstone et al., 2020).

Crevasse-squeeze ridges are diagnostic of surging glaciers (Sharp, 1985). These features form at the termination of a surge through sediment injection into basal crevasses and are preserved by the stagnation and in situ downwasting of the overlying ice (Rea and Evans, 2011). Kettle holes can also form from ice stagnation and meltout during the quiescent phase of tidewater glacier surges (Ottesen et al., 2017).

The presence of landforms indicative of glacier surging within some TVs may imply that they acted as conduits for fast ice flow prior to being fully infilled. The absence of these landforms beyond the TVs indicates that either surges were confined to the TVs due to hydrological factors that initiated the fast flow or that such delicate features were not preserved outside of the TVs. Experiments with silicon models of the subglacial hydrological system demonstrate that TV formation often coincides with surges in the velocity of the model glacier that are triggered by increases in basal water pressure, causing the ice to decouple from its bed (Lelandais et al., 2016). The glacial landforms observed here support a link between the transport of water in TVs and dynamic ice behavior.

The variety of landform assemblages preserved inside the cross-cutting TVs hints at the diverse ice-sheet regimes that formed and filled them. Our HR3-D seismic data suggest that TVs were incised gradually by migrating meltwater channels driven by pulses of meltwater from the ice-sheet surface (Fig. 4A) and were enlarged by ice flow, which induced deformation structures and mass movements along their sides (Prins et al., 2020) (Fig. 4B). Gradual ice retreat from TVs may fill and preserve the most recently active channels as laterally continuous eskers (Fig. 4C), while dynamic ice flow through the TVs, possibly as surges or readvances, followed by stagnation and downwasting, may be indicated by the presence and preservation of crevasse-squeeze ridges and kettle holes (Fig. 4D).

Our study marks the first time that abundant glacial landforms have been convincingly imaged within buried TVs in the North Sea. The presence of eskers and crevasse-squeeze ridges within the mid–upper TV fill packages demonstrates that grounded ice played an active role in TV incision and was present for a substantial time during filling. For these delicate landforms to have been preserved in the geological record, reoccupation of the TVs between different glacial cycles must have been limited, although localized ice readvances may have occurred. This result constrains the formation and infilling of each TV generation to a single glaciation and supports the notion that the multiple generations of TVs present in the North Sea record at least seven glacial advances across northwestern Europe (Stewart and Lonergan, 2011). Greater coverage of HR3-D seismic data on glaciated margins, combined with chronological constraints from shallow drilling, may permit TVs to become a resource to improve understanding of the hydrological systems and dynamics of former ice sheets.

1Supplemental Material. High-resolution 3-D seismic data set resolution considerations and Figures S1 and S2. Please visit to access the supplemental material, and contact with any questions.

We thank bp, Harbour Energy, Equinor Energy AS, Lundin Energy Norway AS, Petoro AS, Aker BP ASA, Total E&P Norge AS, and Petroleum Geo-Services for permission to publish images from the HR3-D seismic data and the Central North Sea MegaSurveyPlus. IHS and Schlumberger provided academic seismic interpretation software licenses. J. Kirkham is supported by Natural Environment Research Council grant NE/L002507/1. This research is supported by the British Antarctic Survey Polar Science for Planet Earth programme. The interpretations made in this paper are the views of the authors and not necessarily those of the license owners. Sean Gulick, Benjamin Bellwald, and Brian Todd are thanked for helpful reviews.

Gold Open Access: This paper is published under the terms of the CC-BY license.