Cool deltas: Sedimentological, geomorphological and geophysical characterization of ice‐contact deltas and implications for their reservoir properties (Salpausselkä, Finland)

Sediments deposited by glacial meltwaters (for example, ice‐contact delta deposits) form permeable packages in the subsurface that can act as reservoirs for both water and hydrocarbons. They are also an important source of aggregate for the construction industry. As reservoirs they are challenging to characterize in terms of their structure, flow and storage properties due to their complex depositional history. In this study, ice‐contact deltas of Salpausselkä I and II end moraines in Southern Finland are studied using a combination of geomorphological mapping, sedimentological studies and near surface geophysical methods. Sedimentary logs from isolated outcrops were correlated to ground penetrating radar (GPR) profiles to unravel the internal structure and depositional history of these ice‐contact deltas. Subsequently, electrical resistivity tomography (ERT) and gravity data were analysed to estimate the depth to bedrock and to model porosity distribution within the sediments. Results of the study suggest that the delta deposits have a broad range of porosities (10 to 42%) with lowest values found in the bottomset beds. The most variable porosities are in the subaqueous ice‐contact–fan zone, and consistently high porosities occur in delta foreset/topset facies. Detailed sedimentary logging linked to the GPR data shows heterogeneities such as mud drapes on foresets and kettle holes which are below the resolution of ERT and gravity methods but significantly affect reservoir properties of the deltas. Moreover, oscillation of the ice‐margin may have introduced larger heterogeneities (for example, buried ice marginal ridges, or eskers) into the sedimentary sequence which are atypical for other Gilbert‐type deltas. Finally, subglacially sculpted, highly variable bedrock topography exerts a major control on sediment distribution within the delta making reservoir volume and quality less predictable. This work has implications for present‐day freshwater aquifers and low enthalpy geothermal energy in southern Finland and other deglaciated regions, as well as hydrocarbon exploration of analogous deposits in the subsurface from Pleistocene and pre‐Pleistocene glaciogenic sequences.

Unconsolidated, Pleistocene and Holocene glaciogenic sediments are especially important as aquifers in regions where repeated episodes of subglacial erosion have removed the older, pre-Pleistocene porous, sedimentary cover, leaving only impermeable crystalline basement across much of the area (Knutsson, 2008;Comte et al., 2012;Ravier & Buoncristiani, 2017). In such places, glaciogenic sands and gravels are often the only viable aquifers for human consumption, agriculture and industry. Understanding their distribution and internal structure is crucial for effective groundwater extraction and also for protection and management (Anderson, 1989;Poeter & Gaylord, 1990;Heinz et al., 2003;Slomka & Eyles, 2013;Best et al., 2015;Erickson et al., 2019;Ó Dochartaigh et al., 2019). Glaciogenic sediments also act as hydrocarbon reservoirs in areas such as the North Sea  and across North Africa (e.g. Osterloff et al., 2004a;Slomka & Eyles, 2013;Kurjanski et al., 2020). The ability to predict physical properties of such sediments largely depends on understanding the depositional processes responsible for their formation (Lønne, 1995;Zieliński & van Loon, 2000;Heinz et al., 2003;Catuneanu, 2006;Hirst, 2012;Slomka & Eyles, 2013;Zecchin et al., 2015;Dietrich et al., 2017a;Winsemann et al., 2018;Lang et al., 2021). This information is typically obtained by detailed sedimentological studies of wells and/or outcrop analogues (Howell et al., 2008). Given that wells and many outcrops are spatially limited and two-dimensional it can be difficult to extrapolate results over a wider area especially given the heterogeneity of glaciogenic sequences (Hirst et al., 2002;Lajeunesse & Allard, 2002;Reinardy & Lukas, 2009;Evans et al., 2012). Shallow geophysical methods can be used to bridge that gap and provide missing information between wells or exposures, but are labour-intensive and costly. Moreover, seismic resolution decreases with depth thus affecting the imaging of deeply buried sediments (Andersen et al., 2012;Ottesen et al., 2012;Rose et al., 2016;Bataller et al., 2019). S-wave seismic with a land-streamer system can be helpful because it enables fast acquisition of data and can provide high-resolution images of the near subsurface Winsemann et al., 2011Winsemann et al., , 2018. On the other hand, aquifer studies often rely heavily on near surface geophysical methods which are characterized either by high resolution and shallow penetration depth (for example, ground penetrating radar) or deeper penetration with lower resolution (for example, electrical resistivity tomography, seismic refraction and gravity), only imaging large contrasts in physical properties and discontinuities. Careful consideration is needed when choosing methods to deploy, ideally allowing identification of key heterogeneities (Gabriel et al., 2003;Andersen et al., 2012;Galazoulas et al., 2015;Paz et al., 2017;Ravier & Buoncristiani, 2017).
The use of modern analogues is a standard practice in evaluations of ancient sedimentary successions (Galloway, 1975;Orton & Reading, 1993;Pye & Lancaster, 1993;Pye, 1993;Hartley et al., 2010;Nyberg & Howell, 2016). This approach is justified when evaluating clastic marine, fluvial or aeolian successions, since their scale and formation processes, with some exceptions, are reasonably similar across geological time. However, it may not be ideal when analysing glaciogenic deposits associated with continental-scale ice sheets, because of the different spatial and temporal scales of the processes involved. Parallels between modern and ancient glacial systems can still be useful, but the differences of scale and magnitude of observed processes have to be taken into consideration (e.g. Spedding & Evans, 2002;Evans, 2006;Evans et al., 2009Evans et al., , 2017Slomka & Eyles, 2015;Lee et al., 2018). The best approach to understanding deep-time, ice-sheet derived sediments is to study glaciogenic deposits remaining on the land surface following the deglaciation of a region (Dietrich et al., 2018;Dowdeswell et al., 2019). Especially valuable study areas are those located where an ice-sheet margin was grounded in water and which are, today, above sea level; of these, the Salpausselkä I and II are excellent examples.
In this study a series of ice-contact deltas that constitute part of the Salpausselkä I and II end moraines in southern Finland, were investigated. The deltas are important aquifers and sources of aggregate for the region. Geomorphological mapping and sedimentological studies were integrated with a suite of near-surface geophysical methods [ground penetrating radar (GPR), electrical resistivity tomography (ERT) and gravity] to identify the sedimentary architecture, key heterogeneities, and provide insights into the formation and reservoir properties. Results of this study are applicable to groundwater reserve estimation and management in glaciofluvial deposits of southern Finland, but also to water and hydrocarbon exploration in glaciogenic sequences in general.
The study area is the Lahti region, Finland, which was at the confluence zone between the Baltic and Finnish Lake District ice lobes during the YD Stadial (Fig. 1). Four outcrops of Salpausselkä I were visited and studied in active gravel pits located in Kukonkangas, 15 km west of Lahti (Fig. 2, locations 1 to 4). Two further Salpausselkä I localities were studied in the area known as Hälvälänkangas, 15 km WNW of Lahti (Fig. 2, locations 5 and 6). Salpausselkä II was studied from sites located near Vierumäki, 20 km NNE of Lahti (Fig. 2,locations 8 to 13). The data were collected from the Syrjälänkangas (GPR, ERT and outcrop), Hyrtiälänkangas (GPR and ERT), Vesivehmaankangas (GPR and ERT) and Mikkolankangas (outcrop and GPR) areas (Fig. 2).

Geology and morphology of the study area
The bedrock geology is dominated by crystalline rocks of the Fennoscandian Shield formed, or metamorphosed, during the Svecofennian Orogeny, 1.92 to 1.77 Ga (Paulamäki et al., 2002;Nironen, 2017). In the Lahti region granite, microcline granite, quartzite, gabbro and paragneiss are dominant (Nironen, 2017). Late Cenozoic glaciations are responsible for the removal of pre-Pleistocene strata and sculpting of the crystalline bedrock surface (Figs 1 to 3).
The sedimentary cover is thin, ranging between 0 to 30 m with the exception of the Salpausselkä I and Salpausselkä II ridges, where thicknesses often exceed 50 m (Fyfe, 1990;Palmu, 1999;Johansson et al., 2011). The Salpausselkä I and Salpausselkä II deposits are stratigraphically higher than the underlying till and sand (Figs 1 and 2).

METHODOLOGY
This study combines outcrop analyses, and acquisition, processing and interpretation of ground penetrating radar (GPR) and electrical resistivity tomography (ERT) data. These were analysed jointly and interpreted in combination with pre-existing Light Detection and Ranging (LiDaR) Digital Elevation Model (DEM) and gravity (Bouger anomaly) data, to provide information on the depositional processes, sedimentary environments, reservoir geometries and properties of the ice-contact deltas. Additional metrics and reservoir properties were calculated.

Geomorphology
A LiDaR DEM (1 m horizontal resolution and 0.3 m vertical resolution) (Fig. 2), provided courtesy of the National Land Survey of Finland and processed by the Geological Survey of Finland, was used to map and analyse the morphology of the study area (Figs 2 and 3). Mapping was conducted using ArcGIS©. The landforms were mapped based on their expression in the DEM, overlain with a 75% transparent multidirectional hill shaded raster, and a 75% transparent slope raster, to best highlight the morphological features. Morphometrics and topographic profiles were also extracted using ArcGIS (Fig. 2).

Sedimentology
Fourteen outcrop sections, ranging between 4.5 to 23 m in height, were logged in detail, up to the decimetre scale, where possible, in active gravel workings. Five sections were located in Salpausselkä I ( Fig. 2A), whilst two further gravel pits in Salpausselkä II were studied (Fig. 2C).

Ground penetrating radar
The GPR profiles were acquired in the Vierumäki region in Salpausselkä II (Fig. 3). Fifteen  80 MHz GPR profiles (14 km) were acquired with a horizontal resolution of 11 cm and a vertical resolution estimated to be 0.2 to 0.6 m, depending on the dielectric properties of the material. The penetration depth is between 5 to 25 m. GPR profiles were collected along existing forest roads, as close as possible to the depositional dip and strike directions (Fig. 3). In some locations, where outcrops were close to the survey line, it was possible to directly correlate GPR results with lithologies and sedimentary structures observed in the outcrop (Fig. 3). For equipment and processing details see Table S1.

Electrical resistivity tomography
The ERT survey (Fig. 3) was undertaken to determine the depth to bedrock and to correlate resistivity zones with the reflection packages identified using the GPR and, in turn, link to the sedimentary facies, to better understand the sedimentary architecture at depth. Three ERT profiles were acquired across the study area ( Fig. 3): ERT1 (1075 m length); ERT2 (720 m); and ERT3 (1875 m), in a north-south orientation (parallel to the depositional dip) and co-linear with GPR profiles GPR 1, GPR 4 and GPR 14, respectively (Fig. 3). Dipole-dipole (DD) and multigradient (mGD) quadripole arrays were used jointly to ensure optimum resolution for imaging the complex arrangement of glacial sediments overlying fractured bedrock (Comte et al., 2012). The depth of investigation (DOI) method (Oldenburg & Li, 1999) was further applied to determine the maximum reliable depth of investigation. For equipment and processing details see Table S1.

Porosity from ERT
The porosity of the individual zones was calculated based on Archie's Law (Eq. 1) which relates pore fluid conductivity (C t ) with the resistivity of the bulk volume (R t ) (Eq. 2) (Archie, 1942): where: ϕ is porosity (0 to 1), R t is resistivity from ERT profiles [Ωm], C w is conductivity of formation waters [S/m], C t is conductivity of bulk rock [S/m], S w is water saturation (0 to 1), n is water saturation exponent and m is cementation factor. Tortuosity factor -α is unknown and was assumed to be unity (1) to enable comparison between results. Combining Eqs 1 and 2 allows porosity (ϕ = 0 to 1) to be calculated. This method was chosen because it is suitable for clay-poor formations of Salpausselakäs (Okko, 1962;Glückert, 1986;Fyfe, 1990;Palmu, 1999). Parameters used for porosity calculation from ERT profiles are listed in Table 3. They were modelled using a Monte Carlo simulation based on a normal distribution of parameters. Input parameters (mean and standard deviation) were calculated based on available literature (see Tables 1 and 2 for details). Water saturation values were determined based on boundary conditions: dry and saturated sediments should have similar porosity values and the resistivity response is mainly governed by water saturation since the sediments within deltaic foresets are well-sorted and uniform. The results were validated by substituting calculated porosity values with average porosities of glaciofluvial sediments sourced from the literature and comparing the resulting resistivity values with resistivity recorded in the field.

Gravity and shallow borehole data
Gravity survey and borehole data, in a series of groundwater investigation reports, were available from the Geological Survey of Finland data repository, (www.hakku.fi; Fig. 3) (Breilin et al., 2005(Breilin et al., , 2006Ahonen et al., 2011). Bouguer anomaly data were sampled approximately every 20 m along profiles (Fig. 3). Final interpretations (after inversion) of gravity profiles as well as published depth-to-basement maps were used to supplement interpretation.

Ice-contact deltas
A series of flat-topped hills are present in the study area (Fig. 2). Their orientation is transverse to the ice flow during the Younger Dryas Stadial (Fig. 1). The ridges can be subdivided into lobate or irregular segments which, in places, amalgamate forming a broad plateau (Fig. 2). Eighty-six segments, 600 m to 19 800 m long (mean 4040 m, median 2830 m) and 190 to 7500 m wide (mean 1490 m, median 1020 m), were identified. Their height varies between the northern, ice-proximal side (55 to 65 m above the neighbouring topography) and the southern, ice-distal, side (45 to 50 m), but can reach up to 100 m near Aurinkovuori (N61.190159, E25.477587). The slope of the ice-proximal side ranges between 8.5°and 21.0°(mean 13.1°, median 12.2°) and is steeper than the ice-distal side (0.1°to 29.0°, mean 6.5°, median 4.8°). Plateau surfaces slope towards the south between 0.1°© Loaded and 2.8°(mean 0.4°, median 0.3°). The distal margin of the plateau is often characterized by a rim of sediment with several major breaks of slope (terraces) below and multiple intermediate, smaller slope breaks. The rim is interpreted as a shoreline delineating the level of the Baltic ice lake to which the delta aggraded, i.e. the emergent delta (top) surface ( Fig. 2) (Glückert, 1995). Subsequent slope breaks are interpreted to represent wave erosion with shorelines forming as a result of the interplay between isostatic rebound and lake-level fluctuations (Jantunen & Donner, 1996;Eronen et al., 2001).

Braided channels
The top surfaces of ice-contact deltas are often undulating, exhibiting a radial, braided pattern of furrows which are largely discontinuous and rarely extend the full breadth of the plateau (Fig. 2). Instead, they cross-cut or merge, creating a radial, dendritic pattern originating from an apex in the elevated, proximal part of the plateau. Multiple dendritic patterns can be observed on larger plateaus. Individual furrows can be traced for 8 to 2150 m (mean 148 m, median 222 m). They are 2 to 110 m wide (mean 50 m, median 35.5 m) and 0.2 to 10 m deep (mean 1.7 m, median 0.9 m) and interpreted to represent a distributary network of braided, channels preserved on emergent delta tops (Donner, 1976).

Kettle holes
Oval and irregular depressions are visible on the delta tops and some eskers (Fig. 2). The area of the depressions ranges between 36 m 2 and 1 042 300 m 2 (median 1428 m 2 ). Individual depressions are 1 to 50 m deep (median 8 m). They are more abundant in the ice-proximal zones. Braided meltwater channels on delta tops can be on both sides of many depressions, which implies that the formation of depressions post-dates the delta formation. The depressions are interpreted as kettle holes formed from melt-out of buried or detached ice blocks (e.g. Wingfield, 1990;Palmu, 1999;Mäkinen & Palmu, 2008;Maries et al., 2017). The size distribution of kettle holes allows division into three categories.
(1) The largest kettle holes occur in the south-east part of Salpausselkä I delta complex, near Lahti ( Fig. 2A). They are positioned in the leeside of a local bedrock high and potentially represent large dead ice blocks that became detached as the ice margin thinned and retreated over the bedrock high.
(2) Large, deep and elongated depressions exhibiting a linear or radial pattern (Fig. 2) are located mainly in the proximal, ice-contact parts of the hills, most likely representing dead ice blocks cut off from the ice margin by meltwater flow.
(3) Small oval and irregular shaped, shallow kettle holes were likely formed by melt out of glaciofluvially deposited ice blocks, based on their distribution in the medial/distal part of the delta (Fig. 2).

Push moraines
Discontinuous, curvilinear or straight ridges are present locally on the proximal slope break. They range between 2 to 7 m (mean 3 m) in height above the top of the plateau and are 20 to 150 m wide (median 42 m, mean 55 m). They are interpreted as push moraines linked to ice-marginal oscillations (Aber et al., 1989;Bennett, 2001;Benn & Evans, 2010). Where such oscillations were more substantial, ice overrode the proximal slope of the deltas and push moraines were on the delta top. (Fig. 2).

Interpretation
Based on the morphologies of the flat-topped hills and associated secondary landforms, the hills were interpreted as ice-contact Gilbert-type deltas deposited in a lake with water depths of 35 to 60 m above the present-day land surface (Gilbert, 1885(Gilbert, , 1890Glückert, 1986;Lønne, 1995;Palmu, 1999;Saarnisto & Saarinen, 2001;Thomas & Chiverrell, 2006;Dietrich et al., 2017b;Kurjanski et al., 2020). Meltwater and sediments were supplied to the ice margin chiefly via subglacial/englacial meltwater channels and deposited under water, proximal to the grounding line as ice-contact, subaqueous fans (Powell, 1990). The sediment supply was sufficiently great that fans aggraded to the contemporaneous water surface resulting in ice-contact delta formation and subsequent progradation. Following delta emergence, the sediments were transported in sub-aerially fluvial channels across the delta-top (topset facies) and were deposited on as sub-aqueous foresets and bottomsets (Gobo et al., 2015;Dietrich et al., 2017b;Lang et al., 2017b). Sediment and meltwater supply were terminated abruptly when the ice margin stepped back, first from the Salpausselkä I and later from the Salpausselkä II moraine. Isostatic rebound and the opening of seaways connecting the Baltic ice lake to the global ocean were responsible for drainage of the lake and the emergence of southern Finland (Bjorck & Digerfeldt, 1984;Glückert, 1995;Stroeven et al., 2016). As a result, the deltas constitute high points in the present-day topography (Fig. 2).

Eskers
Sinuous or curvilinear, often discontinuous sediment ridges with steep sides (Fig. 2) running parallel or sub-parallel to the Younger Dryas Stadial ice-flow direction are observed ( Fig. 1) (Johansson et al., 2011). They are 8 to 40 m high (median 14 m), 16 to 920 m wide (median 78 m). The length of continuous segments ranges between 21 m and 17 400 m (median 174 m). The ridges often climb uphill and do not follow the topography. In places, ridges are flat-topped with numerous circular or irregular depressions. Smaller ridges are often amalgamating and anabranching, exhibiting a radial pattern and merge directly with the Salpausselkä icemarginal ridges (Fig. 2).
The ridges are interpreted as eskers formed subglacially in meltwater conduits (Benn & Evans, 2010;Storrar et al., 2014;Schomacker & Benediktsson, 2017). Meltwater flowed from north to south, driven by the ice-sheet surface slope (Shackleton et al., 2018). Higher and longer eskers represent major meltwater paths draining farther into the interior of the ice sheet, whereas smaller ridges are interpreted as distributary meltwater channels located proximal to the ice margin ( Fig. 2) (e.g. Storrar et al., 2019). Their spatial relationship with ice-contact-deltas implies that the eskers acted as feeders, supplying sediment-rich meltwater to the ice margin.

Ice-contact delta metrics
Eighty-six deltas were identified along Salpausselkä I and Salpausselkä II. Their length, width and outer circumference were measured and a lobosity parameter was calculated following the methodology outlined by Howell et al. (2008) (Fig. 4). The vast majority of deltas (82%) are <5 km wide and <3 km long (Fig. 4). Lobosity varies between 2.2 and 1.09 with the majority of deltas being <1.5 (Fig. 4). Low lobosity values are linked to amalgamated deltas in a 'bajadastyle', where multiple point sourced fans grow to the lake surface forming point sourced deltas, which then coalesce to form one, bench-like delta, for example at Syrjälänkangas and Hyrtiälänkangas (Figs 2 and 4) (Barrett et al., 2020a). Higher-lobosity deltas (>1.7) are clearly point sourced (for example, Vehivesmaankangas - Fig. 2). Two trends can be delineated from the distribution of dimension and lobosity data: (i) amalgamated deltas (length <<< width, with low lobosity); and (ii) lobate deltas (length <= width, with high lobosity) (Fig. 4). Deltas with length <= width but with low lobosity are constrained by topographic highs which prevent them growing laterally. The majority (88%) of deltas with lobosities >1.7 are <2.5 km long and <3.7 km wide ( Fig. 4).

SEDIMENTOLOGY
Sixteen facies were identified, based on the dominant grain-size, sedimentary structures and the general expression in the outcrop (Table 1). Facies were then grouped into eight associations which have been distinguished from logs and outcrops based on the facies successions, dominant grain-size and sorting, sedimentary structures and position with respect to the ice margin at the time of deposition ( Table 2). The latter is interpreted from the morphology of individual deltas (Fig. 2). The distribution of facies associations is described in the context of architectural elements comprising an ice-contact delta and juxtaposed against reflection patterns observed from GPR profiles (Figs 9, 10 and 11).

Architectural elements
Ice-contact slope The ice-contact slope of the delta, where bulldozing and traction by ice are dominant, is characterized by boulder and cobble-sized clasts suspended in a fine-grained, silty and sandy matrix forming subglacial diamicton (FA1) (Table 2, Figs 5, 8H and 8I). Ice-margin advances/oscillations may cannibalize previously deposited, ice-contact fans forming high density debris flows that descended the distal slope. As a result, the ice-contact slope is comprised of subglacial facies association (FA 1) mixed with ice-contact facies association (FA 2) (Figs 9 and 12).
The GPR profiles show that this mid to high amplitude, contorted or subparallel, discontinuous reflection package is present only on the ice-contact slope and in the top 1 to 6 m of the delta (Figs 10 and 11). Frequent hyperbolic reflections indicate large boulders incorporated within the matrix (Fig. 10A). Reflections typically dip in the opposite direction to the delta progradation direction. The depositional limit of this facies in the GPR can be correlated with push moraines visible in the delta morphology linking them to ice overriding (Figs 2 and 12).

Subaqueous ice-contact fan
Highly variable grain size and bedding dips are characteristic of subaqueous ice-contact fans ( Table 2, Figs 5 and 9). Steeply dipping boulder and cobble gravels (FA 2) either pass sharply into, or are interbedded with, planar and low angle sinusoidally stratified and/or trough crossbedded sands and gravels of FA 3 interpreted as cyclic step deposits ( Fig. 8A and D) (Ghienne et al., 2021;Lang et al., 2021). FA 3 passes laterally into planar, ripple and climbing ripple laminated sands with the subordinate silty draping interbeds of FA 5 indicating initial high energy, supercritical flow (jet) deceleration, deposition of coarse fraction at the mouth of the meltwater portal due to flow expansion (FA 2) and transition from supercritical flow (FA 3) to subcritical flow (FA 5) (Figs 8A, 8B, 8C and 12) (Hornung et al., 2007;Winsemann et al., 2009). The facies transition is observed to happen on relatively short distances (tens of metres) laterally, as opposed to hundreds of metres in the down-flow direction (based on GPR data) possibly indicating a very high energy jet and highly supercritical flow conditions (Figs 8A and 12) (Hoyal et al., 2003;Russell & Arnott, 2003;Hornung et al., 2007;Winsemann et al., 2009).
In GPR profiles the subaqueous ice-contact fan is characterized by low to high-amplitude, moderately dipping, semi-continuous reflections ( Fig. 10A and E). Large contrasts in reflection amplitudes are most likely due to large variations in grain size and reflect variable energy of deposition between periods of high discharge (cobbles and boulders) and low discharge/flow FA 1 is interpreted as traction till deposited subglacially . When present on top of the succession, it represents a partial overriding of the delta by the ice sheet. Better sorted lenses/interbeds were either re-deposited as subglacially rafted freeze-on/ripup clasts or deposited by subglacial meltwater at the ice-bed interface (Kessler et al., 2012). Sandy and pebbly matrix together with partially preserved primary sedimentary structures (ripple lamination, low angle lamination) in some deformed beds indicates cannibalization and remobilization of fan/delta sediments FA 2 Ice-contact facies association Cmb, Cx, Gx, Gm, (Sp) This facies association is found in the imminent vicinity of the ice margin. Steeply dipping (5 to 25˚), planar, large-scale cross-bedded and massive, bouldercobble gravels (0.2 to 0.8 m thick) with clasts up to 40 cm in diameter (most <20 cm) (Figs 5, 8A and 8D).
In some beds a fining-upward trend is observed with occasional coarser/armoured top bed surfaces. Sand is either absent or very subordinate and/or confined to lenses in the leeside of larger boulders. Pebble gravel beds exhibit gentle fining-upward trend with more pronounced imbrication towards the top of the bed.
This facies association shows sharp lateral contact with facies association FA3 FA 2 is interpreted to represent direct, ice-contact, fan sediments deposited almost immediately at the terminus of a subglacial conduit from high energy flows (Powell, 1990;Lønne, 1995). Massive or normally graded, planar cobble and boulder gravels may represent local hydraulic jump deposits with base scours larger than the width of the outcrop (ca 20 m) (Russell & Arnott, 2003;Winsemann et al., 2009;Lang et al., 2021). Some of the massive, poorly sorted cobbleboulder gravel beds may represent debris-flow deposits caused by discrete slope failures as the sediment pile passes the slope of repose or is bulldozed by the advancing ice margin. (Hornung et al., 2007).
Coarsening-upward packages with armoured top surfaces were formed as a result of winnowing by high energy jets which removed gravel and pebble fractions (Gobo et al., 2014). Lower energy meltwater pulses deposited medium and coarse gravels with finingupward trend due to waning flow energy. Cross-bedded gravels are interpreted as 2D/3D dunes (Russell & Arnott, 2003). Sand lenses were deposited in the flow shadow behind larger boulders indicating that most sand was winnowed or bypassed this zone. Sharp lateral contact of this facies association with FA 3 is interpreted as deposition in a large-scale (>20 m wide, >10 m high) axial scour of a jet efflux  (Fig. 7H). The succession is often cross-cut by small-scale reverse faults. Gravel beds significantly incise and truncate finer grained strata. The dip of the strata ranges between 2˚and 25˚. Locally, convex lenses of massive cobble to pebble gravels (Cmb, Gm) with subordinate boulders are incised into underlying strata and onlapped by parallel laminated sands and sandy gravels above ( Fig. 8A and C). FA 3 was observed to transition laterally into FA 5 over the distance of tens of metres (30 m in one of the outcrops) ( Fig. 8A and B) FA 3 is interpreted to be deposited subaqueously from fast-flowing glacial meltwater in a proximal subaqueous fan or steep delta slope setting (Lang et al., 2021). Multiple erosional contacts within low angle trough cross-bedded gravels are interpreted as progradation scour fills or chute-and pool deposits (e.g. Lang et al., 2017Lang et al., , 2021Winsemann et al., 2018). Planar and low angle cross-bedded sands and pebbly sands were most likely deposited by breaking and non-breaking antidunes with wavelengths >10 m (Lang et al., 2021).
Horizontal sands and pebbly sands are interpreted as upper-stage plane bed deposits deposited by bottomhugging subcritical flows. It is noteworthy that the large wavelength of some of the bedforms makes it difficult to confidently identify bedforms in outcrop and distinguish between subcritical and supercritical conditions (Lang et al., 2021). Localized lenses of massive or crudely bedded boulder and cobble gravels are interpreted as debris-flow deposits ( (Gobo et al., 2014(Gobo et al., , 2015. Large-scale tangential cross-beds (Stl) are interpreted as aggrading backsets of cyclic steps deposited from supercritical flows filling long-wavelength erosional scours (Winsemann et al., 2009;Dietrich et al., 2016;Lang et al., 2021). Trough cross-bedded sands (Stx) represent 3D dunes and antidunes deposited from transitional (supercriticalsubcritical) flows (Alexander et al., 2001;Dietrich et al., 2016;Winsemann et al., 2018). Sand and gravel-filled scours represent chute-and-pool deposits (Lang et al., 2021). Lenses of clast-supported, massive or crudely bedded conglomerates (Cmb) are interpreted as debris flows associated with high magnitude events transporting the coarse fraction from delta topsets down the delta front (Lin & Bhattacharya, 2020 vary in thickness between 5 to 70 cm. Subordinate, coarser packages of low angle planar and trough crossbedded sands and pebbly sands are locally present. Contorted sands and silts with very prominent soft sediment deformation structures (SSD): load clasts, balland-pillows and complex folds are common ( Fig. 6A and B). In the most distal part of the delta, soft sediment deformation exhibits a characteristic facies succession: medium, well-sorted cross-bedded sands, overlain by contorted sands and silts and capped by draping and planar laminated silts and sands (Fig. 5) This facies association is interpreted as low-energy fan/ delta bottomsets. The fining-upward succession represents waning meltwater pulses transitioning from upper plane through current type A and type B climbing ripples and draping lamination deposited after the flow decelerated below ripple migration threshold (Ashley et al., 1982). The contorted sands with soft sediment deformation structures (SSD) (Fig. 6A and B) are interpreted as slump deposits. From outcrops, it appears that shearing and sliding along a slip surface, possibly gravity driven, occurred within the crossbedded sandy package. As a result, more cohesive silts and fine sands are contorted and folded together with a portion of the underlying coarser sand package, which acted as a 'ball bearing' sub-stratum rather than a clear, planar detachment surface. Following slumping, wavy and planar laminated fine sands and silts were deposited from suspension, capping the sequence. Similar SSDs have often been attributed to earthquakes, resulting from glacial isostatic rebound (Gruszka & van Loon, 2007; Loon & Pisarska-Jamrozy, 2014; van Loon  , 2016). However, slumping in the study area is more likely to be related to autogenic processes of delta progradation, when the slope of the delta exceeds the slope of repose (Pisarska-Jamrozy & Weckwerth, 2013), which may be exacerbated by fluctuating lake-water levels. The finer grain sizes are often characteristic of delta bottomset/prodelta environments (Mastalerz, 1990;Postma, 1990;Nemec & Postma, 1993). However, FA 5 has been also described from a site very close to the ice margin ca 60 m laterally from the mouth of an interpreted meltwater portal ( Fig. 8A (Figs 7D and 8G). Particle diameter is highly variable but rarely decreases below pebble gravels. Individual clasts can reach up to 1 m in diameter in the ice-proximal part of the delta (Fig. 7D). Pebble to boulder imbricated gravel packages can be traced both in the depositional strike and dip directions on the scale of metres to tens of metres (typically <30 m).
Lenticular or U-shaped bodies (length 1 to 8 m, height 0.4 to 1.2 m) with erosional base and filled with massive or crudely bedded sediments which are coarser than neighbouring facies are locally present (Fig. 8G). The coarse sediments are interbedded with subordinate coarse and medium sands (Sp, Spx, Stx) which are either filling voids between larger clasts or form a 'tail' on the lee side of larger boulders or cross-bedded gravel sets The succession has been interpreted to represent a braided delta top, which is reinforced by the morphology of the plateau, with remnants of channel networks still clearly visible (Fig. 2) (e.g. Glückert, 1986;Fyfe, 1990;Postma, 1990;Winsemann et al., 2018). As a result of little or no accommodation and high energy of glaciofluvial depositional environment, multiple erosional internal contacts between facies are present. Lenticular concave up bodies of coarse pebble and cobble gravels are present and interpreted as channelfill deposits. Some isolated pockets of cobble and boulder gravels may also represent supercritical flow deposits associated with local hydraulic jumps formed at obstacles or backsets of antidunes (Russell & Arnott, 2003;Lang et al., 2021). Planar cross-bedded pebble and cobble gravels are interpreted as mid-channel bars, migrating downstream. Bedforms associated with subcritical and supercritical flows are present within the succession indicating variable flow regimes in line with seasonality of glacial meltwater discharge (Marren, 2005). Preservation of sands mainly in flow shadow zones indicate that the meltwater energy was high enough to either erode or prevent sand from depositing,  Fig. 6F).
Boulder and cobble gravel beds are discontinuous with erosional internal contacts (Fig. 6E). The location of the outcrop (sediment log 4 - Fig. 5A) corresponds to a large distributary channel observed in the DEM (Fig. 2 location 4). The channel is cross-cutting some of the distributary channel network on the delta top implying that it was active after the topsets were deposited (Fig. 2). Three similar channels are observed from delta morphology in the study area. They are 900 to 1800 m long, 50 to 250 m wide and incise delta tops between 3 m and 11 m (Fig. 2 -Hyrtiälänkangas and Kukonkangas) topsets either as a response to the Baltic ice lake base level drop or as a consequence of an outburst flood (Brookfield & Martini, 1999;Gomez et al., 2000;Winsemann et al., 2009Winsemann et al., , 2018Gobo et al., 2015;Lang et al., 2021). The channel most likely developed following ice-marginal retreat from Salpausselkä I when an ice-contact lake that developed between the ice sheet and the proximal delta slope of Salpausselkä I overspilled and drained southward into the Baltic ice lake (Fig. 2). High energy of the flow is reflected in sediment grain size within FA 7 which ranges from cobble to boulder-sized gravels and massive or poorly bedded structure of the deposits pointing towards a high-magnitude event ). An outsized clast of +60 cm in diameter is clearly visible in the centre of the section (Fig. 6C). It appears to be disturbing the layers in its vicinity, but evidence for impact structures is absent (e.g. Domack & Lawson, 1985;Bennett et al., 1996). The succession itself is wellexposed in the outcrop (Fig. 7A). In the outcrop this facies association can be observed to form isolated, irregular and laterally thinning fine sediment bodies, onlapping onto the host succession (FA 6, FA 2). Smallscale normal faulting can be observed FA 8 is interpreted to be a kettle-hole lake fill succession that developed during and after a buried ice block was blanketed by delta topset sediments and, subsequently, melted out (e.g. Maizels, 1977;Wingfield, 1990;Benn & Evans, 2010;Götz et al., 2018). Melting of the ice block created localized accommodation, probably a small lake, which was infilled with finegrained sediment. The outsized clast is interpreted as a 'slide-stone', which travelled down along the slope as the walls of the kettle hole (Fig. 7C). Ice melting, leading to the loss of lateral support, is deemed responsible for the presence of small-scale normal faults ( Fig. 7B). Similar normal faults are reported to form due to the reactivation of basement faults but given the morphological context it is unlikely in this case (e.g. Brandes et al., 2011). This association is remarkably similar to the low energy fan/delta bottomsets association (FA 5) making it difficult to interpret it without the context of neighbouring sediments (FA 6) ( Fig. 7A) and the geomorphology revealed by the highresolution DEM   Ice-contact deltas 3075 cessation (fine sands and silts) (Fig. 5B, log 13, and Fig. 12) (Marren, 2005).

Delta foreset
Wedge-shaped, steeply dipping well-sorted, coarse sands and pebble gravels of FA 4 constitute the majority of the delta foresets (Table 2,  In GPR profiles delta foresets are characterized by large, steeply dipping reflections with variable amplitudes (Figs 10 and 11). Some reflections are continuous throughout, and show clear offlap from, delta topsets whereas others are either truncated or show a toplap relationship ( Fig. 10A and C). Variations in reflector architecture are attributed to the changing base level of the Baltic ice lake (Backert et al., 2010;Winsemann et al., 2011Winsemann et al., , 2018. Low reflectivity zones within the delta foresets are interpreted as uniform sediment composition (i.e. grain size) generated under stable discharge conditions ( Fig. 10A and B). Higher reflectivity represents an increase in electrical permittivity contrast, probably resulting from variations in grain size due to flow cessation (lobe switching) and deposition of finer grained sediment drapes from suspension or, discharge increase and deposition of coarser lithologies (boulders or cobbles). In GPR profiles perpendicular to the depositional dip, delta foresets appear as semi-continuous, sub-parallel, sinusoidal reflectors (Figs 10C and 11) prograding at up to ca 90°from the main depositional dip direction (i.e. towards the west or east) (Figs 10C and 11). Variations in the progradation direction often occurred in the vicinity of pronounced bedrock highs, where meltwater flow was locally deflected. Foreset reflections pass distally into more shallow dipping and higher amplitude packages interpreted as delta bottomsets (cf. Delta bottomsets section). In GPR profiles delta bottomsets are identified as downdip continuations of foresets exhibiting a clear change to a shallower dip angle and increase in amplitude (Fig. 10). Reflections are semi-continuous and often onlap onto bedrock topographic highs (Figs 10A, 10C, 10D and 11). Delta bottomsets may pass into, or be Fig. 6. Details of sediments found in outcrops. Log locations indicated by coloured bars and corresponding circled number. (A) and (B) Contorted and fine sediments in the distal delta bottomsets; (C) discordant contact between sands (bottom) and finer sediments above the suspected slump scar in delta foresets, note cut-and-fill structures as well as cross-bedding in the direction opposite to the delta slope possibly indicating backsets cyclic steps and/or upslope-migrating antidunes; (D) interbedded sands with current ripples preserved and gravels with weak imbrication. Note alternating current direction indicators which may point to deposition from 3D dunes or cyclic steps; (E) outcrop of the medial (along the depositional dip) part of the delta. Note large-scale cross-bedding unconformably overlain by massive and crudely stratified cobbles and pebble gravels of a large distributary channel; (F) size of some boulders transported within the channel was in excess of 80 cm in diameter.

Delta bottomsets
interbedded with, prodelta deposits (cf. Prodelta/lake bottom sediment section). They are observed in GPR data as higher amplitude, more continuous reflections attributed to higher water content in the fine-grained sediments (Fig. 10A  and D). The thickness of delta bottomsets is controlled by the underlying bedrock topography and can be extremely variable from 0 m in areas where bedrock crops out to over 60 m in troughs, as confirmed by boreholes (Figs 11 and  12).

Delta topsets
Delta topsets can be observed in the uppermost part of some outcrops (1 to 6 m) as a coarse, bedded package of pebble and cobble gravels with subordinate cobble gravels and sands which truncate the foresets (Figs 5, 7A, 7F, 8F, 9 and 12). The topsets were deposited subaerially by braided glaciofluvial meltwater streams originating at the ice margin (e.g. Fyfe, 1990;Postma, 1990a,b;Thomas & Chiverrell, 2006;Rohais et al., 2008). Limited accommodation and periodic high discharges led to the erosional character of FA 6 (Table 2, Fig. 9). Cobble and boulder gravels were transported mainly as bedload whereas most of the sands and finer fraction bypassed the delta topsets or was deposited during the waning flow stage, and then winnowed by subsequent higher meltwater discharges (Burke et al., 2010;Pisarska-Jamrozy & Zieliński, 2014). Base level fluctuations of the Baltic ice lake led to incision and reworking of some of the topsets (FA 6) and foresets (FA 4) when the lake level fell, versus aggradation of topset sediments, when the lake level rose (Fig. 12)  .
In GPR profiles topsets are characterized by a semi-continuous or continuous, high-amplitude, horizontal reflection package (Fig. 9). The package usually spans the whole length of the delta and terminates at the delta brink. Some topset reflections are observed to offlap and pass into delta foresets showing either a progradation or an aggradation and progradation stacking pattern when the delta brink was responding to rising lake levels ( Fig. 9A and D) (Brookfield & Martini, 1999;Gobo et al., 2015;Winsemann et al., 2018).

Large distributary channel
Several large distributary channels attributed to high magnitude meltwater discharge events are visible in the delta morphology (Figs 2, 10B and 12). In outcrop they are represented by a very coarse, laterally tapering package of crudely bedded boulders and cobbles with multiple reactivation surfaces (Figs 5,6E and 6F).
Several such channels >50 m wide were observed in GPR profiles (Figs 9 and 10C; Appendix S1). They are characterized by concave, mid to high-amplitude basal reflections filled with semi-continuous reflections some of which onlap onto the basal reflector (Fig. 9). Secondary, smaller concave reflectors are present within the channel fill facies and can be correlated to internal erosional surfaces present in the outcrop (Fig. 9). One of the underfilled channels, visible in the delta morphology, can be linked to a semitransparent foreset package interpreted as a product of a large meltwater drainage event (possibly a lake outburst flood) or an incision due to lake level fall (Fig. 10B). The semi-transparent reflection is interpreted as deposited from a sustained flow associated with high-magnitude drainage events resulting in the deposition of a uniform sedimentary package (Clague & Evans, 2000;Winsemann et al., 2011;Westoby et al., 2014).

Kettle holes
Oval or irregular kettle holes cover up to 20% of the delta surface (Fig. 2) (Girard et al., 2015;Maries et al., 2017). Such depressions are often fully or partially filled with fine sediments of FA 8 (Table 2, Figs 5,9 and 12). The sediments are remarkably similar to the delta bottomsets (Fig. 9) and it would be extremely hard to identify a kettle hole fill without the context of neighbouring facies (FA 6 and FA 7) (Fig. 8A, B, C and D) and the geomorphology revealed by the highresolution DEM (Fig. 2). If encountered in well or borehole, without the geomorphological context, Fig. 7. (A) Contact between kettle hole sediments (log 6) and topset facies (log 5). Note the exposed bedrock high. the succession could be misinterpreted as more laterally extensive muds associated with flooding of the delta. (Figs 5 and 12).
In GPR data possible filled kettle holes are identified as high-amplitude convex packages with reflections 'filling' the initial depression (Figs 9 and 10C; Appendix S1). As per bottomset sediments, the elevated amplitude response is linked to increased water retention within the fine-grained sediments.

Esker
Deposits associated with eskers have not been encountered in outcrop and are identified based on their morphological expression linked to the GPR profiles where they exhibit a broad range of responses from medium-amplitude semi-chaotic reflection patterns to almost transparent with some discontinuous, steeply dipping, reflections (Figs 9 and 10D; Appendix S1). Typically, sediments forming eskers range from boulders and cobbles (high discharge) to sands and pea gravels (low/waning discharge) (Fig. 12) (Burke et al., 2015;Knight, 2019). An esker may have a complex geometry due to filling only a part of the subglacial conduit or fully filling a complex subglacial, englacial or supraglacial, or combination thereof, drainage network ( Fig. 12C and D) (Artimo et al., 2003;Mäkinen & Palmu, 2008;Gruszka et al., 2012;Storrar et al., 2019). Higher amplitude reflection packages delineating the top of the esker which are onlapped and or draped by subaqueous ice-contact fan/delta foresets are common diagnostic features (Fig. 12). The higher amplitude is attributed to deposition, from suspension, of overlying fine fractions on top of a submerged esker after the ice margin had retreated. It may not be possible to identify buried eskers/ice marginal sediments solely from the GPR profiles as they are difficult to reliably differentiate from bedrock highs and have similar dielectric permittivity values to neighbouring sediments (Fig. 11).
Prodelta/lake bottom sediment Prodelta sediments were also only interpreted from GPR profiles (Figs 9 and 10A). They are a distal continuation of the deltaic bottomsets ( Fig. 10A; Appendix S1) and are characterized by continuous, subparallel, high-amplitude reflections (Fig. 9). Baltic ice lake fine-grained sediments are known from other localities, to have high silt/clay content, which is consistent with their high reflectivity in the GPR data (Figs 10A and 11) (Donner, 2010;Hyttinen et al., 2011;Zolitschka et al., 2015). Thicker, high reflectivity packages occur in bedrock depressions and in areas down-flow from bedrock knolls ( Fig. 10A and C). In GPR profile 16 (Fig. 11) a higher amplitude layer of silt/clay draping part of the delta is confirmed by wells. It is interpreted to represent deposition from suspension over an inactive, submerged, part of the delta during a period when meltwaters were diverted to the west, possibly due to an ice margin re-advance over the proximal part of the delta (Fig. 2).

Summary of sedimentary environments
From outcrop and GPR data it is apparent that facies distributions within the delta are more complex than typical Gilbert-type deltas. Rapid lateral and proximal-distal facies variations are linked to supercritical/subcritical flow transitions which are characteristic of ice marginal and subaqueous ice-contact fan/delta systems (e.g. Russell & Arnott, 2003;Hornung et al., 2007;Winsemann et al., 2009Winsemann et al., , 2018Dietrich et al., 2016;Ghienne et al., 2021;Lang et al., 2021). Correlation of logs, even at the outcrop scale, is challenging (Figs 5, 7F, 8A and 12). GPR data suggest that subaqueous ice-contact fan deposits extend 200 to 400 m (maximum 600 m) from the former ice margin (Figs 9A and 12; Appendix S1). This distance is a function of the energy of the meltwater jet, angle of repose of the ice-contact fan sediments and the water depth of the Baltic ice lake  (Hoyal et al., 2003;Winsemann et al., 2009;Gobo et al., 2015;Lang et al., 2017Lang et al., , 2021. In ice-contact settings delta and fan facies are genetically linked and exist on a spectrum, i.e. when a subaqueous fan reaches the lake surface or the lake level falls, deposition of a Gilbert-type delta commences. The delta then progrades into the basin until the sediment supply seizes or the lake level rises again above the delta topsets, in which case fan sedimentation takes over and may ultimately overtop the delta (Lønne, 1995;Thomas & Chiverrell, 2006). Deflection of delta foresets around bedrock highs and ponding/preferential settling of finer fractions (bottomsets and prodelta sediments) in bedrock lows can be identified from GPR profiles (Figs 10C and 11). In some cases, older, pre-existing sediments associated with previous ice-margin positions and/or eskers can be buried by delta progradation as observed in GPR data (Figs 9D and 12). Two alternatives, and equally likely, scenarios for the formation of a buried sediment body are considered here (Fig. 12C and D). Irrespective of their genesis, such sediments can be draped with prodelta muds during the time between ice-margin retreat and the progradation of the delta foresets over them (Figs 11 and 12). Subglacial diamicton (traction till) is usually present on the proximal slope and on the delta top, if ice advanced over it, but it may be patchy or absent, if the ice margin was less active (Fig. 12). Underfilled kettle holes are more abundant in the ice-proximal part of the delta (Figs 2 and 12). When kettle holes are filled, typically with finegrained sediments, they introduce further heterogeneity to the deltaic package (for example, Figs 7 A, 7B, 7C and 12). Given the scale of the system, the above detail is unlikely to be imaged in the subsurface by standard seismic exploration methods (e.g. Barrett et al., 2020).

ELECTRICAL RESISTIVITY TOMOGRAPHY
The ERT profiles penetrate the sedimentary succession to depths between 30 m and 70 m ( Fig. 13; Appendix S1). Four distinctive resistivity zones are identified (Fig. 14). The highest resistivity values were recorded in the upper 10 to 20 m of zone 1 (Z1) with values between 12 000 to 120 000 Ωm (mean: 49 700 Ωm, four iterations). In profiles 1 and 3, towards the proximal part of the delta, an intermediate resistivity zone Z2 (8200 to 39 000 Ωm, mean: 14 200) is present. In more distal areas the upper, high resistivity zone makes a sharp contact with either the low resistivity zone (Z3) (300 to 1550 Ωm, mean: 670 Ωm) or with zone Z4, which displays a wider range of resistivity responses (200 to 20 000 Ωm, mean 10 050 Ωm) (Fig. 14).
Outcrops and borehole profiles indicate that clay content is low. This is typical for well-sorted glaciofluvial and deltaic sediments (Glückert, 1986;Palmu, 1999;Zielinski & Loon, 2003;Winsemann et al., 2018;Kurjanski et al., 2020). Based on the geomorphology of the deltas, together with the available outcrops and GPR data, Z1 is interpreted as well-sorted, unsaturated (dry) sands and gravels (Fig. 13). Of note are the extremely high resistivity values, exceeding typical values for glacial sand and gravels (e.g. Palacky, 1988) but consistent within values detected from glacial sediments deposited by the Laurentide ice sheet (Culley et al., 1976) and the Fennoscandian ice sheet in Finland (Sutinen, 1992). Such values are attributed to air-filled pore spaces within the sediments suggesting high permeability and rapid drainage above the water table.
The intermediate zone (Z2) is interpreted as either ice-proximal facies or a heterolithic sandy-silty package deposited in the distal part of the delta (Fig. 13). In comparison to Z1, the lower resistivity values most likely reflect poorer sorting resulting in higher water saturation and locally, lower porosity.
The lowest resistivity values are recorded in the saturated zone Z3, located beneath Z1, but are relatively high when compared to typical freshwater aquifers found in unconsolidated sandy sediments (e.g. Choudhury et al., 2001;Batayneh, 2006;Galazoulas et al., 2015). However, mineralization in water samples shows electrical conductivity values as low as 0.001 S m −1 (mean 0.016 S m −1 ) (Table 3), which is consistent with meteoric water recharging the aquifer and short water retention times. The contrast between Z1 and Z3 is three or four orders of magnitude. In such areas Fig. 9. Architectural elements of ice-contact deltas identified from outcrops and juxtaposed against GPR reflection patterns observed from GPR profiles. Reservoir properties are described qualitatively based on visual inspection of outcrops, and field assessment of grain size, sorting and content of fine fractions. Where no outcrops were available, reservoir properties are assessed based on available literature: (Zolitschka et al., 2015;Maries et al., 2017). the top of the saturated zone, rather than top of the bedrock, is likely to be imaged. It is important to note that average resistivity values for Z3 vary between ERT profiles. In profiles ERT1 and ERT2 the average resistivity is between 1100 to 1500 Ωm (minimum 675 Ωm) whereas in profile ERT3 the average is 801 Ωm (minimum 306 Ωm) (Fig. 14). None of the boreholes along these profiles encountered the water table. Bedrock elevation derived from gravity and borehole data located in the vicinity of the profile ERT 3 data suggests that an aquifer is only locally developed in hollows where the bedrock is over-deepened. This suggests that water saturation is below 100% in most of the zone Z3 as infiltrating meteoric waters are in transit towards the underlying bedrock. Zone Z4 is characterized by a more variable resistivity response and interpreted as bedrock, based on the correlation to exposures (Fig. 7A and E). Table 3 were used to calculate porosity within each zone (Z1 to Z3) (Table 3) (Fig. 14). Only areas above the depth of investigation (DOI) limit were considered. The porosity modelled for Fig. 10. Ground penetrating radar profiles (GPR) with interpretations. For locations see Fig. 3. For GPR profiles 1, 2, 3, 4, 6, 7, 10, 11, 12 and 13, see Appendix S1. Average relative dielectric permittivity (RDP) value of 6.2 was used for time-depth conversion. zone Z1 is between 20% and 30% with extreme values bordering 40% (high) and 15% (low). Such a range is comparable to values reported from unconsolidated glaciofluvial sediments elsewhere (e.g. Olanrewaju & Wong, 1984;Salem, 2001a). It was assumed, based on the interpretation of the GPR profiles and sedimentary logs, that the difference in resistivity between zones Z1 and Z3 is mainly controlled by changes in water saturation rather than porosity. The only expected change is a decrease towards the bottom of the succession as a consequence of the increased fine fraction in the bottomset beds. Assuming the water saturation controls provided constraints for modelling the porosity of zone Z1 and Z3. Modelling of Z3 assumed 95% water saturation to account for a fully saturated aquifer and a transition zone, where saturation increases from close to 0% to 100%. Modelling returned porosity values between 10% and 20% for profiles ERT 1 and ERT 2, whereas in profile ERT 3 porosity of the saturated zone (Z3) ranges between 12.5% and 35%, with most of the values falling between 15% and 25% (Fig. 14). As above, an aquifer is likely to be present only locally in ERT 3, which implies that the average saturation of Z3 is below 95% and likely to be much lower in ERT 1 and ERT 2 (Figs 13 and 14; Appendix S1). Porosity values in ERT 1 and ERT 2 are likely underestimated and only in ERT 3 is the calculated porosity closer to satisfying the conditions outlined above. Porosity falls between 25% and 40% in Z2 with unrealistically high (>48%) and low values (<10%) (Fig. 14). This is attributed to the poorer sorting and larger grain-size variations which are responsible for the broader resistivity range. The unrealistically high porosity values obtained in Z2, >30 to 35%, can be explained by a number of possible (and additive) errors: (i) underestimated resistivity values in Z2 which could be attributed to the inversion producing smooth transitions rather than a well-defined geological unit; (ii) underestimated values for water saturation; and (iii) the presence of clays in Z2, which would result in Archie's model largely overestimating porosity. This is supported by the presence of a fine fraction in the Z2 sediment package, which would favour a higher water saturation due to poorer sorting as seen in outcrop (Figs 6A, 6B, 8E, 8H, 8I and 14).

Depth to bedrock
Estimation for the thickness of the sedimentary sequence was achieved by combining results from gravity measurements and water borehole logs, provided by the Geological Survey of Finland, together with outcrop, GPR and ERT data ( Fig. 13; Appendix S1). The interpretation was performed using a set of rules, outlined in Fig. S2. Areas where at least two, and preferably three, geophysical methods overlap were interpreted with the highest confidence. Uncertainty increased away from boreholes and in places where geophysical methods yield diverging results. Since the gravity profiles were not acquired along accessible routes (roads and paths) it was impossible to acquire co-linear ERT data, which precluded the possibility to invert the ERT and gravity data together (Fig. 3). The depth to bedrock, based on the gravity survey, was used to interpolate the approximate depth to bedrock along all three ERT profiles ( Fig. 13; Appendix S1). The only relatively robust gravity measurements along the ERT profiles are where gravity profiles intersect the resistivity lines (Figs 3 and 13).
The thickness of the delta sediments along the ERT and gravity profiles is highly variable. Less than 5 m of glaciofluvial sediments are present over bedrock highs (Figs 7E,10C and 13). In the central, flat-topped part of the delta, the thickness usually varies between 10 to 45 m, except for deep bedrock troughs, which most likely followed structural lineaments, where borehole data show sediment thicknesses in excess of 90 m (Figs 11 and 13). The geophysical surveys undertaken here and boreholes in the vicinity of ERT 3 indicate a sediment thickness of ca 30 m (Fig. 13).

Integration of geomorphology, sedimentology and geophysics
Sedimentological studies of Quaternary glacial deposits can be supplemented with geophysical methods and geomorphological mapping. Together these offer the advantage of linking 2D exposures and profiles to a 3D geomorphic expression of glacial landforms, which is rarely possible for ancient depositional systems. An iterative, down-stepping approach in terms of penetration depth and resolution, using more than one geophysical method, was found to be most effective when characterizing the icecontact Salpausselkä deltas (Figs 13 and 15). Outcrops were logged and their position within the depositional system was identified using the morphological features observed from the LiDaR DEM (Fig. 2). GPR data was interpreted by linking reflection patterns to corresponding sediments observed in scarce exposures (Fig. 12). ERT and gravity methods were linked/overlapped with structures and sediments identified in the top 10 to 20 m from the GPR and exposures, which allowed extrapolation below the limit of GPR penetration (Figs 3, 11, 13 and 15). This approach provides the most robust interpretation of the sedimentary architecture of the Salpausselkä deltas.

Comparison of geophysical methods
The GPR is invaluable for understanding the subsurface sedimentary architecture and the depth to bedrock. Its utility is limited by the maximum penetration depth, which, in unsaturated sediments, is ca 20 m for a bistatic 80 MHz GPR antenna configuration (Figs 10 and 11). The penetration depth is reduced significantly when the sediments are saturated, and/or clays/silts are present (Fig. 10). Comparison of the ERT and GPR data show that mapping depth to bedrock, based on the loss of GPR signal, leads to an underestimation of the profile crosssectional area by ca 50%, with concomitant effects for calculating sediment volumes (Fig. 15). In bedrock lows, especially in the saturated zone, ERT exhibits a large difference from the gravity data (Fig. 13). This is most likely a result of the ERT picking the top of the aquifer/ saturated zone rather than the top of the bedrock (Fig. 15). This again leads to underestimation of the sediment volumes. An opposite relationship between ERT and gravity (ERT depth to bedrock  Gobo et al. (2015) and Lang et al. (2017). Note the large contrast in grain size, especially in the proximal part, and the role of bedrock and pre-existing sedimentary structures on facies distribution and potential fluid migration pathways. (C) and (D) Two conceptual models for the depositional history of Vesivhmaankangas (study area 2) based on interpretation of GRP and ERT profiles. For detailed sedimentary logs see Fig. 5, for their location see Fig. 2. Fig. 13. Electrical resistivity tomography profiles: (A) profile ERT 3; (B) interpretation and corresponding GPR reflection pattern (Fig. 10D). For profiles ERT 1 and ERT 2 see Appendix S1. For location see Fig. 3.
Fig. 14. Porosity distribution from ERT profiles. Each zone has been calculated using the parameters in Table 3: (A) profile ERT 1; (B) profile ERT 2; (C) profile ERT 3; and (D) frequency distribution of porosity and resistivity within each zone and profile.
< gravity depth to bedrock) can be seen mainly in ice-proximal settings (Fig. 13, Table S3). Borehole GTK111, which directly penetrates the section, shows a close correlation with the ERT result (Fig. 13). Gravity data, showing a shallower depth to bedrock may be attributed to the density of the sediments, used as a parameter during gravity inversion. Proximal, icecontact facies contain more coarse fractions (cobble and boulder) which may have an  (Archie, 1942;Salem, 2001b;Glover, 2016;Byun et al., 2019). † Calculated through Monte Carlo simulation using average resistivity from each ERT profile and boundary condition: porosity Z1~= porosity Z3 and matching mean resistivity value from all the profiles. § Best fit value for 'dry' sediments based on average resistivity within each zone and information from (Glover, 2017). ¶ Standard value for water saturated sands and sandstones (e.g. Archie, 1942;Glover, 2017). ** Calculated based on data from (Olanrewaju & Wong, 1984;Salem, 2001b;Rose et al., 2016). intermediate (as a bulk) density between bedrock (high) and delta sediments (low). In these situations it may be beneficial to implement a four-layer model (dry delta sediments, proximal fan sediments, water saturated sediments and bedrock) instead of the three-layer model used at present (dry sediments, saturated sediments and bedrock) to improve the gravity inversion.

Sedimentary facies distribution and their influence on reservoir properties
Gilbert-type deltas contain significant porous sedimentary bodies that are potential reservoirs for hydrocarbons, repositories for CO 2, sources of low enthalpy geothermal energy and aquifers (Knutsson, 2008;Arola et al., 2014). Their sedimentology has been widely described and they are comprised of predictable facies successions (e.g. Gilbert, 1885;Mastalerz, 1990;Nemec & Postma, 1993;Lønne & Nemec, 2004;Rohais et al., 2008;Colella & Prior, 2009;Gobo et al., 2015;Chavarrías et al., 2018;Winsemann et al., 2018). The Salpausselkä I and II deltas are reliable water reservoirs (Figs 12 and 16) (Fyfe, 1990;Palmu, 1999;Maries et al., 2017;Virtasalo et al., 2019) and the porosity modelling in the current study has suggested that they have porosities within the 20 to 35% range, similar to values reported for other glaciofluvial deposits (Figs 14 and 16) (Olanrewaju & Wong, 1984;Salem, 2001b). Porosity appears to decrease towards the base of the delta (10 to 25% -Figs 14 and 16) due to the smaller grain size of the delta bottomsets/prodelta sediments. The presence of finer grained and/or poorly sorted slump deposits may decrease the reservoir quality further ( Fig. 6A and B). However, the reservoir performance of icecontact deltas can be negatively affected by meso-scale heterogeneities which are directly linked to the depositional processes taking place at the ice margin (e.g. Mitten et al., 2020). Evidence for ice overriding, bulldozing and mixing by ice-margin oscillations are clearly visible in delta morphology (Fig. 2). Eskers and ice marginal ridges observed in the topography on the trend of the deltas indicate that some of the older landforms might have been buried beneath them (Figs 2 and 12). Understanding of such complexities is crucial, especially for deeply buried or low-pressure reservoirs where relatively small changes in porosity/permeability properties can hugely affect the fluid flow.
In ice-proximal environments most of the deposition is associated with highly dynamic point-sourced discharge from subglacial and englacial conduits (Fig. 16). This leads to large, sharp and instantaneous (in geological terms) changes in sediment grain size (boulder gravel to very fine silt) occurring over short distances in all dimensions (x, y, z) and is reflected in the modelled porosity values, ranging from 10% to over 40% in in subaqueous ice-proximal fan facies (Fig. 16) (Rust & Romanelli, 1975;Powell, 1990;Russell & Arnott, 2003). High porosity contrasts can lead to generation of 'thief zones' which are important when a two-phase flow (gas-liquid/water-oil) is considered, since they may lead to early water breakthrough in wells and isolation of gas pockets, effectively making production impossible (Li et al., 2016;Satter & Iqbal, 2016). Glaciotectonic modification due to ice-margin oscillations may effectively decrease sorting of the ice-contact slope and therefore reduce porosity/permeability creating a zone of decreased reservoir potential (Fig. 16). Most subglacial diamictons (traction tills) are poorly sorted (mechanically mixed) with a fine-grained matrix, including substantial amounts of clay/ clay sized material, resulting in low porosity and permeability, making them reasonably good seals (Benn & Evans, 1996;Evans et al., 2006). This is only partially true for the study area since the tills in the region have, dominantly, a silty/fine sand matrix as their provenance is crystalline bedrock. The till is sandy and gravelrich, locally clast-supported as a result of the cannibalization of glaciofluvial deltaic sediments ( Fig. 8H and I). This part of the delta is a waste zone (a lower porosity-permeability, nonproducing zone within a reservoir) rather than a seal capable of retaining a build-up of fluid pressure (Fig. 16).
Typically, deltas prograde into the basin with younger sediments deposited distally away from the efflux point (portal), resulting in proximaldistal younging. This rule does not always hold within ice-contact deltas because, in some cases, backstepping of the ice margin creates accommodation (lake) between the ice margin and previously deposited fan/delta deposits, moraines or eskers (Figs 10D, 12C, 12D and 16) (Ashley, 2002). This space may then be partially, or fully, infilled with younger delta deposits. The older delta may be then re-activated and overtopped, following which progradation further into the basin is re-established (Figs 12 and 16). A lacustrine mud drape covering older deltaic deposits is likely to be deposited in the time between the ice retreat and renewed delta progradation. When such a lake becomes overfilled and/or a significant lake level fall occurs creating a large gradient between two water bodies a large meltwater drainage event may take place. This is observed from the Hyrtiälänkangas Delta, where two large erosional meltwater channels incise the delta topset and feed a small forced regressive delta lobe which can be identified because it lacks the distinctive internal reflectivity in GPR data, that is typical for bottomsets/distal delta settings elsewhere (Figs 2, 10B, 12 and 16). Geomorphology of the delta shows that the drainage channel incised the ice-proximal part of the delta whereas little or no incision occurred in the ice-distal slope. This suggests that at this location fill-and-spill model for lake drainage is more likely than the Baltic ice level fall scenario (Fig. 10A). Such erosional reworking of the delta topsets may deposit good reservoir quality sediments in bottomset or prodelta areas, which typically have moderate/poor reservoir properties (Fig. 16).
Kettle holes are clearly visible on the delta surface ( Fig. 2) (Okko, 1962;Donner, 1976Donner, , 1995Fleisher, 1986;Wingfield, 1990). When infilled and buried, kettle holes may introduce additional heterogeneities into the reservoir. The sediments filling the depressions are likely to be finer grained and of moderate/poor reservoir quality in comparison to surrounding topset, fan or foreset facies (good reservoirs) (Fig. 16). Kettle holes effectively decrease the total volume of the reservoir, their interpretation from well data (for example, based on palynology) may erroneously indicate a thinner reservoir package or a sandy transgressive package deposited on top of the deltas.
The dynamic nature of ice-contact deposition is likely to juxtapose sediments, which would not be in contact in non-glacial deltaic settings (for example, mud drape, esker sediments buried within delta bottomsets, subglacial diamicton above delta topsets and kettle holes) which introduces large contrasts in porosity and permeability to the reservoir (Fig. 16). This can be magnified when the deltas are buried or when reservoir pressures are relatively low, such as in the glaciogenic Pleistocene sediments of the North Sea (Le Heron et al., 2006;Hirst, 2012;Martin et al., 2012;Ottesen et al., 2012;Rose et al., 2016;Kurjanski et al., 2020).
The underlying bedrock topography strongly influences the sediment distribution and reservoir properties of ice-contact deltas by deflecting or focusing meltwater flow which affects facies distribution (Fig. 16). The topography may also subdivide the reservoir, creating disconnected aquifers with independent water tables (buried lakes), or lengthen water migration and retention time reflecting the highly tortuous infiltration path taken to bypass bedrock highs (Figs 15 and  16). This is especially important for reservoir management, and heat transfer in geothermal systems but also for evaluation of deeper reservoirs where the bedrock may not be accurately imaged, leading to errors in volume calculations and overestimations of connectivity.

CONCLUSIONS
In this study the ice-contact deltas of Salpausselka I and II, associated with the Younger Dryas Stadial stillstand of the Fennoscandian ice sheet, were investigated using a combination of geomorphology, sedimentology and geophysical surveys (ERT, GPR and gravity) to provide a full-spectrum evaluation of the ice-contact delta internal structure and reservoir properties and reservoir potential.
High and fluctuating discharge energy from meltwater portals was responsible for the deposition of subaqueous ice-contact fans under highly supercritical flow conditions depositing boulder to cobble sized clasts proximal to the ice margin with large cyclic steps and hydraulic jump deposits downflow. Ice-proximal deposits exhibit the highest grain-size variation distally and laterally. When the fans aggraded to the lake surface, progradation of steep (10 to 35˚) delta foresets commenced, often leading to the coalescence of two or more neighbouring deltas. Deposition on such steep slopes was likely dominated by supercritical flow conditions (formation of cyclic steps and antidunes) which eventually decelerated and deposited finer fractions in the subcritical flow regime as bottomsets. The foresets were capped by subaerially deposited sandur-like coarse glaciofluvial topsets. Short-term fluctuations of the Baltic ice lake level resulted in incision and re-mobilization of parts of the foreset package (lake level fall) supplying coarse sediments to the distal part of the delta or flooding parts of the delta (lake level rise) changing the foreset stacking pattern from progradation to progradation/aggradation.
Overall, ice-contact deltas are excellent reservoirs (porosity 10 to 42%) but, except for the largest of deltas (Figs 2B and 4), they are rarely big enough to constitute a stand-alone exploration target. Alternatively, exploration of multiple, clustered deltas/fans may be economically viable. Their recognition in the subsurface largely depends on a good understanding of the depositional environment and/or illumination of permeable deposits by gas. Their structure is likely to be more complex than that of a typical Gilbert-type delta (Figs 12 and 16). This mainly results from ice-marginal oscillations affecting sedimentation in the ice-proximal part of the delta and the typical temporally varying (diurnal to annual) magnitude of meltwater and sediment flux. Moreover, backstepping of the ice margin (almost instantaneous in geological terms) rearranges sediment input points and provides accommodation space 'behind' or laterally from a previously deposited delta, ice-contact fan or esker (Figs 2 and 12) (Ashley, 2002). This can result in complex, amalgamated deltas with 'hidden' heterogeneities buried within, and abrupt lateral and distal facies changes affecting overall reservoir potential (Figs 2 and 16). The localized presence of highly permeable zones or layers ('thief zones') can lead to preferential fluid flow through those zones and bypass of lower porosity/permeability and early water breakthrough in hydrocarbon producing wells. This may have a significant effect when gas production is considered, especially from shallow gas targets, where formation pressures are low, or in deeper targets, where compaction and cementation may amplify porosity-permeability heterogeneities. Finally, this study shows that sediment distribution is significantly affected by the subglaciallysculpted bedrock topography (Fig. 16). These deltas were deposited over an extremely variable bedrock topography with sediment thickness varying between zero metres, where bedrock highs outcrop at the delta surface, to over 90 m in bedrock troughs. Such an interplay between topography and deposition will likely reroute sediments leading to the development of a complex facies pattern.
Results of this work are relevant for both the aquifer and low enthalpy geothermal energy exploitation in the region as well as similar deposits in the subsurface associated with Pleistocene and pre-Pleistocene glacial periods, where commercially important reserves of water or hydrocarbons may exist. work becomes much more enjoyable when there are good people around to share the burden. We are also grateful for the help and support provided by the team at the Geological Survey of Finland Espoo Office without whom this work would have not been possible. Authors are very grateful for the constructive comments and suggestions offered by two anonymous reviewers which allowed us to greatly improve the manuscript. Least but not last, we would like to thank the editorial board of Sedimentology, Dr Victoria Valdez, Dr Ian Kane and Elaine Richardson for their prompt and efficient handling of the manuscript in times when none of the above should be taken for granted.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. Zielinski, T. and Van Loon, A.J. (2003)

Supporting Information
Additional information may be found in the online version of this article: Figure S1. Sedimentary facies identified in studied sections. See Table 1 for a detailed description. Figure S2. Radar facies identified from GPR profiles. Figure S3. GPR profiles GPR1 GPR2, GPR3, GPR4, GPR6, GPR7, GPR10, GPR11, GPR12, GPR13 with interpretations. Figure S4. Electrical resistivity tomography profile ERT1. Figure S5. Electrical resistivity tomography profile ERT 2. Table S1. Details of geophysical methods. Table S2. Criteria for geophysical data interpretation.