Architecture of tunnel valleys in the southeastern North Sea: new insights from high‐resolution seismic imaging

Tunnel valleys are assumed to form near the margin of ice sheets. Hence, they can be used to reconstruct the dynamics of former ice margins. The detailed formation and infill of tunnel valleys, however, are still not well understood. Here, we present a dense grid of high‐resolution 2D multi‐channel reflection seismic data from the German sector of the southeastern North Sea imaging tunnel valleys in very great detail. Three tunnel valley systems were traced over distances ranging between 11 and 21 km. All tunnel valleys are completely filled and buried. They differ in incision depth, incision width and number of incisions. The tunnel valleys cut 130–380 m deep into Neogene, Palaeogene and Cretaceous sediments; they show a lower V‐shaped and an upper U‐shaped morphology. For individual tunnel valleys, the overall incision direction ranges from east–west to northeast–southwest. Two tunnel valleys intersect at an oblique angle without reuse of the thalweg. These valleys incise into a pre‐existing glaciotectonic complex consisting of thrust sheets in the northwest of the study area. The analysis of the glaciotectonic complex and the tunnel valleys leads us to assume that we identified several marginal positions of (pre‐)Elsterian ice lobes in the southeastern North Sea.


Introduction
Major ice sheet advances beyond the polar limits in the northern hemisphere are well documented for the Quaternary (last 2.6 Ma); the last 900,000 years are assumed to be the period with the most prominent ice sheet advances (Ehlers et al., 2018). Understanding the glacial history of continental shelves is key to reconstructing the configuration of former ice sheets. Fortunately, surface and subsurface landforms serve as an archive that can be used to reconstruct past sedimentation and erosion processes. Identifying these landforms and the corresponding mechanisms of formation can support the reconstruction of ice sheet configurations during the Pleistocene (Aber and Ber, 2007;Dowdeswell et al., 2016). Fluctuations in oxygen isotopes found in deep sea cores have been linked to major cold and warm phases of the planet. The rhythmic pattern of gradual cooling and rapid warming of the climate has been defined as Marine Isotope Stages (MIS). The dating of glaciogenic sediments enabled the scientific community to assign the well-known glaciations around the globe to specific MIS (Ehlers and Gibbard, 2004).
Tunnel valleys are defined as elongated depressions incised into bedrock or unconsolidated sediment (Ó Cofaigh, 1996;van der Vegt et al., 2012;Kehew et al., 2012). Most studies conclude that tunnel valleys result from subglacial palaeo-drainage channels with high flow rates of meltwater eroding valleys into the sediment and bedrock (Ussing and Selskab, 1903;Moores, 1989;Ghienne and Deynoux, 1998;Lonergan et al., 2006). Still, there is no consensus on whether the erosion mechanism is predominantly based on a 'steady-state' or a 'catastrophic' drainage. The 'steadystate' hypothesis assumes the drainage of pressurised groundwater gradually eroding a subglacial channel, which is broadened by fluvial erosion and direct glacial erosion (Boulton and Hindmarsh, 1987;Kehew et al., 2012;van der Vegt et al., 2012). The 'catastrophic' hypothesis assumes subglacial and supraglacial ponding and accumulation of meltwater and subsequent sudden drainage leading to the incision (Wright, 1973;Lewis et al., 2006).
We set out in July 2017 with the expedition AL496 ( Fig. 1 b) as part of the project Nordfriesland-Süd to map the shelf architecture and to identify glacial landforms in order to investigate the glacial history of the survey area. We present high-resolution seismic reflection data, which provide an insight into tunnel valley morphology and infill in very high detail. We aim to image and reconstruct the internal structure and spatial development of tunnel valleys. The dense grid of 2D reflection seismic sections ( Fig. 1 b) serves as the basis for tunnel valley description, classification and spatial extrapolation. This information is used to identify a former ice sheet margin in the German sector of the southeastern North Sea. Additionally, we aim to contribute to the ongoing debate on the formation mechanism of tunnel valleys.

Methods
We acquired 2D reflection seismic and parametric sediment echo sounder data to image subsurface structures from 0 to 800 ms corresponding to approx. 640 m beneath the seafloor (1600 m/s is used for all time-depth conversions in the absence of a precise velocity model [Cotterill et al., 2017;Coughlan et al., 2018]). A total of 1058 km of 2D highresolution multi-channel reflection seismic data were acquired in a closely spaced 2D grid covering approximately 800 km 2 ( Fig. 1 b). We focus on the analysis of the 2D reflection seismic data because these data show sufficient penetration to image the submerged palaeolandscape.

2D multi-channel reflection seismics acquisition
A Sercel Mini GI gun was used in a harmonic mode with a reduced volume of 0.1 L for both the generator and injector chambers for all seismic lines. The main frequency of the source is approximately 300 Hz. A Geometrics GeoEel solid state digital streamer was used to record the seismic response. Our setup included a 35 m tow cable (30 m in water), a 10 m stretch, and seven active sections, each containing eight channels with a channel spacing of 1.56 m. This resulted in 56 active channels and an active streamer length of 87.5 m. The broad frequency spectrum of the source (70-600 Hz) allowed us to image deep sediments as well as to correlate the shallow strata with the simultaneously acquired sediment echo sounder data. The vertical resolution is in the metre range; maximum resolution of about 1 m is achieved in the upper tens of metres below the seafloor. Resolution slightly decreases with increasing penetration due to absorption of high frequencies.

Seismic data processing
This study followed standard reflection seismic processing routines in Schlumberger's Vista Seismic Processing software: (1) Ormsby frequency filtering (55/110 and 600/1200 filter flanks); (2) CMP binning; (3) FK-filtering; (4) normal moveout (NMO) correction using a velocity model based on interactive picking of best-fit velocities; (5) despiking; (6) CMP-stacking; and (7) finite-differences time-migrating these data with a velocity model inferred from an interactive velocity analysis. A shot point distance of 9 m or less allowed the common midpoint (CMP) bin spacing to be set equal to the channel spacing of 1.56 m, resulting in a CMP fold of 6 to 10.
Further special processing involved the removal of multiple reflections, which interfered with primary reflectors due to the very shallow water depth (less than 25 m). Removal of both seafloor and internal reflector multiples was accomplished using a fine-tuned predictive deconvolution with an adapted lag time based on picks of the seafloor and the first multiple. This approach yielded very good results owing to the flat seafloor throughout the study area.

Seismic data interpretation
The fully processed seismic sections were loaded in the IHS Markit Kingdom Software 2017 for interpretation. Kingdom was used to trace horizons and the morphology of subsurface landforms. The horizons were exported for interpolation and plotting with the Generic Mapping Tools open-source software.

Results
The stratigraphic units and subsurface landforms shown in this section were identified based on differences in acoustic properties on the seismic sections. To focus on the description of the tunnel valleys, we refrain from providing a detailed seismostratigraphical analysis of the seismic sections. Instead, we use a more general approach delineating the seven main seismic units (SU) up to 800 ms penetration that have been correlated to the seismic stratigraphy framework of the Federal Agency of Geosciences and Resources of Germany (BGR), which is based on the results of Thöle et al., 2014 andKockel, 2002. Table 1 outlines the correlation of the SU defined here with the stratigraphic framework of the BGR, which includes well control from many deep wells in the German sector of the North Sea.

Stratigraphy of the working area
The seismic sections in the north of the study area reveal the complete sedimentary sequence of the upper Cretaceous to the Quaternary (Fig. 2). Seven SU were identified. Fig. 3 outlines the depth to the major boundaries during late Cretaceous to late Miocene identified as continuous reflectors in the subsurface. These sediments are overlain by Plio-Pleistocene sediments, which build the foundation for the Holocene drape (Thöle et al., 2014). SU1 is the deepest seismic unit characterised by slightly undulating, continuous reflectors of very low amplitude (Fig. 2). Its top is defined by a continuous reflector of significantly higher amplitude than the otherwise very low amplitude reflectors of the thick sedimentary package. This reflector is found at 730 ms two-way-travel-time (TWT; 580 m) (Figs 2 and 3a). Deposits of SU1 represent Cretaceous sediments (Table 1, Thöle et al., 2014).
SU2 is the thickest seismic unit. It is characterised by a transition in seismic facies in the upper third of the unit (at 525 ms, 420 m on the seismic section shown in Fig. 2). Beneath this transition, these data image a thick package of continuous reflectors of very low amplitudes, whereas the reflections show increased amplitudes above the transition. SU2 includes several reflector packages characterised by high amplitudes at their base and low amplitudes at their top (Fig. 2, 525 ms to 475 ms, 420 to 380 m and 440 ms to 400 s, 350 to 320 m). Its shallowest section is characterised by well-stratified low amplitude reflectors. The deposits of SU2 represent Palaeogene sediments (Table 1, Thöle et al., 2014).
SU3 consists of reflectors characterised by an alternating high to medium amplitude pattern from 400 ms (320 m) upwards. The high to medium amplitudes dominate up to 240 ms (190 m), where amplitudes decrease to medium to low. A narrow high amplitude reflector characterises its top, which represents the mid-Miocene unconformity. The mid-Miocene unconformity is truncated by the Quaternary succession (Fig. 3, Lutz et al., 2009) in the central and southern part of the study area and is thus absent from parts of our study area. A high amplitude continuous reflector marks the base of SU3 (Fig. 2), which has been correlated to represent lower to late Miocene sediments (Table 1, Thöle et al., 2014).
SU4 is only found in the northern part of the study area. The abrupt change in seismic facies to very low amplitude reflectors topped by three sharp high amplitude reflectors characterises this unit (Fig. 2). The pattern of alternating high amplitude to low amplitude reflectors, as seen in the underlying units, is absent here (Fig. 2). This clearly delineates SU4 from SU3. The deposits of SU4 most likely represent Pliocene to early Pleistocene sediments because they are constrained by Miocene sediments at their base and the Quaternary succession at their top (Table 1). Due to the lack of data, an unambiguous correlation with Pliocene sediments is not feasible.
SU5 is characterised by alternating strong and weak reflectors in its lower part, which transitions to chaotic strong reflectors in the upper part (Fig. 2). Tracing of reflectors is feasible only in the lower part. The onset of deformation in the lower parts of SU5 delineates this unit from SU4 and links SU5 to the Quaternary succession ( Fig. 2, Table 1). SU5 is thus comprised of mainly Pleistocene sediments that are truncated at its top throughout the study area.
SU6 demarcates the erosional flanks and the fill of tunnel valleys (Fig. 2). These tunnel valley complexes truncate SU1-SU5. SU6 often shows a lower V-shaped unit and an upper U-shaped unit, imaged by a change in flank steepness and a change in seismic facies (Fig. 2). Chaotic low to medium amplitude reflections characterise the lower V-shaped subunit, whereas medium to high amplitude reflectors of good  continuity characterise the upper U-shaped subunit (Fig. 2). Reflectors of the upper U-shaped subunit are sometimes inclined and differ in thickness between individual tunnel valleys. SU6 includes up to four prominent reflectors, which represent subsequent incisions into the tunnel valleys (see 'Tunnel valley morphology, internal structure and spatial extent', below). SU7 is the shallowest seismic unit and its top represents the seafloor (Fig. 2). Its surface is flat throughout the study area. SU7 truncates SU5 and SU6 (Fig. 2). This truncation is represented by an erosional surface that represents the base of SU7. Shallow elongate depressions resembling small channels (up to 300 m width) are often found in SU7 (e.g. Fig. 2 at 25 ms TWT, 20 m, and 2550 to 2800 m along profile). The top sediments of SU7 represent Holocene deposits, whereas the lower sediments of SU7 represent late Pleistocene deposits resting unconformably on SU5 and SU6.

Glaciogenic unconformity
A glaciogenic unconformity has been described before for wide areas of the southern North Sea (Moreau et al., 2012) as well as for the southern Baltic Sea (Al Hseinat and Hübscher, 2014). It formed during the Quaternary succession likely by glacial erosion and is identified based on reflectors lying unconformably on Neogene or older strata. The glaciogenic unconformity has thus been recognised as the onset of the Quaternary, on top of which ice sheets have deposited sediments during the Pleistocene.
The glaciogenic unconformity in the study area is represented by prominent glacial erosional surfaces and the tunnel valley systems, which truncate the lower units SU1-SU4. On top of this unconformity, SU5 and SU6 have been deposited and thus represent the Quaternary succession (Sejrup et al., 1991). The high spatial heterogeneity and often ambiguous seismic representation of the glacial erosional surfaces complicates tracing them reliably over longer distances (Fig. 4). The connection of tunnel valley flanks and glacial erosional surfaces can be used to link both features to the same time of formation.
Tunnel valley morphology, internal structure and spatial extent We mapped three overdeepened buried glacial valley systems consisting of five large tunnel valleys and several smaller valleys (Fig. 5). We focus on characterising the five largest and fully imaged valleys that we traced unambiguously on neighbouring lines. From tracing the landforms over parallel seismic sections, we reconstructed the lateral and vertical extent of the tunnel valleys and a GTC (Fig. 5). Individual plots for TV1-TV3's base reveals an undulating thalweg (Fig. 6). For the sake of clarity, we named the individual main valleys TV#. The fill of the main valleys is denoted by TV#f# and their internal erosional surfaces are denoted as TV#e#.
Although partly surveyed on its edge, the northernmost valley is the broadest and deepest valley with a width of up to 4.5 km and incision depth up to 480 ms (~380 m), (Figs 5 and 7). It is termed tunnel valley 1 (TV1) from hereon. TV1 is characterised by erosional truncations of the surrounding sediments (TV1e1 incised into SU2-SU5) and a subsequent deposition of sediments inside the valley (TV1f1-4 in SU6, Fig. 7). The southern flank of TV1e1 is imaged as a strong reflector and is sharper and steeper (up to 19°) than the northern flank (note that Fig. 7 is roughly eight times vertically exaggerated). The northern flank of TV1e1 is broader, less pronounced (missing a clear seismic reflector) and shallower (approximately 6°). Undisturbed layering of the surrounding strata to the north of the valley (approximately 60 ms to maximum imaging depth at approximately 800 ms) highlights TV1's sharp incision flanks (Fig. 7). Stratification inside this valley is generally not very pronounced and absent in its deeper part; chaotic reflections with low amplitudes and low continuity dominate the deeper part (TV1f1) and preclude the tracing of reflectors. Overall, the valley fill (TV1f1-4) is divided by erosional unconformities marking the different incisions (TV1e2-4), which partly result in sharp seismic reflectors. Two major incisions of 350 ms (280 m) by TV1e1 and 260 ms (210 m) by TV1e2 and a third incision of smaller magnitude up to 100 ms (80 m) by TV1e3 can be identified in TV1 (Fig. 7). The latest incision cuts older incisions at the top (TV1e4). These erosional surfaces delineate several stages in the cut-and-fill structures, which indicate that TV1e1 predates TV1e2, which predates TV1e3, which predates TV1e4. Weak layering is seen in TV1f4,  (Figs 5 and 7), whose lower part connects to TV1e1. A strong reflector at the top of TV1s connects to TV1e3, thus linking the lower fill to TV1f1 and the fill above the reflector to TV1f3.
A GTC has been identified in the northwestern part of the survey area (Fig. 5). Several sheets with a thickness of up to 100 ms (80 m) have been thrust towards the north to northwest (Fig. 7). The GTC initiates below an updipping prominent glacial erosional surface between TV1 and TV2 (Fig. 8) and the thrust sheets dominate the shallow substratum up to 200 ms (160 m) depth in the northwest of the survey area. TV1 cuts the northernmost limit of the GTC (Fig. 7).
Tunnel valley 2 (TV2) cuts approximately 210 ms (170 m) deep and can be traced for about 20 km from northeast to southwest (Figs 5,9,11). The seismic section in Fig. 9 images two prominent reflectors in the strata next to the base of TV2 comprised of north-northwest dipping nearly undisturbed sediments (SU3), clearly incised by TV2e1 and TV2e2 (at 18300 to 18700 m along profile). The average width of TV2 is approximately 1.5 km (Fig. 5). The infill of TV2 is well stratified on most seismic sections, often showing a larger upper U-shaped unit (TV2f2, SU6b) and a comparably small lower V-shaped unit with less internal stratification (TV2f1, SU6a; Figs 9 and 10). The change in seismic facies correlates with the change in flank steepness (Fig. 10). Tracing TV2 from east to west on the perpendicular seismic sections reveals that it originates in the east from well-layered sediments deposited on top of the broad U-shaped base (TV2f2, Fig. 11 at 10400-12500 m along profile). The eastern U-shaped base (TV2e2) shows small depressions in its western rim (Fig. 11 at 12 150 m and 190 ms, 150 m). The upper valley fill stratification (TV2f2) dips slightly towards the south-southeast. Westward of the profile shown in Fig. 11, the tunnel valley fill develops into a pattern characterised by the absence of layering in the lower V-shaped unit, whereas its upper U-shaped part remains stratified (Fig. S2 plots c-f). In the central part of the study area, TV2 has a tributary valley (TV2t1), which joins TV2 about 1.9 km further southwest (Fig. 5).
Tunnel valley 3 (TV3) runs approximately westsouthwest to east-northeast (Fig. 5) and intersects TV2. Tracing the erosional base of both tunnel valleys (TV2e1 and TV3e1) from the easternmost sections towards the west (Figs S2 and S3), it is clear that the valleys initiated as individual valleys in the east. After a few kilometres, they intersect and then separate again (Figs 1 and 11). TV3 clearly does not reuse TV2's thalweg, which can be seen at their intersection (Fig. 11) and parallel seismic sections (Figs S2 and S3). Instead, TV3e1 cut TV2 and progressed towards the west. TV3's width is about 1.5 km, similar to TV2 (Fig. 11). TV3 can be traced for about 18 km. It is incised deeper than TV2 with a mean incision depth of 280 ms (~220 m) and a maximum incision depth of 300 ms (~240 m) into the deposits of SU2 and SU3 (Fig. 11). On the western seismic sections, TV3's lower unit (TV3f1) consists of a chaotic facies but becomes increasingly stratified above (TV3f2) (Fig. S3). Disturbed reflectors in the strata beneath the tunnel valley spatially coincide with the change in flank steepness. The central valleys TV2 and TV3 differ from TV1 in their internal structure and size. They are narrower and incisions into the pre-existing sediments are shallower (Fig. 5). TV2 and TV3 can be separated in plan view based on their mean incision depths alone, which are more or less constant along their thalwegs (Fig. 5).
The southern valley complex consists of two main valleys, TV4 and TV5 (Fig. 12 at 8700-11300 m and 11300-13200 m along profile), and smaller side valleys ( Fig. 12 at 8600-9200 m along profile). These tunnel valley complexes show significantly shallower incision depths of up to 160 ms (approximately 130 m). The flanks connect to a shallow glacial erosional unconformity towards the north (Fig. 12 at 6500-8000 m along profile). Their infill (TV4f1 and TV5f1) is mostly stratified and usually has only a minor lower unit without stratification. TV4e1 and TV5e1 are asymmetric with steeper northern flanks and less inclined southern flanks (Fig. 12). In contrast to the northern valleys (TV1-TV3), TV4 and TV5 are underlain by both inclined Palaeogene (SU2 and SU3) and Cretaceous (SU1) sediments and thus incised into different strata (Fig. 12). The dip of the southern flank of TV4 even matches the dip of the underlying strata (Fig. 12). Fig. 5 illustrates that tracing TV4 and TV5 is ambiguous compared with the northern valleys, despite the dense line spacing. However, the overall incision direction is roughly towards the southwest.

Glaciotectonic structures in the subsurface
Large thrust sheets dominate the shallow part of the subsurface below the Holocene and upper Pleistocene sediments in the northwest (Figs 5, 7, 8, S1) and minor thrust sheets are found just south of TV2 (Figs 9 and S2). The thrust sheets display variable thicknesses between approx. 50 and 100 m and lengths varying between 300 and 800 m depending on their location (shorter in the southeast, longer in the northwest). All thrust sheets are based by undisturbed strata (e.g. Fig. 7) acting as a décollement. The northern extent of thrust sheets is limited by the broad incision of TV1.

Distribution of tunnel valleys
In their deepest part, all tunnel valleys observed in this study incise into well-stratified Neogene and late Palaeogene sediments, which were likely deposited within the Eridanos delta sequence (Overeem et al., 2001;Thöle et al., 2014). Two underlying salt structures, probably reactivated during the Alpine Orogeny (Ziegler, 1990), led to the deformation of Palaeogene and Neogene units SU1-SU4 (Fig. 3): salt pillow 'Karla' updomed the subsurface in the centre and salt pillow 'Rochelsteert' updomed the subsurface in the southeastern part of the study area (Lokhorst, 1998). This uplift resulted in glacial erosional truncation of SU2-SU4 during the Quaternary by massive ice sheets and subsequent deposition of Quaternary sediments. Here, tunnel valleys incised into SU2b and SU3 at the flanks of the salt pillow 'Karla' (Figs 9 and 11).
The subsurface in the southeast shows steeper updoming of sediments (Fig. 3). Tunnel valleys were incised into SU1, SU2 and SU3 in this area (Fig. 12). We consider that the different strata led to asymmetric valley flanks in the southern valley complex (Fig. 11 b), indicating that the difference in underlying strata has influenced the morphology and likely the incision direction of the buried valleys.
However, we are unable to unambiguously validate that the valleys follow a subsurface morphological trend or an associated geological phenomenon, e.g. salt diapirs, as suggested by other authors for other tunnel valleys (Hinsch, 1979;Ehlers and Linke, 1989;Kristensen et al., 2007). Limited by the spatial coverage of our data set, we find that the lateral spacing of the tunnel valleys is more or less constant and roughly coincides with salt structures in the subsurface. Accordingly, the morphology of the tunnel valleys may be primarily controlled by the lithological variations of the substratum as proposed for other areas of the North Sea (Praeg, 1996), but we can neither validate nor rule out substratum controls on the broader drainage pattern.

Internal structure and morphology of tunnel valleys
The different number of subsequent incisions inside TV1-TV3 likely indicates separate phases of cut and fill, a complex ice margin or a combination of both. TV1's incision phases are clearly delineated by individual erosional surfaces that indicate cut-andfill structures. A deep initial incision (TV1e1) into the Miocene sediments (SU3) was followed by minor incisions close to either flank (TV1e2-4). The latest incision (TV1e4) widened the tunnel valley complex to its present lateral extent of approximately 4.5 km (Fig. 7). A similar late widening of tunnel valleys has been described for other tunnel valleys in northwest Europe (Piotrowski, 1994;Jørgensen and Sandersen, 2006). Overall dimensions are consistent with other tunnel valleys, which typically range from less than 100 m up to 4-5 km (Woodland and Woodland, 1970;Le Heron et al., 2004;Jørgensen and Sandersen, 2006;Kristensen et al., 2007;Kehew et al., 2012;van der Vegt et al., 2012;Stewart et al., 2013;Livingstone and Clark, 2016). Common depth-towidth ratios have been postulated to be 1:10 (Ghienne and Deynoux, 1998;Gibling, 2006), which is the same order of magnitude as the tunnel valleys investigated here; namely, ratios between 1:10 and 1:20. Flank angles reported in this article of 6 to 19°fall well within the range reported by Kristensen et al. (2007) for other Pleistocene tunnel valleys in the North Sea. As neither the beginning nor the end of the tunnel valleys are located in the study area, we are unable to provide slopes for these locations to compare with values given in the literature (e.g. Kristensen et al., 2008;Moreau et al., 2012). However, we do observe morphological undulations along the thalweg (Fig. 6), which has been described as characteristic feature for many tunnel valleys (Huuse & Lykke-Andersen, 2000;Kluiving et al., 2003;Kristensen et al., 2007;Stewart et al., 2013). TV2 and TV3 additionally display a pronounced V-shape in their lower parts, which is usually interpreted as inner gorges (Huuse and Lykke-Andersen, 2000;Kluiving et al., 2003;Kristensen et al., 2007;Stewart et al., 2013;Jansen et al., 2014). The southern tunnel valley complexes (TV4 and TV5) show the shallowest incision depths. In addition, TV4 and TV5's internal structure does not show an inclination of reflectors and the stratification of the infill (TV4f1 and TV5f1) is well preserved on most seismic sections. Investigation of the southern end of the study area is limited by strong acoustic blanking due to shallow gas, which precludes tracing TV4 and TV5 for more than a few kilometres. The extensive acoustic blanking in our data renders finding the true extent of the southern tunnel valley complex impossible from these data.

Tunnel valley orientation
The incision direction of the tunnel valleys is roughly east-west to northeast-southwest but there are clear differences in the incision direction of individual tunnel valleys. For example, incision directions of the crossing valleys TV2 and TV3 clearly differ. The incision direction of TV3 is not influenced by the earlier incision of TV2 (Fig. 5). This lack of reuse of TV2's thalweg indicates that TV2 was already completely filled with sediment during TV3's formation. Consequently, both tunnel valleys document independent incisions (TV2e1 and TV3e1) and a time of sedimentation in between. TV3's average incision depth is significantly deeper than TV2's and the incision is oriented more westward (Fig. 5), representing a change both in magnitude and direction of incision. These differences possibly reflect a change in the ice margin's orientation during the formation of TV3 compared with the formation of TV2. Without precise age control, we are unable to attribute the individual incisions to specific ice advances but we consider it likely that individual incisions were formed during separate ice advances (Winsemann et al., 2020). Fig. 12 illustrates how different sediments offer different resistance against incision. The difference in resistance results in steep northern flanks and gentle southern flanks, whose dip compares to the dip of the substrate (Ehlers and Linke, 1989;Jørgensen and Sandersen, 2006;Stackebrandt, 2009). The laterally variable substratum in the south influences the incision direction and depth, which is reflected in the complex incision pattern and dip of the flanks of the southern valley complexes TV4 and TV5 (Fig. 5). In contrast, TV1-TV3 incised into mainly Neogene strata with little lateral heterogeneity and overall reduced erosion resistance, which led to higher incision depths, wider valleys and better traceability (Fig. 5) assuming similar magnitudes of drainage and duration.
Comparison with Lutz et al. (2009) The spatial interpolation of tunnel valleys based on our dense seismic grid led to a 3D model for TV1-TV5, which we used to define valley flanks (Fig. 5). Lutz et al. (2009) presented the most recent compilation of the distribution of tunnel valleys in the German sector of the North Sea based on a number of surveys. Fig. 5 provides a comparison with the compilation of Lutz et al. (2009). We find that former interpretations of incision directions interpreted before deviate significantly from our results; in parts, incision directions differ by up to 90°, which has implications for the configuration of the former ice sheet margin covering the study area. Furthermore, the width of tunnel valleys differs significantly from earlier studies (Lutz et al., 2009), and tributary tunnel valleys are now visible. Our data reveal greater sinuosity in the west, whereas sinuosity is low or absent in the east (Fig. 5). The comparison also reveals that a GTC has been previously classified as tunnel valleys (see grey-shaded area in Fig. 5). Here, Lutz et al. (2009) postulated a broad tunnel valley, whereas our data clearly show glaciotectonic disturbance of the strata (Figs 7 and 8). Only our new densely spaced data set allowed us to define the true distribution and directions of the tunnel valleys. The mapping of Lutz et al. (2009) was based on a broad line spacing and assumed mainly southward ice advances (Figs 1 and 5). This comparison shows that a very close line spacing or 3D seismic data sets are vital to examine the spatial characteristics of tunnel valleys in detail.

Tunnel valley formation
The glaciogenic unconformity terminates towards the west of the study area represented by an upward dipping strong reflector, which spatially coincides with the onset of glaciotectonic thrusting (Fig. 8). The geomorphological characteristics of this unconformity as well as the spatial correlation with glaciotectonic thrusting suggests that this unconformity represents the onset of ice-marginal thrusting, thus reflecting a temporary termination zone of a previous glaciation. To the west, no pre-Holocene glacial erosional unconformities are identified.
Towards the east, multiple glacial erosional unconformities between late Miocene and Holocene strata indicate several phases of glacial erosion. These possibly reflect periods of sea level low stands during glaciations. Based on the assumption that the Elsterian glaciation was the first one to cover vast areas of the North Sea (Ehlers et al., 2011;Winsemann et al., 2020), we suggest that abrasion by an extensive ice lobe during this glaciation led to the deepest glacial erosional unconformity, which connects to TV1-TV3 flanks (Fig. 4). Thus, the deep tunnel valleys incised into Miocene strata and truncating the glacial erosional unconformity likely formed during an Elsterian ice advance.
The separation into V-and U-shapes of the tunnel valley base is unambiguous for TV1-TV3. This morphology and the undulating thalwegs (Fig. 6) would be unexpected to result from fluvial erosion. From the subsurface data alone, it is difficult to identify whether a pre-glacial depression existed that was occupied by an ice stream, and whether subsequent subglacial drainage led to the meltwater erosional surfaces. Any morphological markers of this depression would be eroded by the subsequent incision of tunnel valleys. However, taking into account the undulating thalwegs, steep-sided V-shape inner gorges and overdeep meltwater erosion into significantly older undisturbed strata, we consider it unlikely that the tunnel valleys initially formed in a fluvial regime (Huuse and Lykke-Andersen, 2000;Jansen et al., 2014). Based on the morphological evidence, we therefore opt for a subglacial formation during high-pressure drainage conditions. TV2 does not show indications of multiple incisions. Consequently, we assume that TV2 may include the complete and undisturbed sedimentation cycle of a filled and buried tunnel valley, assuming that no subsequent drainage event has completely removed a previous infill. The clear separation of a lower V-shaped (inner gorge) and an upper U-shaped unit with the absence of an erosional unconformity (Fig. 11) suggests a single ice advance for its formation.
The southern shallow and more complex tunnel valley systems (TV4 and TV5) are separated from their deeper counterparts (TV1-TV3) by at least one glacial erosional unconformity throughout the study area. This separation suggests that an ice sheet retreat and subsequent ice sheet advance separate their formation (Moores, 1989;Praeg and Long, 1997;Kristensen et al., 2007;Sandersen et al., 2009;Stewart and Lonergan, 2011;Moreau et al., 2012). If we consider the partial lack of internal erosional stages in TV2 and TV3 but assume an Elsterian origin (Lutz et al., 2009), then the shallow glacial erosional unconformities and TV4 and TV5 may have formed during a later glaciation-deglaciation cycle. Previous studies have shown that the Saalian glaciation probably reached the study area (Lambeck et al., 2006), but the Weichselian ice sheets did not (Ehlers, 1990). Therefore, we suggest that the shallow glacial erosional unconformities and connecting valleys (TV4 and TV5) have likely formed during a post-Elsterian glaciation (Winsemann et al., 2020).
Time control for the formation of individual tunnel valleys can only be achieved by direct sampling and dating, but such data are not available. The relative timing of incision does not exclude multiple incisions during a single glaciation. Assuming the formation of all tunnel valleys during a single glaciation-deglaciation cycle would imply that the Saalian glaciation did not produce a prominent glacial erosional unconformity. However, we consider it unlikely that the Saalian ice sheet reached the area without producing a glacial erosional unconformity.

Orientation of former ice margins
The distribution and direction of tunnel valleys provide indications on the configuration of ice sheet margins. Woldstedt (1922) and Huuse & Lykke-Andersen (2000) have established that tunnel valleys usually form within ice sheet margins and that their orientation is usually parallel to ice flow (see van der Vegt et al. (2012) and references therein). Following this assumption, the incision directions of tunnel valleys TV1-TV5 indicate ice sheet advances oriented from north to south and northeast to southwest. However, these derived directions conflict with the north-northwest thrust direction of the GTC between TV1 and TV2, which suggests either a change in the orientation of the ice margin between the formation of the GTC and the tunnel valleys, or that the tunnel valleys did not flow perpendicular to the ice sheet margin.
Assuming the formation of tunnel valleys parallel to ice flow during a single stagnating ice advance, the distribution of tunnel valleys indicates an arcuate ice margin (Fig. 13), in analogy to the findings of Benediktsson et al. (2015). In this scenario, the GTC formed at the northern to western margin of the ice lobe. The thrusting of sheets towards the northwest indicates that glaciotectonic stress acted along an arcuate ice margin resulting in a change of thrust direction from north to west (Fig. 5). The subsequent multi-phase formation of TV1 may be attributed to drainage of meltwater beneath dynamic ice sheet margins (Fig. 13).
This scenario supports the assumption of an arcuate ice lobe advancing towards the west or southwest. However, it is likely that the GTC and the tunnel valleys formed during a number of ice advances, which matches well both the incision direction of the valleys and the thrust directions of the GTC (Fig. 13). Without dating, we are unable to assign an exact age to the GTC. However, based on the incision of valleys, it is clear that the GTC formed before the tunnel valleys. Assuming the bigger and deeper tunnel valleys (TV1-TV3) formed during the Elsterian glaciation (Lutz et al., 2009;Lang et al., 2012), we suggest an early Elsterian or pre-Elsterian age for the GTC in accordance with Winsemann et al. (2020). The morphological features thus likely document a number of ice advances into the southeastern North Sea with different orientations (Fig. 13).

Conclusions
This study imaged tunnel valleys in the German sector of the southeastern North Sea in very high detail. These data shed light on the complex formation of tunnel valleys and an associated GTC. Based on our results and previous research in the area, we assume that the deep tunnel valleys in the study area (TV1-TV3) are older (formed during the Elsterian glaciation), whereas the shallower tunnel valleys (TV4-TV5) are younger (possibly formed during the Saalian complex). However, in the absence of age control, we are unable to assign definitive absolute ages to the structures apart from their relative time of formation.
We conclude that pressurised subglacial meltwater was responsible for their formation. The complex fill and incision directions indicate that a number of ice advances covered the study area. The direction of incision as well as the direction of thrusting in the GTC advocate for different orientations of several ice advances during a previous glaciation of the North Sea. An arcuate shape of the ice margin is likely the reason for the thrust directions transitioning from north to west.
The study shows that the classic model of few continentalscale advances during former glaciations does not hold true. Rather it seems that ice sheets have been highly dynamic and their margin changed shape and direction over short time scales, resulting in complex incision patterns of buried tunnel valleys and glaciotectonic deformation.

Supporting information
Additional supporting information may be found in the online version of this article at the publisher's web-site. Figure S1. Compilation of seismic sections imaging tunnel valley 1 (TV1) from east to west on multiple parallel lines. TV1 is kept in the middle of the plot. Inset gives location of the seismic section. Figure S2. Compilation of seismic sections imaging tunnel valley 2 (TV2) from east to west on multiple parallel lines. TV2 is kept in the middle of the plot. Inset gives location of the seismic section. Figure S3. Compilation of seismic sections imaging tunnel valley 3 (TV3) from east to west on multiple parallel lines. TV3 is kept in the middle of the plot. Inset gives location of the seismic section.
streamlining of the article and provided additional information on the study area. SK and KS led the data acquisition, which was aided by AL and DU. All authors have contributed with valuable comments, discussions and guidance for the improvement of the article.
Conflicts of interest-The authors declare no conflicts of interest.