A Late Pleistocene channelized subglacial meltwater system on the Atlantic continental shelf south of Ireland

The study of palaeo‐glacial landforms and sediments can give insights into the nature and dynamics of ice sheets. This is particularly the case with regards to the subglacial record, which is challenging to observe in contemporary glaciated settings and hence remains only partially understood. The subglacial hydrological system is an essential component of ice dynamics, where increased water pressure enhances ice motion and sediment deformation, thus reducing ice‐bed contact. Tunnel valleys are large, sinuous, steep‐sided incisions that, together with smaller scale meltwater channels, indicate subglacial meltwater discharge beneath large ice sheets. Through the use of high‐resolution marine geophysical data, a system of buried and exposed tunnel valleys, possible subglacial or proglacial meltwater channels and palaeo‐fluvial valleys have been identified across the shelf of the Celtic Sea between Ireland and Britain. The presence of steep‐sided and overdeepened tunnel valleys is indicative of a large channelized meltwater drainage system beneath the former Irish Sea Ice Stream, the most extensive ice stream to drain the last British–Irish Ice Sheet. After the rapid ice expansion across the Celtic Sea shelf around 28–26 ka, the tunnel valleys were carved into both bedrock and glacigenic sediments and are associated with rapid ice stream retreat northwards into the Irish Sea Basin between 25.6 and 24.3 ka. The presence of a major subglacial meltwater system on the relatively shallow shelf suggests that significant erosive meltwater discharge occurred during the last deglaciation and highlights the important contribution of meltwater to the retreat of the British–Irish Ice Sheet on the continental shelf.

During the last glaciation, the Irish and Celtic Seas ( Fig. 1) were occupied by the British-Irish Ice Sheet (BIIS), the westernmost ice mass of the Eurasian Ice Sheet (Hughes et al. 2016). With an extensive marineterminating margin, the BIIS attained a shelf-edge position along the western Atlantic margin around 27-25 ka and the Irish Sea Ice Stream (ISIS) was its largest outlet (Knutz et al. 2001;Scourse et al. 2009Scourse et al. , 2019Chiverrell et al. 2013;Smedley et al. 2017b;Lockhart 2019;Fig. 1). Following advance to its maximum extent at 25.6AE0.5 ka, the ISIS retreated from the shelf edge of the Celtic Sea back into the Irish Sea Basin (ISB; Fig. 1) by 24.3AE0.2 ka (Small et al. 2018). Here, the pace of retreat decelerated due to topographic controls, a fall in relative sea level and reduced palaeotidal amplitude (Ward et al. 2016;Scourse et al. 2019;Coughlan et al. 2020;Van Landeghem & Chiverrell 2020). Although the ISIS occupied a large area (likely over 50 000 km 2 ), its geomorphological record on the sea floor remains poorly defined. This, in part, is due to the presence of a series of mega-ridges dominating the present-day sea floor of the Celtic Sea shelf. These may relate, at least in part, to postglacial mega-tidal processes (Scourse et al. 2009;Lockhart et al. 2018) and hence overprint the glacial landform record.
Elsewhere in the Irish and Celtic Seas, several large, overdeepened valleys have previously been identified and interpreted as tunnel valleys. Wingfield (1989) undertook an early geomorphological characterization of these channels and classified them into three generations spanning up to four glacial phases (Late Elsterian, Saalian and Weichselian glaciations; <500 000 years ago). This work attracted a critical response, especially with regard to adjacent incisions mapped in the English Channel, where sedimentological and other evidence demonstrated a fluvial origin of the channels (Hamblin et al. 1990). More recent studies have demonstrated that the overdeepened Irish Sea valleys were formed as timetransgressive erosional features during the last glacial period (Callaway et al. 2011;Coughlan et al. 2020;Van Landeghem & Chiverrell 2020).
Tunnel valleys are large, overdeepened channels of variable dimensionslengths greater than 100 km, widths of up to 10 km and depths of 500 mincised into bedrock or soft sediments (e.g. patterns, can have overhanging tributaries, start and terminate abruptly, and are often associated with subglacial meltwater landforms likely formed close to the ice margin, such as eskers ( O Cofaigh 1996;Stackebrandt 2009;van der Vegt et al. 2012;Huuse & Kristensen 2016;Montelli et al. 2020;Ottesen et al. 2020). While there are several theories about the formation of tunnel valleys, subglacial meltwater is considered the primary erosive agent in their evolution, either as a single or multiple catastrophic events or under more steady-state conditions ( O Cofaigh 1996;Kehew et al. 2012;van der Vegt et al. 2012;Livingstone & Clark 2016;Lewington et al. 2020). They provide drainage routeways for large volumes of water and sediments at the base of an ice sheet ( O Cofaigh 1996;Greenwood et al. 2016;Livingstone & Clark 2016).
As the subglacial hydrological system is an essential component of ice-sheet drainage, it plays an important role in ice-sheet dynamics (Zwally et al. 2002;Price et al. 2008;Doyle et al. 2014). Subglacial water flow is intimately related to ice-sheet motion (e.g. enhancing glacier sliding) and responsible for subglacial sediment deformation. The investigation of former ice-sheet beds through their geomorphological record is, therefore, an important contribution to the understanding and modelling of the wider significance of subglacial hydrology in the Earth system (Price et al. 2008;Dowdeswell et al. 2016;Livingstone & Clark 2016).
Recent data from the Celtic and Irish Seas have contributed to an enhanced understanding of the regional reconstruction and ensuing modelling of the former BIIS dynamics (Smedley et al. 2017a, b;Chiverrell et al. 2018Chiverrell et al. , 2020Lockhart et al. 2018;Small et al. 2018;Scourse et al. 2019Scourse et al. , 2021Van Landeghem & Chiverrell 2020). A recent synthesis of geochronological data (Scourse et al. 2021) demonstrated that the pace of deglaciation within the basin was mainly influenced by internal controls such as topography rather than external forcing.
The aim of this paper is to investigate the formation mechanisms and evolution of a large network of subglacial tunnel valleys across the floor of the Celtic Sea through the spatial integration of new seabed and sub-seabed acoustic data. First, we describe and interpret the geomorphology of the submarine landforms in the context of the regional shallow-acoustic stratigraphy. We then discuss the wider implications of this geomorphology for the dynamics of the former BIIS. The presence of subglacial valleys on the Celtic Sea shelf implies abundant meltwater production and erosion underneath an extensive area (~25 000 km 2 ) of the former ice sheet, perhaps influencing the last retreat of ice from the shelf.
Regional stratigraphy of the Celtic Sea The regional Quaternary stratigraphy of the Celtic Sea shelf above the acoustic basement is characterized by two vertically stacked units: the glacigenic Upper Little Sole Formation (ULSFm) and the glacifluvial or tidal Melville Formation (MFm). These units are mainly located in the mid-and outer shelf and are draped by modern sediment. Together, the ULSFm, MFm and modern deposits record the transition from glacial to glaciomarine and postglacial sedimentation (Pantin & Evans 1984;Scourse et al. 1990Scourse et al. , 2018Scourse et al. , 2019Praeg et al. 2015;Ward et al. 2016;Lockhart et al. 2018). On the basis of a single foraminifera assemblage collected from a core on the outer shelf, Pantin & Evans (1984) previously interpreted the ULSFm as Late Pliocene/Early Pleistocene. Whilst the ULSFm unit is mainly exposed on the lowermost section of the mega-ridges, the topmost seismic unit of the ridges has been correlated with the MFm (Pantin & Evans 1984). Unfortunately, no direct geochronological data exist for the MFm.
New geophysical and sediment core data across the Irish-UK Celtic Sea shelf redefined the regional stratigraphical units both spatially and temporally (Lockhart et al. 2018), based on their acoustic character, stratigraphical organization of glacigenic sediments and their association with the mega-ridges on the mid-and outer shelf ( Fig. 1; ibid.). The revised Celtic Sea stratigraphy (Lockhart et al. 2018) suggests that the glacigenic sediments, attributed as the ULSFm, date to between 27-24.3 ka (Praeg et al. 2015;Scourse et al. 2019). These often stiff and deformed glacigenic deposits of the ULSFm record the oscillatory nature of the ice sheet during its Late Pleistocene retreat across the shelf (Scourse et al. 2021). The glacigenic sediments were later partially eroded in the course of the ensuing marine transgression, during which the overlying sand ridges formed on the shelf. The sandy core of the largest ridges is considered to be composed of the MFm, which has an inferred age of 24-13 ka based on dates from underlying (ULSFm) and overlying (modern sediment) units. The uppermost depositional unit is discontinuous and of variable thickness and indicates marine sediment deposition during a continued rising postglacial sea level (Pantin & Evans 1984;Scourse et al. 2009). Although these postglacial marine sediments have been chronologically constrained to between 13 and 4 ka on the inner shelf, they lack dating on the mid-to outer shelf, where they are presumed to become older towards the shelf break (Furze et al. 2014;Lockhart et al. 2018). This postglacial-uppermost unit, overlaying ULSFm and MFm, corresponds to the seismic facies SF1 identified by Lockhart et al. (2018).
As a result of this new interpretation of the regional Celtic Sea stratigraphy (Lockhart et al. 2018), the revised age of 27-24.3 ka for the ULSFm refutes the previous interpretation of a Late Pliocene/Early Pleistocene age for this unit (Pantin & Evans 1984). In turn, the age of the Melville Formation is now hypothesized to be between 24-14 ka, with the uppermost deposits younger than 13.9 ka (Furze et al. 2014;Lockhart et al. 2018). As the data used in this paper cross some of the seismic lines presented by Lockhart et al. (2018) and Lockhart (2019), we follow their revised stratigraphy to corroborate the interpretation for the sedimentary units foundwithin our data set.

Material and methods
This paper uses low-and high-resolution bathymetry, backscatter and sub-bottom profiler data. Regional multibeam echo-sounder (MBES) bathymetric datawere compiled from the EMODnet portal at 115-m spatial resolution (EMODnet Bathymetry Consortium 2018). New, high-resolution MBES and shallow seismic data were collected by the Integrated Mapping for the Sustainable Development of Ireland's Marine Resource (INFOMAR) programme between 2016 and 2019, onboard the Marine Institute of Ireland vessels the RV 'Celtic Explorer' and RV 'Celtic Voyager'. The surveys were conducted using a Kongsberg EM2040 system following hydrographic standard practice IHO S-44 1a (IHO 2020). Bathymetric and backscatter data were delivered at 5-m resolution after processing in CARIS Hips & Sips and QPS FMGT Geocoder, respectively. Shallow seismic data (3.5 kHz) were acquired onboard the RV 'Celtic Explorer' using an IXBlue Echo 3000 (2017 and 2018 data sets) and Knudsen 3260 chirp source, and onboard the RV 'Celtic Voyager' using a Sonar Equipment Services pinger source. Additional seismic data were collected during the BRITICE-CHRONO cruise JC106 of the RRS 'James Cook' in 2014 and CV14007 cruise of the RV 'Celtic Voyager' using a Kongsberg SBP-120 chirp system (2.5-6.5 kHz) and a multi-tip Geo-Source 200-400 sparker, respectively. The JC106 data set (Lockhart et al. 2018;Lockhart 2019) was used in this work to extend the mapping of the features, in addition to the INFOMAR and CV14007 sub-bottom data (Fig. 1A).
Global Mapper and ArcGIS were used for geomorphological mapping of spatial data, whereas seismic profiles were visualized and interpreted using IHS Kingdom software. The stratigraphical correlations are based on Lockhart et al. (2018; location of their seismic lines is shown in Fig. 1A). Sediment velocity for calculation of depths below the sea floor was estimated to bẽ 1600 m s À1 after Lockhart et al. (2018). The geomorphological interpretation was carried out using a combination of terrain shading and first terrain derivatives (slope and aspect) and, for the purpose of this paper, was limited to channel-like features. Measurements of individual channel morphology (e.g. width, depth and length) were carried out systematically in ArcGIS and ArcGIS Pro, using the digitized sea floor and buried features. Locations and depths (xyz) of the latter were used to perform accurate georeferenced calculation within ArcGIS. Width (distance between digitized valley shoulders) and absolute depth (i.e. from base to valley shoulder) were computed using semiautomatic tools (i.e. Generate Transect Along Lines, Intersect and Stack Profile) within the GIS suite. The measurements were retrieved from perpendicular segments at 1-km distance along the reach of the channel (centreline), using the minimum and maximum values as measured in the bathymetric data at each transect. The depth values of the sea floor valleys have been calculated based on the MBES bathymetry; the buried valley depths have been measured from base to valley shoulder within the seismic profiles.
The channel lengths were measured on each channel centreline, representing the polygon main axis derived from the geomorphological analysis. Where two or more minor-branching channels (i.e. tributaries) coalesce, the longest section was used to determine the length.
The orientation measurements of the channels were carried out semi-automatically by dividing the centreline (the main axis of the valleys and ridges) in equal segments (100-m-long units) to capture variations in internal direction along the reach of the channels. The orientation measurements were performed separately for each group of interpreted landforms. Since the ISIS palaeo-ice flow is known across this region to follow a NE-SW trough axis (e.g. Small et al. 2018) these orientation measurements are presented as monodirectional rose diagrams. The orientation of the interpreted palaeo-fluvial valleys shows instead the predominant down-slope direction.

Results and interpretation
Analysis of bathymetric and sub-bottom data across the Celtic Sea shelf has revealed awealth of channels forming a large drainage network on the inner and mid-shelf. In the following sections, the characteristics and distribution of these channels are presented, and their interpretation discussed.

Sea floor drainage network of the inner and mid-shelf
Across the inner and mid-Celtic Sea shelf, around 80 channels are visible on the sea floor on both bathymetric and sub-bottom data. The exposed channels are predominantly SSW oriented, parallel to the southern Irish coastline (Figs 2, 3). The channels are between 2.4-76.1 km long, between 71-2566 m wide, with the tendency to widen towards their termini, and display straight to low sinuosity profiles (Figs 3A, D, S1A, S2A). Their depth, compared to the surrounding sea floor, is on average 3.4 m and the longer channels tend to be slightly deeper reaching up to a maximum of 29.3 m (Figs 3A, B, D, E, S1A). They start and terminate abruptly and show both U-shaped and V-shaped cross-profiles at their termini (Fig. 3A, B and related topographic profiles). They have overdeepened basins along the thalwegs, as shown by both bathymetric and seismic profiles (profiles X-X', Y-Y' and Z-Z', Figs 3B, D, 4B). Several tributary tunnel valleys form a branching/joining geometry pattern linked to the larger tunnel valley (Fig. 3A); whereas some tributaries appear to be overhanging valleys positioned above the main channel (Fig. 3B, D).
The gradient of the channel sides is generally steep (7°-9°) with a subtle terracing on the sides ( . The localized bathymetric lows are better viewed in the seismic data where the basement has been carved to form an irregular surface along the reach of the channel (Fig. 4). The thalweg has been often infilled by sediments, probably related to more recent sediment transport across the region (Figs 3A, D, 4). The backscatter data display low intensity values (≤À25 dB) within the channels and moderate to high intensity (≥À25 dB) in the surrounding areas (Fig. 3C, F). Localized zones of relatively high backscatter returns are also visible along the channel sides and terraces. The intensity of the backscatter data, coupled with the visual interpretation of the bathymetric data, suggest that high backscatter returns are most likely areas of harder substrate, including bedrock (McCullagh et al. 2020;Toole et al. 2020). Low intensity backscatter was groundtruthed as fine-grained muddy sediments (Doyle et al. 2019;Toole et al. 2020), which seem to partly fill the thalwegs of most of the channels, as is also visible in the seismic data (Fig. 4).
In the sub-bottom data, the base of the channels corresponds to a high amplitude reflector. This reflector forms the upper boundary of the bedrock/acoustic basement that at times also crops out at the sea floor (Fig. 4B). The sedimentary unit within the channels is discontinuous, as it infills depressions, is of low amplitude and generally around 3 m thick, but up to 6-7 m thick in places (Fig. 4A). Its character is consistent with that of the superficial drape or uppermost sediment deposit identified across the inner and outer Celtic Sea shelf. This deposit represents postglacial sediment deposition during rising sea level (Pantin & Evans 1984;Furze et al. 2014;unit SF1 of Lockhart et al. 2018). The presence of a sedimentary infill unit in the channel means that the automatically extracted depth of the exposed channels (as derived from the bathymetric data) is an underestimation of their true depth. Based on the seismic data, the total depth may be up to 3 m deeper (e.g. Fig. 4). In the central part of the Celtic Sea, to the east of the largest channel, several, discontinuous, narrow (on average~400 m wide), winding and comparatively small-scale ridges (~1-4 m elevation) were identified (Figs 5, S3). They are found as small individual segments (~1 km long) that, together with the longer ones (up to 10.5 km in length), are arranged in network. Adjacent ridges centrelines are mostly parallel, WSW-oriented and spaced between 0.8 and 2.5 km (Figs 2B,5). These ridges are also subparallel to the exposed and buried channels.
The exposed tunnel valleys are also associated with comparatively small-scale ridges, here interpreted as eskers based on their long and sinuous morphology and dimensions within their network (Figs 5, S3). The general geometry of eskers varies between singlestraight to sinuous-undulating morphologies that can be arranged as individual segments or in network (up to several km; Perkins et al. 2016;M€ akinen et al. 2017;Storrar et al. 2020). Although it is common to find fragmented eskers buried within tunnel valley infills (e.g. Boulton et al. 2001;Storrar et al. 2020), this has not been observed in our data. The general orientation of the eskers is subparallel to the tunnel valleys, which suggests a common genesis and/or that they were formed as part of the same flow regime. Ice-sheet reconstruction across Fig. 5. Small-scale, winding and discontinuous ridges arranged in network mainly WSW-oriented and found east of the northern tunnel valleys (general location in Fig. 2). Black arrows indicate some of the ridges and are oriented following the ridges' long-axes, as also shown by the rose diagram in Fig. 2B. BOREAS the shelf (e.g. Small et al. 2018) indicates a NE-SW aligned ISIS trough axis direction. Hence, the orientation of both tunnel valleys and eskers appears consistent with these ISIS models and parallel to the palaeo-ice flow. This may indicate a time-transgressive formation of the eskers during stages of ice-sheet retreat over a rigid bed (e.g. crystalline bedrock; Arnold 2019). A more rapid and instantaneous meltwater release would have resulted instead in a larger range of directions (Hewitt 2011), which is not the case here.

Buried channels on the mid-shelf
With increasing distance across the mid-shelf, sea floor channels are no longer visible on the shelf and instead are replaced by a series of buried channels, as imaged on the sub-bottom data (RV 'Celtic Explorer' cruise CE1701; Figs 1A, 2B, 6).
These buried channels are predominantly WSW-and SSW-oriented, and their widths and depths range between 165 and~2000 m and up to 20 m, respectively (Figs 2B, S1B). Their lengths are between 3-28 km based on the interpolated system lengths (L i ; Table S1).
The buried channels are cut into and infilled by a low to medium amplitude acoustic unit with subhorizontal reflections (light blue in Fig. 6), which in turn is draped by an uppermost, low amplitude acoustic unit (orange unit in Fig. 6). Given the regional seismic units recognized in the area (Lockhart et al. 2018), the channels are positionedwithin and towards the top of the Upper Little Sole Formation (dated at 27-24 ka; ibid.). The superficial drape (orange unit in Fig. 6) can be correlated to the postglacial marine sediment found both on the inner and outer shelf and formed during rising postglacial sea level (Pantin & Evans 1984;Furze et al. 2014;unit SF1 of Lockhart et al. 2018). This unit also infills the inner and mid-shelf tunnel valleys cut into the acoustic basement/ bedrock (Fig. 4), as mentioned in the previous section, and therefore we assume that it is the same uppermost sediment drape extending from the inner to mid-and outer shelf.
The buried channels show almost identical orientations, dimensions and geomorphological characteristics to their sea floor counterparts (Figs 2B, S1B, Table S1). In fact, the closely spaced seismic lines of CE1701 clearly show undulating thalwegs and abrupt starts and Fig. 6. Fence diagram of the buried tunnel valleys including seismic profiles from the CE1701 data set ( Fig. 2A). The continuity aswell as the abrupt start and terminations of the tunnel valleys is highlighted by the grey areas. The labelled tunnel valleys (i, j, k and l) correspond to those indicated by the arrows in CE1701 -080 shown on the bathymetry in Figs 1C, 2A. The stratigraphical interpretation is based on Lockhart et al. (2018) and the uncertainty in the spatial continuity of the seismic facies is indicated by question marks. The horizontal scale is the same for the entire data set. terminations (Fig. 6). Geographically they also appear to be the continuation of the sea floor channels based on their position and orientation (Figs 2B, 6). We therefore interpret these buried channels as tunnel valleys and as the continuation of the meltwater drainage system revealed at the sea floor. We hypothesize that, combined, both the buried and exposed channels were formed during the continuous headward erosion of the ice sheet during its final retreat (cf. Livingstone & Clark 2016).
Within our data set, the examination of the data further offshore, towards the outer shelf, has shown that similar geomorphological features as described above are no longer visible in the seismic data.

Channels close to Irish and UK coastlines
Closer to the Irish and UK coastlines (Fig. 2), on both the inner and mid-shelf, several channels were also mapped ( Fig. 7). On the Irish side, the channels to the west (Fig. 7A) differ from those to the east (Fig. 7B) based on their morphology and orientation. The westernmost channels (Fig. 7A) have a consistent direction parallel to the coastline, trending WSW (shelf-break directed; Fig. 2B) with minor/rare branching and very low sinuosity. The channels have gentle slopes, U-shaped profiles, are on average 600 m wide, have absolute depth averaged 2.5 m along their length (2.5-29 km; Fig. S1C, Table S1) and show an irregular topography along the thalweg (profile X-X'; Fig. 7A). They start and terminate abruptly, and cut perpendicularly through a sequence of arcuate ridges (white dotted areas in Fig. 7A), which have been interpreted as possible moraines or grounding-zone wedges based on their geomorphological similarity with landforms previously mapped around the western Atlantic Irish continental shelf (Benetti et al. 2010;Peters et al. 2015;Callard et al. 2018).
The Irish easternmost and the UK channels are on average 7.3 km long, they run perpendicular to the coastline, with a prominent SSE and NNW downstream profile, respectively, along the thalwegs and directed The absolute depth of these channels can reach up tõ 12 m but is on average 3.4 m (Fig. S1D). They are, on average~1 km wide and gently V-shaped along the reach of the channels (Fig. 7B, C). The dendritic drainage pattern together with size, orientation and sloping towards the Celtic Sea shelf, are all characteristics that show the similarity between the eastern Irish and UK inshore channels.
Given these morphological characteristics, two possible explanations arise for the formation of the two sets of channels described above. They could be: (i) meltwater conduits, either subglacial or proglacial, or (ii) palaeo-fluvial valleys. The abrupt starts and termini, localized lows and the absence of a constant downstream profile of the westernmost Irish shelf channels (Fig. 7A) are all inconsistent with the pattern of channels formed by subaerial water flow (Grau Galofre et al. 2017). These westernmost channels might therefore represent subglacial or proglacial meltwater conduits. In this case, we refer to meltwater channels as mesoscale features forming in a variety of environments (i.e. subglacial, proglacial, lateral; Glasser & Bennett 2004) (Fig. 7A). This is typical of meltwater conduits mainly driven by the ice-surface slope and topography beneath the ice. Thus, the resulting drainage pattern can contrast with the present-day topography (Bennett & Glasser 2009), as is the case here. Additionally, the orientation of the westernmost channels is subparallel to the tunnel valleys and eskers (Fig. 2B). As the latter represent the main glacial drainage pattern, these marginal meltwater features also follow the direction of the ice flow.
Although some of the meltwater channels can reach up to 28 km in length, the depths appear shallower and the widths narrower compared to those of the tunnel valleys (Fig. S1A, B, C, Table S1). Additionally, the westernmost channels appear to cross-cut the arcuate ridges at the sea floor (white dotted areas in Fig. 7A), suggestive of an event post-dating an ice stillstand. The channels could therefore represent former proglacial meltwater conduits, explaining the partial erosion of the palaeoglacial features preserved at the sea floor (Fig. 7A).
The easternmost set of channels offshore Ireland, on the other hand, and those offshore the UK might have a fluvial origin, based on their morphology and orienta-tion with respect to the coastlines: SSE for the former and NNW for the latter (Figs 2B, 7B, C). The channels might have formed in the postglacial environment prior to the Holocene marine transgression. Relative sea-level (RSL) model predictions from glacial isostatic adjustment (GIA) simulations for 20 ka within the study area (Brooks et al. 2008;Bradley et al. 2011;Scourse et al. 2018) suggest values between À80 and À90 m for both the SW of Ireland and SW of England. This is lower than the channel depths both to the east of inshore Ireland (~-70 m; Fig. 7B) and northwest of north Cornwall (Fig. 7C), indicating emergence of the shelf and potential for fluvial incision of the exposed surface.
These channels seem to be aligned, but not necessarily clearly connected, with onshore river systems of Cork and Waterford in Ireland across the inner shelf and aligned with rivers flowing from northwestern Cornwall (including, and north of, the river Camel; Fig. 1C) in the UK. Other N-S oriented channels were previously mapped offshore Co. Waterford (SE Ireland; Fig. 1C), where they were interpreted as palaeo-river channels potentially connected to the River Barrow-Suir drainage system during sea level lowstand (T oth et al. 2020). Therefore, it is possible that the valleys mapped in this paper were also used as preferential pathways by fluvial processes after the final disappearance of ice from the area. Interestingly, the potential fluvial systems cannot be traced upslope or northward to their submerged channels (across Devonian and Carboniferous rocks) from the data available. As also visible from Fig. 7B, C, there is no actual connection between current onshore river systems and those mapped offshore. In fact, during previous surveys (e.g. NERC project NE/H024301/1, 'Late Glacial Sea Level Minima'), no extensions of channels near the coastline were found. This, therefore, results in a strip of land, between the present-day coastline and the beginning of the offshore palaeovalleys, where no channels appear to have been preserved. This may be partially due to postglacial modification following the supposed age-formation of the channels. However, the available data do not allow us to fully test the different hypotheses for the formation of the easternmost Irish and UK channels and our interpretation remains only tentative.

Discussion
A system of tunnel valleys has been mapped at the seabed and sub-surface across the Celtic Sea (Fig. 8) in the first comprehensive investigation of these features at a regional scale. The tunnel valleys form an extensive and complex meltwater drainage system in this sector of the former BIIS, potentially linked to tunnel valleys recognized in the Irish Sea (e.g. Coughlan et al. 2020). The sets of channels closest to the Irish and British coastlines have been given two different interpretations. Those to the west of southern Ireland, parallel to the coastline, are either subglacial or proglacial meltwater channels (Fig. 7A); whereas the other sets of channels, perpendicular to the southeastern Irish and UK coastlines, might be palaeo-fluvial valleys (Fig. 7B, C). Collectively, they represent different stages in the evolution of the drainage system linked to the history of the Irish Sea Ice Stream (ISIS) and were presumably formed during or after its deglaciation. The following discussion focusses on the formation of the tunnel valleys and their implications for subglacial processes of the former ISIS.
The tunnel valleys documented in this study provide geomorphological evidence for a large-scale (~25 000 km 2 ) subglacial drainage system associated with a grounded ISIS on the Celtic Sea shelf. They imply an increased meltwater influence for this sector of the BIIS compared to the Malin Sea and Scottish shelves, located further north, where evidence of tunnel valleys has not been observed (Bradwell et al. 2008;Arosio et al. 2018;Callard et al. 2018;Tarlati et al. 2020). Although there is no direct chronological constraint on the age of the tunnel valleys, the correlation with the regional stratigraphy suggests formation of the channelized system during the last glaciation (Praeg et al. 2015;Lockhart et al. 2018;Scourse et al. 2019Scourse et al. , 2021Chiverrell et al. 2020;Van Landeghem & Chiverrell 2020). The tunnel valleys appear to be cut into both the basement (Fig. 4) and the glacigenic sediments of the Upper Little Sole Formation (ULSFm; Fig. 6). This indicates that the erosion occurred towards the end of or after the deposition of the ULSFm (27-24.3 ka; Praeg et al. 2015;Scourse et al. 2019) during a period of ice retreat, as the tunnel valleys are all located towards the top of the formation. While the channels cut into the bedrock are widely distributed across the shelf and visible at the sea floor, those that cut into glacigenic sediment are mostly found across the mid-shelf as buried channels (Figs 6, 8).
The tunnel valleys extend up to~76 km in length, although the majority are <20 km long and the average is 13.8 km (Figs S1A, B, S2A, Table S1). Although some of them were found to reach a width of~2.5 km, the majority are between 200-1200 m with an average of 811 m (Fig. S1A, B). In comparison, reported dimensions of tunnel valley length from the North Sea range from a few km (short tunnel valleys~13.5 km long) to around 80-100 km with a few exceeding >100 km (e.g. Stewart et al. 2013;Ottesen et al. 2020) and those from the southern Laurentide Ice Sheet extend for up to 65 km with a median length of 6.4 km (e.g. Livingstone & Clark 2016). While the lengths of the tunnel valleys in the Celtic Sea are comparable to those elsewhere, their widths appear relatively narrower compared to tunnel valley widths reported from Europe (e.g. North Sea; 300 m up to 10 km; Stewart et al. 2013;Huuse & Kristensen 2016;Ottesen et al. 2020) and North America (500-3000 m; e.g. Kehew et al. 2012;Livingstone & Clark 2016). In Britain and Ireland, large subglacial drainage features Fig. 8. Geomorphological map showing the coastlines and related key locations for the interpreted channelized features, both exposed and buried. Offshore the Irish and UK coastlines the channels are drawn with dotted lines due to their uncertain interpretation. The isochrones (thick coloured lines) indicate the ice retreating from the shelf-edge (green line at 25.6AE0.5 ka) to a northward direction (modified from Scourse et al. 2021). The red dashed line represents the boundary between the inner and mid-shelf. TVs = tunnel valleys; MCs = meltwater channels.
(tunnel valleys up to 3 km in width and 200 m in depth) cut into lithified bedrock have been reported to be narrower than those formed in areas with unconsolidated sediments (e.g. Knight 2002;Evans et al. 2005). The fact that most of the Celtic Sea tunnel valleys are cut into lithified bedrock (i.e. Palaeogene and Upper Cretaceous chalk; Fig. 9) can, thus, possibly explain their narrower width compared to the North Sea tunnel valleys, which are eroded into the Pleistocene sediment substrate. Several other factors may contribute to variations in tunnel valley widths (e.g. Sjogren et al. 2002;Livingstone & Clark 2016): thermal or basal conditions over space coupled with the amount of meltwater discharged over time; the rate of ice retreat if the widest segments developed during slow or stable ice withdrawal; migrating of stream resulting in branching and terracing of the valley floor. Largely, the Celtic Sea tunnel valleys display relatively narrow and shallow (max. 29.3 m) U-and V-shaped cross-sectional profiles, with some bifurcation and branching geometry (Fig. 3). The geographical distribution of the mapped tunnel valleys compared with the reconstructed isochrones ( Fig. 8) lead us to exclude that width is related to time of valley formation (e.g. stillstands). Rather it suggests that the tunnel valleys formed and filled relatively quickly, thus preserving comparatively steep tunnel valley walls (Figs 4, 6). This kind of geometry has been reported from similar palaeo-glaciated contexts, for which time-transgressive channel formation has been suggested (e.g. Huuse & Lykke-Andersen 2000;Praeg 2003;Stewart et al. 2013).
As geology seems to have an impact on the shape and dimensions of tunnel valleys (e.g. Knight 2002;Evans et al. 2005;Sandersen & Jørgensen 2012), their distribution is displayed in Fig. 9 over the underlying geological units and fault systems. This investigation enabled us to qualitatively examine if the geological factors may have influenced the ice-bed conditions and the occurrence, in terms of location and orientation, of tunnel valleys in the Celtic Sea. The buried tunnel valleys were eroded into the glacigenic sediments of the ULSFm; thus, these valleys' underlain geology (predominantly Upper Cretaceous and partly Oligo-Miocene mudstone and siltstone; Fig. 9) has poor relevance in their formation.
The exposed tunnel valleys are incised into bedrock, mostly Palaeogene rocks (siliciclastic, argillaceous and sandstone deposits), at times Upper Cretaceous chalk, and more rarely in Oligo-Miocene deposits (mudstone and siltstone, more widespread towards the outer shelf; Figs 1B, 9) (BGS 2013). In terms of erodibility of the geological units, the chalk is likely to be more compacted and cemented (i.e. lower efficiency of glacial erosion) than the Palaeogene and Miocene deposits (Moreau et al. 2012). As almost 70% of all tunnel valleys erode into Palaeogene rocks, there seems to be a possible link between the location of the tunnel valleys in the inner/ mid-shelf and the underlying geology. However, as all tunnel valleys exposed on the Celtic Sea shelf are largely of similar morphology and dimensions across the area, it appears that the underlying basement geology does not clearly affect their morphology. In other settings, as observed in the North Sea Basin, tunnel valleys incised in chalk are narrower and more V-shaped than those formed in unconsolidated sediment or less resistant strata (Praeg 2003;Moreau et al. 2012). In general, also, the valleys cross from chalk (Upper Cretaceous) to Palaeogene rocks with no change in orientation and only one isolated tunnel valley, located in the SE corner, appears aligned with a Permo-Triassic unit (Fig. 9). NE-SW trending faults are present within the Celtic Seabasin (Figs 1B,9). These faults date back to a period of active rifting during the Triassic-Jurassic and Cretaceous, which was followed by a Tertiary inversion (e.g. Roberts 1989;Shannon 1995;Rodr ıguez-Salgado et al. 2019 and references therein). Recent data from the western Celtic Sea (Rodr ıguez-Salgado et al. 2019) have shown that the major NE-SW trending faults may not extend through the Upper Cretaceous, except where inverted. The tunnel valleys are located close to NE-SWoriented faults in the central Celtic Sea and/or near geological boundaries (Figs 1B, 9;GSI 2012;BGS 2013). The close association with such geological structures might indicate erosive exploitation of pre-existing weaknesses into the substrate (Lloyd 2015). Our investigation, however, did not show major faults involving the valleys in correspondence of the Upper Cretaceous chalk and Palaeogene rocks; thus, a structural control on the position of the valleys cannot be asserted with certainty.
Overall, it appears that the tunnel valley locations coincide with the main ISIS axial trough and the distribution of Palaeogene rocks (Figs 1C,9). This might indicate that the position of the ice stream coupled with the underlying lithology may have governed the locus of meltwater accumulation and erosion at the ice-bed on the Celtic Sea shelf. It is possible that the basal topography of the basin might have also influenced the large-scale spatial structure and dynamics of the ice margin (Ign eczi et al. 2018). In turn, this could have facilitated the erosion, meltwater accumulation and eventually influenced the location of the tunnel valleys on the shelf.
The ISIS advance across the Celtic Sea shelf during the last glaciation was short-lived and characterized by a surge-type event over soft marine sediments (Scourse et al. 1990;Scourse & Furze 2001). The advance likely started after 28.3 ka and by 25.6 ka the ice was grounded close to the shelf edge (Praeg et al. 2015;Smedley et al. 2017a;Lockhart et al. 2018;Scourse et al. 2019). The ISIS advance from the narrow St George's Channel occurred at a rate of~350 m a À1 (Smedley et al. 2017a), with a sudden release of ice flowing through the channel (Scourse et al. 2021). During this phase of rapid ice advancing onto clay-rich marine sediments (e.g. Scourse et al. 2019), the hydraulic regime at the ice-bed was possibly not efficient enough to drain all the basal meltwater through the substrate as groundwater flow. An increase in basal water pressure would have led to a separation of the ice from the bed (Piotrowski 1999). This would have enhanced the basal sliding and surging, resulting in ice stream instability, increase in erosive potential and the rapid advance of the ISIS at this stage (Van Landeghem & Chiverrell 2020;Scourse et al. 2021). It is therefore unlikely that the tunnel valley network observed in the Celtic Sea developed under ice surging onto the continental shelf, as its presence would act to reduce ice-bed separation and consequently to stabilize the ice margin (e.g. Kamb 1987).
The ISIS grounding phase at the shelf-edge of the Celtic Sea was relatively short-lived (Scourse et al. 2019). After the initial ice pullback by 25.2 ka, the ice margin receded rapidly towards a mid-shelf position by 25 ka (Fig. 8; Scourse et al. 2021). During the overall retreat, the tunnel valleys likely cut into the ULSFm as timetransgressive features, similar to the carving of the Wicklow Trough at a later stage of retreat in the Irish Sea (Coughlan et al. 2020). It is also possible that the sediment transported by meltwater contributed as a major source for the MFm sediments that form the megaridges. The presence of eskers to the east of the large tunnel valley in the central Celtic Sea (Figs 5, 8) strengthens the interpretation of a time-transgressive formation of the tunnel valleys. Despite the fact that the eskers are restricted to a relatively small area (~275 km 2 ; Figs 5, 8), their overall distribution, together with the presence of the tunnel valleys, reflects the subglacial drainage organization beneath the ice sheet. In fact, the eskers' long-axes tend to be perpendicular to the ice margin (Boulton et al. 2009;Greenwood et al. 2016;Lewington et al. 2020). The absence of eskers in other parts of the Celtic Sea may be related to preservation issues but also to changes in the direction of the ice retreat from a NE-SW flow to a more northward flow (Van Landeghem & Chiverrell 2020).
The presence of tunnel valleys incised into the bedrock on the inner and part of the mid-shelf and infilled by modern sediment (Fig. 4) suggests that their downcutting also occurred during the rapid retreat of the grounded ice. It is probable that discharges of hydrostatically pressurized meltwater were significantly erosive. Erosional processes at the ISIS bed could also account for the limited sediment deposition over bedrock in the inner part of the shelf (Figs 3, 4), compared to the midshelf (Fig. 6). The contrast in sediment thicknesses is perhaps similar to previous advance phases, during which erosion and advection of subglacial sediment to the ice margin occurred (cf. Boulton et al. 1996).
The tunnel valleys have not been observed in the innermost part of the shelf, but they reappear in the Irish Sea marking a deceleration of the grounded ice margin, from a rapid to slow pace of retreat (Coughlan et al. 2020;Van Landeghem & Chiverrell 2020). Different factors could have contributed to this: variations in meltwater availability, changes in bedrock geology and increased rate of retreat as the ISIS was leaving the Celtic Sea (Scourse et al. 2021). It has also been inferred that the reverse slope of the Celtic Deep might have favoured such rapid retreat (between 230 and 650 m a À1 ; ibid.) limiting the time for headward erosion or producing a segmented tunnel valley (Livingstone & Clark 2016). Recent analyses from the southern sector of the former Laurentide Ice Sheet demonstrated that the formation and morphology of tunnel valleys can be tightly coupled to variations in basal thermal regime, bed topography, time of erosion and climate (Livingstone & Clark 2016). It is not possible to identify a main driver, but it is probable that a combination of the above factors resulted in the distribution of tunnel valleys observed in the Celtic Sea compared to those identified further north, in the Irish Sea.
The presence in the Celtic Sea shelf of a channelized subglacial drainage network highlights the key role of meltwater erosion during ice evacuation from the continental shelf. Our study also shows that the large subglacial system might have formed within a restricted time frame of~1-2 ka (Scourse et al. 2021) by the meltwater delivered to the bed as the ice sheet deglaciated (Fig. 8).

Conclusions
• Using high-resolution bathymetric and seismic data, an extensive system of exposed and buried tunnel valleys has been mapped in the central part of the Celtic Sea shelf (Fig. 8). • The network of tunnel valleys consists of around 80 steep-sided incisions, up to 76 km long, largely located on the inner and mid-shelf. Recent geochronological and stratigraphical reconstructions in the same area suggest that the subglacial meltwater drainage system may be linked to the last Irish Sea Ice Stream (ISIS) retreat phase across the shelf, between 25.6 and 24.3 ka (Fig. 8). • Bedrock lithology and structure in the form of preexisting faults have been qualitatively investigated to link a possible control on the channels' formation with meltwater exploiting pre-existing weaknesses to erode the tunnel valleys. Observations show that there is no obvious control of the basement lithology and structural geology on tunnel valley distribution and orientation (Fig. 9). Rather, their location largely coincides with the main ISIS axial path, indicating that the position of the ice stream may have controlled the location of the greatest meltwater erosion rather than the inherited structural setting. • The overall retreat occurred within~1-2 ka, with the formation of overdeepened and steep-sided tunnel valleys and highlights the erosive power of subglacial meltwater at this stage of fast ice withdrawal from the shelf. The presence of highly erosive meltwater discharge on the relatively shallow Celtic Sea shelf further emphasizes the contribution of meltwater to this area of the BIIS, in contrast with its sectors located further north (e.g. Malin Sea shelf). • Tunnel valleys are absent from the inner part of the shelf, but they reappear further north in the Irish Sea Basin (Wicklow Trough), likely as part of a segmented system of tunnel valleys through the Irish Sea Basin and onto the Celtic Sea shelf.

Supporting Information
Additional Supporting Information may be found in the online version of this article at http://www.boreas.dk.  S2. The histograms show the frequency distribution of the exposed channels (central Celtic Sea shelf) and the easternmost Irish and UK channel (perpendicular to the coastlines) lengths, and related basic statistical measures. Due to the small total count (n ≤ 10), the lengths of the other landforms have been included in Table S1.  Table S1. Total length values for the landforms described as channels. For the buried channels, Li corresponds to the 'interpolated system length' used to calculate the system continuity. Li is coupledwith letters that refer to the labelled buried tunnel valleys shown in Fig. 6. For the exposed tunnel valleys and palaeo-fluvial valleys histograms of length are included in Fig. S2. Nevertheless, their length measurements have also been added in the table below for thoroughness.