The Heligoland Glacitectonic Complex in the southeastern North Sea: indicators of a pre‐ or early‐Elsterian ice margin

Glacial landforms are abundant in the North Sea basin and are often used to reconstruct the impact and dynamics of ice sheets during the Pleistocene. Geophysical methods have allowed the mapping and structural analysis of glacial landforms at the surface and in the subsurface to estimate the position of former ice margins in the North Sea. However, the glacial history of the southeastern North Sea remains underexplored. In this study, we present a structural analysis of Late Pliocene to Late Pleistocene sediments based on a dense grid of 2D high‐resolution multi‐channel reflection seismic data from the German sector of the southeastern North Sea. We show that the Heligoland Glacitectonic Complex (HGC) is larger than previously assumed (700 km2, 32×22 km) and characterized by three distinct zones of thrusting and deformation on two décollements. The kinematic restoration of seismic cross‐sections and dip measurements of thrust faults demonstrate that the HGC was formed by an ice lobe advancing from the southeast. To explain the origin of the HGC, we provide alternative models for its formation during a single ice advance or two ice advances in the study area. Furthermore, we validate the early or pre‐Elsterian age of the HGC based on nearby Elsterian tunnel valleys, and conclude that salt structures in the subsurface may have influenced its location.

Northwest Europe holds spectacular examples of landforms formed during the Pleistocene glaciations where geomorphological records document landscapes shaped by highly dynamic ice sheets and meltwater. Numerous studies have demonstrated the significance of glacial meltwater in the formation of associated glacial landforms in this region (Piotrowski 1994(Piotrowski , 1997Huuse & Lykke-Andersen 2000a;Ehlers & Gibbard 2004;Kristensen & Huuse 2012). The intense analysis of outcrops and seismic data have described a multitude of glacial landforms, some of which show an orientation that can be used to deduce their formation mechanism and the direction of ice movement, such as drumlins, mega-scale glacial lineations, eskers, grounding-zone wedges, tunnel valleys, and glacitectonic complexes (GTCs) ( O Cofaigh 1996;Huuse & Lykke-Andersen 2000a;Boulton et al. 2007;Dowdeswell et al. 2008;Bennett & Glasser 2009;Stewart et al. 2013;Ottesen et al. 2014). GTCs, in particular, may be used as a direct indicator of ice movement if these structures have not been overprinted by subsequent ice sheets; thrust sheets separated by thrust faults are oriented perpendicular to the former ice margin and may, therefore, be used to derive the direction of ice movement based on kineto-stratigraphy (Berthelsen 1978).
GTCs consist of a number of thrust sheets intersected by reverse faults, for which the following three main driving mechanisms have been described: gravity sliding, gravity spreading, and 'push from the rear' (Elliott 1976;Chapple 1978;Siddans 1979;Pedersen 1987). Gravity sliding is defined as the deformation of a sedimentary body by down-slope movement under its own weight, which is typically considered a mass-wasting process (Pedersen 1987). Whilst reverse faults have been recognized at the toe of such mass-wasting deposits, most authors conclude that this scenario is unlikely to result in large glacitectonic complexes (Aber & Ber 2007). Gravity spreading is commonly defined as the formation of listric reverse faults pointing towards the foreland, that is induced by a static load acting on the sediments. This static loading results in a lateral pressure gradient and lateral stress and, subsequently, thrusts can exert additional static load, which leads to thrust fault development limited to a distance in the foreland at which the shear strength of the sediments is higher than the lateral stress (Bucher 1956;Pedersen 1987;Aber & Ber 2007). In contrast, the 'push from the rear' mechanism, or 'bulldozing', assumes a moving sheet that exerts lateral stress on the foreland, resulting in the formation of thrusts by simple displacement on a weak basal layer (Chapple 1978;van der Wateren 1985). As there are few glacigenic examples that can be attributed to gravity sliding, many studies conclude that either gravity spreading or 'push from the rear' mechanisms-and often in combination-were the dominant drivers of glacitectonic disturbance (Aber & Ber 2007).
The analysis of GTCs differs for onshore and offshore realms. Furthermore, the structural elements of GTCs are easily accessible in onshore outcrops and may be measured in situ (Schwarzer 1985;van der Wateren 1995;Lee & Phillips 2013;Phillips et al. 2013;Pedersen 2014;Livingstone & Clark 2016;Gehrmann et al. 2019) whereas offshore analyses rely on seismic data. Highresolution 3D seismic data would allow a detailed structural analysis of offshore GTCs, but data of sufficient resolution are rarely available (Rafaelsen et al. 2007). Therefore, 2D reflection seismic data are typically used to study GTCs. This approach is an approximation based on the available profile orientations and the limitation of measurements of apparent dips and strikes. Consequently, a high number of crosssections is necessary to estimate the direction of ice movement (Pedersen & Boldreel 2017). Increasing the density of 2D seismic grids facilitates the reconstruction of GTCs, which in turn can facilitate the interpretation of former ice margins (Cotterill et al. 2017;Vaughan-Hirsch & Phillips 2017;Phillips et al. 2018).
Glacitectonic thrusts have been compared to thinskinned tectonics (Croot 1987;Aber et al. 1989) and mimicked in sandbox experiments (K€ oster 1959;Schlische et al. 2014) to understand the mechanics of deformation in detail. Based on the mechanics of glacitectonic deformation, the concept of restoring balanced crosssections has been introduced to validate interpretations of the strata before deformation (Croot 1987;Woodward et al. 1989;Pedersen 2005Pedersen , 2012Pedersen , 2014Gehrmann et al. 2019;Winsemann et al. 2020). This restoration of crosssections helps to quantify the displacement of the source material, which in turn provides information about the geometry of the ice margin and thrust direction (van der Wateren 1995; Pedersen & Boldreel 2017;Winsemann et al. 2020). On this basis, glacitectonic landforms may be used as a direct indicator of ice-sheet dynamics and combined with the mapping of till provinces that have traditionally been carried out to determine the extent of former ice sheets (Berthelsen 1978;Winsemann et al. 2020).
Our study aims to investigate the structural evolution of the Heligoland Glacitectonic Complex (HGC) northwest of Heligoland. We focus on (i) the true spatial extent of the HGC, (ii) the geomorphic parameters of the HGC, and (iii) the reconstruction of ice motion and ice margins during the formation of the HGC. To address these issues, we present high-resolution seismic reflection data that provide new insights into the HGC in unprecedented detail.

Geological setting
The German sector of the North Sea is part of the Central European Basin System (CEBS) that stretches from the North Sea to Poland and from Norway to the German midlands (Maystrenko et al. 2008). The sedimentation of evaporites during the Zechstein period in the Late Permian resulted in widespread halokinesis across the CEBS, which produced deep rim synclines and local subsidence (Maystrenko et al. 2005b). The formation of the North Sea dates back to the Late Jurassic to Early Cretaceous, when the uplift of its margins was forced by the tectonic regime due to the opening of the northeast Atlantic Ocean and the Alpine Orogeny (Ziegler 1990(Ziegler , 1992. Subsequently, up to 3000 m of Oligocene to Holocene sediments were deposited, which have been influenced locally by salt structures (Caston 1979;Overeem et al. 2001;Maystrenko et al. 2005a;Th€ ole et al. 2014).
During the Middle Miocene to Early Pleistocene, the Eridanos River connected the North Sea Basin via the CEBS and supplied sediments to the Eridanos Delta (Gibbard & Lewin 2016). The southeastern North Sea, in particular, is an area of deposition strongly influenced by the Eridanos River discharge, with sandy, silty and clayey material dominating the units (Overeem et al. 2001;Th€ ole et al. 2014;Gibbard & Lewin 2016), and which form the basis of the glacigenic sediments deposited during the Pleistocene (Huuse 2002;Gibbard & Lewin 2016). The influence of local salt structures on sedimentation patterns and structures formed during the Pleistocene is still under debate (van der Wateren 1995; Sirocko et al. 2008;Lang et al. 2014;Wenau & Alves 2020).
An increasing body of evidence shows that ice sheets reached far into the North Sea and controlled erosion and sedimentation processes during the Pleistocene (Cameron et al. 1987;Wingfield 1990;Ehlers & Wingfield 1991;Praeg 2003;Moreau et al. 2012). Recent studies have led to an ever-finer discrimination of ice advances and suggested highly dynamic ice sheets in the North Sea (Lonergan et al. 2006;Graham et al. 2007;Stewart & Lonergan 2011;Ottesen et al. 2020). The North Sea, in particular the German sector of the southeastern North Sea between Amrum and Heligoland, was not covered by ice sheets during the Last Glacial Maximum, yet glacial deposits and glacigenic structures are widespread in the subsurface of the North Sea as a result of earlier glaciations (Ehlers et al. 2011;Batchelor et al. 2019). Compared to the Miocene and Pliocene deposits, ice sheets during these earlier glaciations deposited a wider range of grain sizes in the form of shallow glaciomarine and glaciolacustrine sediments across vast areas of the North Sea that were locally disturbed by subglacial erosion and proglacial deformation (Huuse & Lykke-Andersen 2000a, b;Moreau et al. 2012;Ottesen et al. 2020;Wenau & Alves 2020).

Material and methods
We acquired 2D reflection seismic and parametric sediment echo-sounder data to image subsurface structures and landforms from 0 to 800 ms (approximately 600 m) beneath the sea floor. A total of 1058 km of 2D high-resolution multi-channel reflection seismic data were acquired in a closely spaced 2D grid covering BOREAS approximately 800 km 2 with a mean profile spacing of 800 m (Fig. 1B).

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

Seismic data processing and interpretation
The data processing followed standard reflection seismic processing routines in Schlumberger's Vista Seismic Processing software: (i) Ormsby frequency filtering (55/ 110 and 600/1200 Hz filter flanks), (ii) CMP binning, (iii) FK-filtering, (iv) normal moveout (NMO) correction using avelocity model derived from the seismic data, (v) de-spiking, (vi) CMP-stacking, and (vii) finitedifference migration of the data in the time domain with a velocity model inferred from an interactive velocity analysis. A shot point distance of 9 m or less allowed the common midpoint (CMP) bin spacing to be set equal to the channel spacing of 1.56 m, resulting in a CMP fold of 6 to 10.
Special processing involved the removal of multiple reflections that interfered with the primary reflectors owing to the shallow water depth (<25 m). Removal of both sea floor and internal reflector multiples was accomplished using fine-tuned predictive deconvolution with an adapted lag time based on picking of the sea floor and the first multiple. This approach yielded acceptable results owing to the flat sea floor throughout the study area. The fully processed seismic sections were loaded into IHS Markit Kingdom Software (2018) for interpretation. The Kingdom software was used to trace the horizons and the morphology of the subsurface landforms. The horizons were then exported for interpolation and plotting using the Generic Mapping Tools (GMT) open-source software.
Petex' MOVE software was used for the kinematic restoration of the cross-sections and to measure the apparent strike and dip of faults (note that the true dip likely deviates). MOVE's module '2D Move-on-Fault' using the Trishear method produced the most plausible results when restoring both thrusts and folds to the predeformation state based on key reflectors. As a result of the restoration, we calculated the shortening of the section by comparing its deformed length with the undeformed length.

Results
The stratigraphical units and subsurface landforms shown in this section were identified based on differences in the acoustic properties of the corresponding seismic sections and correlated with the stratigraphical framework of the Federal Agency of Geosciences and Resources of Germany (BGR). We omit a detailed seismo-stratigraphical analysis of the seismic sections as our focus is on glacial landforms. Rather, we provide a summary of the stratigraphy of the study area (Table 1); a more detailed description is provided by Lohrberg et al. (2020).

Stratigraphy of the study area
The seismic sections image the sedimentary sequence of the upper Miocene to the present (Fig. 2, Table 1). The lowermost part of the seismic sections consists of reflectors characterized by an alternating medium to low amplitude pattern, which dominates the complete unit SU1 (Fig. 2). A narrow high-amplitude reflector characterizes the top of this unit, which represents the mid-Miocene unconformity ( Fig. 2; Huuse & Clausen 2001;. The mid-Miocene unconformity is truncated by the Quaternary succession (Lutz et al. 2009) in the central and southern parts of the study area, and is, therefore, absent in some parts. The reflectors below the mid-Miocene unconformity ( Fig. 2; approximately 0.45 s and deeper) represent lower to mid-Miocene delta deposits (Th€ ole et al. 2014).
The mid-Miocene unconformity marks an abrupt change in seismic facies from high-amplitude to lowamplitude reflectors in SU2 (Fig. 2); the characteristic pattern of alternating high to low amplitude reflectors, as seen in SU1, is absent here. Reflectors have low amplitudes but are highly continuous and likely reflect sandy to silty sediments.
The deposits of SU3 (up to 0.35 s) most likely represent Pliocene to Early Pleistocene sediments because they are constrained by late Miocene sediments at their base and Pleistocene sediments at their top ( Fig. 2). Due to the limitations of the stratigraphical framework, an unambiguous correlation with Pliocene sediments was not possible. The d ecollement reflectors D1 and D2 constrain SU3; D2 is found at the top of the unit, whereas D1 is found at the bottom of the unit (Fig. 2).
The deposits of SU4 range from approximately 0.25-0.05 s (Fig. 2). This part of the seismic section is characterized by alternating strong and weak reflectors in its lower part, which transitions to increasingly chaotic strong reflectors in the upper part ( Fig. 2; 0.15-0.25 vs. 0.05-0.15 s). Unambiguous tracing of reflectors is feasible only in the lower part of SU4. The deformation of the strata starts from 0.25 s and continues upwards. Reflectors R1-R4 were introduced for this unit as they formed the basis for the reconstruction of deformation (see below). SU4 mainly comprises Pleistocene glaciomarine and glaciolacustrine sediments and is truncated at its top throughout the study area. The shallowest seismic unit is SU5; its top represents the sea floor and its surface is flat throughout the study area. The SU5 base truncates the structures beneath. This truncation is represented by an erosional surface marking the base of the unit. The reflectors on the top of this erosional surface represent Holocene deposits, consisting of shallow marine sands.

Glacitectonic deformation
The seismic sections in the western part of the study area show strong layering in the subsurface and significant dips of reflectors in SU4 (Fig. 3). The full profiles show that the dipping reflectors are underlain by non-dipping undisturbed reflectors (SU1, SU2 and SU3) and that the dipping reflectors are bound to the southeast by a tunnel valley and a large filled depression.
A sequence of four main reflectors (R1-R4) was identified across the sections due to their alternating pattern of low and high amplitudes on top of a prominent, near-horizontal reflector D2 (Fig. 4A). R1-R4 are disturbed and separated into single sheets by steep and often curved fault planes every few hundred metres. They have been truncated at their top, and discordant horizontal reflectors overlie the erosional truncation. At the southeastern end of the individual sheets, R1-R4 are almost horizontal, whereas their dip increases toward the northwestern end of the individual sheets. Seismic amplitudes decrease from southeast to northwest, correlating with the increasing dips of the reflectors.
Considering the geometry of the faults as well as the dips of the reflectors and the similarity of the seismic facies (SU4 with R1-R4), we identified the individual sheets on top of D2 to be comprised of the same material but separated by reverse faults (Fig. 4A). Therefore, the sheets are termed 'thrust sheets' from here on, and D1/D2 are referred to as d ecollements (detachments). Despite the alternating pattern of reflectors, their overall amplitudes decrease from the top to the bottom inside each thrust sheet, and the amplitudes are generally inversely correlated with the dip of the reflectors (Fig. 4A). Amplitudes are lowest where the strata are highly deformed, potentially on a sub-resolution scale (Fig. 4B, C).
Similar to the individual sheets imaged on level D2, we identified a section along which individual sheets developed on a deeper level D1 (Fig. 5). On top of D1, reflectors show overall lower amplitudes, and the tracing of reflectors is ambiguous. Therefore, it is not possible to trace key reflectors and estimate the shortening. However, the structures are very similar to the sheets on level D2, which qualifies these structures as thrust sheets at a deeper level.

Location and extent of the HGC
A series of thrust sheets and a large-scale deformation over several kilometres qualifies this set of structures as a glacitectonic complex or hill-hole pair. Because of its proximity to the island of Heligoland, it has been named the HGC (Winsemann et al. 2020). The extent of the HGC was identified based on thrusts and folds in the northwestern part of the study area and extends for approximately 170 km² in the available seismic sections; however, the correlation with the results of Winsemann et al. (2020) reveals that the thrust-fold-belt stretches far  (2020). BOREAS beyond our study area. Combining the extent of glacitectonically disturbed strata in both studies and the associated source depression, the extent of the HGC was found to be exceptionally large (700 km²), stretching for 32 km across and 22 km in the thrust direction (Fig. 6A). The coverage of seismic sections is good, yet the density of seismic sections decreases towards the west where the HGC is best imaged (Fig. 6A). A total of 163 faults on two levels (D1 and D2) were identified in the 70 seismic sections (Fig. 6B). The true number of faults in the study area is likely lower as the number given here is based on tracing the faults in single seismic profiles, and most faults are likely to extend across several seismic profiles. Therefore, the true number of faults can only be determined by correlating and tracing individual faults across multiple seismic profiles, which was not feasible with the 800-m line spacing of our data. Indeed, the correlation between individual fault planes across sections was only possible at the intersections of these seismic profiles.
Below the thrusts and folds, we identified an upper intact d ecollement horizon (D2) from which the thrusts arise at steep angles, and along which the sediment packages have been sheared (Fig. 5). D2 is located approximately 240 m below the sea floor and is identified on the western boundary of the study area. Below D2, we identified a second d ecollement horizon (D1) in the centre of the study area, which can be traced 3.8 km towards north-northwest and ends abruptly (Fig. 5). Thus, the thrusts and folds emerging from D1 are limited in lateral extent to a few kilometres. The thrusts and folds connected to D1 are smaller than those connected to D2 (note that the tops of most of the D2 thrusts have been eroded) and seismic amplitudes are generally lower (Fig. 5). The tops of the D1 folds have partly been modified by thrusts sheared from D2, producing an erosional truncation of most of the D1 folds and thrusts. However, D2 undulates where D1 is present beneath, while D2 is entirely flat beyond the extent of D1. The most western D1-fold fades out and leaves undisturbed strata to the west, similar to the thrusts and folds on D2. The thrusts and folds on both d ecollements D1 and D2 are limited to the east by a large depression, which is separated by a prominent reflector M1 (Fig. 5).

Restoration of deformation
On average, the thrust faults dip towards the southeast (Fig. 6B), indicating that thrusting was directed toward the northwest. The dip angles are nearly identical for the faults emerging from D1 and D2 (Fig. 6B). Note that the dip angles are apparent dips owing to the limitations of 2D seismic sections. A comparison of the restored cross-section with the present state reveals a net shortening of approximately 5.1 km, which translates to 33% (Fig. 7).

The Heligoland Glacitectonic Complex, SE North Sea
The translation of material left a source depression towards the southeast of the glacitectonic complex, which has been filled with undifferentiated sediments similar to a glacial basin (Fig. 8) and is large enough to supply the thrusted strata (Fig. 3). While its eastern boundary is not covered by seismic data, its western and northwestern boundaries are well imaged on seismic sections in the centre of the study area. In these sections, a high-amplitude reflector (M1) was identified, which was dipping towards the southeast and separated the large source depression from the thrust-and-fold belt in the northwest (Figs 3, 8). Below M1 towards the west, the amplitudes and continuity of the reflections are low. Farther to the west, the first thrust sheets become identifiable, characterized by increased continuity and higher amplitudes of the well-laminated thrusted strata. Above M1 towards the east, the continuity and amplitude of the reflections increase. The fill of this source depression is characterized by thin, high-amplitude, and weakly stratified reflectors compared to well-stratified thrusted strata. Depressions of different scales are found in the source depression infill, which shows no morphological link to M1 and the thrusted strata below (Fig. 8).

Zones of deformation
The data show that the HGC can be subdivided into (i) a proximal ice-contact zone with steep faults (up to 35°), (ii) an inner proglacial zone characterized by lower angle thrusts and folds, and (iii) an outer proglacial zone dominated by long-wavelength folds and decreasing thrust angles, which flattened out and resulted in the absence of glacitectonic disturbance (Figs 3, 9).
Zone I: ice proximal. -The zone of steepest thrusts extends over approximately 3.7 km north-northwest. The apparent dip of the faults ranges from 21°to 51°and all thrusts are based on D2 (Fig. 9A). Due to the steepness of these thrusts, their tops were significantly higher and, therefore, prone to erosion. Accordingly, their tops are not preserved in Zone I, which precludes the identification of hanging-wall anticlines and associated structures. Zone I shows the highest density of faults along D2, which led to strong compression of the strata in individual thrust sheets (up to 21% compared to the thrust flat). In rare cases, a clear footwall ramp is visible, and the base of the thrust sheet can be traced down to D2.

BOREAS
In most cases, however, the thrust sheets appear to emerge directly from D2 and it is difficult to assess where the thrusts are initiated.
Zone II: transition zone. -The transition zone is characterized by a lower dip of thrusts and hanging-wall anticlines at a lower vertical level such that they are partly preserved in the data (Fig. 9B). The dip of thrusts decreases toward the west, whereas the lateral length of individual thrust sheets increases (Fig. 9B). In most cases, the top of the hanging-wall anticlines is eroded such that no full anticline is preserved in the transition zone. The compression of the strata in the deformed part is considerably lower than in Zone I.
Zone III: ice distal. -The ice distal part of the data records complete folds in the subsurface (Fig. 9C). The folds often show flat undisturbed strata towards the former ice front that develop into steep folds directed away from the former ice front (Fig. 9C). Synclines that develop at the end of hanging-wall anticlines or between two hanging-wall anticlines are preserved in this zone. The overall compression of the strata is low in Zone III except for the zones of compression in the centre of the folds. Patches of chaotic low-amplitude reflections often coincide with the hanging-wall anticlines above.

Proximity to an ice-sheet margin
The distinct morphological characteristics of Zones I to III of the thrusts on D2 represent different distances from the former ice margin. Such zonation is a characteristic feature of propagating thrust stacks Phillips 2018). Zone I is closest to the ice margin, which is reflected in the steepest thrust angles, the highest spatial density of thrusts, and the longest displacements along the faults (Fig. 9A). Here, the forward movement, geometry, and weight of the ice margin bulldozed the initial pile of sediments. This led to subsequent failure and additional thrusts along the d ecollement due to the lateral stress introduced by gravitational forces of the ice margin and the initial pile of sediments. As a result, a sequence of thrusts with a decreasing dip of fault planes formed towards the foreland. Farther from the ice margin, the lateral stress decreased and the deformation changed from thrust-dominated to internal deformation, which is reflected in the overturning of the hangingwall anticlines (Fig. 9B). This transitional zone is termed Zone II. Deformation most distal to the ice margin further transitioned to the fold-dominated type, where fault planes were less frequent. This part is referred to as Zone III. Here, internal deformation allowed for such strong overturning of hanging-wall anticlines that reflectors in the zones of strongest deformation are imaged with significantly lower amplitudes compared to the rest of the seismic sections due to transmission loss (Fig. 9C).
The location of the ice margin may be partly reflected by the position of M1, the thrust directions, and the extent of the thrust-and-fold belt. M1 directly overlies the initial thrusts imaged in Zone I, thereby marking the westernmost limit of the source depression (Fig. 8).
Considering the transition of the apparent dips of the faults (Fig. 9), an arcuate ice lobe extending towards the northeast best fits the observations. Accordingly, M1 may represent the grounding line of the ice sheet that carved out the glacial basin that was filled during deglaciation. Extensive blanking by shallow gas in the south of the study area precludes tracing M1 towards the south and, thus, inhibits the recognition of the ice margin in this area.
Two d ecollements -How many phases of thrusting? D1 and D2 show different extents and stratigraphical levels. The lack of undisturbed sediments between the two d ecollements indicates a contemporaneous formation of thrusts and folds at both levels. However, none of the thrusts are linked, which indicates some kind of separation. Although D2 thrusts may have displaced or eroded sediments formed between the formation of D1 and D2, we could not identify any disturbance of the sedimentary sequence north of the thrusts. Here, the full sedimentary sequence and the absence of truncations strongly indicate undisturbed sedimentation, which indicates the contemporaneous formation of D1 and D2. Given the lackof conclusive evidence, we suggest two possible scenarios for the formation of the observed d ecollement levels, both of which are illustrated in Figs 10 and 11 and discussed in the following sections.
Scenario I -Two phases of thrusting. - Figure 10 illustrates a possible mode of formation for two separate thrust phases and resulting thrusts and folds. In this Fig. 9. A composite seismic section showing the zones of thrusting and folding along a north-northwest to south-southeast transect in the west of the study area. The location of the profile is shown in Figs 1 and 6. Note that measured dips are apparent dips due to the fixed orientation of the seismic section. BOREAS scenario, the sediments were deposited prior to their deformation on either d ecollement, implying that there must have been some time for sedimentation between the phases of deformation.
In Scenario I, an ice sheet reached the area and its static load translated to lateral stress close to the ice margin as a result of its morphological profile (Aber et al. 1989;Andersen et al. 2005). Excess meltwater introduced into the subsurface under high hydraulic heads led to overpressures in distinct subsurface layers in front of the advancing ice margin, which activated faulting in the foreland along zones of weakness due to decreased normal pressures (Rubey & Hubbert 1959;Mathews & Mackay 1960;Croot 1987;Hart et al. 1990;Boulton & Caban 1995;Andersen et al. 2005). Once the initial fault was established, the ice margin pushed the sediment pile (or thrust sheet) up and along the developing thrust fault. At this stage, the thrusted sediment pile forced lateral stresses due to the static load and, thus, stimulated the activation of the next thrust fault in its foreland bygravity spreading (Andersen et al. 2005). The lateral stress imposed by the static load faded out towards the foreland (possibly during a temporary stillstand) of the ice margin, which led to decreasing compressive stresses and decreasing displacement in the foreland. Farther into the foreland, lateral stresses were compensated by internal deformation rather than displacement, which resulted in folding (Fig. 10C). This mechanism continued during the advance of the ice margin and led to the displacement of Pliocene and Early Pleistocene strata over several kilometres toward the west and northwest up to depths of 150 m. At some point, the flow of ice outpaced the displacement during the deformation, possibly following a temporary standstill. As a result, the ice sheet advanced over the thrusts and folds, thereby eroding their tops and preserving only their lower and middle parts.
During the retreat of the ice sheet, diamicton and glaciolimnic/glaciomarine sediments were deposited on top of the erosional surface and filled the source depression from which the sediment piles originated. The dipping erosional truncation M1 remains a relic of the temporarily grounded ice margin and serves to identify its position during the initiation of thrusting (Figs 8, 10E).
A re-advance of the ice margin then led to the displacement and thrusting of the glacial sediments deposited on top of the erosional surface during the first phase, using the erosional surface as the upper d ecollement level D2. Displacement along this erosional surface led to increased abrasion by the sliding sediment piles, thereby increasing the smoothness of the erosional surface. During the melting of the ice sheet, glacial sediments were likely again deposited on top of the thrusts and folds as well as in the source depression.
Truncation of D1 thrusts and folds as well as a second source depression would be indicative of Scenario 1; however, this model fails to explain the intact hangingwall anticline in Fig. 5 at the westernmost end of D1 deformation, and only one source depression was identified in the data.
Scenario II -One phase of thrusting. - Figure 11 illustrates an alternative mode of formation to produce the same subsurface structures by a single ice advance. In this scenario, all sediments were deposited prior to their deformation on both d ecollements. This would only be feasible in a parallel mode of formation for thrusts and folds on both d ecollement levels.
The formation of two d ecollements during one ice advance follows the principles illustrated in Fig. 10. In contrast to Scenario I, a second suitable d ecollement horizon already existed in the subsurface before an ice sheet reached the area. The static load of the ice margin was so high that lateral stresses were imposed at great depths, thereby activating both d ecollements at the same time. Similar to Scenario I, excess meltwater facilitated faulting, and gravity spreading led to the progression of thrust sheets. At some point, the internal friction must have been sufficiently high to cease deformation and propagation along the lower d ecollement (D1). The continuation of deformation along the upper d ecollement (D2) suggests a change in the environmental setting during the formation of the glacitectonic complex. It is possible that the farther propagation along D2 may be associated with a temporary stillstand of the ice margin and a subsequent re-advance mediated by D2.
The intact hanging-wall anticline in Fig. 5 at the westernmost end of D1 deformation and the lack of truncation of D1 thrusts and folds are key elements in the interpretation of the formation of the glacitectonic complex during a single thrust phase. Still, we suggest that a temporary stillstand of the ice margin is likely. The truncation of the D2 thrusts and folds indicates that the ice margin progressed over the glacitectonic complex after its formation.

Time of formation in the context of other glacigenic structures
Several tunnel valleys found in the study area cut the glacitectonic complex ( Fig. 6; Lohrberg et al. 2020). This unambiguously shows that the glacitectonic complex is older than the tunnel valleys. Without exact time control, we can only rely on the presumed Elsterian age for the tunnel valleys (Lutz et al. 2009;Lohrberg et al. 2020). Consequently, the glacitectonic complex formed during the Elsterian glaciation or earlier. Batchelor et al. (2019) compiled evidence for ice sheets dating back to MIS 24, and showed that the Elsterian ice sheet (MIS 12) predominantly advanced over Central Europe from a northerly direction. Modelling also shows that there was an ice advance into the North Sea from the east during MIS 16 (Toucanne et al. 2009;Batchelor et al. 2019). This would pre-date the presumed MIS 12 incision of the tunnel valleys and, thus, reconcile both the relative time of formation and the thrust and incision directions. Inevitably, evidence for pre-MIS 12 glaciation is sparse Fig. 11. Model of formation for Scenario II with one ice advance. The sequential schematic plots illustrate the formation of thrusts and folds on both d ecollement levels D1 and D2 by the same ice margin, resulting in thrusting in the form of a duplex stack. BOREAS due to the erosive nature of subsequent ice sheets. Nevertheless, considering that there could have been a balance between deposition and erosion during glaciations, it may be possible that only the top of the HGC was eroded during subsequent glaciations such that its middle and lower parts have been preserved since MIS 16. However, due to the lack of accurate dating of interstadial deposits, this remains speculative. Individual ice lobes were agile and did not have to move parallel to the overall long-term forward direction of the ice sheet. Nevertheless, it seems peculiar that an ice lobe exceeding 30 km in width (Figs 6A, 12B) progressed in the opposite direction of the main advance direction. This is especially true in the absence of obvious geomorphological barriers that could result in ice-sheet splitting. Bendixen et al. (2018) have shown glacitectonics dating back to the MIS 16. Their data document an ice advance in the central North Sea that produced spectacular glacitectonic disturbances of the strata. Owing to the availability of dated materials, the formation of the glacitectonic complex can be unambiguously inferred to pre-date the formation of tunnel valleys in the same area. Despite the subsequent major glaciations (MIS 12, MIS 6-10 and MIS 2-5), the glacitectonic complex has been preserved as an MIS 16 geomorphological record. Although the glacitectonic complex identified by Bendixen et al. (2018) lies deeper than the HGC (216-288 vs. 30-100 m b.s.f.), a similar setup may still be observed. The location and depth of the features initially preclude a direct comparison based on the assumption of similar sedimentation rates andwithout considering subsidence; however, sedimentation rates in proximity to the landmasses of northern Germany and southern Denmark are highly variable, and the influence of tidal erosion is difficult to estimate. Accordingly, the possibility that most of the Quaternary deposits are missing in the southeastern North Sea should be considered because of the interplay of erosive forces, subsidence of the North Sea basin, and Neogene uplift of Scandinavia (Huuse 2002). These considerations are in accordance with the regional stratigraphy of our data, which shows generally shallow Miocene and Palaeocene reflector packages close to the shore (Lohrberg et al. 2020). Thus, we suggest that our data record mainly the older glaciations of the North Sea (MIS 16 to MIS 12) with a minor influence from MIS 6-10 (up to 15 m b.s.f.) and a weak influence from MIS 2-5 in the form of subaerial drainage systems (up to 8 m b.s.f.).

Location of the HGC and nearby salt structures
Two salt pillows underlie the Cenozoic strata in the study area (Lokhorst 1998). No direct spatial correlation has been established between glacial landforms with subsurface salt structures in the past, even though spatial correlations seem likely (van der Wateren 1995; Wenau & Alves 2020). Considering the widespread permafrost at the time of glaciation, the increased heat flux above salt Fig. 12. A southeasterly ice advance from the Baltic region is the most likely scenario for the formation of the Heligoland Glacitectonic Complex (HGC). A. Overview map for regional context. B. Detailed map illustrating the advance of an ice margin into the study area from the southeast. The 'hill' represents the areawhere the strata have been thrusted and deformed, whereas the 'hole' represents the source depression left behind after the displacement of the material by the ice margin. pillows may have been a slight push on the tipping point, leading to the formation of these structures (Sirocko et al. 2008;Lang et al. 2014). Proximal to salt structures, the permafrost depth was likely reduced so that glacial meltwater could penetrate into highly permeable Neogene aquifers confined by clay layers. Fault zones above the salt structures may serve as a means by which meltwater recharged the subsurface and facilitated the build-up of subsurface overpressures (Wenau & Alves 2020). Distal to the salt structures, the deep permafrost may have acted as abarrier for meltwater drainage, which inevitably led to high overpressures in the aquifers. The resulting overpressure generated a weak layer proximal to the salt structure. The overlying ice sheet excavated a glacial basin and bulldozed the material, thereby inducing a lateral stress regime (Figs 10, 11;Vaughan-Hirsch & Phillips 2017). With the weak layer preconditioned, large thrust sheets were able to move along the slippery d ecollement layer in the form of a propagating thrust stack (Croot 1987).
The gentle dip of the d ecollement towards the west indicates that a morphological barrier is not needed for bulldozing to take effect. Rather, it seems that an overpressured aquifer, possibly linked to a rapidly advancing ice sheet, may be the key prerequisite for the formation of a glacitectonic fold-and-thrust belt (Sigf usd ottir et al. 2020). Increased heat flow over salt structures in the subsurface controls the depth of permafrost across vast areas and, thus, the hydraulic connection to aquifers (Sirocko et al. 2008). Therefore, the distribution of salt structures may influence the zones of overpressure, thereby indirectly controlling the spatial distribution of glacitectonic complexes in the North Sea. To further investigate this relationship, spatially extensive high-resolution data sets covering wide parts of the North Sea are needed for a statistical evaluation.

Conclusions
Our high-resolution 2D reflection seismic data set images the HGC in unprecedented detail. A comparison with earlier studies shows that the HGC extends farther east than previously thought; its areal extent (approximately 700 km 2 ) makes the HGC one of the largest submarine glacitectonic complexes known worldwide.
Structural analyses of seismic sections showed that the HGC may be subdivided into three distinct zones characterized by different dip angles and densities of thrust faults as well as different degrees of internal deformation. These zones outline an arcuate ice margin that advanced from the southeast. The two levels of thrusting and their different lateral extents may be explained by thrusting during a single or, alternatively, two separate ice advances. The intact hanging-wall anticline at the westernmost end of the lower level of thrusting, the lack of erosional truncation of lower thrusts and folds, and the presence of only one source depression favours the model of thrusting during a single ice advance along two pre-conditioned d ecollements.
Age control is highly limited owing to the lack of deep wells. The incision of large Elsterian tunnel valleys through parts of the HGC implies that the HGC is of early Elsterian age (MIS 12) or could be the result of a pre-Elsterian glaciation in the North Sea.
Although large salt structures led to the updoming of Cenozoic strata below the HGC, a direct spatial correlation with the glacial landforms remains uncertain. It seems likely, however, that the salt structures influenced the thermal regime beneath the ice sheets and may have contributed to laterally uneven meltwater supply below the ice sheet. This possibly led to increased pore-water pressure, which may have facilitated the formation of the HGC.