Evidence for a grounded ice sheet in the central North Sea during the early Middle Pleistocene Donian Glaciation

The central North Sea has been covered by extensive ice sheets, sourced in the surrounding landmasses, several times during the Pleistocene, leaving behind a wide array of preserved landforms including glaciotectonic thrust structures (Ehlers et al. 1984; Cameron et al. 1987; Huuse & Lykke-Andersen 2000a; Andersen 2004; Carr 2004; Buckley 2012) and tunnel valleys (Cameron et al. 1987; Wingfield 1989, 1990; Praeg 2003; Kristensen et al. 2007; Stewart & Lonergan 2011; Buckley 2012; van der Vegt et al. 2012; Stewart et al. 2013). The rise in quality and coverage of offshore seismic reflection data, particularly 3D seismic data, has increased the number of identified glacial features as well as improved our understanding of the full 3D geometry of these preserved landforms (e.g. Huuse & Lykke-Andersen 2000a; Praeg 2003; Andersen et al. 2005; Buckley 2012). It is therefore timely to begin to integrate these interpretations with the more traditional investigation of glacial features, particularly glaciotectonics, from onshore outcrop studies (e.g. Aber & Lundqvist 1988; Hart & Boulton 1991; Harris et al. 1997; Bakker & van der Meer 2003; Astakhov 2004; Pedersen 2005; Lee 2009). In this study, interpretation and analysis of 2D and 3D seismic data from the central North Sea reveal thin-skinned subsurface structural features linked with glaciation (folds, faults and glaciotectonic fabrics) of crucial importance in reconstructing palaeo-environments and understanding ice sheet advance and retreat patterns during the Pleistocene. In addition, the structures provide important information on glacial deformation events (e.g. bed deformation, meltwater channels and depositional landforms) and thus the origin, nature and evolution of lowland glaciations in NW Europe.

During the last few decades, much of the research focused on the NW European region by Quaternary scientists has been on the major glacial periods of the Late Pleistocene; that is, the Elsterian (Marine Isotope Stage (MIS) 12), Saalian (MIS 10 – 6) and particularly the Weichselian (MIS 4 – 2). This is the case in both onshore and offshore regions in relation to the extent of ice sheets, morphology, glaciodynamics and glaciotectonics (Huuse & Lykke-Andersen 2000a, b; Bennett 2001; Sejrup et al. 2003; Andersen 2004; Svendsen et al. 2004; Pedersen 2006; Phillips et al. 2008; Lee 2009; Lee et al. 2012). Researchers use the Marine Isotope Stages, deduced from δ¹⁸O analysis of benthic foraminifera from globally distributed deep-sea core samples, as a proxy for the global ice volume, as a function of global sea-level change, and temperature (Raymo & Ruddiman 1992; Lear et al. 2000; Zachos et al. 2001; Miller et al. 2005, 2011; Kleman et al. 2008). The Pleistocene is characterized in the oxygen isotope curve by cycles between high δ¹⁸O values (4 – 5‰; Lisiecki & Raymo 2005) corresponding to glacial periods and low δ¹⁸O values (<3‰; Lisiecki & Raymo 2005) corresponding to interglacial periods.

These glacial–interglacial cycles are seen from the onset of the Quaternary at 2.58 Ma (MIS 103); however, in NW Europe glaciogenic records from prior to the Elsterian glaciation are very scarce. This is in part due to the lack of evidence preserved onshore, with much controversy over any potential pre-Elsterian deposits, owing to the overall erosional nature of the ice sheets in upland areas such that younger ice sheets have removed much of the evidence of any older ones that may have existed. The central North Sea, on the other hand, acted as a sediment trap for much of the Quaternary, accumulating 1.1 km of Quaternary deposits in the deepest parts of the basin, creating an extensive sedimentary archive sufficient to preserve glaciogenic records against erosion by later ice sheets (Caston 1977; Holmes 1977; Cameron et al. 1987; Fyfe et al. 2003; Rasmussen et al. 2005; Nielsen et al. 2008; Knox et al. 2010; Anell et al. 2012; Goledowski et al. 2012; Ottesen et al. 2014; Lamb et al. 2017b). However, the comparative lack of research into the Middle Pleistocene history of the central North Sea basin since the initial investigations during the 1980s and 1990s, despite the rapid increase in quality and quantity of the seismic data, has also limited our understanding of any pre-Elsterian glaciation.

Many of the seismic data available for study of the North Sea Pleistocene were designed and acquired for use in the hydrocarbon industry, but it is possible to extract very detailed near-surface information from the data to study the Quaternary deposits in unprecedented spatial and temporal detail, even compared with dedicated high-resolution seismic surveys acquired for shallow hazard site surveys (see Praeg 2003; Bulat 2005).

This study takes advantage of the excellent coverage of the central North Sea by continuous 3D seismic data to fully investigate a glaciotectonic complex identified on high-resolution seismic data (Fig. 1) as probable pre-Elsterian in age. The paper discusses in detail the structure and morphologies of the overall complex, identifying a series of smaller complexes and glaciotectonic units within the larger feature. The complex is interpreted with respect to a regional chronostratigraphic framework, based on palaeomagnetic correlation, and discussion of the nature of the ice front and the ice flow direction is used to reconstruct a probable series of events during the formation of the glaciotectonic complex.

The North Sea is an epeiric sea bordered by Norway to the NE, Denmark to the east, Germany and the Netherlands to the south and the UK to the west (Fig. 1). It is generally characterized by water depths less than 100 m, with limited accommodation space left to fill (e.g. Huuse & Lykke-Andersen 2000a). The structural configuration of the North Sea is largely a result of the Late Jurassic–Early Cretaceous rifting followed by subsequent thermal cooling and subsidence combined with basin infill (Ziegler 1990; Nielsen et al. 2009). Subsidence continued during the Pliocene and Quaternary (Caston 1977; Holmes 1977; Cameron et al. 1987; Rasmussen et al. 2005; Ottesen et al. 2014).

In the Northern North Sea, the Pleistocene succession consists of shelf-progradational clinoforms, which mark an increase in sediment supply from southern Norway concurrent with the first major Northern Hemisphere glaciations (Gregersen et al. 1997; Huuse 2002; Ottesen et al. 2012). In the central North Sea, Lower and early Middle Pleistocene deposits consist of a single seismo-stratigraphic unit, the Aberdeen Ground Formation. On seismic sections, the unit is characterized by laterally continuous, high-amplitude reflections (Cameron et al. 1987; Buckley 2012) and shallow cores suggest that the sediment consists of predominantly delta-front or open shallow marine sediments evolving to glaciomarine towards the top of the Formation, as evidenced in both sediments and microfossils (Cameron et al. 1987; Sejrup et al. 1987; Buckley 2012). Age-control on the Aberdeen Ground Formation is limited, with most of the available dating tied to changes in the Earth's magnetic field (Stoker et al. 1983). The Jaramillo Palaeo-magnetic Event (1.07 Ma; Fig. 2; MIS 31; Channell & Kleiven 2000; Lisiecki & Raymo 2005; Channell et al. 2009) and Brunhes–Matuyama Palaeo-magnetic Reversal, used as the marker of transition from the Lower to Middle Pleistocene (780 ka ± 5 ka BP; Fig. 2; MIS 19; Stoker et al. 1983; Lisiecki & Raymo 2005; Channell et al. 2009), are identified and dated in borehole 77/02 (cored by the British Geological Survey (BGS); Stoker et al. 1983; Sejrup et al. 1991; Gatliff et al. 1994). The top of the Aberdeen Ground Formation is defined by a regional glacial unconformity and dissection by a wide network of subglacial tunnel valleys, formed by the action of sub-glacial drainage, known as the Ling Bank or Swarte Bank unconformity (Cameron et al. 1987; Wingfield 1989; Sejrup et al. 1991; Gatliff et al. 1994; Huuse & Lykke-Andersen 2000b; Kluiving et al. 2003; Praeg 2003; Kristensen et al. 2008; Stewart & Lonergan 2011; Stewart et al. 2012, 2013; van der Vegt et al. 2012). The age of this unconformity and the tunnel valleys that form is poorly constrained owing to the erosional nature of the unconformity. However, it has been approximated from comparison with tunnel valleys of onshore mainland Europe where the valley infill and associated facies have been dated to the Holsteinian interglacial (MIS 11; Ehlers et al. 1984; Scourse et al. 1998; Huuse & Lykke-Andersen 2000b; Kluiving et al. 2003; Stewart & Lonergan 2011) giving a minimum age for formation of during the Elsterian glaciation.

Our degree of understanding of the North Sea glaciations is directly correlated to their age, with the Weichselian and Late Saalian complex geological history now well understood (MIS 6–MIS 2; e.g. Graham et al. 2011; Lee et al. 2012; Stewart et al. 2012). The early Saalian and Elsterian (MIS 12–MIS 10) are less well constrained chronologically; however, it is widely accepted that the North Sea was repeatedly covered by large ice sheets from Fennoscandia and Britain during the Middle and Late Pleistocene (Graham et al. 2011; Lee et al. 2012; Moreau et al. 2012; Stewart et al. 2012). The sediments preserved in the North Sea for this period, generally between 100 and 300 m thick depending on tunnel valley fill, are dominated by glacial erosion and depositional processes building an archive of glacial activity (Graham et al. 2011). Although the three-glaciation model for the period between MIS 12 and MIS 2 remains in use there is some debate on how this model fits in with recent studies, which have identified a minimum of seven generations of ice advances between MIS 12 and MIS 2, based on crosscutting relationships between sub-glacial tunnel valleys, and whether these represent individual glaciations or multiple ice advances within a single glaciation (Kristensen et al. 2007; Stewart & Lonergan 2011; Stewart et al. 2013).

On- and offshore evidence of pre-Elsterian glaciations throughout NW Europe is limited (Graham et al. 2011; Lee et al. 2012), with the onshore regions being particularly susceptible to glacial erosion by younger ice sheets whereas offshore regions remain under-investigated. Part of the cause of this situation has been a largely unspoken bias towards the period following the ‘Mid-Pleistocene Transition’ (1.2 – 0.8 Ma) during which time the frequency and magnitude of the glacial–interglacial cycles switched from a 41 kyr cycle to a longer and higher amplitude 100 kyr cycle, presumed to be caused by Milankovitch forcing (Raymo 1994; Raymo & Nisancioglu 2003; Lisiecki & Raymo 2005, 2007; Pedersen 2012). Shorter and weaker glacial cycles are considered prohibitive to the formation and preservation of glacial features. The earliest evidence of North Sea glaciation is Menapian (c. 1.2 Ma) according to Gibbard (1988) and Carr (2004), due to the presence of subglacial diamicton in a borehole located in the British sector of the northern North Sea (Sejrup et al. 1987) and the Troll core in the Norwegian Channel (Sejrup et al. 1995, 2003); however, this evidence has rarely been taken to imply large-scale ice sheet formation over the North Sea basin itself. In comparison studies of the Middle Pleistocene, glaciations of Eastern Europe and Northern Russia yielded evidence in the form of subglacial diamicton termed the Donian till (Astakhov 2004; Velichko et al. 2004), which has subsequently been dated and used to propose an extensive Cromerian/Donian glaciation, which has been widely correlated to MIS 16. Bijlsma (1981) found evidence of Fennoscandian cobbles and boulders in Menapian deposits of the Baltic River system, and inferred significant inland ice over Fennoscandia during this period, but made no suggestion of extension of ice beyond the highland regions, in agreement with work done by Zagwijn (1986). In the UK onshore discussion over deformed glacial till deposits in Norfolk first suggested an MIS 16 glaciation (Lee et al. 2002, 2004); however, recent studies concluded that the Happisburgh Till must relate to a younger glaciation as the underlying interglacial deposits correlate to substages within the Cromerian (Lee 2009; Preece et al. 2009). The only evidence published thus far for extensive grounded pre-Elsterian ice in the central North Sea has been the identification of glaciotectonic thrust features in the upper Aberdeen Ground Formation in the central North Sea (Buckley 2012, 2017). The thrusts are observed below the Ling Bank Unconformity, thus implying the existence of a pre-Elsterian grounded ice sheet in the central North Sea (Buckley 2012).

A variety of morphological features associated with grounded glaciations have been observed for the Late Pleistocene throughout NW Europe and the North Sea in outcrop, bathymetry data and seismic data. Common features include mega-scale glacial lineations, tunnel valleys, glaciotectonic fabrics and iceberg ploughmarks (e.g. Sættem 1990; Huuse & Lykke-Andersen 2000a; Sejrup et al. 2003; Aber & Ber 2007; Graham et al. 2007; Ottesen et al. 2008; Stewart et al. 2012; Høyer et al. 2013; Hughes et al. 2014). Glaciotectonism arises when an advancing ice sheet subjects the substrate to different stresses horizontally and vertically at the ice front, primarily caused by the weight of the ice as well as the forward driving motion. Primary mechanisms include, but are not limited to, bulldozing, gravity spreading and subglacial deformation, which act in response to specific glacial and subglacial conditions to form thrust sheets (Aber & Lundqvist 1988; Huuse & Lykke-Andersen 2000b; Pedersen 2000, 2014; Andersen et al. 2005; Buckley 2012).

Permafrost is often associated as a significant mechanism in the form of glaciotectonic complexes, as permafrost tends to stiffen the substrate such that stresses applied by the ice sheet are transmitted over a larger area (Etzelmüller et al. 1996; Boulton et al. 1999; Bennett 2001; Madsen & Piotrowski 2012). The relationship between permafrost and deformation is complex, with evidence that it is favourable to proglacial deformation as it forms an impermeable cap to allow the build-up of groundwater pressure (e.g. Etzelmüller et al. 1996; Bennett 2001; Madsen & Piotrowski 2012) whereas other evidence suggests that push moraine formation is precluded where permafrost is too strong (Etzelmüller et al. 1996; Bennett 2001). However, it has been emphasized by various researchers (e.g. Croot 1987; Van der Wateren 1994) that permafrost is not the crucial factor in deformation of large thrust faults, and large glaciotectonic complexes may form without the presence of permafrost.

Glaciotectonic fabrics include a variety of thrust complexes as well as specific features such as hill–hole pairs. The term ‘push moraine’ is often used to describe all forms of glaciotectonic complexes, ice-marginal or sub-marginal moraines (e.g. Bennett 2001; Pedersen 2005) and is used in this paper to encompass glaciotectonic structures and landforms (e.g. hills, ridges and plains). Throughout the last two decades, glaciotectonic complexes from the Late Pleistocene have been identified and investigated in both on- and offshore regions of northern Europe (Vangkilde-Pedersen et al. 1993; Van der Wateren 1994; Huuse & Lykke-Andersen 2000a; Pedersen 2000, 2005, 2014; Andersen 2004; Lee 2009; Buckley 2012).

Hill–hole pairs are distinctive glaciotectonic landforms that consist of a discrete hill of ice-shoved material located immediately down-glacier from a source depression. Typically they are found to be between 30 and 200 m in height and cover an area between 1 and 100 km² and form at or close to the ice margin (Aber et al. 1989; Sættem 1990; Bennett 2001). Push moraines are characterized by the presence of thrust faults and imbricate stacks (Pedersen 2005) and are identified as thin-skinned deformation on seismic sections (Vangkilde-Pedersen et al. 1993; Harris et al. 1997; Huuse & Lykke-Andersen 2000a; Pedersen & Boldreel 2015). Push moraines do not occur at every ice margin, but reflect the interplay of specific glacial conditions, or conditions favourable for deformation in the glacial foreland, or a combination of both (Bennett 2001).

For example, gravity spreading is the result of the weight of the ice causing a change in the gravitational forces that act on a once stable sediment package (Aber & Ber 2007), and the thrusts of push moraines formed by gravity spreading often die out with increasing distance from the load. Push from the rear, on the other hand, causes glaciotectonic deformation owing to the forward motion of a glacier as it advances and piles up material at the ice margin during an advance, forming a push moraine (Van der Wateren 1994; Huuse & Lykke-Andersen 2000a; Bennett 2001; Andersen et al. 2005; Aber & Ber 2007). The formation of push moraines may be conceptually related to factors such as friction on the décollement surface, foreland geology and rheology, and glacier stress (Bennett 2001). These considerations are used to estimate the mechanisms responsible for the thrust complexes reported herein. In doing so, the different glacial morphologies contribute clues to the reconstruction of glacial dynamics, maximum extent and flow patterns of the ice sheet.

The study uses two 3D seismic surveys shot for the hydrocarbon industry: the high-resolution Broadseis^TM and the CNS MegaSurvey datasets. The Broadseis^TM dataset is a high-resolution 3D seismic reflection dataset with a vertical resolution of 7 – 8 m and a bin spacing of 12.5 m (Fig. 1). The CNS MegaSurvey has a vertical resolution of 10 – 15 m with a bin spacing of 50 m (Fig. 1). Three-dimensional seismic reflection volumes allow mapping with dense spatial coverage across large areas using semi-automated interpretation techniques, the use of arbitrary lines and time-slices that allow horizontal slicing of the stratigraphy to fully visualize the planform expression of stratigraphic and structural features (Brown 1999; Praeg 2003; Cartwright & Huuse 2005).

In this study the glaciotectonic complex investigated was identified first on a number of time slices, to identify the areal extent and narrow down the study area, before it was examined using a series of arbitrary seismic lines oriented perpendicular to the thrust lineaments. This allowed, in the first instance, the overall character of the complex to be identified and compared with analogues to confirm the interpretation as well as identify any potential errors owing to seismic processing. A north–south-oriented seismic acquisition footprint is present in the upper parts of the high-resolution dataset, but the noise is systematic and mainly seen on crosslines and time slices. However, careful analysis allowed detailed near-surface mapping to be carried out (Figs 3a, b and 4; Bulat 2005).

Mapping of the complex was then undertaken primarily using the in-line seismic sections, as these were least affected by the acquisition footprint, with both x-line and arbitrary lines used to ensure consistency of mapping. Three pronounced reflections representing the bounding of the complex and a regionally significant and continuous horizon were identified and mapped to form the basis of the stratigraphical framework for the complex. The strong and continuous reflection characteristics of these key surfaces allowed for consistent auto-tracking across the study area. A final surface, identifying the top of the area of stratigraphic interest, was mapped manually using a grid pattern at the base of the regional tunnel valley unconformity.

The thrust units were mapped using a fault mapping tool to connect the thrust lineaments along strike more easily, whereas truncations and unconformities related to the thrust were mapped as individual horizons. The mapping was broadly carried out in a west-to-east direction at regular intervals of 50 m, with smaller intervals used where the structure was more complex, and all glaciotectonic structures were mapped before the complex was placed into relative chronological order. Following mapping the complex was divided into a number of units, subareas (referred to as complexes) and fault types. Units represent significant stratigraphical packages bounded by the four regional surfaces and correlated to the wider regional framework and chronology. The subareas are smaller complexes that subdivide the full glaciotectonic complex by ice flow direction and relative position within the stratigraphy, representing significant dynamic changes along the ice front. Finally the fault types are principally divided by either localized truncation surfaces or a change in the nature of the deformation and are primarily for descriptive purposes, but also allow for the order of events to be more clearly discussed.

Correlation to the regional framework was undertaken by extending the mapping of the two regionally significant reflections into the CNS MegaSurvey, primarily northwards to the location of BGS borehole 77/02 (58°29.50'N, 00°30.30'E). Core 77/02 was drilled into the Quaternary sediments, penetrating to a depth of 217 m below seafloor at a water depth of 147 m, and palaeomagnetic dating was conducted on samples (Stoker et al. 1983; Sejrup et al. 1991). The magneto-stratigraphic dates were converted from metres to two-way travel time (TWT) assuming a standard seismic velocity of 1800 – 2000 m s⁻¹ as a representative value for the relatively shallow subsurface (Praeg 1996, 2003; Jørgensen et al. 2003). Two palaeomagnetic dates were found to strongly correlate to reflections mapped in the study area, the base of the Jaramillo Palaeomagnetic Event (c. 1.07 Ma; Fig. 2; Stoker et al. 1983; Sejrup et al. 1991; Lisiecki & Raymo 2005) and the Brunhes–Matuyama Reversal (c. 0.78 Ma; Fig. 2; Stoker et al. 1983; Sejrup et al. 1991; Lisiecki & Raymo 2005).

A previously unknown thrust complex has been observed on the high-resolution BroadSeis^TM dataset covering an area of 700 km² in the central North Sea, 57°16'N to 57°27'N and 0°55'E to 1°18'E (Figs 3 and 4). The thrust complex lies at a depth of c. 240 – 320 ms TWT (c. 216 – 288 m; Figs 5 – 9) between the Brunhes–Matuyama Reversal (at >320 ms TWT, 288 m depth) and the Ling Bank Unconformity (between 180 and 350 ms, c. 162 – 315 m) (Fig. 5). The stratigraphic package is split into three units of pre-thrusting, thrust complex and post-thrusting to fully describe the setting and evolution.

The stratigraphy below the thrust complex is divided by three pronounced reflections r₁–r₃, where r₁ is the oldest and r₃ is the youngest. These deposits consist of substrate consistent with the regional stratigraphy and are considered here for correlation to the regional chronostratigraphy. The two lowermost reflections r₁ and r₂ have been correlated to core 77/02 north of the study area. The reflections define two units: unit 1 bounded by r₁ and r₂ and unit 2 bounded by r₂ and r₃.

The lower boundary of unit 1 is located at a depth of c. 475 – 485 ms TWT (427 – 437 m; Figs 5 and 6) and is defined by reflection r₁. r₁ is a smooth, continuous high-amplitude and nearly horizontal reflection. The upper boundary of unit 1 is defined by reflection r₂, a nearly horizontal high-amplitude reflection at a depth of 395 – 410 ms TWT (369 – 356 m; Figs 5, 6 and 9). The internal seismic reflection pattern of unit 1 consists of horizontal reflections of medium amplitude and is interpreted to have a thickness of 75 – 107 ms (68 – 96 m) thinning towards the north and NE and downlapping onto r₁ in the area covered by the CNS MegaSurvey (Fig. 5).

Unit 2 is bounded by r₂ at its base and by r₃ at the top. r₃ is a high-amplitude reflection and is interpreted as the uppermost undeformed layer just below the thrust complexes. It is located at a depth of c. 310 – 320 ms TWT (279 – 288 m; Figs 5 – 9). The internal reflection pattern shows parallel reflections of medium to high amplitude. The thickness of the unit is 65 – 100 ms TWT (59 – 90 m) in the area that is covered by the BroadSeis^TM dataset. The unit thins towards the north and the internal reflections downlap onto r₂ towards the north in the area where the CNS MegaSurvey is located (Fig. 5).

The thrust complex has been divided into three subareas, identified as thrust complexes 1–3, based on changes in the orientation of the thrust blocks observed on the 294 ms TWT time slice (Fig. 3) and changes in dip direction seen in vertical sections (Figs 6 – 9). The entire thrust complex is visible on time slices between 315 and 265 ms TWT and the interpretation of the 294 ms time slice is illustrated in Figure 3b and c. Figure 4 shows a closer view of the planform texture of each of the thrust complexes. The thrust complexes are characterized by a very different reflection pattern compared with the near-horizontal reflections of the background stratigraphy (Figs 5 – 9).

Thrust complex 1 is located in the eastern part of the dataset (Figs 3 and 6) and covers an area of 30 km² with a maximum axial length of 9 km (Fig. 3). Its planform is that of an elongated ellipsoid (Figs 3 and 4b) with thrust slices aligned along a NW–SE orientation (Figs 3 and 4b). A SW–NE seismic cross-section through the complex shows the décollement surface (r₃) with overlying reflections forming anticlinal and thrust folding on top of each other stretching over a distance of 3.6 km to the SW of the profile (Fig. 6). To the NE the section reveals a depression of <30 ms TWT depth (27 m) infilled by onlapping sediments (Fig. 6), which is observable even without vertical exaggeration on the cross-section (Fig. 6c). The combination of SW-moving thrust-cored anticlines and a depression to the NE is similar to a hill–hole pair (Fig. 6). The average thickness of the anticline is c. 70 – 100 m. The hill–hole pair is oriented in a SW–NE direction and suggests that the ice-push originated from the NE.

Thrust complex 2 is split into a northern and a southern section (Figs 3, 4a, 7 and 8). In time slice the thrust complex appears as two bands of lineations, representing the crests of the folded beds, bending in an arc towards north (Fig. 3). The southern complex stretches over an area of 100 km² with an axial width of 25 km from east to west, and the northern complex stretches over an area of 27 km² with an axial width of 13.5 km east to west. They are separated by a basin-like structure (Figs 3 and 8).

The southern complex is thought to consist of two types of thrusts, A and B, because of observed changes in seismic reflection pattern (Fig. 7), and the distribution of these areas is illustrated in Figure 3c. The thrusts terminate downwards at the décollement surface (r₃), which dips slightly towards the south, at a depth of c. 310 ms TWT (280 m) in the proximal area and at a depth of c. 322 ms TWT (290 m) in the distal part of the complex (Fig. 7). Figure 7a shows a cross-section, with a vertical exaggeration of five, through the complex where the southernmost reflections are characterized by medium amplitude, located above the décollement surface and are interpreted as symmetrical folds (anticlines; A; Fig. 7b). In total, seven folded anticlines are interpreted with a truncated upper boundary (Fig. 7). North of these anticlines, a series of steeply SW-dipping reflections are stacked obliquely in an imbricated manner as a series of monoclines (Fig. 7). In total, 18 imbricates are identified in area B in the western part of the southern complex, and c. 55 imbricated reflections are identified to the east. The upper boundary of the complex shows truncation of the imbricates (Fig. 7). The thickness of each thrust and fold is c. 40 – 50 m. Along the eastern margin of the complex, reflections are influenced by a separate thrust complex (thrust complex 3; below) and thrusts are seen to dip in a more southerly direction (Fig. 3b and c; annotated by overlapping blue and green colours). The imbricated reflections die out along the décollement surface towards the south.

The northern complex is also divided into two distinct areas of thrusting; C and D (Figs 3c and 8). The décollement surface (r₃) is interpreted at a depth of c. 310 ms (280 m; Fig. 8). Area C consists of an anticline followed by one dipping high-amplitude reflection on the proximal side ending with slightly folded reflections. Area D consists of an anticline with five dipping reflections stacked obliquely in an imbricated manner on the proximal side of the anticline. A syncline is located at the southern termination of area D, draping onto the northern flank of the anticline of area C. The extent of area D is limited by the anticline. The thicknesses of the imbricated thrusts vary within the two areas. In area C, the thickness of each structure is c. 100 – 110 m and no truncation of reflections is observed. In area D, the thickness of each structure is c. 70 – 80 m and reflections are truncated in the northernmost section. Around 380 – 450 ms TWT, seismic multiples of the thrust complex can be observed cross-cutting the pre-thrusting stratigraphy (Fig. 8a).

Separating the northern and southern parts of thrust complex 2 is a basin structure (area E; Fig. 8). The lower boundary of the basin correlates to the southern extent of the northern thrust complex and the northern extent of the southern thrust complex. The infill consists of parallel reflectors onlapping the sides of the basin (Fig. 8a) and small thrusts are observed along the northern limit of the basin, deforming the lower boundary and the northernmost infill (Fig. 8).

Thrust complex 3 is found in the eastern part of the study area bounded by reflections r₃ and r₄ (Figs 3 and 9). On the time slice the complex is seen as lineations oriented NW–SE, bending slightly towards the NE at the northernmost limit of the complex (Fig. 3). The complex covers a total area of 187 km²; however, this area is truncated by tunnel valleys splitting the total complex into three smaller regions, all of which are interpreted to be the same complex because of the similarity of the thrusts observed (Fig. 3). Within each of these regions two distinct types of thrusts/folds are observed, described as areas F and G (Figs 3c and 9). In area F, the thrusts/folds consist of northeastward-dipping, high-amplitude reflections (Fig. 9), which weaken towards the SW along the décollement surface. Adjacent to this, area G consists of medium-amplitude reflections dipping NE that are stacked obliquely in an imbricated manner (Fig. 9), similar to the thrusts observed in area B of the southern thrust complex. The upper boundary of the complex is truncated by the overlying reflections, interpreted as reflection r₄. In area F, the thickness of each unit between reflectors is c. 70 – 90 m and in area G c. 50 – 80 m (Fig. 9). In the NW of the third thrust complex, a more chaotic thrust pattern is observed on the time slice, probably owing to the interaction with thrust complex 2 (Figs 3 and 4).

Above the glacial thrust complex, the undeformed post-thrusting stratigraphy consists of medium-amplitude, semi-horizontal reflections with a lower boundary seen as reflections r₄ (Figs 5 – 9). These reflections are truncated by an extensive erosion surface (Figs 5 – 9), giving the unit a variable thickness of c. 50 – 90 m dependent on the depth of the erosive surface above. The major erosion surface truncates all three thrust complexes across the study area (Figs 3 and 8). The geometry of this surface in cross-section and time slice suggests a correlation to the regional Ling Bank unconformity formed by the large network of Middle to Late Pleistocene tunnel valleys.

The overall stratigraphic framework of the thrust complex, correlated to the regional geographical setting, places the thrusts between the Brunhes–Matuyama Palaeomagnetic Reversal (r₂, 780 ka ± 5 ka BP, MIS 19) and the Ling Bank Uncomformity (r₄, c. 450 ka, MIS 12), in agreement with the chronostratigraphic calibration of the subsurface succession to core data. The thickness of stratigraphy between the Brunhes–Matuyama and the thrust complex (c. 45 – 72 m) and between the thrust complex and the Ling Bank unconformity (c. 50 – 90 m) strongly suggests that the thrust complex, representing a large-scale glaciotectonic complex, pre-dates the Elsterian glaciation (MIS 12), with the most likely interpretation at MIS 16. The thrust complex has been subdivided into three areas representing a minimum of three ice sheet advances and retreats originating from the north and NE based on thrust orientation and cross-cutting behaviours.

The pre-thrust units 1 and 2 correspond to the Aberdeen Ground Formation of the central North Sea seismic stratigraphy as described by Stoker & Bent (1985), Cameron et al. (1987, 1992), Sejrup et al. (1991), Gatliff et al. (1994) and Stoker et al. (2011). The Aberdeen Ground Formation is described as hard, heavily over-consolidated grey clay with sporadic shell fragments and plant remains, occasionally containing lenses and laminae of silt and fine-grained sand (Gatliff et al. 1994; Graham et al. 2011; Stoker et al. 2011). It is interpreted mainly to correspond to a deltaic to shallow glaciomarine environment during a time period when sea-level lowstands became increasingly deeper and the North Sea accommodation space decreased (Stoker & Bent 1985; Cameron et al. 1987; Sejrup et al. 1987, 1991; Knudsen & Sejrup 1993).

Reflections r₁ and r₂ bound unit 1 and correspond to the base of the Jaramillo Palaeomagnetic Event (c. 1.07 Ma; Figs 5 and 6; Stoker et al. 1983; Channell & Kleiven 2000; Lisiecki & Raymo 2005) and Brunhes–Matuyama Palaeomagnetic Reversal (780 ka ± 5 ka BP; Figs 4 – 6; Stoker et al. 1983) respectively. Reflection r₃, representing the top surface of unit 2, is interpreted as the décollement surface within the Aberdeen Ground Formation and dips slightly towards the SSW (Figs 6 – 9). Décollement surfaces are basal detachment surfaces associated with compressional settings such as folding and thrusting. The décollement surface was probably formed when ice sheet loading induced sufficient differential lateral surface stresses to weaken already over-pressured clay, causing deformation along the horizon (Huuse & Lykke-Andersen 2000a; Andersen et al. 2005).

The thrust complex comprises a system of folds and thrust faults that slid on the basal décollement surface (Aber & Ber 2007) within the Aberdeen Ground Formation. It is located above the Brunhes–Matuyama reflection (r₂) and below the Ling Bank Unconformity, thus representing deposits formed sometime during MIS 19 – 13. Given the palaeogeographical location of the thrusts, and the lack of evidence for large-scale structural events, it is most probable that the thrusts were part of a large glaciotectonic complex formed as a result of the deformation of sediment beneath or in front of an ice sheet. The sediments themselves are described as shallow glaciomarine (Gatliff et al. 1994; Graham et al. 2011; Stoker et al. 2011), which lends further support to the interpretation, as it is probable that a shallow marine basin may become terrestrial during glacial lowstands and a large ice sheet more easily explains extensive thrusting in such sediments where extension owing to slope failure is more likely than compression. The complex shows a resemblance to glaciotectonic structures as observed elsewhere in NW Europe by Klint & Pedersen (1995), Huuse & Lykke-Andersen (2000a), Andersen et al. (2005), Buckley (2012) and Pedersen (2014). Notably, the morphology of the thrusts in the distal region of the complex appears to resemble that of thrusts that have previously been associated with gravity spreading (area A and F; Figs 3, 7 and 9; e.g. Huuse & Lykke-Andersen 2000a; Pedersen 2000; Bennett 2001; Andersen et al. 2005) whereas the more proximal imbricated stacks resemble push moraines (area B, C, D and G; Figs 3 and 7 – 9; e.g. Aber & Lundqvist 1988; Klint & Pedersen 1995; Huuse & Lykke-Andersen 2000a). A discussion of each of the complexes follows below.

Thrust Complex 1 is located in the northeastern part of the study area (Fig. 6) and consists of multiple stacked anticlines forming a mounded feature towards the SW with a depression directly adjacent towards the NE. This is observed to be very similar to the typical morphology of a hill–hole pair. Hill–hole pairs consist of a ridged hill and depression at or close to an ice margin and can be associated with glaciotectonic deformation in frozen bed conditions (Figs 3, 4b and 6; Aber et al. 1989; Sættem 1990; Bennett 2001; Aber & Ber 2007). Material is carved out of the depression by the advancing ice sheet and then thrust up to form a hill. The preservation of the overlapping anticlines towards the SW of the thrust complex 1 (Fig. 6), and the lack of truncation of these features, suggests that the ice front did not advance over the entire structure at the time of formation, nor did any later ice advance disturb or remove the hill-feature. Based on the onlapping of horizontal reflections onto the sides of the hole-feature (Fig. 6) it appears that the depression was infilled post-formation of the complex, although it cannot be ruled out that ice may have partly re-covered the depression post-formation. The orientation of the thrust sheets forming the hill (strike NW–SE; main raft dip NE) suggests an initial ice flow from the NE creating the hill–hole pair.

Thrust Complex 2 is located in the northern part of the study area (Fig. 7) and comprises a southern and northern complex separated by a basin-like structure. The reflection pattern of the southern complex is very different from that of the northern complex, indicating differences in the ice force-mechanisms creating the thrusts, such as push–squeeze, thrusting under the ice load and/or ice load causing gravity spreading. Based on the evidence discussed below, thrust complex 2 is interpreted to have been formed during at least three ice advances.

The thrust pattern of the southern part of thrust complex 2 is subdivided into areas A and B based on a difference in the seismic reflection pattern (areas A and B; Fig. 7). The area A thrusts, as described above, consist of a series of anticlines that are backed by the steeply dipping, imbricated thrusts of area B. These characteristic reflections closely resemble detached thrust structures observed elsewhere in the North Sea (Huuse & Lykke-Andersen 2000a; Andersen et al. 2005) in terms of geometry, size and stacking pattern (Fig. 7). The difference in thrust pattern observed between area A and area B is interpreted to have been caused by different force mechanisms acting ahead of a single ice advance: a combined gravitational spreading created the thrusts in area A (Fig. 7) and pushed from the rear (i.e. classical push moraines) in area B (Fig. 7).

To visualize the mechanisms creating the imbricated push moraines, we refer to the schematic model and classifications of push moraines produced by Bennett (2001). Imbricated push moraines are classified by their distinct morphological diversity and formation and schematic models and aspects (Bennett 2001). Specifically, the imbricated thrusts of area B shows significant similarities to a multi-crested push moraine transverse to the direction of the ice flow. Multi-crested push moraines are imbricated stacks of steeply inclined thrust sheets with each thrust rising from a basal décollement surface, and further classification indicates that the imbricated thrusts resemble fold–thrust-dominated moraines (Bennett 2001). Under this model of ice push from the rear, friction along the décollement plane must be relatively low for the creation of the wide thrust belt (Bennett 2001).

Evidence to support this is found in the eastern part of the southern complex, where south-dipping thrusts have been identified. These southward-dipping thrusts resemble the thrusts formed at Møns Klint, Denmark (Pedersen 2000), where superimposed thrusting was caused by a readvance from a different direction. The south-dipping thrusts were originally imbricated thrusts that have been reorganized as a result of a later ice advance from the NE, most probably the advance that created thrust complex 3 (see below). Similarly, the truncation of the upper boundary of the thrusts within the southern complex may represent glacial truncation by an additional ice sheet overriding the complex.

The northern part of thrust complex 2 is subdivided into two areas, C and D, related to the identified change in the seismic reflection pattern (Fig. 8). The pattern of reflections is similar to that observed in the southern part of the complex, suggesting a similar formation of push from the rear; however, the typically higher angle of normal thrusts in area D, with additional reverse thrusts containing diapiric folds, is suggestive of greater friction at the base of the ice that may suggest formation during permafrost conditions (Bennett 2001). The southern flank of the area D thrusts drapes over the northern tip of the area C folds (Fig. 8), although there is no direct truncation. This indicates that the thrusts and folds were formed during two separate ice sheet advances. Additionally, the preservation of the crests of the thrust sheets indicates that the ice did not advance over the thrust belt after deformation. Given the evidence for overriding ice in the southern complex and the probable ice advance direction (north to south), the northern part of the thrust complex (areas C and D; Fig. 8) is probably younger than the southern area (areas A and B; Fig. 7), perhaps created during a later ice advance.

The basin feature within area E is interpreted as a relatively low area separating the two parts of the thrust complex. Infilling of the basin most probably took place during deglaciation or interglacial times prior to erosion of tunnel valleys that occurred during a subsequent glaciation, leading to a lack of structural and stratigraphic continuity (Figs 3c and 8). The formation of this basin between the two thrust sheets of complex 2 is key evidence for proglacial deformation.

The curved lineaments seen on the time slice (Fig. 3) represent the stratigraphic reflections from the thrust blocks and can be used as a proxy for thrust strike direction. The thrust lineaments in thrust complex 2 have an arcuate form from east to west (Figs 3 and 4a) and cross-sections (Figs 5 – 9) indicate that the thrusts show a dominant north-to-south compression direction. Thus, the complex was most probably formed by ice flow from a northern direction.

Thrust Complex 3 is located in the northeastern part of the study area (Fig. 3) and consists of imbricated thrust structures (area G; Fig. 9) and folded anticlines that weaken along the axis of deformation (area F; Fig. 9). Similar to thrust complex 2, the change from folded anticlines to steeply dipping imbricated thrust sheets is interpreted to be evidence for differential forces ahead of an advancing ice sheet. The folded anticlines were formed by gravity spreading as a function of the ice load (area F; Fig. 9b), whereas the steeply dipping thrust sheets were formed by push–squeeze processes ahead of the advancing ice (area G; Fig. 9). The imbricated structure of thrusts (area G; Fig. 9) resembles multi-crested push moraines oriented transverse to the direction of the ice flow, much like the imbricated thrust in thrust complex 2. The thrusts originated from a basal décollement surface with a low to medium degree of friction, creating a wide thrust belt, and deformation may have occurred under permafrost conditions.

The separation of thrust complex 3 from thrust complex 2 is identified by differences in the dip direction of the thrusts. The lineation of the thrust blocks in thrust complex 3 on timeslices, from NW to SE, as well as the dip of thrusts in cross-section towards the NE (Figs 3, 4c and 9) indicate that these features were formed by a northeastern ice advance, rather than a northerly one. The upper parts of the thrusts in thrust complex 3 are observed to be truncated by overlying reflections linked to possible overriding by the ice that formed the thrust complex.

A minimum of three ice advances are interpreted based on the dip direction of the thrusts and on truncations of the complex by overriding sediments. Because of the resolution of the dating available it is unclear whether these ice advances represent full retreat–advance episodes or fluctuations along a dynamic, lobate ice front. Without further information either interpretation is possible; however, the details of each ice advance can be considered relative to each other. The initial phase of the ice sheet advance into the study area was from the NE, forming the hill–hole pair in the eastern part of the study area. This advance did not override the hill–hole feature at this time. Following this event, ice advanced from the north and the southern area forming thrust complex 2, which was subsequently overridden by the ice sheet. As the ice sheet retreated and advanced again from the north, the northern part of thrust complex 2 was created. The timing of formation of thrust complex 3 may have been simultaneous with the formation of either the southern or the northern part of complex 2. Although thrust complex 3 is overridden by ice in the same manner as the southern portion of complex 2, as evidenced by the truncation of thrusts by overlying sediments, there is evidence of superimposion of complex 3 on the eastern margin of the southern complex. This is suggestive of a later event, perhaps concurrent with the northern complex. However, this needs more research to be further constrained. It is known that ice lobes can lead to thrusting in different directions (Aber et al. 1989; Pedersen 2000; Andersen 2004), which may be how the thrusting was created with ice flow directions from both north and NE (Fig. 10). This was probably a dynamic ice front.

Undeformed sediments, which are seen as parallel reflections on seismic cross-sections, are deposited on top of the thrust complex after deglaciation (Figs 5 – 9). Various studies using 3D seismic data from the North Sea concur with the interpretation of extensive tunnel valley erosion, and in total at least seven glacial cycles have been identified (Graham et al. 2011; Stewart & Lonergan 2011; Stewart et al. 2012). Stewart & Lonergan (2011) estimated that the age of the oldest buried tunnel valleys was probably from the Elsterian glaciation.

This 50 – 90 m thick unit overlies the thrust complex, and is truncated by large erosional features. In some places, these features cross-cut into the glacial thrust complex itself. The erosional features are interpreted as the widely mapped network of subglacial tunnel valleys, of Elsterian to Weichselian age, that cross the North Sea (Huuse & Lykke-Andersen 2000b; van der Vegt et al. 2012; Stewart et al. 2013). The presence of a significant thickness of undeformed sediments between the thrust complex and the tunnel valleys indicates that the thrusts were formed during an earlier ice advance, allowing a long period of time of undisturbed deposition. Considering the substantial thickness of the intervening interval, it seems likely that the glacial tectonic episode belongs to a significant pre-Elsterian glaciation (i.e. in MIS 16).

The formation of glaciotectonic thrusts leads to horizontal shortening of the stratigraphy owing to lateral movement of sediments by the ice. To better understand the development of the glacial deformation outlined in this paper, we calculated the extent of glacial tectonic shortening. By measuring the length of the individual thrusts, adding them together and measuring the distance the thrusts cover, the shortening can be estimated along cross-sections perpendicular to thrust strike. For the southern part of thrust complex 2, the total length of all the thrusts combined is c. 11 km whereas the thrust complex itself is c. 6 km in length. This gives a shortening of 5 km or c. 40%. In the northern part of thrust complex 2 the total length of the thrusts sheets is c. 5 km in a complex 2.5 km long giving a shortening of 2.5 km or 50%. For thrust complex 3, the cumulative length of the imbricated thrust blocks is c. 7 km in a complex 4.5 km in length, giving a shortening of 2.5 km or 35%.

Various researchers have calculated horizontal shortening in different thrust complexes; Huuse & Lykke-Andersen (2000a) suggested horizontal shortening of c. 40% of glaciotectonic thrust structures of Elsterian and possible Saalian age in the eastern Danish North Sea, Pedersen (2005) calculated c. 50% shortening at the Rubjerg Knude Glaciotectonic complex in Denmark, and Croot (1987) showed a total shortening of 54% in present-day glaciotectonic structures found in Iceland. The calculated horizontal shortening in the present study shows values of 35 – 50% and is thus in the range of the studies mentioned above. This indicates that the deposits were moved a significant distance through glacial processes. The calculated shortening of the thrust complex may underestimate the total shortening owing to the upper truncation of the imbricated thrusts in the southern part of complex 2 and in complex 3.

The maximum and minimum age of the glaciotectonic thrust complex is constrained by the Brunhes–Matuyama Reversal below and the tunnel valley unconformity above (Stewart & Lonergan 2011; Stewart et al. 2013). Tunnel valleys have been identified from the Elsterian glaciation to the Last Glacial Maximum (c. 21 ka BP); however, dating of the erosional features is difficult and is restricted to the absolute dating of infrequent interglacial deposits found within the tunnel valley fill and relative dating based on cross-cutting relationship.

Further constraining the age of the glaciation is made difficult by the lack of more detailed chronostratigraphic information. What few cores have penetrated below the Ling Bank unconformity have had little detailed work done on them to establish glacial and interglacial events preserved in the sediments. The thickness of the sediment package between the thrusts and the Ling Bank unconformity (r₄) strongly implies that the thrusts were formed significantly prior to the tunnel valleys, and hence the Elsterian, and the same can be said of the period between the Brunhes–Matuyama Reversal (r₂) and the décollement surface (r₃). However, this leaves a significant period from MIS 18 to MIS 13 in which to narrow down the probable glacial period.

The use of the marine oxygen isotope record as a proxy for the glacial–interglacial cycles is considered to provide age estimates. The global record of δ¹⁸O measured from benthic foraminifera is considered to be a proxy for global ice volumes, as ice preferentially takes up ¹⁶O compared with ¹⁸O when water freezes, enriching seawater and encouraging preferential uptake by benthic foraminifera leading to an increase in the δ¹⁸O ratio (Raymo 1994; Miller et al. 2011). Attempts have been made to correlate the δ¹⁸O curve to the degree of glacial ice extent. For example, Kleman & Stroeven (1997) built a framework for the Fennoscandian glaciation under the assumption that the Fennoscandian glaciation varied in approximate alternation and proportion with the global ice-volume signal. However, this model is fraught with caveats because the significant local variations in precipitation, as well as the influence of temperature, can have a large effect on the correlation of a local record to the global averaged record.

In this study, we use the δ¹⁸O record primarily to identify periods of glacial activity globally, by identifying Marine Isotope Stages (MIS) representing the alternating glacial and interglacial cycles. In terms of the thrust complex studied, there are four interglacial stages (δ¹⁸O <3.7‰; MIS 19, 17, 15 and 13) and three glacial stages (δ¹⁸O >4.5‰); MIS 14 (δ¹⁸O ∼4.5‰), MIS 16 (δ¹⁸O <5‰) and MIS 18 (δ¹⁸O ∼4.65‰) (Fig. 2; Lisiecki & Raymo 2005; Channell et al. 2009) during the period between the Brunhes–Matuyama Reversal and the best-estimate date of MIS 12 for the Ling Bank Unconformity. Any one of the three glacial stages could have resulted in a significant Fennoscandian Ice Sheet forming the glaciotectonic complex observed in this study. It is, however, considered that of these three the most likely stage for this glaciation is MIS 16. As well as being the largest and longest sustained of the three glacial stages, evidence for an extensive MIS 16 glaciation corresponding to the Donian Stage has been established in Russia (Astakhov 2004; Velichko et al. 2004).

Estimates of accumulation rates during this period, assuming an MIS 16 age, seem to agree with this interpretation. Up to c. 135 m of undeformed sediment accumulated after the Brunhes–Matuyama Reversal and prior to thrusting (calculated from the undeformed sediment package immediately adjacent to the thrust complex) whereas c. 90 m of undeformed sediment accumulated post-thrusting and prior to the Ling Bank Unconformity. This results in relatively stable accumulation rates of 0.75 and 0.6 m ka⁻¹ for an MIS 16 date, whereas using either an MIS 18 or MIS 14 date creates a significant bias towards one of the two time periods. These calculations, however, are estimates only and do not take into account the significant glacial erosion that probably occurred during this period owing to the overriding ice sheets.

The current understanding of the palaeogeography of the pre-Elsterian central North Sea is limited by the available data and the shallow focus of most previous studies. Shallow cores indicate that the top of the Aberdeen Ground Formation is pro-deltaic to glaciomarine in nature (Cameron et al. 1987; Sejrup et al. 1987; Buckley 2012) whereas seismic stratigraphy suggests that the glaciotectonic complex discussed sits on the topset of a small clinoform (c. 100 m in height, Fig. 5), significantly smaller than those observed further south in the earlier part of the Quaternary and Pliocene (Cameron et al. 1987; Rasmussen et al. 2005; Lamb et al. 2017a,b). Clinoforms are commonly used to estimate water depth within the basin and for clinoforms of 100 m in height water depths are likely to fall into the region of 100 – 120 m. This suggests that during the period of glaciotectonic deformation the North Sea was probably a shallow and broadly flat sea during interglacial periods; however, whether the subaqueous conditions continued during glacial periods is more difficult to discern and would have been the result of complex interactions of ice extent and thickness, flexural isostacy, eustacy, sediment and meltwater flux, and the presence or absence of dammed meltwater lakes.

During glaciations in the youngest 1 myr, under the influence of 100 kyr cycles, the North Sea was broadly terrestrial during glacial periods and large ice-dammed lakes have been hypothesized to have existed in the North Sea as a result (Gupta et al. 2007; Murton & Murton 2012; Cohen et al. 2014; Sejrup et al. 2016). Evidence for grounded ice comes from the presence of subglacial features such as tunnel valleys (Kristensen et al. 2007; Stewart & Lonergan 2011; Stewart et al. 2012, 2013), mega-scale glacial lineations (Sejrup et al. 2003; Graham et al. 2007) and large numbers of moraine ridges preserved on the sea floor (Bradwell et al. 2008), although the presence of glacial lakes is mainly inferred from erosional sequences and depositional structures close to the ice front (Murton & Murton 2012; Cohen et al. 2014). The shift from terrestrial to shallow marine between glacial and interglacial cycles is a reflection of high-amplitude changes in global sea level, glacial isostatic adjustments and dynamic topography, with the Weichselian to Holocene sea-level change being of the order of 120 m. During the Middle Pleistocene, eustatic changes were substantially lower amplitude (typically <60 m), and the basin is thus less likely to have been subaerial during glacial lowstands.

The global sea-level curve for the Quaternary as produced by Miller et al. (2011) gives approximate estimates for sea-level stands during the glacial–interglacial cycles. Although care must be taken to consider the global curve in light of patterns rather than exact figures, owing to regional and local influences on sea level, it provides a robust dataset for suggesting testable patterns. On the sea-level curve of Miller et al. (2011), MIS 16, the estimated age of the glaciotectonic complex, shows a sea-level change closer to that of the later 100 kyr cycles (>100 m change) than that of the early 41 kyr cycles (<60 m) (Miller et al. 2011). The probably shallow water depths of 100 – 120 m coupled with the glacial lowstand suggests that it is not unreasonable to assume that during the glacial period during which the thrust complex formed, the shallow marine setting became terrestrial. Even if the sea-level changes during the glacial–interglacial cycle were not sufficient to expose the sediments to a terrestrial environment, the marine basin would have been sufficiently shallow (<50 m) that it could not support the significant weight of floating ice.

The glaciotectonic complex indicates that the ice sheet that formed the complex had to be grounded in the deformed region. Hill–hole pairs, such as the one in thrust complex 1 (Fig. 6), are generally formed at the ice margin (Sættem 1990; Bennett 2001) whereas thrust complexes 2 and 3 are interpreted to be ice-proximal push moraines, with truncation of the imbricate thrusts observed in both complexes, which can be taken as suggestive of overriding of the complexes by an ice sheet to a certain degree (Figs 7 and 9). In grounded ice sheets the interplay between permafrost and meltwater controls the deformation of the sediments below the ice sheets (Etzelmüller et al. 1996; Boulton et al. 1999; Bennett 2001; Madsen & Piotrowski 2012). However, it is difficult to draw conclusions based on the presence of thrust structures alone. The fill of the basin-like feature within area E (Figs 3c and 8) may be an indication of meltwater features within the study area. However, the limited extent of the high-resolution seismic data, and the influence of the tunnel valleys from the regional unconformity above the succession, may prevent the observation of these features more clearly. The presence of permafrost has low preservation potential and probably similarly is uncertain, as textures and features that may act as indicators for past permafrost are likely to be below seismic resolution. To fully assess the meltwater–permafrost setting a close study of a series of cores would be needed to combine with the seismic data (e.g. Harris et al. 1997), although numerical modelling suggests that large-scale glaciotectonic deformation can happen without permafrost (Andersen 2004).

The length and height of thrust blocks provide clues to the glacial load and the presence of a décollement layer at depth, and the degree of deformation was also influenced by the thrust sediments and the porewater available. Length and height data for the thrust blocks from this study and others have been plotted in Figure 11 to build a database of geometric data that can be used to evaluate similarities and differences in glaciotectonic architecture during different glaciations. The thrust complexes originate both on- and off-shore and from different glaciations: Weichselian (Klint & Pedersen 1995; Harris et al. 1997; Pedersen 2000; Jakobsen & Overgaard 2002; Rafaelsen et al. 2007; Thomas & Chiverrel 2011; Pedersen & Boldreel 2015), Saalian (Andersen 2004), Elsterian (Huuse & Lykke-Andersen 2000a) and pre-Elsterian glaciation (Buckley 2012).

The plot shows a wide size range of glaciotectonic complexes throughout northern Europe (Fig. 11). The largest thrusts are found offshore whereas the onshore thrusts are of a smaller magnitude. However, no pronounced pattern is observed, suggesting that there are no significant differences between glacial thrust complexes at present located onshore compared with those located offshore. The size of the thrusts in the present study is well within the range of the sizes interpreted in other studies. The thrust complexes found offshore may be better preserved owing to greater offshore subsidence and sedimentation rates offsetting erosion during ice advances.

The glaciotectonic complex described in the present study therefore can be interpreted to have formed at an oscillating ice front from a large, grounded ice sheet during the MIS 16 glaciation, at which time the North Sea was probably largely terrestrial, as in later glaciations. The ice flow direction was from the NE, suggestive of a large Fennoscandian Ice Sheet source. The implications of the North Sea record on a larger scale are dependent on the size of the British Ice Sheet during the same period; a damming of the northern end of the basin would result in the formation of a large lake south of the study area. The study area lacks any evidence of the British Ice Sheet; however, it is not unreasonable to assume that a Fennoscandian Ice Sheet of this size would have met a correspondingly large British Ice Sheet, although the notion of confluence between the two currently remains conjectural. Further research into the area SE of the thrust complex–ice-front using similar high-resolution datasets may reveal evidence of a proglacial lake or free-flowing drainage conduit at the MIS 16 glacial limit. This would give further evidence for the MIS 16 glaciation reaching the central North Sea and any potential for ice sheet confluence during this period.

Interpretation of high-resolution 3D seismic data from the central North Sea has led to the discovery and subsequent characterization of a well-preserved glaciotectonic complex. The thrusting that occurred affected the Aberdeen Ground Formation and adds to the complexity of the Pleistocene stratigraphy.

The initial stage of glaciation of the study area is interpreted to be in the form of a hill–hole pair, suggesting a frozen bed condition. This is followed by multiple advance and retreat phases of the ice sheet forming two thrust complexes during a minimum of three glacial advances. Analysis of the thrust complex structure suggests that the two primary processes of formation were gravity spreading and push from the rear (push moraines). The morphology of the thrust complex resembles that of thrusts associated with gravity spreading (areas A and F; Figs 3, 7 and 9) in the distal parts and imbricated push moraines in the more proximal part (areas B, C, D and G; Figs 3 and 7 – 9). During deformation of thrust complex 2, a basin was formed between the thrust sheets.

The interpreted imbricated thrusts are classified as fold–thrust-dominated moraines of considerable width seen on the time slice (Figs 3 and 4) as thrust blocks and may have been formed under permafrost conditions; however, that is not a requirement for deformation. The calculated horizontal shortening of the thrust complex is 35 – 50% with individual thrust blocks being some 3 – 500 m long by 30 – 70 m thick. These parameters are comparable with those for thrust complexes of Elsterian, Saalian and Weichselian age of on- and offshore NW Europe.

Chronostratigraphic and stratigraphic evidence suggests that this complex is from an extensive pre-Elsterian glaciation (pre-MIS 12). We suggest that thrusting from the grounded ice sheet most probably occurred during the Donian MIS 16 (630 ka; Fig. 2), some 200 kyr earlier than previous ice sheet models predicted. This ice sheet originated from the north and NE, probably from southern Norway (Fig. 10). This discovery highlights the fragmented knowledge of the glaciation history of NW Europe and the under-explored record of glaciation within the North Sea Basin and other offshore sediment archives.

We would like to thank CGG for providing the Broadseis^TM dataset, and PGS Ltd for providing the MegaSurvey^TM dataset. We are grateful for their permission to publish the examples shown herein. Schlumberger is thanked for provision of the Petrel license at the School of Earth, Atmospheric and Environmental Sciences, University of Manchester, UK and Department of Geosciences and Natural Resource Management, Geology section, University of Copenhagen, Denmark. We would like to thank the reviewers for their constructive comments, which have led to a greatly improved paper.