Deep-water sedimentation processes on a glaciated margin: The Foula Wedge trough mouth fan, West of Shetland

Trough Mouth Fans (TMF) are sedimentary depocenters located at glaciated continental margins and consist predominantly of glacigenic debris flow deposits. The Foula wedge is a Pleistocene TMF accumulated offshore West of Shetland over the Northeast Atlantic margin. This study presents an analysis of a 3D seismic reflection dataset imaging the distal Foula wedge basin fan deposits between 1010 and 1100 m water depth, directly downslope from a gully system which was active untill the end of the last deglaciation. Results reveal, in un- precedented detail, the basal surface of this fan system and its internal complex architecture. Features typical of both debris flow deposits and turbidites are identified, including a basin channel network with linear and diverging erosional features forming distinctive terminal lobes, stacked and backstepping. The study links the seafloor morphology of the basin fan with its subsurface geomorphology, showing connection with the down- slope gully system to the east. It presents evidence for a complex distal depositional system on glaciated margins, characterised by heterogeneous sediment delivery processes and deposits. A conceptual evolution model is proposed, with a glacigenic debris flow-dominated TMF at the LGM, subsequently influenced by meltwater discharges, with deposition occurring as a function of the shelf margin and slope paleo-morphology, slope substrate composition, interaction of downslope and along slope processes and ice-margin dynamics.

This study focuses on the distal part of the Foula wedge TMF, offshore West of Shetland in the north-eastern Atlantic Ocean (Figs. 1 and 2). It uses 3D seismic reflection data to characterise the architecture of the basin fan deposits, located downslope from a gully system and provides new insights regarding the style of formation of this deep-water fan system. The aim is to understand the development of the fan system within the context of the glaciated margin but also to identify the main sediment transport processes. Building upon existing literature (e.g. Stoker and Varming, 2011;Stewart and Long, 2012), this work further investigates the role of the substrate lithology on the evolution and preservation of the downslope gully and fan systems. By integrating the results with available chronological information, an updated conceptual  Fig. 2. B) Location of study area (red box; 3D globe map from EMODnet, 2020). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 2. A) Regional seafloor shaded relief map illuminated from the north-east encompassing the West Shetland slope (WoS) and Faroe-Shetland Channel. The study area (3D seismic data extent) is delineated by the blue polygon. Bathymetry is derived from the 3D seismic, AFEN, 2000 and EMODnet, 2020 dataset combined with bathymetry maps from Long et al., 2004 andLong, 2016, geo-referenced andcomposed in PetrosysPRO. This shows seafloor morphology associated with the LGM within the study area and its vicinity. Shaded colours refer to the regional geological framework (after Stoker, 1995, Masson et al., 2003Knutz and Cartwright, 2003;Knutz and Cartwright, 2004;Davison, 2004;Hohbein and Cartwright, 2006;Stoker and Varming, 2011;Stewart and Long, 2012). BGS borehole and sub bottom profiler line locations are also shown for reference (full details in Table 1); B) A representative downslope geoseismic section (for location see 2A) showing shallow seismic stratigraphy across the Foula wedge encompassing the study area (after Stoker, 1995;Davison, 2004;Smallwood, 2004;Stoker and Varming, 2011). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) model for the evolution of this sector of the margin is proposed, suggesting an increased role for meltwater processes soon after the Last Glacial Maximum (LGM).

Geological setting
The study area is located on the eastern side of the Faroe-Shetland Channel (Fig. 1), which is part of the Faroe-Shetland sedimentary basin. A major structural high, known as the Cambo high, represents the pre-Cretaceous basement overlain by a thick post-rift succession that includes a Cenozoic sequence of marine deposits (Stoker and Varming, 2011;Ellis and Stoker, 2014;Stoker, 2016;Hardman et al., 2019). This sequence was impacted by extensive Paleocene/lower Eocene igneous activity (Hardman et al., 2019) and successively characterised by shallow to deep-water sedimentation (Lamers and Carmichael, 1999;Stoker and Varming, 2011;Stoker et al., 2013). The late Neogene was characterised by onshore uplift resulting in tilting of the margin with accelerated offshore subsidence. This, combined with increased sediment supply due to Pleistocene glaciations, changed the configuration of the Faroe-Shetland basin and lead to the narrower present-day basin architecture (Stoker and Varming, 2011).

Plio-Pleistocene to present
The Late Pliocene to Pleistocene depositional history of the margin is controlled by a combination of processes resulting from glacialinterglacial cycles related to ice margin variations and relative sea level fluctuations (Kurjanski et al., 2020), and bottom currents (Knutz and Cartwright, 2003), leading to intermittent phases of erosion, deposition, and reworking (Stoker and Varming, 2011 and references therein;Stewart and Long, 2012;Bradwell et al., 2019). On the east side of the Faroe Shetland Channel, over the continental slope offshore West of Shetland, the Pleistocene sequence consists of glacial sediments deposited by an unusually dynamic ice sheet. This sequence forms coalescing sedimentary wedges hundreds of metres in thickness with a Table 1 Additional available data within the study area and key regional geological background data sources.

Data type and source
Value for this study Data: Shallow gravity cores logs. Format: PDF. Source: Siccar Point Energy. Reference: Fugro, 2011 (Confidential). Additional information on shallow sediments in Stewart and Long, 2012. Reported information and general seafloor sediment descriptions were used to support seafloor sediment interpretations as part of the seafloor morphology provinces mapping. Average locations from Fugro, 2011 images plotted on Fig. 4A. Stewart and Long (2012) study extent shown in Fig. 2A.
Data: seafloor photography still images around 204/5a-1 exploration well. Format: jpg. Source: Siccar Point Energy. Reference: Fugro, 2011 (Confidential) Data: Seafloor photography still images around 204/10a-5 exploration well. Format: Jpg. Source: Siccar Point Energy. Reference: Fugro, 2017 (Confidential) Locations were plotted on seafloor maps and related images were used to ground truth seafloor sediment interpretations and supporting seafloor morphology province classification. Average locations plotted on Fig. 4A To provide greater regional context. Support morphology correlation with slope/shelf regional area to understand deep-water sediment delivery processes. All data were imported/converted and combined using PetrosysPRO to then compose the various bathymetry maps ( Figs. 1,2A,12,13 and 14) and derive the profile for Fig. 2B. Edited RGB raster were used to infill gaps or extent maps to provide greater regional context. seaward migration of the shelf-break up to 50 km (Stoker and Varming, 2011;Stewart and Long, 2016;Ballantyne and Small, 2019;Bradwell et al., 2019). Along the south-east of the Faroe Shetland Channel this includes the prograding glacial Rona and Foula wedges and, to the far north-northeast, the North Sea Fan (Stoker, 1995;Stoker and Varming, 2011;Figs. 1 and 2). Directly to the northeast of the study area, the stratigraphy also comprises an up to 225 m thick contouritic deposit, part of the Plio-Pleistocene West of Shetland Drift (WSD-slope) -which pinches out to less than 50 m in thickness towards the southwest Cartwright, 2003, Knutz andHohbein and Cartwright, 2006; Fig. 1).
The study area is located in water depths ranging from 700 m to 1100 m (Figs. 1 and 2), covering the southernmost part of the distal Foula wedge basin floor fans Long, 2012, 2016), and the southern end of the West Shetland Drift (Knutz and Cartwright, 2003;Knutz and Cartwright, 2004). These deposits are defined by Stoker and Varming (2011) as part of the Faroe Shetland Neogene 1 (FSN-1) megasequence. On the West of Shetland slope, although the boundary remains poorly defined, pre-glacial Pliocene to mid Pleistocene Morrison 1 sequence and the glacigenic middle to upper Pleistocene sequence known as Morrison 2 are overlain by the late to post-glacial MacAulay sequence (Fig. 2B, C). Borehole data from the shelf break 82/11 (Rona wedge), slope 85/01 (Foula wedge) and the base of slope 99/3 (distal Rona wedge) (Table 1), together with legacy British Geological Survey (BGS) sub-bottom profiler lines (Fig. 2), show that the Morrison 2 sequence consists predominantly of glacigenic sediments. These include massive clay diamicton debris flow packages interbedded with thin glaciomarine and contouritic sandy-clay and sand layers underlying late glacial glaciomarine sediments deposited from suspension and postglacial hemipelagic clay with sandy-clay of the MacAulay sequence (Stoker, 1995;Stoker and Varming, 2011). The Morrison 2 sequence thins downslope, and it is locally absent, passing laterally into the lower slope/basin areas comprising the undifferentiated Plio-Pleistocene Faroe-Shetland Channel sequence (including the WSD, as part of the FSN-1 megasequence, Stoker and Varming, 2011).
Sedimentation rates experienced during the Holocene are typically less than 1 cm/kyr and strong bottom currents, with estimated velocities up to 1.0 m/s in deep-water (Masson, 2001), control modern sediment deposition and reworking of the Pleistocene deposits (Masson, 2001;Masson et al., 2003Masson et al., , 2004Stoker and Varming, 2011;Stewart and Long, 2012). The present-day seafloor morphology (Fig. 2) is therefore relict from the LGM (Stoker and Varming, 2011) and can be divided into: 1) an upper slope iceberg ploughmark zone; 2) an upper/mid slope comprising GDF lobes (as part of the Foula wedge); 3) a downslope gully system and triangular shaped glacigenic basin floor fans which extend up to 15 km into the Faroe-Shetland Channel associated with the Foula wedge (Stoker and Varming, 2011;Long, 2012, 2016). Based on Stoker and Varming (2011), the basin floor fans form the upper part of the Morrison 2 sequence and are overlain by less than 5 m-thick deposits of the MacAulay sequence (Fig. 2).  (Fugro, 2017) showing seafloor sediment types within SMP2, and a well-defined comet mark developed around a large boulder; the coarser sediments are associated with glaciomarine sedimentation e.g., iceberg dropstones. C) Transverse bathymetry profile (P1) showing seafloor incision cutting across SMP2. D) Downslope 3D seismic section (P2) encompassing SMP1 and SMP2 displaying the seismic character of the underlying vertical section (unit 1 and unit 2a) and showing the large bathymetric hollow present over the west edge of SMP1.

Glacial history
While a generally accepted model for the glaciation of this region exists (Clark et al., 2018;Hall et al., 2019;Ballantyne and Small, 2019;Bradwell et al., 2019 and references therein,), it lacks detailed and/or spatially continuous chronological constraints, especially in the marine realm. Extensive shelf wide glaciation around the British Isles probably occurred throughout the Pleistocene (Rea et al., 2018), but detailed information regarding the extent of the ice masses and sediment input onto the West of Shetland slope before the LGM are lacking (Ballantyne and Small, 2019;Bradwell et al., 2019). During the last glaciation, marine-terminating ice expanded across the continental shelves within the period 35-32 k yr BP (Ballantyne and Small, 2019).
The topography and bathymetry appear to have directed westward flowing ice between the Orkney and Shetland islands to the Rona and Foula wedges Varming, 2011, Stewart andLong, 2016;Ballantyne and Small, 2019;Bradwell et al., 2019;Hall et al., 2019;Fig. 1). The sediment sequences within the wedges, together with moraine evidence, indicate ice-proximal deposition associated with the arrival of a combined British Irish ice sheet (BIIS) and the Fennoscandian ice sheet (FIS) at, or near, the West of Shetland shelf-break between 26 and 25 k yr BP (Stoker and Varming, 2011;Clark et al., 2018;Hall et al., 2019;Ballantyne and Small, 2019;Bradwell et al., 2019). By 23 k yr BP, following separation of the BIIS and FIS an independent ice sheet/ice cap (Shetland ice cap) may have dominated over the Orkney-Shetland platform (Bradwell et al., 2019). The Shetland ice cap was highly dynamic between ~21 to 18 k yr BP, including multiple ice lobe re-advances in varying flow directions associated with high calving fluxes in a marine-influenced sector (Bradwell et al., 2019) directly to the east of the study area. The marine terminating ice sheet margin retreated to near the present-day Shetland coastline between 17.0 and 16.5 k yr BP and deglaciation of the land was completed by 15 k yr BP (Bradwell et al., 2019).

Available data
The 546 km 2 dataset used for this study is a single azimuth 3D seismic reflection full-stack time-migrated volume, acquired in 2012 and provided by Siccar Point Energy (Fig. 2). The data are zero phased and the polarity convention is such that an increase in acoustic impedance is a negative number, displayed as a black trough in the seismic sections ( Fig. 3-D). The sample interval is 4 milliseconds (ms), and the in-lines and cross-lines are at 12.5 m spacing, which can be assumed equal to the horizontal resolution. The dominant frequency in the section down to 100 ms below seafloor is approximately 50 Hz yielding a theoretical vertical resolution of approximately 8 m, assuming an average formation velocity of 1600 m/s and using the ¼ wavelength Rayleigh method (Kallweit and Wood, 1982;Simm and Bacon, 2014). The bathymetric map is derived from the two-way travel time (TWT) to depth converted seafloor pick, using a reference sound velocity in water of 1467 m/s. The regional bathymetry maps are derived from multiple datasets in Petrosys PRO which includes depth raster images, gridded bathymetry data (AFEN, 2000;GEBCO, 2014;EMODnet, 2020) and bathymetric map images extracted from published studies (e.g. Long et al., 2004;Long, 2012, 2016) geo-referenced for this study (a fuller description wit data source and references is presented in Table 1.) Additional data, listed and described in Table 1, were used to assist the interpretation and integration with the regional geological framework (Figs. 1, 2, 4A, 12, 13 and 14).

Seismic data interpretation
Interpretation of the 3D seismic data was performed using both IHS Kingdom and Open DTect software to maximise attribute analysis. Resulting maps were compiled in Petrosys PRO software.
The seafloor seismic reflection and three other significant subsurface reflections were traced by a combination of both manually and 3D auto tracking functions. A single average formation velocity of 1600 m/s was used to convert subsurface reflections from time to depth domain.
Seismic facies analysis follows techniques described by Mitchum et al. (1977).
To enhance interpretation and visualisation of seismic attributes, including root mean square (RMS) amplitude (to accentuate patterns and give indication of sediment type), dip of maximum similarity (to enhance discontinuities and edges), relative acoustic impedance (to identify changes in physical properties), shaded relief (to display digital elevation data and allow a 3D visualisation), and frequency decomposition (to highlight thickness variation within stratigraphic units and lateral continuity), were all calculated. Following the approach proposed by Posamentier et al. (2007), multiple 3D seismic interpretation techniques (e.g., surface flattening, slicing, 3D viewer and combined maps) were used and attribute extractions were performed either within constant windows, defined by the relative picked reflections or, flat windows between a specific start and end time interval. Spectral decomposition attribute analysis was carried using OpenDtect software. The frequency range within the seismic volume was screened using the 'evaluate attribute' tool on several time slices and key surfaces and enabled the quick identification of optimal frequency bands to enhance internal structure. Red-green-blue (RGB) colour blended maps were then used to delineate seismic geomorphology supporting the interpretation of depositional processes (Posamentier et al., 2007).
This integrated approach aided the detailed mapping of individual features and discontinuities observed within the 3D volume (Brown, 1996;Chopra and Marfurt, 2005). The morphological features were mapped in Petrosys PRO and organised within individual geospatial layers and were subsequently used to perform geostatistical analysis. Linear features were exported as text files and segment coordinates were used to calculate length and bearing in MS Excel. Polygon feature statistics and dimension analysis were undertaken predominantly in Petrosys PRO. Rose plots were generated using GeoRose software.

Seafloor geomorphological characterisation
Based on a range of attribute characteristics (Fig. 3), the seafloor morphology was subdivided in three Seafloor Morphology Provinces (SMP): SMP1 and SMP2 over the basin, and SMP3 over the lower continental slope ( Fig. 4-A). The 'hollow' of Bulat and Long (2005) is observed in the south-west (Fig. 4). Downslope profiles show a depression up to 10 m below the surrounding seafloor. This represents an almost totally infilled relict erosional scarp, part of the Westray hollow complex, which has been interpreted to be formed by intense bottom-current activity following the Paleogene/Neogene compressional tectonic phases (Smallwood, 2004).
A low relief channel-like feature, 8.7 km in length and 625 m in width, is observed to the east of this depression. This feature is better imaged by the seafloor dip ( Fig. 3B) and RMS (Fig. 3C and D). No direct camera/sample data are available over SMP1, but immediately to the north-east, Stewart and Long (2012) indicate generally clay-gravel sediments at the seafloor.

Seafloor morphology province 2
SMP2 (Fig. 4A) covers approximately 342.5 km 2 and is characterised by a seafloor that is generally smooth with gradients less than 1 • . Sediment waves are identified to the north-west and west of the hollow (SMP1) with local seafloor gradients up to 3 • (Fig. 4B). Along the eastern edge of this province, seafloor gradients are locally up to 12 • , associated with a linear southeast/northwest oriented incision which terminates at the boundary between SMP2 and SMP1 (P1 Fig. 4). This incision, 18 m deeper than the surrounding seafloor, has a symmetrical profile in cross section and is up to 140 m in width (Fig. 4C). It also exhibits a higher seismic RMS amplitude response than the surrounding SMP2 area (Fig. 3C).
Photographs and video from 2011 and 2017 seafloor surveys (Table 1) imaged gravelly sandy-clay with coarse gravel and occasional pebbles, cobbles, and few sparse boulders/drop stones over the region. These photographs also reveal comet marks and smaller scale comet tails composed of coarse sandy sediments developed around single boulders (Fig. 4B).

Seafloor morphology province 3
SMP3 ( Fig. 4A) covers approximately 66.5 km 2 and is characterised by an undulating seafloor with elongated linear bedforms parallel to the continental slope with local dip up to 6 • (Fig. 3B). Stewart and Long (2012) indicate seafloor sediments comprising predominantly sandygravels in this area.

Seismic stratigraphy and geomorphology
Seismic stratigraphic units 1 and 2a and 2b are identified and correlate with the seafloor morphology provinces, where unit 1 underlies SMP1, unit 2a underlies SMP2, and unit 2b underlies SMP3 (Figs. 4D and 5). The study focuses on unit 1, while unit 2 is discussed to clarify spatial geometric relationships, seismic attributes, geological context, and substrate type.

Unit 1
Unit 1 is observed in the basin area extending approximately 20 km laterally and 12 km into the Faroe Shetland Channel. It is likely that it extends farther basinward where it thins below data resolution.
The top of this unit corresponds to the seafloor reflection while the bottom is marked by a high amplitude irregular reflection, marker 2 (Figs. 4 to 7). This represents an unconformity, and the strong phase reversal (opposite to the seafloor reflection) indicates a sharp decrease in acoustic impedance in the underlaying unit 2 ( Fig. 7D and E). In seismic sections, unit 1 exhibits a chaotic to transparent seismic facies with discontinuous low to medium amplitude internal reflections (Fig. 7). Only one additional reflection (marker 1) can be traced through the unit, which divides it into an upper and a lower part (Figs. 6 and 7). The upper part reaches a thickness of 32 m (ca. 40 ms TWT) in localised zones within the central section of the study area (Fig. 6A), while the lower part reaches a maximum thickness of about 75 m (ca. 95 ms TWT), towards the south-west edge where the unit infills a topographic depression (Fig. 6B). Seismic sections reveal that this depression corresponds to the large bathymetric hollow observed at the seafloor within SMP1, where it is wider and extends to the south and west beyond unit 1 and is defined here as a paleo-hollow (Figs. 6C and 7C). The edges of the paleo-hollow and the undulating seafloor observed to the west are related to either unit 2 or, a thin (below seismic resolution) unit 1 The combined unit 1 thickness map shows a generally thinning downslope trend, with increased thickness over the central-east edge of the area (Fig. 6C).
Seismic attribute analysis of unit 1 reveals a more complex basal morphology than that observed at the seafloor (e.g. Fig. 8A versus 8C and 8D versus 8F) and highlights an articulated network of channels and lobe/fan features developed in two predominant directions, westsouthwest and north-northwest (Fig. 8C). The architecture of these features is best appreciated using dip of maximum similarity seismic attribute time slices (either flattened or un-flattened), calculated at regular time intervals within unit 1 (Fig. 8). These time slices show a mottled to blocky seismic texture, predominant in the lower part of unit 1 (Fig. 8G). RGB spectral frequency decomposition time, and surface slices, highlight internal details and connectivity of the channels, but also give an indication of unit thickness (Figs. 6D and 8D to G).
Morphological features have been mapped and subdivided into three categories (Fig. 9): 1) Type-1 linear features, 2) Type-2 linear features and 3) lobe-shaped features, defined as terminal fans (TF), located at the end of Type-1 features.   These time slices show a combined seismic response due to the vertical resolution resulting from overlaying information imaged from above and beneath the base of unit 1. The colour variation enhances the articulated basin channel network and terminal fans revealing texture and edge details and improving understanding of the spatial relationship and thickness variation in response to frequencies; this seismic attribute clearly distinguishes blocky textures as well as shallow diverging Type 2 incisions. A1) shows the blocky seismic character of unit 1 where it infills the paleo-hollow (TF-7); A2) shows a younger Type 1 channel (CH-6) clearly cutting an underlying terminal fan (TF-8) and shows terminal rims at the end of the Type 2 incisions; B1) shows the characteristic fan-like pattern of the north-westernmost TFs (Group 2). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) include ten channel-like incisions (blue lines on Figs. 9 and 12) ranging between 180 m and 810 m in width and between 1500 m and 9000 m in length. These show two predominant directions; towards the westsouthwest, Group 1 and north-northwest, Group 2 (Fig. 9). Group 1 channels are characterised by knick-points (CH-4 and CH-5) and significant changes in directions (CH-6) which are not seen in Group 2, but this might simply reflect limited data coverage over the slope. Channel CH-6 appears to resemble the shape of the shallow relief channel observed at seafloor within SMP 1 (Fig. 3).
The southernmost channels (CH-1 to CH-5) show a relatively lower RMS amplitude response and a smoother cross-section compared to the other channels (Fig. 10).
Channel CH-6, the most representative example of this type of feature, shows an overall downslope shallowing and a general basinward reduction in thickness of unit 1 except for the paleo-hollow area (Fig. 7). Cross-sections highlight the irregular surface at the base of channel CH-6 and its terminal fan TF-7, with paleo-topography varying vertically (< 25 m) due to the presence of linear incisions (Type 2) running along the thalweg (Fig. 7B). The channel abruptly terminates against the paleo-hollow ( Fig. 7C and E). The clear arcuate paleo-hollow edge on its east flank, observed on all seismic attribute maps, relates to the surface of basal unit 1 cutting directly into older Eocene sediments (Smallwood, 2004), rather than unit 2 (Fig. 10E).

Type 2 linear features -Basin linear scouring.
Type 2 features consist of a series of sub-parallel, linear incisions of low to medium relief developed within the Type 1 channels roughly following the direction of the main channel thalweg (Fig. 9). These incisions are absent, or below resolution, in the southernmost channels of Group 1, while downstream they become narrower and divergent forming a fan-like pattern (terminal fans, Fig. 9), which are described in section 4.2.1.3.
Type 2 features have an average length of 1355 m, ranging from 30 m to 7410 m, with the shortest features formed on the terminal fans and the longest ones found along the main channels. The width varies between approximately 15 m and 65 m, with the minimum detectable width limited by the 3D data resolution.
The orientation of the Type 2 incisions is predominantly westsouthwest / east-northeast (Fig. 9B) dipping downstream, consistent with the main channels (Type 1), with several features steered to the north-northwest (northern part of Group 2) and others to the westsouthwest (Group 1 and southern portion of Group 2).
Cross-sections were extracted at regular intervals of approximately 500 m along the main channel feature CH-6 and its terminal fan to show how the morphology varies downstream (see Fig. 7A for profile locations). Type 2 incisions generally have a symmetrical V-shaped crosssection, up to 12 m deeper than the adjacent surface and become more asymmetrical down the channel and shallower across the terminal fans. Within the terminal fans the RMS amplitude is generally lower than adjacent positive relief areas while it is higher within the main channels (Fig. 10).
An isolated linear incision up to 50 m wide and 10 m deep with a northwest/southeast orientation is observed on the southeast edge of Group 1 (Fig. 9). This feature is deeper than the other Type 2 incisions, it does not follow any Type 1 channel and is characterised by a low seismic amplitude response. It appears to be cut by at least three younger eastwest oriented Type 1 channels (CH1, CH2 and CH5, Fig. 11). This narrow incision can be traced up slope where it becomes shallower and visible at the present-day seafloor as part of the extensive gully system observed up slope from the study area (see section 5.2).

Distal terminal fan/lobes features.
At least 16 terminal fans/ lobes (TFs) were identified at the end of Type 1 channels, with TF-1 to TF-9 related to channel Group 1 and TF-10 to TF-16 related to channel Group 2 (Fig. 9). Some of these terminal fans/lobes are partially masked by Type 1 channels, while others may represent multiple systems as their internal morphology shows a variation in Type 2 lineation directions and relief which may indicate the presence of separate TFs (e.g., younger systems overlying older ones, section 5.2.1). However, due to limited vertical resolution, they have been mapped as single terminal fans/ lobes.
Mapped terminal fans/lobes have an aerial extent ranging between 1.2 km 2 to 8.5 km 2 , while the cross-slope width varies between 0.7 km and 2.3 km. The larger ones, such as TF-8 and TF-12, extend further downslope and exhibit the greatest number of internal Type 2 incisions. RGB colour blended frequency maps (Fig. 11) reveal the fan-like pattern of the north-westernmost TFs, which comprises numerous low-relief lineations diverging in various directions associated with systems further downstream (Fig. 11B and B1). Type 2 lineations generally end in a sub-circular depression, up to 6 m deep, delimited by a small rim (Fig. 9, orange dots, and Fig. 11A2).

Unit 2
Unit 2 is visible across the whole study area and is bounded at the top by the seafloor reflection (where unit 1 is missing), or marker 2 and, at its base, by marker 3. Unit 2 comprises unit 2a over the basin and unit 2b over the lower slope (Fig. 5). Unit 2a is partially eroded where the overlying unit 1 is present. Unit 2a consists of moderate amplitude welllayered reflections up to approximately 80 m (100 ms TWT) in thickness (Figs. 4 and 5) and unconformably overlies an older Eocene sequence (Smallwood, 2004;Stoker and Varming, 2011). Unit 2a infills the paleohollow and where it is not eroded at the top by the unit 1 surface, to the west. The well-layered acoustic character is locally variable, exhibiting basal downlap and onlap, lateral accretion, thinning and undulating subparallel reflections (Figs. 4 and 5). Over the lower slope area, unit 2b reaches up to approximately 55 m (70 ms TWT) and the bottom section is characterised by a high amplitude negative seismic ('hard') event, exhibiting in both strike and dip sections and an undulating morphology with regular spaced crests and troughs (Fig. 5). The geometry and orientation of these features form a sediment wave field, comparable to the one identified by Hohbein (2005), on high resolution seismic 3D data, approximately 126 km to the north-east (Fig. 1), within the West of Shetland Drift (Knutz and Cartwright, 2003), which is of early Pliocene age (Hohbein and Cartwright, 2006).

Regional framework integration
SMP1 and associated seismic unit 1 are referred here as the glacigenic basin floor fan system -unit 1 (GBFFS-unit 1) fed by a gully system above, detailed in Stewart and Long (2012). Together they represent the middle to lower continental slope and distal basin part of the Foula wedge TMF (Stoker et al., 2006;Stoker and Varming, 2011;Long, 2012, 2016). SMP2, SMP3, and their associated seismic units 2a (basin) and 2b (lower slope), are interpreted to represent the undifferentiated Plio-Pleistocene FSN-1 megasequence, including the West of Shetland Drift contouritic sequence (Hohbein and Cartwright, 2006;Stewart and Long, 2012;Batchelor et al., 2021;Fig. 2). Within this sequence, Hohbein (2005) identified a 1700 km 2 sediment wave field of Early Pleistocene age (Hohbein and Cartwright, 2006) extending along the continental slope. Considering the extent of the WSD-slope (Fig. 1), the comparable sediment wave level observed over the study area within unit 2b (refer to Fig. 5) likely represents the southern tail of this wave field. Depths below seabed of this level may vary due to the different sedimentation rates and Faroe-Shetland channel configuration.
The late glacial to post-glacial MacAulay sequence should overlie units 1 and 2 but is thin and below 3D seismic resolution (Stoker and Varming, 2011;Fig. 2).

Glacigenic basin floor fan system -Unit 1
The lower continental slope gullies and the basin floor fan system have been numbered for simplicity from northeast to southwest and GBFFS-unit 1 includes fans F8 and F9 (Fig. 12A). Fan numbers corresponds to gully numbers that appears to have fed them last. The seafloor incision observed over the eastern edge of the study area, terminating at the GBFFS-unit 1, represents the end section of downslope gully G-10 (Fig. 12A). Basin floor fans 1-7 were delineated based on the regional bathymetry alone, while 3D seismic attribute analysis (Fig. 12C   Fig. 12. Regional shaded bathymetry image derived from the AFEN (2000) dataset combined with bathymetry maps from Long et al. (2004), geo-referenced and composed in PetrosysPRO (see Fig. 2 for location) versus the basal seismic morphology showing the basin development of the downslope gully system (blue polygon marks limit of the 3D seismic extent). A) Large scale seafloor basin floor fan deposits and subdivision in relation to visible lower slope gullies, numbered from northeast to southwest. Fan numbers corresponds to gully numbers that appears to have fed them last. Fans F8 and F9 can be further subdivided based on 3D seismic attribute analysis. Relationship with the underlying seismic geomorphology is shown in insets C) and D). B) Unit 1 basal surface seismic geomorphology interpretation (refer to Fig. 11) with the overlying bathymetry showing the relationship with the seafloor morphology. Dashed blue arrows represent flow lines, associated with the main Type 1 channels originating from gully G-8, while dashed purple arrows show those originating from gully G-9; C) Frequency spectral decomposition calculated at the seafloor overlying shaded relief bathymetry. The dashed yellow polygon marks the area interpreted to represent the latest stage of the unit 1 deposition. D) Dip of maximum similarity seismic attribute at unit 1 basal surface overlaying the shaded bathymetry. The large scale basin fan interpretation also shown. Light pink and light red shaded polygons indicate areas where deposits generated from G-9 and G-8 respectively are dominant. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) and D) allowed differentiation of basin fan F9 from F8.
The fan boundaries show a general lateral overlapping trend to the southwest with F1, developed from gully G-1, and partially covered by younger systems emanating from gullies G-2, G-3, G-4, and G-6 (fans shaded in orange in Fig. 12A). Gully G-6 bifurcates at the slope-break and feeds two basin fans, F6' and F6", with the former being the older (Fig. 12A). The oldest basin fans were likely fed from gullies G-5, G-7 and G-10, now masked by the relatively younger fans (F4, F6 and F9), while the youngest basin fan (F8) is developed from gully G-8 (Fig. 12A).

Correlation between seafloor and unit 1 basal morphology: Sedimentation patterns
The bathymetry shows poorly defined shallow channels cutting the crests of the basin floor fans (Figs. 8A, D and 12C), including GBFFS-unit 1, giving the appearance of levée-style construction (Davison, 2004;Stewart and Long, 2012). Results from this study reveal for the first time the basal geomorphology and internal architecture of the GBFFS-unit 1 which, combined with frequency spectral decomposition calculated at the seafloor, show a far more detailed and complex structure than previously described (Figs. 8 and 12). This is characterised by an articulated channel network (Type 1), sublinear incisions (Type 2), and numerous terminal fans/lobes (TFs).
Type 1 channels, CH1 to CH5, observed within GBFFS-unit 1 relate to the basin fan F9 developed from gully G-9 and crosscut the basin extension of gully G-10; while Type 1 channels CH6 to CH10 appear to develop from gully G-8 (Fig. 12B) and form basin fan F8.
Seismic sections across the paleo-hollow show a multi-phase infill by contouritic sediments of unit 2a (Smallwood, 2004). Subsequently channel CH-6, the youngest feeder channel, eroded and partially infilled the paleo-hollow (Fig. 10). The abrupt change in paleo-topography (Fig. 7) likely confined and accelerated the flow within the paleohollow, favouring basal erosion of the contouritic deposits (unit 2a) which was subsequently followed by accumulation and infilling. The blocky seismic texture observed on time slices (Fig. 11A1) could be associated with the presence of mud clasts, generated during the basal erosion (Haughton et al., 2003;Talling, 2013;Peakall et al., 2020) or due to the heterogeneous nature of the material transported by the gravity flow (CH-6) i.e. iceberg rafted debris, including boulders/dropstones as observed at present-day seafloor (Fig. 4A).
The degree of basal erosion is also strongly related to the sediment composition of the substrate (Posamentier and Kolla, 2003;Peakall et al., 2020), which is reflected on the amplitude seismic response within the GBFFS-unit 1 (Fig. 10A). A high RMS response is interpreted to be associated with coarser infill or the presence of the more competent substrate of the underlaying Eocene sequence (Smallwood, 2004;Stoker and Varming, 2011), east of the paleo-hollow (Fig. 10E). The basal surface marks a decrease in relative acoustic impedance (Fig. 7E), suggesting a transition to a more coherent and cohesive clay/fine grained unit underneath, which is in agreement with the contouritic nature of unit 2a (Stoker and Varming, 2011). Alternatively and in areas where it LGM; the last glacially derived morphology is preserved at present-day seafloor (1 and 2) due to the low post glacial sedimentation rate; earlier phases are recorded within the GBFFS-unit 1 and at its basal surface (2a, 2b) and are interpreted based on internal architecture, seismic geomorphology, and preservation of the gully system (1). Blue circled numbers identify the same features on different views. 1) Downslope gully system -Top right: 3D view of the regional seafloor morphology and 3D seismic reflection subsurface section, top left: 3D view of the regional seafloor gradient. Bathymetry is derived from the 3D seismic, AFEN (2000) and EMODnet (2020) dataset combined with bathymetry maps images from Stewart and Long, 2012;Stewart and Long, 2016; 2) Basin fan system including GBFFS-unit 1; 2a) and 2b) seismic geomorphology details within GBFFS-unit 1: examples of frequency spectrum decomposition showing complex basin fan architecture with distributary channels/terminal fans. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) cuts directly the Eocene sequence, it may indicate over compaction at the basal surface within the overlying GBFFS-unit 1, as seen on masstransport deposits (Shipp et al., 2004;Maselli and Kneller, 2018), and/or presence of coarser grained basal lag deposits as sampled at the base of the Pleistocene section at borehole 99/3, located ~51 km south west over the Rona wedge debris flow lobes (Davison, 2004; BGS log - Table 1 and Fig. 2 for location). The latter explanation would account for the imprint of type 2 incisions which resemble groove-like features widely observed at the base of sand/sandstone layers in deep-water successions (Dzułynski and Sanders, 1962;Ricci Lucchi, 1969;Smith, 1972).
The basal reflection is therefore interpreted as a diachronous surface which shows distinct and well-defined acoustic patterns related to multiple phases of cut and fill. The apparent lack of well-defined basal infill deposits (either fine-or coarse-grained sediments) directly above this prominent reflection is probably due to the limited vertical resolution, and consequent tuning effects, of the 3D seismic dataset (Fig. 10E). The absence or reduced thickness of fine-grained basal deposits could be also attributed to winnowing of finer sediments by bottom currents enhanced during interglacial or warmer periods (Stoker, 1995;Weaver et al., 2000;Knutz and Cartwright, 2003;Knutz and Cartwright, 2004;Long et al., 2004;Davison, 2004;Hohbein and Cartwright, 2006;Long et al., 2011). In the Faroe Shetland Channel bottom currents have flowed south westerly since the Early Pliocene (Hohbein and Cartwright, 2006) and may have also promoted the southwest deflection observed in some of the Type 1 channels (green arrows in Fig. 8 C).

Glacigenic basin floor fan system -Unit 1 internal architecture
The Faroe Shetland Channel basin floor fans are part of a system which has migrated at both local and regional scales over time. Within the study area, the two main fans F8 and F9, forming GBFFS-unit 1 (Fig. 12A), developed from two (G-8 and G-9) of the ten downslope gullies observed over the lower continental slope (Fig. 12B). Gully G-8 likely represents the last active conduit feeding fan F8 (Fig. 12C and D). In particular the area of F8 is marked by a distinct frequency spectrum decomposition response at the seabed (Fig. 12C), which indicates an increased thickness compared to the adjacent areas as part of the upper part of GBFFS-unit 1.
Cross-cutting (Type 1) channels indicate that deposits derived from the westernmost gully G-10 are masked by the imprint of the more recent events associated with gullies G-8 and G-9.
The architecture of the GBFFS-unit 1 is reminiscent of a turbidite channel-lobe system (Fig. 13). For example, the terminal fans/lobes geometry is similar in appearance to the distal fringes of the Mississippi fan hybrid dendritic channel system (Talling et al., 2010(Talling et al., , 2012 and in the northeast they show multiple generations of diverging incisions (Type 2) with at least four back-stepping phases (Figs. 11B and 13). This retrograding stacking pattern is common to distal turbidite basin lobes like the Niger delta submarine fan (Zhang et al., 2016). The Type 2 incisions, however, also resemble the feather-like pattern identified within the North Sea Fan glacigenic debris flow deposits (Gafeira, 2009), recently re-interpreted as channel-levee deposits within a glacigenic turbidite system (Bellwald et al., 2020). Comparable linear erosional features interpreted as divergent internal furrows (here Type 2) and a series of stacked frond-likes lobes (here TFs) are observed on the mid-Norwegian paleo-slope within the mid/late Pleistocene glacigenic sequence of the Naust formation (Montelli et al., 2018). Montelli et al. (2018) interpreted the frond-shaped features as debris flow deposits developed from turbidity currents and proposed that the variation in shape and location of these deposits is related to sediment mechanical properties and local changes in bathymetry.
The GBFFS-unit 1 also has features diagnostic of other gravity-flow deposits such as a mass transport complex (Bull et al., 2009;Gafeira, 2009) or low sinuosity deep-water channelised debris flows (Posamentier and Kolla, 2003). This poses the question of what were the sediment delivery processes responsible for such complex architecture.
The narrow incisions (Type 2) observed within the main distributary channels (Type 1) are interpreted as groove scours generated by individual large clasts being transported in high energy clast-rich gravity flows (Dzułynski and Sanders, 1962;Ricci Lucchi, 1969;Allen, 1984). This appears to have been the dominant process active at that point in the development of the deep-water GBFFS-unit 1. As the linear incisions diverge downstream forming terminal fans/lobes (Fig. 11) and become progressively less incised, it is possible to infer a basinward transition towards partially unconfined flows, as proposed for diverging grooves observed at the base of debris flow-deposits offshore Indonesia (Posamentier and Kolla, 2003).

Downslope sedimentation processes and deep-water deposits
Basin deposits, like the GBFFS-unit 1, are rarely investigated and are generally defined as stacked debris flow fan/lobes or channel-levee deposits. They exhibit a chaotic seismic character on sub-bottom profiler data and a hummocky top surface (Nygård et al., 2007;Stoker and Varming, 2011;Long, 2012, 2016;Fransner et al., 2018). Using high quality 3D seismic data, similar to other studies in the Barents Sea (Laberg et al., 2010) and Norwegian regions (Gafeira, 2009;Montelli et al., 2018), this study has shown that these deposits may have an internal structure which can be used to infer the depositional processes.
A multi-generation debris floor fan system was discussed by Long (2012, 2016), and this study shows that at least two of the downslope Gullies, G-8 and G-9, diverged into large cross-cutting distributary channels (Type 1) observed at the base of GBFFS-unit 1 and were previously unidentified ( Fig. 12B and D). Each gully has acted as the primary downslope conduit multiple times before the bulk of the sediment delivery switched to another gully. The GBFFS-unit 1 basal surface, revealed here in unprecedented detail, is interpreted as the result of the erosive action by multiple gravity-driven flows, not necessarily debris flows or turbidites, and represents a diachronous surface resulting from a combination of several events. Diamictic glacigenic debris-flow deposits (like those encountered at borehole 85/1 [BGS, 1991, borehole log BGS 1985) which formed the upper slope of the Foula wedge TMF were initially incised by gravity-driven flows. Downslope, these flows remained erosive due to the relatively steeper lower slope and less resistant substrate, consisting of cohesive fine grained contouritic sediments of the West of Shetland Drift -slope deposits (Figs. 2 and 13;Hohbein and Cartwright, 2006;Stewart and Long, 2012). This interpretation is supported by the downslope gully morphology: low sinuosity, with U-shaped low angle flanks in the upper/mid slope, becoming narrower and steeper, V-shaped, in the lower slope (Stewart and Long, 2012). A number of gullies converge downslope rather than branching out (G4 and G6). This is opposite to observations on other high latitude TMFs where gullies usually disappear/branch out in the mid slope (Gales et al., 2013a(Gales et al., , 2013b or transform into large, low energy turbiditic channels distally (Noormets et al., 2009;Pedrosa et al., 2011;Gales et al., 2013b;Ó Cofaigh et al., 2018;Fransner et al., 2018;Newton et al., 2021). It also contrasts the classic models of debris flows decelerating downslope forming lobes and transforming into lower density turbidity currents across the basin. The substrate lithology and the relatively high gradient of the lower slope have controlled development of the gully system (Fig. 13), the dynamics of the subsequent gravity-driven flows and the formation of the whole basin floor fan system. As these confined high energy density flows reached the basin, they still had the capacity for widespread channelised erosion (Type 1), before eventually decelerating (within the terminal fan/lobes) and gradually transforming into more dilute flows.
The deposits forming the GBFFS-unit 1 show a high degree of vertical and lateral variability which is inferred to be similar across the rest of  Stoker, 1995, Davison, 2004Hohbein and Cartwright, 2006;Stoker and Varming, 2011;Stewart and Long, 2012;Ballantyne and Small, 2019;Bradwell et al., 2019). Blue circled letters identify the same features on section and map views. Stage 0: Pre LGM setting dominated by along-slope currents (contour-currents) with unknown glacigenic input; Stage 1: classic GDF dominated TMF building over the upper slope; Stage 2: significant influence of meltwater pulse discharges soon after the LGM, associated with frequent ice margin fluctuations and the unusual dynamics of this margin alternating to increase/decrease the GDF and contourcurrent activity. This results in a complex downslope gully and basin fan system controlled by paleo-morphology and substrate composition; Stage 3: The last ice recession is dominated by meltwater generated sediment-rich gravity flows and turbid-surface plumes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the basin fan system (F1 to F6, Fig. 12A). This variability is a function of the gravity flows also reworking sediments deposited by ice-rafting, bottom currents and late suspension settling. This interpretation agrees with the lithofacies analysis from well BH 99/3 on lower slope lobes of the neighbouring Rona wedge (Davison, 2004, BGS, 1999. In summary, GBFFS-unit 1 is interpreted to consist of: 1) relatively older lag deposits/debrites generated by either a) multiple independent long runout gravity-driven flows; b) multiple debritic pulses similar to linked debrites in fluvial systems (Haughton et al., 2003(Haughton et al., , 2009Talling et al., 2012;Peakall et al., 2020); 2) turbiditic channelised/lobe deposits derived from longitudinal gravity flow transformation, as described in many classic Pleistocene TMFs (e.g. Pedrosa et al., 2011;Laberg et al., 2018;O', Cofaigh et al., 2013 and the ancient Bolla Bollana TMF in South Australia (Le Heron et al., 2014); 3) the youngest deposits are from lower energy turbidites/plume suspension settling forming the upper part (F8 in Fig. 12C). While short offset reprocessing of the 3D dataset could improve resolution and imaging details, thus allowing further refinement of this interpretation, borehole/core samples will be required for ground truthing.

Evolution of the Foula Wedge basin floor fan system and implications for glaciated margin development
Existing geological information combined with the results from this study, on sediment delivery processes that influenced the development of the Foula wedge basin floor fan, can be summarised as follows (Fig. 13): a) an upper slope comprising stacked lobate glacigenic debris flows; b) a mid slope sediment bypass zone; c) a lower slope erosional zone where gravity-driven flows incise low sinuosity gullies up to 25 km in length; d) a transition to the basin floor with sufficient gradient to sustain high-density flows carrying clast-rich sediments to the fan system in the basin.
The next sections tentatively link these observations to the existing chronology, wider paleoglacial reconstructions and regional geological framework (Figs. 13 and 14).

Stage 0 -Pre last glacial maximum
Existing reconstructions indicate that extensive shelf-wide glaciation by the British Irish ice sheet may have started at near the beginning of the Pleistocene (Rea et al., 2018). Glacially derived sediments formed prograding wedges 50-200 m thick repeatedly during subsequent glacial cycles (Stoker, 1995;Hohbein and Cartwright, 2006;Stoker and Varming, 2011).
The pre-LGM slope and basin sequence is interpreted to consist of predominately contouritic and glaciomarine sediments of the Plio / Pleistocene Morrison 1 and lower Morrison 2 in the mid/upper slope (Stoker and Varming, 2011) and the West of Shetland Drift in the lower slope and basin (Hohbein and Cartwright, 2006;Stoker and Varming, 2011;Batchelor et al., 2021;Fig. 14).

Stage 1 -Last glacial maximum
This portion of the margin is traditionally interpreted as the confluence zone between the British-Irish and the Fennoscandian ice sheets which reached the shelf edge approximately 26-25 k yr BP during the LGM (Bendixen et al., 2018;Clark et al., 2018;Rea et al., 2018;Hall et al., 2019;Ballantyne and Small, 2019;Bradwell et al., 2019). At this time, the upper to mid slope was probably dominated by the emplacement of massive diamicton deposits from glacigenic debris flows forming classic lobate structures over the upper to mid slope of the Foula wedge (Stoker, 1995;Stoker and Varming, 2011;Long, 2012, 2016;Bradwell et al., 2019). This constructional phase of shelf-margin development is interpreted as a GDF-dominated TMF (Fig. 14).

Stage 2 -Post last glacial maximum
Following the shelf-wide LGM glaciation, ice mass reduction led to disconnection of the BIIS and FIS between 23 and 21 k yr BP. This promoted the development of the dynamic Shetland ice cap, which experienced numerous ice lobe advances and ice divide reconfigurations (cross-cutting flow directions), between 23 and 20 k yr BP (Bradwell et al. (2019). The most extensive advances that, neared the shelf break, produced GDFs on the upper continental slope which alternated with meltwater-driven sedimentation farther down the continental slope originating when the ice margin was more distal to the shelf break (Fig. 14). Subsequently, the continental slope became dominated by erosional rather than depositional processes, with high energy gravitydriven flows constrained within gullies enabling them to reach the foot of the slope and travel across the basin. Flows within the gullies comprised heterogenous sediments, derived from the upper slope diamicton (GDF), with material from iceberg rafting and suspension settling and possibly mud clasts derived from basal erosion of both lower slope and basin contouritic sequences (unit 2a and 2b).

Stage 3 -Final ice-sheet/cap recession
A distinct meltwater pulse at ~19 k yr BP has been proposed which ties in with a rapid retreat of the marine terminating margin of the Shetland ice cap (Sejrup et al., 2016;Clark et al., 2018;Hjelstuen et al., 2018;Ballantyne and Small, 2019;Bradwell et al., 2019). Due to the paleo-topography, this deglaciation phase was possibly accelerated at the Foula Bight, directly up slope from the study area (Davison, 2004;Stoker and Varming, 2011, Fig. 2). The overdeepened Papa Basin linked to the Foula Bight embayment, a zone of deeper water deflecting the shelf break approximately 10 km landward (Davison, 2004, Fig. 2), would have facilitated higher ice flux and a draw-down of the ice surface thus directing the sediment-rich subglacial meltwater discharge. Periodically these were of sufficient density to initiate the high energy gravity flows responsible for the development of the youngest downslope gullies, which, as indicated by Stewart and Long (2012), can be traced back into GDF lobes deposits of the upper slope Foula Wedge. These gravity flows together with iceberg rafting became the dominant sediment delivery process to the basin (Fig. 14 stage 3).
Ultimately, both slope and basin deposits were covered by the ~5-mthick McAulay sequence consisting of ice-rafted and finer sediments derived by suspension settling from meltwater plumes and hemipelagic rainout. These were continuously reworked by the action of strong bottom currents (Masson, 2001;Stoker and Varming, 2011;Stewart and Long, 2012) as part of the late glacial to post-glacial deposition.

Wider implications
According to the model presented here, at the LGM, the Foula Wedge is initially a classic GDF-dominated TMF system (Fig. 14, stage 1), building out the glacier fed sediment depocenter beyond the shelf edge (Vorren et al., 1989, Vorren andLaberg, 1997, Type 1 TMF from Pope et al., 2018). Subsequently, it becomes influenced by meltwater driven density flows associated with frequent ice margin fluctuations (Fig. 14, stage 2), which entrain sediments as they move down the continental slope and favoured long run out gravity-driven flows, similar to the Type 3 TMF (Pope et al., 2018) including the Kveithola TMF, Barents Sea (Lucchi et al., 2013).
Stage 2 marks a transfer of sediment deposition to the basin with a non-depositional/sediment by-pass zone in the mid slope (Fig. 14, stages  2 and 3). This change in erosional-depositional style may be related to the density of flows discharged across the shelf-break, being a function of meltwater and sediment volumes which varied in response to the ice sheet/ice cap dynamics (Laberg et al., 2018;Pope et al., 2018). It may also reflect more localised conditions, for example an ice margin undergoing changes due to increased calving rates, fluctuating ice margins, meltwater discharge and ice divide migrations (Stoker, 1995;Bjarnadóttir et al., 2014;Dowdeswell et al., 2016) which would fit with the ice sheet/cap fluctuations in the region (Bradwell et al., 2019).
Seafloor morphological differences with the relatively larger adjacent Rona wedge system (Fig. 2), also suggest that the local paleo-slope morphology, substrate composition and the shelf configuration may have had a significant impact on subglacial water-flow paths (Clark et al., 2012;Dowdeswell et al., 2016;Stewart and Long, 2016), and consequently on the development of the gully and basin floor fan systems.
During Stage 3, when the ice margin is more distal to the shelf break, the Foula bight embayment focuses ice and meltwater flux allowing the gravity flows, which can still erode the slope and reach the basin floor, to become the dominant sediment delivery process. The presence of cohesive contouritic deposits over the lower slope, combined with the slightly higher slope gradient, has favoured the erosion of the narrow gullies and extension of the Type 1 channels down to the basin. Based on the features identified within GBFFS-unit 1, it appears that each gully may have acted as the main downslope pathway multiple times. This doesn't exclude the possibility that a gully could have continued to be an open conduit for lower energy flows while new gullies were initiated by higher energy flows elsewhere. This work demonstrates the heterogeneity and dynamic nature of TMFs and highlights that local factors, such as substrate composition, paleo-slope morphology, and bottom current strength, can significantly influence the behaviour of gravity flows and their deposits.
Evidence from the basin floor fans suggests a transition in sediment delivery processes controlled by meltwater-driven sedimentation, which impact the large scale and long-term evolution of TMFs. Borehole data (85/01 logs in Stoker, 1995) also indicated that gravity flow deposits are separated by units associated to along slope sedimentation, reflecting reduced glacial sediment supply to the slope and interpreted here to relate to dynamics of the ice margin, for example retreat or ice divide migration leading to reduced ice flux/iceberg calving or redirected subglacial meltwater discharge. The important role of meltwater inputs identified here agrees with others TMF studies in both mid and high latitudes, for example on the Laurentian fan (Piper et al., 2007), the Disko fan (Ó, , and the North Sea fan (Bellwald et al., 2020).

Conclusions
1) 3D seismic reflection data have allowed an improved understanding of a glacigenic basin floor fan system (GBFFS-unit 1) that, together with the associated downslope gully system, forms part of the distal Foula wedge TMF, offshore West of Shetland. Seismic geomorphology reveals a subsurface distributary channel network and terminal fans resembling a deep-water turbidite lobe/channel-levee system. However, the presence of features, such as groove scours, and a distinctive blocky seismic texture, suggests that long runout debris flow/cohesive gravity flow processes operated over the basin area of the fans. This indicates a more complex depositional system than previously thought. 2) Integration with regional information shows that this portion of the basin fan system was fed by two of the main mid/lower-slope gullies. Both the basin fans and gullies are part of the Foula wedge TMF. The period, soon after the LGM was characterised by numerous ice margin fluctuations including rapid retreats and re-advances. It is suggested that gravity driven flows were initiated by sediment charged meltwater, which varied as function of the dynamics of the combined British-Irish and Fennoscandian ice sheets and the subsequent independent Shetland ice cap, along with the local shelf configuration and paleo-slope gradient and composition. 3) It is proposed that sedimentation was dominated by GDFs during the LGM. Subsequently, frequent meltwater discharge pulses generated gravity driven flows down the slope, which eroded and transported sediment. Iceberg rafting and variations in bottom currents added to the sedimentary processes dominating the Foula wedge TMF evolution soon after the LGM. The relatively high gradient of the lower slope and presence of pre-LGM contouritic deposits have facilitated clast-rich gravity flows to remain erosive and confined within deeply incised gullies to the bottom of the slope and across the basin. 4) The GBFFS-unit 1 clearly demonstrates the heterogeneity of basin floor fan deposits and provides new insights on the deep-water sedimentation processes and patterns operating on the paleo-slope/ basin environment associated with a mid-latitude TMF.

Data availability
The 3D seismic reflection data used in this study have been provided by Siccar Point Energy under a non-disclosure agreement between Siccar Point Energy and the University of Aberdeen. This dataset is not publicly accessible, so please contact the Authors to obtain information on how to access the data. Seafloor photograph presented in Fig. 4B is courtesy of Siccar Point Energy as part of the survey data package provided by Fugro (2017). AFEN (2000) bathymetry data for this study were obtained via Fugro GB (North) Marine Limited in 2019. EMODnet Bathymetry is Open Source at: http://www.emodnet-bathymetry.eu. Boreholes and sub-bottom profiles from the British Geological Survey (BGS) are available under Open Government Licence at GeoIndex Offshore | BGS, Interactive map at http://mapapps2.bgs.ac. uk/geoindex_offshore/home.html?_ga =2.227456346.671481726.1639228033-183223392.1589786690#.
Regional topography presented in Fig. 1 is provided by EMODnet Bathymetry Consortium 2020. Bathymetric data presented in Figs. 2A, 12, 13 and 14 are derived from a combination of multiple sources, which have been geo-referenced and composed in PetrosysPRO: 1) the seafloor horizon extracted from the 3D seismic reflection data provided by Siccar Point Energy, 2) AFEN (2000) dataset, 3) EMODnet (2020) dataset, 4) Bathymetry map images extracted from Long et al. (2004) and Stewart and Long (2016).
For all data listed above, excluding the 3D seismic dataset, refer to Table 1 for further details on data type and source links.

Declaration of Competing Interest
All the authors declare no conflict of interest.