Buried late Quaternary channel systems in the Danish North Sea – Genesis and geological evolution

Quaternary Science Reviews 223 (2019) 105943

Contents lists available at ScienceDirect

Quaternary Science Reviews journal homepage: www.elsevier.com/locate/quascirev

Buried late Quaternary channel systems in the Danish North Sea e Genesis and geological evolution
sik Prins*, Katrine Juul Andresen Lasse Te SeisLab Aarhus, Department of Geoscience, Aarhus University, Høegh-Guldbergs Gade 2, 8000, Aarhus C, Denmark

a r t i c l e i n f o

a b s t r a c t

Article history: Received 1 May 2019 Received in revised form 10 September 2019 Accepted 11 September 2019 Available online 19 September 2019

The Danish North Sea has been subject to extensive seismic surveying (both 2D and 3D) during the past four decades. The 3D seismic data from the area has mainly been used for hydrocarbon exploration, but in the last decade, it has also been utilized for investigating the Quaternary sediments and especially buried incised valleys. The erosional and depositional processes involved in forming such valleys are important in paleolandscape reconstruction, but are often difficult to determine, since sediment deposited by flowing water, displays ambiguous geometric structures. In this study, we combine the analysis of infill geometries, from high-resolution 2D seismic data, with plan-view morphology analysis, from 3D seismic data, and lithological data from borehole reports, in order to constrain the origin of buried channels in the Danish North Sea. This data integration makes differentiations of a subaerial vs. subglacial origin more robust. The study focusses on analyzing a proglacial to fluvial channel system, which was initiated as a subglacial tunnel valley system during the Last Glacial Maximum, and subsequently evolved as a fluvial system after deglaciation. The channel system was active in a marsh environment during a period of subaerial exposure in the late Weichselian to early Holocene. During the Holocene transgression, the system was gradually flooded leaving estuarine tidal deposits in the underfilled channels. Once the area was completely flooded, the channel system was buried under marine sands. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Quaternary Paleogeography Geomorphology North sea River valleys Tunnel valleys Seismic Doggerland

1. Introduction Buried valleys are a very striking morphology in the Quaternary succession of the North Sea (Fig. 1) (e.g. Van der Vegt et al., 2012). They develop under a series of different geological conditions and analysing such valleys is thus important in order to reveal the Quaternary geological evolution in the North Sea. The lack of lithological samples and age dating from the buried valley infill has however generally made this analysis difficult. This study combines available data (3D seismic data, 2D seismic data and geotechnical data) from the Danish North Sea to describe a buried channel system from seismic facies analysis and morphological analysis of seismic time slices. During the Quaternary, the North Sea area was characterized by cyclic changes between glacial and interglacial conditions (Knudsen and Sejrup, 1993). Three major glaciations have

* Corresponding author. E-mail addresses: [email protected] (L.T. Prins), [email protected]. dk (K.J. Andresen). https://doi.org/10.1016/j.quascirev.2019.105943 0277-3791/© 2019 Elsevier Ltd. All rights reserved.

traditionally been accepted for the mid to late Quaternary: the Elsterian (Marine Isotope Stage (MIS) 12, ca. 480-410 ka), the Saalian (MIS 10-6, ca. 370-135 ka) and the Weichselian (MIS 5d-2, ca. 109-12 ka) (Ehlers et al., 2011; Houmark-Nielsen, 2011). The central parts of the North Sea Basin, which are the focus of this study, are characterized by Quaternary sediments reaching more than 1000 m in thickness (Ottesen et al., 2014). Subsidence along the basin axis and uplift at the margins can be observed on regional seismic profiles (Huuse and Lykke-Andersen, 2000) and the Quaternary subsidence generally follows the thermal subsidence pattern established during the Cenozoic in response to the Jurassic rifting episode of the central North Sea (Ahmadi et al., 2003; Westaway, 2017). The rapidly oscillating climate in the Quaternary (illustrated by variations in the d18O ratio from deep sea cores Fig. 2) resulted in different depositional environments over a relatively short period. In the early Quaternary (from around 2.75 Ma) the ice sheets were mainly restricted to the British and Norwegian landmasses (Ottesen et al., 2018) and glacial erosion delivered more than 1000 m of sediment into the Central North Sea Basin from the

2

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

Fig. 1. Bathymetry (EMODNET) and overview of the study area in the central North Sea. Black squares highlight the focus of the study area (1) as well as other studies from adjacent areas (2 e Hepp et al., 2017 and Coughlan et al., 2018, 3 e Müther et al., 2012, 4 e Cotterill et al., 2017). Mapped Quaternary valleys, are outlined in different colours, the system that we focus on is outlined in black. Cross cutting relationships suggest that other valleys in the area are older than the system described in this study. Interpretations of LGM ice extent variations are also indicated, as well as tilting lines and apparent position of the LGM forebulge as presented by Kiden et al. (2002). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

surrounding land areas. In the southern North Sea Basin, sediment input came mainly from the large Baltic/European rivers (Overeem et al., 2001). Conditions changed when grounded ice sheets moved into the North Sea Basin in the early Middle Pleistocene (0.6 Ma) (Bendixen et al., 2018), however Stewart (2009) reported tunnel valley formation in the witch ground basin dating back to before 0.78 Ma, as the oldest morphological evidence in the central North Sea. Ice sheets covered most of the North Sea during the past three major glaciations. During the Weichselian at the Last Glacial Maximum (LGM) (ca. 25-19 ka) however, a large part of the southern North Sea was ice-free (Carr et al., 2006; HoumarkNielsen, 2011; Hughes et al., 2016). There is a growing consensus that the tripartite ice sheet model is too simple (Graham et al., € se et al., 2012) and that 2010; Stewart and Lonergan, 2011; Bo there were several episodes of ice advance and retreat during each glaciation, thus making it uncertain whether there were three or more glaciations in the North Sea. In the late Weichselian and early Holocene, glacial lakes and marshes dominated the landscape, and there was exposed land

between Britain and the mainland Europe (the former Doggerland), that potentially hosted pre historic settlements (Coles, 1998; Fitch et al., 2005; Cotterill et al., 2017; Coughlan et al., 2018). Pronounced erosional channels in the upper 400 m of the Quaternary succession are primarily interpreted as subglacial tunnel valleys, although some have also been identified as subaerial river systems (Gaffney et al., 2007; Cohen et al., 2014; Van Heteren et al., 2014; Cotterill et al., 2017). Tunnel valleys are assumed to drain towards the ice margin, and are therefore used to interpret maximum ice sheet extent during glaciations (Wingfield, 1989; Moreau et al., 2012; Phillips et al., 2017). Discriminating subglacial valleys from fluvial or submarine channels based on seismic data alone includes analysis of channel dimensions such as width/depth ratio (Gibling, 2006), plan-view morphologies, thalweg variations - subglacial valleys typically show an undulating channel thalweg due to formation by pressurized meltwater (O’Cofaigh, 1996; Kristensen et al., 2007; Van der Vegt et al., 2012), and infill geometries such as up-ice dipping clinoforms that have described from subglacial tunnel valleys (Praeg, 2003; Kristensen et al., 2008). Many other infill geometries are

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

3

Fig. 2. Overview of the late Quaternary stratigraphy in the study area in relation to Marine Isotope Stages (Johnson et al., 1993). Thickness of individual units from the reference borehole (BH 550513.5). The depth of channels indicated in ‘Depositional environment’ is therefore not to scale. ‘Depositional environment’ for BH 550513.5 serves as a good example of the three defined regional geology units (Unit1-3). Lithostratigraphic units from adjacent areas (Cotterill et al., 2017; Coughlan et al., 2018) are presented as well as the foraminiferal stratigraphy by Knudsen (1985). NZGT: Nieuw Zeeland Gronden Terschellinger Bank Mb.

however ambiguous and do not indicate a specific environment, but rather a variety of different depositional environments. In addition, the infill of buried valleys is not necessarily related to the process that eroded the valleys (e.g. glacio-genic deposition) but may have formed much later. Some plan view morphologies such as scroll bars and oxbow lakes are indicative of fluvial processes and would not be expected in a subglacial setting. Subglacial valleys are often characterized by straight or slightly sinuous anastomosing plan-view morphologies (O’Cofaigh, 1996; Huuse and Lykke-Andersen, 2000; Stewart et al., 2013; Stewart, 2016) but can however also display a dendritic plan view pattern as observed by Praeg (1996). The presence of past ice sheets continue to influence the North Sea region today through glacio-isostatic adjustments (Steffen and Wu, 2011). There are no glacial rebound data from the central North Sea, but models based on data from Fennoscandia collected from 1892 to 1991 (Ekman, 1996) show that relative subsidence due to glacial rebound is expected in the central North Sea. Furthermore, models generally predict uplift close to the land areas in both the Norwegian, British and Belgian sectors, and predict subsidence in

the central North Sea (Emery and Aubrey, 1985; Ekman, 1996; Bradley et al., 2009). These models calculate uplift rates for the land areas based on onshore data points. A degree of uncertainty on the isostatic adjustment would therefore be expected for the central North Sea. The climatically induced sea level rise that formed the Yoldia Sea (ca. 15 ka), was not sufficient to submerge the entire area between Europe and Britain and parts of the area were subaerially exposed until the Holocene transgression at around 7 ka (Fitch et al., 2005; Jensen et al., 2005; Cotterill et al., 2017). Age constraints for the Quaternary in the North Sea are limited, particularly for the late Quaternary (Lonergan et al., 2006; Stewart and Lonergan, 2011; Hughes et al., 2016). For the Norwegian sector, Reinardy et al. (2018) have introduced amino zones based on amino acid and strontium measurements which, to some extent has improved the dating there, but reworking has resulted in some uncertainty on some of these dates. Amino stratigraphy does have some limitations in relation to reworking, temporal resolution and temperature history that can have a significant effect on the age estimates (Penkman, 2010). For the Danish sector, the absolute dating of channels is

4

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

not very precise, and relative dating based on seismostratigraphy, cross cutting relationships and correlations to the sparse well dates that are available (Knudsen, 1985; Knudsen and Sejrup, 1993) are the best tools available. The present study describes the late Quaternary geological evolution of a channel system located in the vicinity of the maximum ice extent during the LGM in the central Danish North Sea using extensive 2D and 3D seismic data and available borehole data interpretations (Fig. 1). The channel system is positioned close to the margin of the LGM ice extent and the uncertainty of the exact position of the LGM extent enables the system to be interpreted as either a subglacial system, a subaerial system or a combination of the two. The objective of this study is thus to analyse the channel systems-dimensions, geometries and morphology in order to interpret whether the system was formed in a subglacial or subaerial environment. 2. Data The seismic data from the Danish North Sea used in this study includes a 3D seismic time cube called the Contiguous-MigrationArea (from here referred to as Cube 1) publicly available (GEUS e Frisbee, N.D), a depth-converted 3D seismic cube called the Danish-Regional-Merge (from here referred to as Cube2) (provided by Total), six 2D regional sparker lines from the NS10 survey, publicly available (GEUS e MARTA, N.D), six high-resolution 2D seismic tie lines (collected in 2014 by Gardline Geosurvey ltd., provided by Total), and three 2D seismic site surveys (SS1-3) acquired in relation to platform constructions at the Adda, Tyra and Broder Tuck fields (Table 1 and Fig. 3) (Adda and Tyra data was provided by Total, Broder Tuck was provided by GEUS (The Geological survey of Denmark and Greenland)). The merged nature of the 3D seismic cubes are particularly visible on shallow time slices where the acquisition footprint show the different acquisition directions. Cube1 consists 11 merged 3D seismic surveys with a total of 4840 inlines and 4390 crosslines, with an inline and crossline spacing of 12.5 m. The data is zero phase (wavelet is symmetrical about zero time) with a downward increase in acoustic impedance represented by negative amplitudes (European polarity). It reaches a depth of 2000 ms TWT (two-way-travel time). For Cube1, the frequency range is between 6 and 120 Hz with two dominant frequency populations at 30 Hz and 55 Hz resulting in a vertical resolution, here defined as a quarter of the dominating wavelength (l/4) (Brown, 1999), of z16 m and 8.5 m and a horizontal resolution defined here as half the dominating wavelength (l/2) (Brown, 1999) of z32 m and 17 m, using a velocity of 1905 m/s (cf. Kristensen and Huuse, 2012). Cube2 consists of 13 merged 3D seismic surveys, with a total of 10465 inlines and 9487 crosslines, it has an inline and cross-line spacing of 12.5 m Cube2 is similarly zero phased with European polarity. It penetrates down to 1000 m. The dominating wavenumber is 0.03 m-1 corresponding to a wavelength of l ¼ 33 m and a vertical and horizontal resolution of 8 m and 17 m respectively. The 2D seismic data (regional lines, tie lines and site surveys) have a much better vertical resolution down to ca. 1 m for the

regional lines and 4e6 m for the site surveys (see Table 1). Geotechnical borehole reports are publically available (GEUS e MARTA, N.D) with an archive number (BHxxxxxx.x). Although the majority of the geotechnical boreholes in the study area are very shallow (<5 m), eight “deep” boreholes (22e80 m long) are available and they are the primary source for lithology and age information (Figs. 3 and 4). 3. Methods 3.1. Interpretative workflow Our methodology workflow (Fig. 5) includes three main steps. It combines the benefits of the good coverage from the conventional 3D seismic data with the high resolution from the 2D seismic data. Due to the generally poor resolution of the 3D seismic cubes and very poor resolution in the near surface parts, it is difficult to distinguish geological features from noise in 3D seismic crosssections except from the presence of very large channels (Step 1). The lateral changes in amplitude are however strong enough for smaller features to be visible in time slices. By combining 3D seismic time slices and 2D seismic sections, these smaller channels can be distinguished from the noise and mapped (Step 1). Water layer multiples are usually problematic when interpreting the shallow subsurface. Here, we have actively used the multiples to improve the interpretation of the lateral channel positions. Since the multiples have a larger stacking fold due to their apparent deeper stratigraphic position, the seismic image from the multiple has a better signal to noise ratio with less acquisition footprints than the primary reflections (Fig. 6). This improves the certainty of the mapping of the lateral extent of channels on 3D seismic time slices. The first and second multiples express a higher amplitude contrast across the time slice than the primary reflection (compare Fig. 6a, b and c) resulting in a much higher certainty of the position of the channels, but a higher uncertainty of the depth to the channel. In some cases this method enables the identification of internal features of the channels (such as bars) and in turn channel morphology (Fig. 7). Precautions should be made when interpreting geology using multiples as they can obscure the relative age interpretation. If channels are present in the same stratigraphic levels as the multiples, there is a risk of classifying the ‘true’ channels and the ‘multiple channels’ into the same age category, when they may be of completely different ages (see Figs. 6 and 7 for examples). The benefits of using 3D seismic data compared to 2D seismic datasets for mapping of geological structures such as channels, were highlighted, by Cartwright and Huuse (2005) and Kristensen et al. (2007), among many others - a dense 2D seismic grid supplemented with a good understanding of the geological evolution will, at best, allow for a simple interpolation of channel paths between seismic lines when mapping meandering channel systems. As a supplement to the spatial mapping from the 3D seismic time slices, the higher resolution 2D seismic data have been used to resolve infill geometries, seismic facies, dimensions and morphology of the channels (Step 2, Fig. 5). In Step 3, these

Table 1 Seismic frequencies and resolution for the different 2D seismic data (SS1-SS3 ¼ site surveys). See Fig. 3 for location of surveys. Data type

Seismic source

FSpectrum

Fdominant

Vertical resolution

Horizontal resolution

SS1 SS2 SS3 Regional lines Tie lines

Sparker Sparker Sparker Sparker Sparker

4e350 Hz 15e260 Hz 6e320 Hz 100e1000 Hz 50-600 HZ

100 Hz 75 Hz 120 Hz 500 Hz 150 Hz

5m 6m 4m 1m 3m

9m 12 m 8m 2m 6m

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

5

Fig. 3. Overview of the mapped channel system. The system is divided into four segments: the northern channels (green), the central channel (red), the southern channels (blue) and the eastern segment (turquoise). Measurement locations are indicated by m1-m20 (see Fig. 8 for values). Angle of confluence between channel segments are drawn (only for angles different from 90 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

interpretations are compared to geotechnical borehole reports and logs and, when necessary updated (Fig. 5), ensuring the most confident mapping of the buried channel systems. We associate the seismic facies with depositional environments in order to determine the origin of the system. This is done using the seismic geomorphological principles (Posamentier et al., 2007), comparison with previous studies describing similar channels from the southern North Sea (e.g. Hepp et al., 2017; Coughlan et al., 2018) and by using the available borehole data. 3.2. Depth conversion For a qualitative seismic-to-well-tie and simple depth conversion we have used an interval velocity (Vint) between 1600 m/s and 1650 m/s (cf. Hepp et al., 2017). We can then approximate the depth through the linear relation:

Depth½m ¼ Vint ½m=s  0:5  TWT½sz0:813  TWT½ms

(1)

3.3. Paleo-flow direction In order to interpret paleo-flow direction in the channels, several observations were used. These include: i) mapped planview morphology (Figs. 3 and 7); ii) measured angle of confluence between channel segments (Fig. 3); iii) measured width of channels (Fig. 8); iv) measured channel relief (depths) defined as

the relief from the shoulder to the bottom of the channel (Fig. 8); v) measured relative paleo-base level defined by depth from present day mean sea level to the base of the channels (Fig. 8). Present day sea surface was chosen for the reference datum because relative sea level changes can be regarded as a vertical parallel shift of a horizontal datum line, which means that the relative distances between the sea surface and the channel shoulder, and the sea surface and the channel base, will be constant over time. The present day sea surface as a reference datum is thus independent of deposition, erosion and tectonic movement, and is in our opinion, the most robust reference datum for measuring paleo-base level in buried valleys. The maximum difference in channel relief of the first order channels is 8 ms TWT (first order channels receive water from second order channels. Second order channels deliver water to first order channels and receive water from third order channels, and so on). This corresponds to a depth of 6.5 m, which is close to, but still more than the vertical resolution of the 2D lines included in this study (Table 1). The low resolution of the 3D seismic data has not allowed for the mapping of the channel thalweg. Instead, measurements from high-resolution 2D seismic lines and site surveys were made in order to determine the depth to the base of the channels (Figs. 3 and 8). This also means that an arbitrary line displaying the thalweg along its course could not be produced, which would have been a strong tool in assessing the conditions under which the system was formed.

6

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

Fig. 4. Regional geology as defined from regional 2D seismic lines correlated with borehole 550517.11 (a), SS3 borehole (b) and borehole 550513.5 (c). In d), a borehole transect across the study area showing lithologies is presented in (d) (for position of seismic lines and boreholes see Fig. 3 and 6). The regional geology is subdivided into three units (Unit 1e3). Please note that BH550513.5 is positioned 900 m from the seismic line.

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

7

Fig. 5. Three-stepped workflow chart describing the iterative process of integrating 3D seismic data with high-resolution 2D seismic data and borehole data.

4. Observations and results

conditions like today.

Several buried channel systems of varying relative burial depth and age have been identified in the study area (Fig. 1). The focus of this study is however on a single large coherent system buried (ca. 3e13 m) below the seafloor. In order to understand the geological context of the channel system we also describe the regional geology (Fig. 4).

4.1.2. Unit 2 The top of the Unit 2 equalling the base of Unit 1 is erosive, while the base of Unit 2 is an undulating surface (Fig. 4). Unit 2 comprises subparallel reflections of medium amplitudes (Figs. 4 and 9) and is dominated by clay in the available boreholes (Figs. 2 and 4c). The thickness of the unit varies between 8 and 16 m but some of the boreholes terminate within the unit (e.g. BH550407.27, Fig. 4c) indicating that the unit is thicker in some areas. According to the geotechnical report BH550513.5, this unit is inferred to be of Weichselian age. The unit corresponds to unit S2 of Knudsen (1985), where it is described as being non-marine, presumably of Eemian to Weichselian age, based on relative dating.

4.1. Regional geology Based on the high-resolution 2D seismic data the shallow sequence (topmost 100 m) has been subdivided into three seismic units (Unit 1e3). This subdivision is based on seismic appearance and lithology from the available boreholes (Fig. 4). Seismic-to-welltie is based on borehole (BH) 550517.11 (Fig. 4a), SS3 borehole (Fig. 4b) and BH550513.5 (Fig. 4c). BH550513.5 is situated 950 m from the seismic line and although this is a relatively large tie-in distance, the depth conversion of the TWT for the main seismic surfaces (using Eq. (1)) matches well to the depths of the lithological boundaries in the borehole, hence justifying the correlation. In Fig. 4b as well as time slice analyses (Fig. 6) a deeper and larger channel system is visible. Spatially it partly overlaps with the shallow channel system (see Fig. 1). This deeper system as well as parts of the system described in this study has previously been described by Müther et al. (2012). They interpreted the deep system to represent a subglacial tunnel valley of Saalian age, and further interpreted the channel system described in this study to be directly linked to the deep system. Our observations however clearly show that the two systems are stratigraphically separated (see Figs. 4b and 6) and that the shallow channel system therefore must be younger than the Saalian. Knudsen (1985) defined three stratigraphic zones (S1-3) based on foraminifera for the Skjold area, which is the same area in which BH550416.7 is located (Fig. 6). This stratigraphy corresponds to the regional geology units defined in the following (Fig. 2). 4.1.1. Unit 1 Unit 1 comprises the uppermost 3e20 ms TWT (2.5e16 m) below the seafloor and is generally characterized by highamplitude and continuous parallel sub-horizontal reflections. It is thickest close to the central channel. In the boreholes (Fig. 4d) Unit 1 varies from 3 to 13 m in thickness and consists of marine sands with some silty intervals and peat layers towards the base (BH550513.5) (Figs. 2 and 4c and d). Channels are cut into the erosional surface that represents the base of Unit 1. The lower part of Unit 1 thus is characterized by infilling of the eroded channels (Fig. 4a). The infill will be described further in sections 4.3 and 5.1. This unit corresponds to the S1 unit of Knudsen (1985) (Fig. 2), which is of Holocene (Flandrian) age, and represents a transition from non-marine over shallow brackish conditions to marine

4.1.3. Unit 3 Unit 3 exhibits a more chaotic reflection pattern compared to Unit 1 and 2 (Fig. 4a) and is in its lower parts characterized by large scale incisions (e.g. Fig. 4b). None of the available boreholes reach the bottom of the unit, and the base of Unit 3 is therefore not defined in the boreholes. The upper part of Unit 3 is predominantly sandy (Figs. 2 and 4c). According to the geotechnical report BH550513.5, the unit is expected to be of Eemian age but may be older. Through the correlation between the general geology units and the stratigraphy (S1eS3) presented by Knudsen (1985) we find that the thickness and lithology of Unit 1 corresponds to S1, which is a transitional unit marking the change from non-marine to marine conditions, this unit is definitely of Holocene (Flandrian) age. S2 that consists of Eemian - Weichselian deposits at a depth of 38e46 m, comprised of a silt layer followed by a sandy layer. This corresponds to the lower part of Unit 2 in BH550416.7 (Fig. 9) and is associated with the Eemian interglacial. This puts the channel in a stratigraphic window between Eemian and the Holocene, entailing a Weichselian origin of the system if it was initiated subglacially and an early Holocene age if it was initiated subaerially. 4.2. Morphology of the channel system The Unit 1 is thickest around the central parts of the channel system, meaning that there was a topographic low in that area when the system was active. This has had an influence on the morphology of the system, and has enabled the subdivision of the system into four segments based mainly on catchment but also on morphology (Figs. 3 and 4): i) A central channel (red) that is more than 40 km long and is comprised of several semi-straight segments with sinuousities ranging from 1.03 to 1.3. Measurements of the depth to the bottom of the channel show variations along its course, but a dominant dip direction towards the west (Fig. 8). ii) Several large, up to 18 km long, relatively straight NNE-SSW

8

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

Fig. 6. Time slices at different TWT in 3D seismic cube C1 showing the same channel system (see Fig. 1 for position). a) Displays the most correct vertical position of the system at 60 ms TWT (ca. 56 m below the seafloor), but b) and c) (first and second multiple) resolve the lateral extent of the system in greater detail.

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

9

Fig. 7. Time slice image at 96 ms TWT for the eastern channel segment displaying different plan view morphologies typical for fluvial systems (black outline). The apparent crosscutting relationship with an older channel is a result of the water layer multiple that images the channel system we describe. For location see Fig. 6.

striking northern channels (green), that act as large tributaries to the central channel, with a catchment situated north of the central channel. Based on only a few measurement points these channels have a thalweg dip towards the central channel (Fig. 8). iii) Two southern channels (blue) that also act as tributaries to the central channel, and both originate from one or two catchments south of the central channel. The southern channels show varying thalweg dip directions (Fig. 8). iv) An eastern segment comprised by smaller channel segments (turquoise). The plan-view morphologies of the entire channel system are best characterized as dendritic with three hierarchal orders of channels above seismic resolution (Figs. 3 and 8). The northern channels appear more straight and parallel than the rest within all orders of channels. The angle of incident of branching channels varies between 31 and 110 , with most angles falling into three populations one around 35e55 , one around 70 , and finally a large population around 90 degrees (Fig. 3). The geometric variations are presented in Fig. 8. For the central channel, depth to the bottom of the channel varies between 83 ms and 94 ms TWT (67.5e76 m) with the deepest parts towards the west (m1-m6 Figs. 3 and 8). The channel relief (erosional depths) within the central channel varies from 13.5 to 33 ms TWT (11e27 m) also deepening towards the west (m1-m5 Figs. 3 and 8). The central channel is up to 300 m wide, whereas the northern channels typically are only 150 m wide with bottom depths of 77e84 ms (63e68 m) and erosional depths of 11e21 ms (9e17 m). The width/depth ratio range is from c. 9.5e31.7 for the central and southern channels, with higher values (20.2e40.1) for the northern channels averaging 26.6 (Fig. 8). Towards the east, the central channel connects to the eastern segment, where a series of planview features resembling scrollbars, meander bends and cut-offs are present. Oxbow lakes and abandoned channels can also be recognized (Fig. 7). 4.3. Infill characteristics The infill of the channel system is described in order to identify possible diagnostic fluvial or subglacial geometries.

4.3.1. Central channel The infill variations of the central channel described from 2D seismic lines is shown in Figs. 9b, 10b and 11. Towards the east (Fig. 10b) the central channel splits into two smaller channels infilled by two different seismic facies. The infill of the northern part of the central channel is characterized by very high-amplitude dipping reflectors that fill the entire channel (Fig. 10b). Where present, these high-amplitude reflections define the base of Unit 1. Further west (Fig. 9b), the central channel changes character to a more well-defined erosional channel with distinctive infill reflections (partly due to better quality of the seismic data). Here, a clearly defined erosional base, consisting of three troughs, progressively deepening towards the east, is observed (Fig. 9b). In the shallowest trough, there is a series of imbricated reflections dipping towards the centre of the channel. The two other troughs contain concave reflection sets with transparent infill. The infill of all three troughs is truncated by a high-amplitude erosional surface (E1). A second erosional surface (E2) defines the top of the channel infill. Furthest to the west (Fig. 11), the bottom of the central channel becomes more symmetric and the infill is characterized by semihorizontal subparallel reflections (Fig. 11a). There is a characteristic bulge of approximately 5 ms TWT in height (4.5 m), which changes from being asymmetric in Fig. 11a to more symmetric in Fig. 11b. The distance along the central channel from Fig. 11a to b is ca. 4 km indicating the minimum length of the bulge, assuming it is coherently continuous between the two seismic lines. The overburden appears undisturbed by the bulge. A minor velocity push down (2e3 ms TWT) can be observed immediately below the channel in Fig. 11a. 4.3.2. Southern channels Transects of the southern channels are shown in Fig. 9c, d and 10b. The channel in Fig. 9c displays an erosional asymmetric bottom, steeper on the eastern flank, which consists of multiple troughs similar to the central channel (Fig. 9b). In Fig. 9c imbricated reflections dip towards the channel centre on the western flank, while the central part is filled by concave reflections. The entire channel is truncated by a high-amplitude erosional surface which correlates to the erosional surface of the central channel (E1). Three distinct amplitude anomalies occur within the erosional surface.

10

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

Fig. 8. Schematic drawing highlighting the geometrical variations of the channel system based on 20 measurement points (m1-m20). See Fig. 3 for position of the measurement points.

These are likely to represent polarity reversals possibly caused by gas accumulations in the sediments. The channel in Fig. 9d is much narrower (150 m wide) compared to the central channel (z300 m) (Fig. 9b). The infill reflections all dip towards the east and are truncated by the same erosional surface as observed elsewhere. Above the erosional surface downlaps are observed. These downlaps are most likely related to reworking of the sand unit above. Towards the east, the southern channels become wider (Fig. 10b) and appear to consist of an upper unit that is uniform and comprised of low-amplitude subparallel reflections overlying a lower unit displaying high amplitude dipping reflections. 4.3.3. Northern channels The northern channels are described using two seismic lines that intersect the channel. One intersection towards the north (Fig. 10a) and one in the south (Fig. 10b), where the southern and central channels converge. The cross-sectional expression of the northern channel ranges from asymmetrical with a steep western flank (Fig. 10a) to smaller

individual symmetrical channels (Fig. 10b). The infill is subdivided into an upper and a lower sequence. The lower sequence comprises high-amplitude semi-chaotic reflections, with imbricated reflections along the flanks of the channels (Fig. 10a and b). The upper sequence consists of low-amplitude horizontal reflections onlapping the flanks of the channels. The top of the infill is truncated by an erosional surface in a similar to the surface described for the central channel and the southern channels. 4.3.4. Eastern segments Only one 2D seismic line intersects the eastern channel segments. The quality of the seismic line in this area is poor, but the infill facies generally resembles the infill facies visible from Fig. 10. Despite the poor quality, the base reflection of Unit 1 could be traced from the central channel area to the area of the eastern channel segments confirming that the eastern channel segments belong to this system. As the entire system is incised into the same horizon, we assign the same age interval (Weichselian-early Holocene) to all the tributaries. However, we are well aware, that channel systems like

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

11

Fig. 9. a) High-resolution 2D seismic tie line (see Fig. 3 for position) intersecting the channel system at three places, revealing distinct geometries and seismic facies of the infill. The segment in b) shows lateral migration or multiple generations of internal channels. Similar morphologies can be seen in c) and d). e) Displays a synthetic log of the borehole 550416.7 located 300 m south of the seismic line.

this are time transgressive, and that the different parts of the system may not have been active simultaneously. Better age constraints are needed in order to determine the exact temporal evolution of the system.

the channels above the erosional surface.

4.3.5. Seismic facies of the infill We have classified the infill of the channel system into three seismic facies based on the above described characteristic differences in reflection patterns, amplitude and continuity (Figs. 9e11). Facies 1 (green in Figs. 9e11) is characterized by narrow bands of very high-amplitude subparallel reflections. This facies is typically observed between- and at the base of the channels in the eastern and northern parts of the system. Facies 2 (yellow) is defined by medium to high-amplitude semi-chaotic imbricated to sub-parallel reflection patterns. Internal erosional or troughshaped reflections also occur. Facies 3 (purple) comprises transparent to low-amplitude parallel to subparallel reflection patterns that onlap the channel sides and typically occur in upper parts of

5. Discussion In the light of the known geological development of the central and southern North Sea region during the latest Quaternary (Huuse and Lykke-Andersen, 2000; Ottesen et al., 2014, 2018; Hughes et al., 2016), the late Weichselian to early Holocene age of the channel system invites for a number of origins to be considered. These include subglacial channels and subaerial channels formed in either proglacial, fluvial or tidal environments. Key observations related to i) infill facies and depositional environments, ii) planview morphologies iii) channel geometries, iv) stratigraphic (relative age) and spatial position (subglacial or subaerial) are discussed in the following.

12

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

Fig. 10. Regional sparker lines crossing the channel system in the east and north (see Fig. 3 for position). a) Intersection, without and with interpretations, of one of the northern tributaries, where channel infill is subdivided into two seismic facies. Note the lateral variations in channel bottom depth indicating some degree of lateral migration. b) Seismic line crossing the area where the central, northern and southern channels converge. See text for further descriptions.

5.1. Depositional environments from infill facies To address the origin of the channel system we associate depositional environments for the interpreted seismic facies and regional geology units. Of the available boreholes only one (BH550416.7, from the Skjold area) intersects the channel system (Fig. 3). Unfortunately, there is no direct tie to a high-resolution 2D seismic line, but the plan-view morphologies (Fig. 3) show that BH550416.7 is situated centrally in one of the southern channels approximately 300 m south of the section shown in Fig. 9a and c. BH550416.7 comprises four units (Fig. 9e). i) A sandy upwards coarsening unit in the top (0e13.4 m below seafloor (mbsf)), which corresponds to Unit 1 comprising the marine sand unit deposited during and after the Holocene transgression (Knudsen, 1985; Hepp et al., 2017; Coughlan et al., 2018). ii) A clay unit (13.4e33.5 mbsf) consisting of 4 m of heavy (plastic) clay followed by 16 m of horizontally laminated clay with silt streaks. This unit comprises the channel infill iii) an alternating silty sand and clay succession (33.5e51.0 mbsf), that represents Unit 2. iv) A sandy unit in the base (33.5e80.1 mbsf) representing Unit 3.

5.1.1. Seismic infill facies 1 Seismic Facies 1 is prevalent in the bottom of the northern channels and in areas between closely spaced channels. In some boreholes (550513.5 and SS3 borehole) peat layers are found at the base of Unit 1, as a transition between the sand dominated Unit 1 and the clay dominated Unit 2 (Fig. 4). Hepp et al. (2017) and Coughlan et al. (2018) describe similar bundles of high-amplitude reflections at the bottom of the widest parts of their channels to be sequences of silt and peat deposited in a marsh environment. On this basis we interpret Facies 1 to represent interchanging peat and clay layers created mainly in a marsh environment. 5.1.2. Seismic infill facies 2 Seismic Facies 2 is interpreted on the basis of both infill geometries and the shape of the channels. The presence of imbricated reflections can be explained by slumping of channel edges, lateral migration of the channel, multiple generations of cross-cutting channels (i.e. caused by avulsion), or point bars in relation to meandering fluvial channels. The most distinctive imbricated reflections are observed at curves in the channel systems that display an asymmetric channel profile (e.g. Fig. 9b and c). Straighter segments of the channel system typically display subparallel semihorizontal facies (e.g. Fig. 11a). Based on these observations we

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

13

Fig. 11. High-resolution 2D seismic tie lines in the western part of the central channel (see Fig. 3 for position) showing a distinct bulge within the infill. The overlying succession appear undisturbed by the bulge.

find it most likely that the imbricated reflections represent point bars, deposited in relation to changing flow velocities at curves in the channel system, in turn implying a fluvial origin. The distinctive channel profile with multiple troughs observed in Fig. 9b and c as well as the general asymmetry of many of the channels suggest a certain amount of lateral migration during erosion of the channels. This is a typical characteristic feature of subaerial meandering channel systems (e.g. Lewin and Macklin, 2003). Concerning the infill lithologies, the 300 m offset between the only borehole BH550416.7 intersecting a channel and a highresolution 2D line (Fig. 9a, b and e) makes it unclear if the 16 m thick succession of horizontally laminated clay with silt streaks relates to Facies 2 or 3. Gibbard and Lewin (2002) found that interglacial river systems in lowland Britain sometimes adopted a stable fine-sediment-dominated meandering mode, and a silty clay succession could thus be of fluvial as well as of tidal-marine origin.

No final conclusions can be made based on the lithologies, but the existence of the seismic diagnostic fluvial features listed above, allows for Facies 2 to be interpreted as fluvial deposits, and the silty clay interval in BH550416.7 is attributed to this facies. 5.1.3. Seismic infill facies 3 The low-amplitude subparallel reflection pattern of Facies 3 indicates deposition in a relatively low energy environment, such as estuarine or marine environments. The change to such environments can be explained by a gradual inundation and transition from fluvial to tidal and later marine conditions. Similar inundation systems are described by Hepp et al. (2017) in the southern North Sea that contain tidal or marine silty clay deposits overlying basal peat layers. These deposits show similarities with the laminated clay in BH550416.7, but contain a basal peat layer absent in this borehole. The concave shape of the erosional surface (E1) observed

14

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

between Facies 2 and 3 (e.g. Fig. 9) indicates generally that the channels were underfilled before deposition of Facies 3, allowing submarine and/or estuarine tidal currents to flow into the channel during a period of rising sea level. 5.1.4. Elongated bulge The pronounced bulge in the western part of the central channel (Fig. 11) defines the top of Facies 2 and is co-located with the highamplitude erosional surface (E1) at the interface between Facies 2 and 3. The undisturbed overburden, proximity to the seafloor and a lack of physical evidence of gas above the channel exclude velocity effects as an explanation for the bulge. Elongated bulge-shaped depositional features include eskers (subglacial) (e.g. Fiore et al. (2011) and contourites (submarine). Contourites, normally have a res well-defined bottom reflection and internal architecture (Fauge et al., 1999) which cannot be observed on the seismic lines despite of the very good seismic resolution (<1 m) and considerable height of the bulge (4.5 m). However, since the top of the bulge is associated with the distinct erosional surface it is more likely that the bulge represents an erosional remnant similar to a very long longitudinal bar. The present day oceanographic circulation pattern in the North Sea, particularly in the Wadden Sea and the English Channel is completely dominated by tidal currents. If the conditions allowed it, such currents may similarly have created a set of tidal bars (i.e. the bulge) in the channel system with inlet and outlet on each side. Based on the above discussion we suggest an erosional origin to the bulge, with tidal channels on its sides. We interpret the low-amplitude parallel reflections of Facies 3 to correspond to a marine or tidal/estuarine environment. However, in order to constrain the exact origin of Facies 3 and Facies 2, more borehole data tied to high-resolution 2D seismic lines are needed. The variations in seismic facies in this channel system resemble the energy changes expected from the non-marine to marine transition described in the foraminiferal stratigraphy from the Skjold field by Knudsen (1985) as well as the evolution described by Hepp et al. (2017). Therefore we find that interpreting a similar evolution is reasonable. However similar seismic facies are also reported from tunnel valley infill, usually with a more pronounced chaotic seismic expression at the base of systems (O’Cofaigh, 1996; Lutz et al., 2009; Kristensen and Huuse, 2012). 5.2. Plan view morphologies The dendritic plan view morphology of the investigated channel system is a typical characteristic feature for subaerial river systems. Similar dendritic patterns have however also been described for some subglacial tunnel valley systems (Praeg, 1996) and the overall dendritic pattern is hence not a conclusive argument for either a subaerial or subglacial formation. The plan view morphologies observed for the eastern channel segments however include a number of distinct fluvial features, such as scroll bars, splays, oxbow lakes, large meander bends and abandoned channels (Fig. 7). Looking at the angle of confluence between the branching channels, these are mostly around 40 e55 , 70 or 90 , which is to be expected for a dendritic river system (De Serres and Roy, 1990; Hackney and Carling, 2011), since the angle of confluence is controlled by the shape of the basin, in areas without significant geological structures constraining the water flow (Hackney and Carling, 2011). In a subglacial setting, however, the water flow is controlled mainly by the hydraulic gradient, which in turn is determined by

the ice surface (O’Cofaigh, 1996). A dynamic ice cover would expectedly lead to variations in the hydraulic gradient with time, which in turn would force tributary confluence angles to vary according to the change in hydraulic gradient. Such a mechanism was hypothesized by Moreau et al. (2012) as an argument for multiple ice front positions related to confluence points of large tunnel valleys in the southern Dutch and British sectors of the North Sea Basin. The tributaries of the dendritic channel system in this study are however morphologically different than the tributaries described by Moreau et al. (2012), as they act as smaller tributaries feeding a larger system (i.e. the northern and southern channels feed into the central channel) and not large converging tunnel valleys. Furthermore, the observed angle of confluence for the channel system varies less than expected for subglacial drainage systems (cf. Moreau et al., 2012), favouring a subaerial origin. Determining a subglacial or subaerial origin solely on the basis of confluence angles is however speculative and needs to be combined with other evidence. 5.3. Channel geometries The width to depth ratio (9e40) of the channel system matches both a subglacial and subaerial origin according to Gibling (2006)'s “valley fills within subglacial or proglacial settings” classification. That definition is however broad and comprises valleys with widthto-depth ratios in the interval of 2.5e42. If subaerial, the dimensions of the system puts it in the category of “fixed river systems”. These overlaps between definitions suggest that W/D ratio alone cannot be used to determine whether the system was formed in a subglacial environment or not. Assuming that the system would not erode below sea level, if it was a subaerial system, a comparison between the paleo sea level curves and the erosional depths could determine whether it was possible that the system was subaerial. The maximum depth to the base of the channel is 94 ms (76 m). But, even with the modelled LGM sea level of 130 m (Lambeck et al., 2014), that measurement is not indicative of the depositional environment, as the age interval for the system is too wide compared to the sea level changes that occurred from the LGM to the early Holocene. Complete reoccupation of subglacial channels by subsequent fluvial systems are rare and have previously been questioned as a mechanism for this channel system by Müther et al. (2012). They describe the channel system as a subglacial landform, but argue that the dendritic nature of the system is unlikely for a tunnel valley catchment. However Praeg (1996) describe dendritic patterns in tunnel valleys in the British sector, and even though complete fluvial reoccupation is rare it is still possible that the system was occupied by proglacial rivers upon deglaciation, and continued to evolve as a fluvial system after the ice had retreated. 5.4. Subglacial, subaerial or both The following are arguments for either the system being either subaerial, subglacial or a combination of the two. 5.4.1. Subaerial formation In a subaerial model of formation, the paleo transport direction should be expressed, or explained unambiguously in both the plan view morphology and the thalweg dip direction. The dendritic plan-view geometries and the confluence angle of tributaries towards the central channel suggests a SW to NE directed flow within a tributary system. The regional paleogeography, with the main depocenter of the North Sea situated north of this study area (Coles, 1998; Ottesen et al., 2018) and no obvious

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

depocenter to the south also favours a tributary system with a NE transport direction. Distributary systems are usually associated with deltas or alluvial fans, neither of which there is evidence for in our data. Thus, we do not find any reason to believe that this was a distributary system. The large tributaries from the south and the north (southern and northern channels, Fig. 3) strike along the dip of the surrounding topography, also pointing towards a tributary system. Coffey and Shaw (2017) and Devauchelle et al. (2012) describe how both tributary and distributary systems dominated by groundwater flow will form angles around 72 . The angle distribution in this system suggests that it might, to some degree, be controlled by groundwater flow. The combination of angle populations around 40 e55 , 70 and 90 further indicate that this tributary system had an overall paleo flow direction from southwest to northeast, with the largest tributaries entering the system from the NNW (Fig. 3). However, an opposite dominating, but undulating, flow direction from NE to SW is indicated by the apparent thalweg dip of the central and southern channels (Fig. 8). The simplest approach explaining this apparent inverse paleoflow in a subaerial model of formation is the effect of the ice sheet forebulge and subsequent isostatic response to it upon deglaciation. Theoretical calculations of the forebulge and isostatic compensation puts the tilting line near this study area (Fjeldskaar, 1994) and the approximate position of the glacial forebulge coincides with the position of the channel system (Fig. 1). This explains the apparent inversion of flow direction by local variations in isostatic compensation. In the “southern channels” part of the system (Fig. 3), third order tributaries show an apparent flow direction towards larger channels as expected (Fig. 8). When considering the glacio-isostatic variations it would be expected that channel segments perpendicular to the tilting lines (i.e. the first order and some second order channels) are more susceptible to apparent flow inversion than channel segments parallel to the tilting lines (i.e. 3rd order channels). This is the case for the southern tributaries (Figs. 3 and 8). However, some of the second order channels have been inverted like the first order channels. All these observations suggest that it is unlikely that the system was active as a subaerial fluvial system only.

5.4.2. Subglacial formation In a subglacial model of formation, the thalweg profile is expected to undulate due to the high water pressure. This is indeed the case for the central channel system. However most of the second order tributaries show a thalweg dip direction towards the central system, along the general dip direction of the surrounding topography. Assuming that tunnel valleys drain meltwater towards the ice margin, tributaries entering a central channel from both upice and down-ice directions this close to the LGM ice sheet margin is unlikely. This is emphasized by the dendritic nature of the system, as well as expected channel convergence towards the ice margin (Moreau et al., 2012), which is not observed. This suggests several episodes of tunnel valley formation in order to form the plan view pattern seen in this system. If the channel system was initiated subglacially it must have been during the LGM, as it is the only time during the Weichselian that ice could have extended far enough into the Danish North Sea (Boulton et al., 1985; Ehlers et al., 2011; Hughes et al., 2016). From the directionality of the southern and northern channels and the position so close to the LGM ice margin, we find it unlikely that the system was active as a subglacial tunnel valley system only.

15

5.4.3. Subglacial initiation, fluvial reoccupation In a combined subaerial and subglacial model of formation, the undulating thalweg profile of the central channel would have been created by pressurized meltwater and upon deglaciation, fluvial reoccupation of the system continued its morphological evolution into a large dendritic tributary fluvial system with both a northern and a southern catchment and a depocenter to the northeast. Hughes et al. (2016) have placed their maximum confidence LGM extent line (Fig. 1) without morphological evidence in their database or reference list (DATED-1) in this study area. We therefore expect that the maximum confidence line of Hughes et al. (2016) was positioned there through interpolation between measurement points in Britain (Clark et al., 2010) and possibly Denmark, as well as morphologies in the Dogger bank area that were correlated through the area. If this system was initiated in a subglacial setting, it would also provide morphological evidence for the DATED-1 database. We find it most likely that the system was initiated as a subglacial system, which evolved into a fluvial system upon deglaciation. 6. Conclusions and implications The study has taken advantage of 3D conventional seismic data and high-resolution site surveys acquired for locating Hydrocarbon related infrastructure in order to map Weichselian-Holocene channels and their infill in the central North Sea region. The study shows that:  The channel system was initiated as a subglacial tunnel valley system, formed during the LGM (Fig. 12a). This is evident from measurements of the depth to the base of the channel showing an inverted transport direction compared to the results from the time slice analysis. The channel system was subsequently reoccupied by rivers, and evolved into a fluvial system draining both a northern and a southern catchment towards the northeast, during the late Weichselian-early Holocene (Fig. 12b). During the early Holocene, marine incursion created estuarine tidal conditions in the channels (Fig. 12c) before the system was completely drowned and buried during the late Holocene (Fig. 12d).  The seismic facies from this study confirm the evolution from non-marine to marine conditions that occurred in this area during the Holocene as presented by Knudsen (1985).  The system could have evolved entirely as a fluvial system, with post burial inversion of the apparent flow pattern, caused by glacio-isostatic effects. However, we find that a subglacial origin is more likely due to a very deep erosion and variations in the depth to the channel bottom along its course, indicating formation by pressurized water.  We can confirm that the maximum confidence line of ice sheets advance during the LGM was at least as far south as presented by Hughes et al. (2016), and this study provides morphological evidence for that model and represents the first mapping of a Weichselian tunnel valley system this far south in the Danish North Sea. With better constraints on infill lithology and age, the system could potentially provide new information about the deglaciation history in the area after the LGM. The study has also shown that the combination of 3D time slice analysis with high-resolution 2D seismic data and borehole data enables mapping of subtle shallow features such as several generations of buried valleys. Differentiating between subglacial and subaerial valleys based solely on seismic data has, proven difficult, and in order to be certain of a subglacial or subaerial origin, a high-

16

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

Fig. 12. Conceptual diagram showing the different depositional environments at different times. a) Generation of tunnel valleys during the LGM. b) Rivers from the north and the south-southwest converge and transport sediment across the marsh land towards north-northeast. Silt and peat along with fluvial sediments are deposited (Facies 1 and 2). c) The ice sheets are no longer present in the North Sea and marine incursion has initiated from the northeast. Estuaries are formed in the former rivers that are now under tidal influence and Facies 3 is deposited. d) The transgression is complete and the marine sands of Unit 1 are deposited now overlaying and burying the channel system.

resolution 3D image of the valley thalweg must be constructed along with more core data and dating of infill and substrate. The mapping has, however improved our understanding of the relative age of crosscutting channels as well as the interpretation of the channel infill and let us conclude that the system is a subaerial proglacial or fluvial system initiated in a subglacial environment. Acknowledgements Extensive thanks go to the Danish Hydrocarbon Research and
sik Prins during Technology Centre (DHRTC) who is funding Lasse Te his PhD project. DHRTC, Total and GEUS are furthermore thanked for providing seismic and well data for the study and for allowing us to publish our results. Schlumberger, Elllis and IHS Kingdom are thanked for granting academic licenses for their software to SeisLab Aarhus without which we could not have performed the study. Ole Rønø Clausen is thanked for constructive feedback on the manuscript. Appendix A. Supplementary data Supplementary data to this article can be found online at

https://doi.org/10.1016/j.quascirev.2019.105943.

References Ahmadi, Z.M., Sawyers, M., Kenyon-Roberts, S., Stanworth, C.W., Kugler, K.A., Kristensen, J., Fugelli, E.M.G., 2003. Paleocene. In: Evans, D., Graham, C., Armour, A.A., Bathurst, P. (Eds.), The Millennium Atlas: Petroleum Geology of the Central and Northern North Sea, vol. 2003. The Geological Society of London, pp. 235e259. Bendixen, C., Lamb, R.M., Huuse, M., Boldreel, L.O., Jensen, J.B., Clausen, O.R., 2018. Evidence for a grounded ice sheet in the Central North Sea during the early middle Pleistocene donian glaciation. J. Geol. Soc. 175, 291e307. https://doi.org/ 10.1144/jgs2017-073. €se, M., Lüthgens, C., Lee, J.R., Rose, J., 2012. Quaternary glaciations of northern Bo Europe. Quat. Sci. Rev. 44, 1e25. https://doi.org/10.1016/j.quascirev.2012.04.017. Boulton, G.S., Smith, G.D., Jones, A.S., Newsome, J., 1985. Glacial geology and glaciology of the last mid-latitude ice sheets. J. Geol. Soc. 142, 447e474. https:// doi.org/10.1144/gsjgs.142.3.0447. Bradley, S.L., Milne, G.A., Teferle, F.N., Bingley, R.M., Orliac, E.J., 2009. Glacial isostatic adjustment of the British Isles: new constraints from GPS measurements of crustal motion. Geophys. J. Int. 178 (1), 14e22. https://doi.org/10.1111/j.1365246X.2008.04033.x. Brown, A.R., 1999. In: Interpretation of Three-Dimensional Seismic Data, fifth ed., Ch. 1. AAPG. 13: 978-0891813521. Carr, S.J., Holmes, R., Van der Meer, J.J.M., Rose, J., 2006. The last glacial maximum in the north Sea Basin: micromorphological evidence of extensive glaciation. J. Quat. Sci. 21, 131e153. https://doi.org/10.1002/jqs.950. Cartwright, J., Huuse, M., 2005. 3D seismic technology: the geological “ Hubble ”’.

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943 Basin Res. 17, 1e20. https://doi.org/10.1111/j.1365-2117.2005.00252.x. Clark, C.D., Hughes, A.L.C., Greenwood, S.L., Jordan, C., Sejrup, H.P., 2010. Pattern and timing of retreat of the last British-Irish Ice Sheet. Quat. Sci. Rev. 44, 112e146. https://doi.org/10.1016/j.quascirev.2010.07.019. Coffey, T.S., Shaw, J.B., 2017. Congruent bifurcation angles in river delta and tributary channel networks. Geophys. Res. Lett. 44, 11427e11436. https://doi.org/10.1002/ 2017GL074873. Cohen, K.M., Gibbard, P.L., Weerts, H.J.T., 2014. North Sea paleogeographical reconstructions for the last 1 Ma. Neth. J. Geosci. 93, 7e29. https://doi.org/10. 1017/njg.2014.12. Coles, B.J., 1998. Doggerland: a speculative survey. Proc. Prehist. Soc. 64, 45e81. https://doi.org/10.1017/S0079497X00002176. Cotterill, C.J., Phillips, E., James, L., Forsberg, C.F., Tjelta, T.I., Carter, G., Dove, D., 2017. The evolution of the Dogger Bank, North Sea: a complex history of terrestrial, glacial and marine environmental change. Quat. Sci. Rev. 171, 136e153. https:// doi.org/10.1016/j.quascirev.2017.07.006. Coughlan, M., Fleischer, M., Wheeler, A.J., Hepp, D.A., Hebbeln, D., 2018. A revised stratigraphical framework for the Quaternary deposits of the German North Sea sector: a geological-geotechnical approach. Boreas 47, 80e105. https://doi.org/ 10.1111/bor.12253. De Serres, B., Roy, A.G., 1990. Flow direction and branching geometry at junctions in dendritic river networks. Prof. Geogr. 42, 194e201. Devauchelle, O., Petroff, A.P., Seybold, H.F., Rothman, D.H., 2012. Ramification of stream networks. Proc. Natl. Acad. Sci. 109 (51), 20832e20836. https://doi.org/ 10.1073/pnas.1215218109. Ehlers, J., Gibbard, P.L., Hughes, P.D., 2011. Introduction. In: Ehlers, J., Gibbard, P.L., Hughes, P.D. (Eds.), Quaternary Glaciations - Extent and Chronology, A Closer Look. Elsevier, pp. 1e14. Ekman, M., 1996. A consistent map of the postglacial uplift of Fennoscandia. Terra. Nova 8, 158e165. Emery, K.O., Aubrey, D.G., 1985. Glacial rebound and relative sea levels in Europe from tide-gauge records. Tectonophysics 120, 239e255. https://doi.org/10.1016/ 0040-1951(85)90053-8. res, J., Stow, D.A.V., Imbert, P., Viana, A., 1999. Seismic features diagnostic of Fauge contourites drifts. Mar. Geol. 162, 1e38. Fiore, J., Girardclos, S., Pugin, A., Gorin, G., Wildi, W., 2011. Würmian deglaciation of western Lake Geneva (Switzerland) based on seismic stratigraphy. Quat. Sci. Rev. 30, 377e393. https://doi.org/10.1016/j.quascirev.2010.11.018. Fitch, S., Thomson, K., Gaffney, V., 2005. Late Pleistocene and Holocene depositional systems and the palaeogeography of the Dogger bank, North Sea. Quat. Res. 64, 185e196. https://doi.org/10.1016/j.yqres.2005.03.007. Fjeldskaar, W., 1994. The amplitude and decay of the glacial forebulge in Fennoscandia. Nor. Geol. Tidsskr. 74, 2e8. Gaffney, V., Thomson, K., Fitch, S., 2007. Mapping Doggerland, the mesolithic landscapes of the southern North sea. Archaeopress 131, 1e10. GEUS e Frisbee, N.D. (https://frisbee.geus.dk/frisbee). GEUS - MARTA (N.D), (https://eng.geus.dk/products-services-facilities/data-andmaps/marine-raw-material-database-marta/). Gibbard, P.L., Lewin, J., 2002. Climate and related controls on interglacial fluvial sedimentation in lowland Britain. Sediment. Geol. 151, 187e210. Gibling, M.R., 2006. Width and thickness of fluvial channel bodies and valley fills in the geological record: a literature compilation and classification. J. Sediment. Res. 76, 731e770. https://doi.org/10.2110/jsr.2006.060. Graham, A.G.C., Lonergan, L., Stoker, M.S., 2010. Depositional environments and chronology of Late Weichselian glaciation and deglaciation in the central North Sea. Boreas 39, 471e491. https://doi.org/10.1111/j.1502-3885.2010.00144.x. Hackney, C., Carling, P., 2011. The occurrence of obtuse junction angles and changes in channel width below tributaries along the Mekong River, south - east Asia. Earth Surf. Process. Landforms 36, 1563e1576. https://doi.org/10.1002/esp.2165. €rz, T., 2017. Tributaries of the elbe palaeHepp, D.A., Warnke, U., Hebbeln, D., Mo ovalley: features of a hidden palaeolandscape in the German bight, North Sea. In: Under the Sea: Archaeology and Palaeolandscapes of the Continental Shelf, Ch. 14. https://doi.org/10.1007/978-3-319-53160-1. Houmark-Nielsen, M., 2011. Pleistocene glaciations in Denmark: a closer look at chronology, ice dynamics and landforms. In: Ehlers, J., Gibbard, P.L., Hughes, P.D. (Eds.), Quaternary Glaciations - Extent and Chronology, A Closer Look. Elsevier, pp. 47e58. https://doi.org/10.1016/B978-0-444-53447.00005-2. Hughes, A.L.C., Gyllencreutz, R., Lohne, Ø.Y.S., Mangerud, J.A.N., Inge, J., 2016. The last Eurasian ice sheets e a chronological database and time-slice reconstruction, DATED-1. Boreas 45, 1e45. https://doi.org/10.1111/bor.12142. Huuse, M., Lykke-Andersen, H., 2000. Overdeepened quaternary valleys in the eastern Danish North Sea: morphology and origin. Quat. Sci. Rev. 19, 1233e1253. https://doi.org/10.1016/S0277-3791(99)00103-1. Jensen, J.B., Bennike, O., Lemke, W., Kuijpers, A., 2005. The storebælt gateway to the baltic the storebælt gateway to the baltic. Geol. Soc. Den. Greenl. Bul. 7, 45e48. Johnson, H., Richards, P.C., Long, D., Graham, C.C., 1993. United Kingdom Offshore Regional Report: the Geology of the Northern North Sea. HMSO, London, 111pp. Kiden, P., Denis, P., Johnston, P., 2002. Late Quaternary sea-level change and isostatic and tectonic land movements along the Belgian-Dutch North Sea coast: geological data and model results. J. Quat. Sci. 17, 535e546. https://doi.org/10. 1002/jqs.709. Knudsen, K.L., 1985. Foraminiferal stratigraphy of quaternary deposits in the roar, Skjold and dan fields, Central North sea. Boreas 14, 311e324. ISSN 0300e9483. Knudsen, K.L., Sejrup, H.P., 1993. Pleistocene stratigraphy in the Devils Hole area, central North Sea: foraminiferal and amino-acid evidence. J. Quat. Sci. 8, 1e14.

17

Kristensen, T.B., Huuse, M., 2012. Multistage erosion and infill of buried Pleistocene tunnel valleys and associated seismic velocity effects. In: Huuse, M., Redfern, J., LeHeron, D.P., Dixon, R.J., Moscariello, A., Craig, J. (Eds.), 2012. Glaciogenic Reservoirs and Hydrocarbon Systems, vol. 368. Geological Society, London, Special Publications, pp. 159e172. https://doi.org/10.1144/SP368.15. Kristensen, T.B., Huuse, M., Piotrowski, J., Clausen, O.R., 2007. A morphometric analysis of tunnel valleys in the eastern North Sea based on 3D seismic data. J. Quat. Sci. 22, 801e815. https://doi.org/10.1002/jqs. Kristensen, T.B., Piotrowski, J.A., Huuse, M., Clausen, O.R., Hamberg, L., 2008. Timetransgressive tunnel valley formation indicated by infill sediment structure , North Sea e the role of glaciohydraulic supercooling. Earth Surf. Process. Landforms 33, 546e559. https://doi.org/10.1002/esp. Lambeck, K., Rouby, H., Purcell, A., Sun, Y., Sambridge, M., 2014. sea level and global ice volumes from the last glacial maximum to the Holocene. In: Proceedings of the National Academy of Sciences Oct 2014, vol. 111, pp. 15296e15303. https:// doi.org/10.1073/pnas.1411762111. Lewin, J., Macklin, M.G., 2003. Preservation potential for late quaternary river alluvium. J. Quat. Sci. 18, 107e120. https://doi.org/10.1002/jqs.738. Lonergan, L., Maidment, C.R., Collier, J.S., 2006. Pleistocene subglacial tunnel valleys in the central North Sea basin: 3-D morphology and evolution. J. Quat. Sci. 21, 891e903. ISSN 0267e8179. Lutz, R., Kalka, S., Gaedicke, C., Lutz, R., Winsemann, J., 2009. Pleistocene tunnel valleys in the German North Sea: spatial distribution and morphology. Z. Dtsch. Gol. Ges. 160, 225e235. https://doi.org/10.1127/1860-1804/2009/0160-0225. Moreau, j., Huuse, M., Janszen, A., Van Der Vegt, P., 2012. The glaciogenic unconformity of the southern North Sea. In: Huuse, M., Redfern, J., LeHeron, D.P., Dixon, R.J., Moscariello, A., Craig, J. (Eds.), 2012. Glaciogenic Reservoirs and Hydrocarbon Systems, vol. 368. Geological Society, London, Special Publications, pp. 99e110. https://doi.org/10.1144/SP368.5. Müther, D., Back, S., Reuning, L., Kukla, P., Lehmkuhl, F., 2012. Middle Pleistocene landforms in the Danish Sector of the southern North Sea imaged on 3D seismic data. In: Huuse, M., Redfern, J., LeHeron, D.P., Dixon, R.J., Moscariello, A., Craig, J. (Eds.), 2012. Glaciogenic Reservoirs and Hydrocarbon Systems, vol. 368. Geological Society, London, Special Publications, pp. 111e127. https://doi.org/10. 1144/SP368.7.
Cofaigh, C., 1996. Tunnel valley genesis. Prog. Phys. Geogr. 20, 1e19. O Ottesen, D., Dowdeswell, J.A., Bugge, T., 2014. Morphology, sedimentary infill and depositional environments of the early quaternary North Sea basin (56-62N). Mar. Pet. Geol. 56, 123e146. https://doi.org/10.1016/j.marpetgeo.2014.04.007. Ottesen, D., Batchelor, C.L., Dowdeswell, J.A., Løseth, H., 2018. Morphology and pattern of quaternary sedimentation in the north Sea Basin. Mar. Pet. Geol. 98, 836e859. https://doi.org/10.1016/j.marpetgeo.2018.08.022. Overeem, I., Weltje, G.,J., Bishop-Kay, C., Kroonenberg, S.,B., 2001. The late cenozoic eridanos delta system in the southern North Sea Basin: a climate signal in sediment supply? Basin Res. 13, 293e312. https://doi.org/10.1046/j.1365-2117. 2001.00151.x. Penkman, K., 2010. Amino acid geochronology: its impact on our understanding of the Quaternary stratigraphy of the British Isles. J. Quat. Sci. 25, 501e514. ISSN 0267e8179. Phillips, E., Hodgson, D.M., Emery, A.R., 2017. The quaternary geology of the north Sea Basin. J. Quat. Sci. 32, 117e126. https://doi.org/10.1002/jqs.2932. Posamentier, H.W., Davies, R.J., Cartwright, J.A., Wood, L., 2007. Seismic geomorphology e an overview. In: Posamentier, H.W., Wood, L.J., Cartwright, J.A. (Eds.), Seismic Geomorphology: Applications to Hydrocarbon Exploration and Production, vol. 277. Geol. Soc. London, Special publications, pp. 1e14. Praeg, D., 1996. Morphology, Stratigraphy and Genesis of Buried Mid-pleistocene Tunnel-Valleys in the Southern North Sea Basin. PhD Thesis. University of Edinburgh. Praeg, D., 2003. Seismic imaging of mid-Pleistocene tunnel-valleys in the North Sea Basin-high resolution from low frequencies. J. Appl. Geophys. 53, 273e298. https://doi.org/10.1016/j.jappgeo.2003.08.001. Reinardy, B.T.I., Petter, H., Hjelstuen, B.O., King, E., Augedal, H., 2018. A Quaternary aminostratigraphy constraining chronology of depositional environments in the North Sea Basin. Mar. Geol. 402, 139e152. https://doi.org/10.1016/j.margeo. 2017.11.004. Steffen, H., Wu, P., 2011. Glacial isostatic adjustment in Fennoscandia e a review of data and modelling. J. Geodyn. 52, 169e204. https://doi.org/10.1016/j.jog.2011. 03.002. Stewart, M.A., 2009. 3D Seismic Analysis of Pleistocene Tunnel Valleys in the Central North Sea. Ph. D thesis. University of London, p. 319. Stewart, M.A., 2016. Assemblage of buried and seabed tunnel valleys in the central North Sea: from morphology to ice-sheet dynamics. Geol. Soc. London, Memoirs 46, 317e320. https://doi.org/10.1144/M46.140. Stewart, M.A., Lonergan, L., 2011. Seven glacial cycles in the middle-late Pleistocene of northwest Europe: geomorphic evidence from buried tunnel valleys. Geology 39, 283e286. https://doi.org/10.1130/G31631.1. Stewart, M.A., Lonergan, L., Hampson, G., 2013. 3D seismic analysis of buried tunnel valleys in the central North Sea: morphology, cross-cutting generations and glacial history. Quat. Sci. Rev. 72, 1e17. https://doi.org/10.1016/j.quascirev.2013. 03.016. Van der Vegt, P., Janszen, A., Moscariello, A., 2012. Tunnel valleys: current knowledge and future perspectives. In: Huuse, M., Redfern, J., LeHeron, D.P., Dixon, R.J., Moscariello, A., Craig, J. (Eds.), 2012. Glaciogenic Reservoirs and Hydrocarbon Systems, vol. 368. Geological Society, London, Special Publications, pp. 111e127 (1). https://doi.org/10.1144/SP368.13.

18

L.T. Prins, K.J. Andresen / Quaternary Science Reviews 223 (2019) 105943

Van Heteren, S., Meekes, J.A.C., Bakker, M.A.J., Gaffney, V., Fitch, S., Gearey, B.R., Paap, B.F., 2014. Reconstructing North Sea paleolandscapes from 3D and highdensity 2D seismic data: an overview. Neth. J. Geosci. 93, 31e42. https://doi. org/10.1017/njg.2014.4. Westaway, R., 2017. Isostatic compensation of Quaternary vertical crustal motions:

coupling between uplift of Britain and subsidence beneath the North Sea. J. Quat. Sci. 32, 169e182. https://doi.org/10.1002/jqs.2832. Wingfield, R.T.R., 1989. Glacial incisions indicating middle and upper Pleistocene ice limits off Britain. Terra. Nova 1, 538e548.