Geophysical reconstruction of the late Holocene proximal proglacial landsystem at Skeiðarársandur, southeast Iceland

Sandur plains are extensive sedimentary bodies formed in proglacial settings by glaciofluvial processes. Icelandic sandar have been hypothesised to be comprised of thick alluvial successions that provide critical insight into the processes that contributed to their formation and evolution. However, to‐date, most sandar research has focused on the analysis of sedimentary successions associated to topographically‐confined and small‐scale systems. These, however, do not capture the variety or scale of processes that influence sandar architecture. Therefore, detailed subsurface analysis of large‐scale and unconfined sandar is vital to understand how these systems respond to fundamental drivers, such as: (i) glacier oscillations, and (ii) episodic sediment flux from glacier outburst floods (aka. jökulhlaups). We report an extensive, low‐frequency (40 & 100 MHz) ground‐penetrating radar (GPR) survey of the ice‐proximal component of a large‐scale ( ~ 1300 km 2 $\unicode{x0007E}1300\,{\text{km}}^{2}$ ) active sandur in southeast Iceland. A bright and continuous reflection (PR1) is identified within all radargrams and provides a boundary between pre‐LIA and LIA to present‐day sedimentation. GPR data reveals a ~50 m thick ice‐proximal sediment wedge that is attributed to jökulhlaup and surge‐related glaciofluvial activity during the Little Ice Age (LIA). An approximate rate of deposition of 0.2–0.65 m a− 1 ${{\rm{a}}}^{-1}$ has been calculated for the sediment wedge above PR1. Furthermore, we propose an extensive, sandur‐wide, pre‐LIA ice marginal limit of Skeiðarárökull, southeast Iceland, based on observations reported here and in previous work.


Introduction
Sandur plains (plurally: sandar) are extensive sedimentary bodies formed in proglacial settings by sediment deposition from meltwater systems emanating from the margin of a glacier or ice-sheet. Sandar exist in both topographically confined (e.g. valley sandur) and unconfined (e.g. coastal sandur) proglacial settings (e.g. Maizels, 1979). Confined sandur plains are dominated by coarse-grained glaciofluvial sediment deposited by a single braided fluvial system originating from the ice-margin (Krigström, 1962;Bluck, 1974). These systems often display limited variation in their down-stream fining pattern alongside a high degree of lateral variability (Bluck, 1974;Aitken, 1998). Sandar developed within topographically confined settings, such as valley sandur, can lead to the development of thick successions of sediment infill that provide a record of deglacial sedimentation and high magnitude glaciofluvial events (e.g. Bluck, 1974;Haraldsson and Palm, 1980;Haraldsson, 1981;Carrivick et al., 2004a, b;Russell, 2007Russell, , 2009aStorms et al., 2012). Where glacier outburst floods are common the stratigraphy of confined sandur settings can be entirely dominated by depositional and erosional features (Carrivick et al., 2004a, b;Russell, 2007Russell, , 2009a. Unconfined sandur have often been considered as a genetic equivalent to alluvial fans (Boothroyd and Ashley, 1975;Boothroyd and Nummedal, 1978). However, these original models of unconfined sandur environments (e.g. Boothroyd and Ashley, 1975;Boothroyd and Nummedal, 1978) only provide a snapshot of sediments near the sandur surface and fail to characterise both the spatial and vertical variability of the sandur stratigraphy, which has been suggested to be highly complex as a result of the interplay of glacier fluctuations and variations in the magnitude and frequency of glaciofluvial deposition (Olsen and Andreasen, 1995;Aitken, 1998;Russell and Marren, 1999;Marren, 2002aMarren, , b, 2005Russell, 2009b;Salamon and Zieliński, 2010;Pisarska-Jamroży and Zieliński, 2014;Pisarska-Jamroży, 2015;Flindt et al., 2018;Mleczak and Pisarska-Jamroży, 2021;Mleczak et al., 2021). Over geological timescales the proximal source area, and thus the make-up of the distal sedimentary component, of sandur systems can vary as a result of oscillations of the glacier margin and changes in meltwater outlet location, adding complexity to the stratigraphy of the outwash system (Marren, 2004;Mleczak and Pisarska-Jamroży, 2021). Improved understanding of sandar sediment architecture and stratigraphy is vital to understand how the proximal-to-distal components of these vast depositional systems respond to glacier margin variations, and changes in sediment supply and glaciofluvial activity (Marren, 2004).
Sandur systems and their sediment architecture are assumed to be primarily influenced by the position of the ice-margin at the time of glaciofluvial deposition and the magnitude and frequency regime of meltwater systems that drain the glacier (Marren, 2002b(Marren, , 2005. This complexity is not represented effectively in early sandur models. For example, the 'coarse braiding' identified as an important component of proximal sandar is a response to ice-margin retreat-influenced incision and terrace formation (Galon, 1973a, b;Klimek, 1973;Gomez et al., 2000;Marren, 2004;Toomath, 2013, 2014). The development of entrenched channels and proximal depressions with ice-margin parallel drainage networks are characteristic of episodes of ice-margin retreat (Price, 1969;Galon, 1973a, b;Klimek, 1973;Gomez et al., 2000;Marren and Toomath, 2013), but are not necessarily representative of sandur sediment architecture. During periods of ice-margin advance and stability, aggradation of the sandur system becomes dominant, as numerous braided networks emanate from outlets across the ice-margin delivering sediment to the proglacial environment (Maizels, 1979;Marren, 2002b). Changes in the magnitude and frequency regime during this period can result in variations in the aggraded sediment succession (e.g. Marren, 2005). For example, high-magnitude low-frequency glaciofluvial events (e.g. jökulhlaups) can lead to the formation of thick flood successions interbedded with thin, stacked bar core units deposited by normal ablationcontrolled meltwater activity (Marren, 2005).
This study aims to: (i) identify major bounding surfaces and/ or unconformities within the proximal subsurface of an active sandur in SE Iceland; (ii) establish the impact of glacier fluctuations and episodic sediment flux from high magnitude outburst flood events on the 'proximal' sediment architecture of Skeiðarársandur; and (iii) reconstruct ice-proximal processes in the late Holocene. We utilise extensive low-frequency (40 & 100 MHz) ground-penetrating radar (GPR) surveys of an iceproximal sandur in southeast Iceland to achieve these aims.

Study Area: Skeiðarársandur and Skeiðarárökull
Skeiðarársandur, south-east Iceland, is the largest (~1300 km 2 ) active proglacial outwash system in the world (Fig. 1). The sandur is fed by meltwater systems emanating from the margin of the Skeiðarárökull outlet glacier (Fig. 1), a 1400 km 2 surgetype (Table 1) outlet of the Vatnajökull ice-cap (Björnsson et al., 2003;Roberts, 2005;Magnússon et al., 2011). The sandur plain has a shallow gradient and currently extends approx. 25 km from the terminus of Skeiðarárökull to the coast.
In western Skeiðarársandur, a large and partially buried icemarginal landsystem related to the Sandgígjur (alternatively known as 'Sandgígur' in historic maps and previous literature) moraines (Fig. 1) has been identified (Harrison et al., 2022). This moraine and associated ice-contact/end-moraine fan are hypothesised to represent a pre-LIA terminus position of Skeiðarárökull (Harrison et al., 2022). Surface and subsurface components of the Sandgígjur moraine indicate a maximum height of 67 m (based upon a radiowave velocity of 0.13 m − ns 1 ) (Harrison et al., 2022). A total of approx. 1 km 3 of glaciofluvial sediment, associated to jökulhlaup events, has buried the moraine following glacier retreat (Harrison et al., 2022). It is highly likely that the Sandgígjur moraines continue eastwards across the sandar. Therefore, further geophysical investigations of the proximal regions of Skeiðarársandur are crucial to identify the continuation of the Sandgígjur moraines to the east.
It is theorised that jökulhlaup events are critical in the sedimentation and formation of the vast outwash systems that dominate the south Iceland lowlands (Þórarinsson, 1974;Maizels, 1993b;Gudmundsson et al., 2002;Marren, 2005). Numerical modelling predictions suggest that the stratigraphic succession of Skeiðarársandur is composed of up to 94% of jökulhlaup-related deposits (Maizels, 1993b). However, these numerically modelled derived predictions are only constrained by historical observations (e.g. Þórarinsson, 1974;Ives, 2007) and recent observations of jökulhlaup sediment flux at Skeiðarársandur (e.g. Magilligan et al., 2002;Smith et al., , 2006Roberts, 2005;Roberts et al., 2005). In recent Holocene history (e.g. last 1000 years) jökulhlaups have been a common occurence at Skeiðarársandur, with a marked increase in recorded jökulhlaup events following the end of the Little ice-age (LIA) (Fig. 2).
The earliest recorded jökulhlaup is thought to have inundated Skeiðarársandur as early as 1201 as a result of drainage from the Graenalón ice-dammed lake (Þórarinsson, 1939;Roberts et al., 2005;Ives, 2007, Table 2). However, the majority of the jökulhlaups that have drained from Graenalón have predominately inundated the western regions of Skeiðarársandur routing primarily through the Súla meltwater system   Table 2). Volcanically generated jökulhlaups at Skeiðarársandur are believed to have increased in frequency following 1783 when the Laki eruption likely reactivated the main calderas within the Grimsvötn subglacial volcanic system (Guðmundsson, 1992). Jökulhlaups generated from eruptions of Grimsvötn have predominately routed through the Skeiðará meltwater system and meltwater outlets in the central portions of the sandur (e.g. Háöldukvísl, Fig. 1) (Þórarinsson, 1974, Table 2). However, these floods can be so large that they inundate the entire sandur (e.g. Stórahlaup in 1861Þórarinsson, 1974. Volcanogenic jökulhlaups at Skeiðarársandur can transport large ice-blocks (up to 20 m in diameter) and large pieces of detached glacier margin (~1 km in length and 150 m in height) onto the sandur, as well as delivering sediment fluxes up to × m 1.8 10 8 3 (Þórarinsson, 1974;Snorrason et al., 2002;Russell et al., 2006;Blauvelt et al., 2020, Table 2).
The large ice-block depressions located in the more central portions of the proximal sandur (e.g. Fig. 1 B) have previously been geophysically characterised (Blauvelt et al., 2020). Their formation has been suggested to be a product of transportation and deposition of a large segment of glacier ice (~1 km in length and 15 m thick) coincident with numerous ice-blocks during the 1903 jökulhlaup (Blauvelt et al., 2020, Table 2). Aggradation of the sandur and burial of the ice is believed to be a result of further jökulhlaup activity (1922)(1923)(1924)(1925)(1926)(1927)(1928)(1929)(1930)(1931)(1932)(1933)(1934)(1935)(1936)(1937)(1938) in the Háöldukvísl region (Blauvelt et al., 2020, Table 2). However, wider subsurface investigations around these features have not been conducted, limiting what can be inferred about the extent of buried ice deposition and post-depositional modification within the region. Historical accounts indicate that approximately fourteen jökulhlaups have inundated the central portions of Skeiðarársandur from as early as 1598until 1940(e.g. Þórarinsson, 1974Björnsson, 1997;Ives, 2007;Blauvelt et al., 2020, Table 2), where ice-margin recession resulted in the formation of incised channels and the rerouting of meltwater into a large ice-marginal depression ( Fig. 1) (Blauvelt et al., 2020). The impact of this persistent jökulhlaup deposition during the LIA has resulted in the aggradation of central portions of Skeiðarársandur to~125 m above sea-level (a.s.l.), compared to only 90 m a.s.l. in western portions of the sandur (Blauvelt, 2013;Blauvelt et al., 2020).
Skeiðarársandur is a prime target for investigating the largescale sediment architecture of sandur environments due to limited vegetation cover, the lack of soil development and the GPR-friendly properties of outwash sediments. The proximal non-active terraces of Skeiðarársandur are primarily formed of unsaturated sands and gravels (e.g. Robinson et al., 2008) which are an ideal material for the propagation of radiowaves through sediment (Neal, 2004;Reynolds, 2011). Furthermore, historical observations and accounts combined with the abundance of surface geomorphological evidence provide independent constraint against geophysical profiles (e.g. Tables 1, 2). The majority of past GPR surveys at Skeiðarársandur have been limited to specific landforms and surveys targeting recently emplaced jökulhlaup deposits (e.g. Cassidy et al., 2003;Burke et al., 2008Burke et al., , 2010aBlauvelt et al., 2020). The lack of large sedimentary exposures south of the iceproximal depression, combined with the sparse nature of past geophysical surveys provides a new opportunity to characterise the large-scale sediment architecture of abandoned and preserved proximal-sandur terraces by utilising extensive 'reconnaissance' low-frequency GPR surveys.

Súla 1977c
Graenalón ice-dammed lake drainage. Oct 1997. Súla 1978 Graenalón ice-dammed lake drainage. Súla 1979 Graenalón ice-dammed lake drainage. Súla 1980 Graenalón ice-dammed lake drainage. Súla 1981 Graenalón ice-dammed lake drainage. Súla 1982 Medium hlaup -1983a Very minor hlaup -1983b Graenalón ice-dammed lake drainage. Súla 1984 Graenalón ice-dammed lake drainage. Súla 1986a Minor hlaup -1986b Graenalón ice-dammed lake drainage. Súla 1987 Graenalón ice-dammed lake drainage. Súla 1988 Graenalón ice-dammed lake drainage. Súla 1989 Graenalón ice-dammed lake drainage. Súla 1990 Graenalón ice-dammed lake drainage. Súla 1991 Skeiðarárjökull surge. Minor hlaup -1992 Graenalón ice-dammed lake drainage. Súla 1994 Graenalón ice-dammed lake drainage. Súla 1995 Graenalón ice-dammed lake drainage. Súla 1996a Minor hlaup -1996b Massive  resolution and complementary sub-surface information. Positional data were acquired with a handheld GPS streamed into the GPR recording software during radar surveys. A minimum velocity of 0.13 m − ns 1 , attained through WARR (wide angle reflection and refraction) surveys at the Sandgígjur moraines (Supp. Fig. S.1), was used for time-depth conversion of the radar profiles. Our velocity value corresponds to the radiowave velocity expected in unsaturated sands and gravels (Neal, 2004;Reynolds, 2011). This is supported by the observation of the spring line intersecting with the sandur surface in more medialto-distal reaches of the sandur plain at an elevation of~37 m a.s.l. This is well below our survey elevation of 60-120 m a.s.l. Based upon our determined velocity value and the centrefrequency of the 40 MHz GV7 we achieve a theoretical vertical resolution of 0.81-1.63 m. This limits our ability to resolve sediment structures and/or units that are less than 1.63 m in thickness. However, the low centre-frequency allows the 40 MHz GV7 to achieve penetration depths of up to 100 m (e.g. Harrison et al., 2022). With the 100 MHz GV7, vertical resolution is improved, though penetration range is limited to a maximum of 40 m. Radargrams were processed using typical processing steps for low-frequency GPR data (e.g. Harrison et al., 2022) and interpreted in Reflexw geophysical software (http://www.sandmeier-geo.de/reflexw.html). Topographic information used for elevation correction of GPR profiles was extracted, using the XY trace position from the handheld GPS, from metres above sea-level (m a.s.l) converted ArcticDEM tiles (Porter et al., 2018) to radargram trace headers. OpendTect seismic interpretation software was used for 3D visualisation of processed data (https://www.opendtect.org/). Descriptions of radar reflections used for interpretation of radar data are taken from Neal (2004).

Radar facies
In order to facilitate the analysis and interpretation of radargrams collected during this study a number of radar facies were identified, based upon analogous information present within the literature (e.g. Tables 3, 4, 5, 6), and constraining geomorphological evidence (e.g. Fig. 1). Radar facies identified within this study are grouped into four categories: (A) glaciofluvial (Table 3); (B) morainic (Table 4); (C) buried ice meltout (Table 5); and (D) subaqueous and basin-fill deposits (Table 6). Detailed description of the structures and facies of each radargram collected in the proximal reaches of Skeiðarársandur are presented below, followed by their interpretation, and geomorphological significance for the reconstruction of the proximal regions of Skeiðarársandur.

Radargrams SKCP
Radargram SKCP is an across-sandur (E-W) profile that follows an old access road extending eastwards roughly 12.5 km from the Gígjukvísl meltwater system (Fig. 1 A). A moderatelycontinuous high-amplitude reflection, hereafter known as the bounding surface PR1, is present throughout the majority of the profile (Fig. 3). We can be confident the PR1 surface is not the water table as it truncates primary reflections within the radargram, rather than cross-cutting the primary structures (Neal, 2004). The amplitude of PR1 varies and the reflection has an undulating morphology, with an elevation range of 20 m (between 37-57 m a.s.l.). Average thickness of material   Rist (pers. comm.) in Nummedal et al. (1987).~37 jökulhlaups, between 1951 and 2003, have drained from Graenalón with an ever decreasing discharge i.e. peak discharges only -− m s 100 300 3 1 in 2001 . Boundaries of thick sediment units associated to large aggradational jökulhlaup events.
Relevant sources: Beres et al. (1995Beres et al. ( , 1999; Jakobsen and Overgaard (2002) above the reflection is~33 m, with the thickest overlying materials located 4 km and 7.5 km along the radargram (Fig. 3). The majority of reflections in the radargram are above PR1. However, during the first 4 km of the radargram reflections can be seen below PR1 up to a depth of~60 m (an elevation of 10 m a.s.l.). These reflections are predominately characterised by a set of 5-6 moderately continuous, low-amplitude, sub-parallel reflections 0.75-1.75 km along the radargram (Fig. 3 C). Uniform deposit thicknesses of~3.5 m or~6 m can be observed between these reflections. An oblique eastward dipping reflection (slope of ∘ 2. 3 ) can be identified in the first 500 m of the radargram, originating prior to the start of the PR1 bounding surface, and covering an elevation range of 20 m (20-40 m a.s.l.). Reflections in the first 3 km of the radar survey are dominated by a 40 m thick unit of moderately continuous, sub-parallel reflections, including PR1 (Fig. 3 C). Another strong continuous reflection (PR2) is apparent above PR1 (Fig. 3 C). However, PR2 is limited to only the first 2.75 km of SKCP and is primarily located adjacent to the Gígjukvisl meltwater channel and in front of the large LIA push-moraines (Fig. 3 C). PR2 delineates the upper boundary of a~10 m thick unit of planar sub-parallel reflections (Fig. 3 C). Above this unit of sub-parallel reflections there are two obvious packages (0.8-2.5 km and 3.1-4.5 km along the radargram) of highly chaotic reflections with high levels of hyperbolic diffractions (Fig. 3 C). These reflection packages occur in obvious lows in the PR1 reflection and where the sandur surface morphology is characterised by an elevation high. It is worth noting that this portion of the radargram is directly in front of the large LIA push moraines and traverses two obvious outwash fans (Fig. 1 A). The morphology of this reflection package, including PR1, shows morphological similarity to the sandur surface.
Bounding surface PR1 is not identifiable over the entirety of the radargram, with an obvious break in the reflections 5-7.4 km along the radargram. Coincident with this gap in PR1 are large meltout structures in the surface geomorphology (e.g. Fig. 1 B) (Blauvelt et al., 2020), as well as the presence of highly chaotic, trough-shaped and oblique reflections (    Potential indication of buried ice present within the sub-surface.  Basin fill sediment architecture possibly linked to the presence and infilling of a proglacial lake system. ( Fig. 3 D). Discontinuous and sub-horizontal reflections can be identified in between the trough-shaped and very oblique reflections. Reflections in the eastern-most parts of the radargram (i.e. 8-12.5 km along the profile) become difficult to identify (Fig. 3). The returned power in this section of the radargram is of a much lower amplitude, and reflection characteristics are increasingly chaotic with a high number of hyperbolic diffractions. PR1 becomes visible again during this portion of the radargram; however, it is more discontinuous than in the more western portions of SKCP (Fig. 3). We can be confident that this is a continuity of the PR1 reflection due to its stratigraphic position, similar radar properties, and its presence in multiple intersecting radargrams. 100 MHz radargrams coincident with the 40 MHz data (Fig. 4) provide added context to the long across sandur profile (e.g. Fig. 3). The 100 MHz radargrams are broken into two parts (Fig. 4 & Supp. Fig. S.2). PR1 is identifiable within the first 2 km of SKCP-100-A (Fig. 4), at an elevation range of 41-49 m, and demarcates the lowermost reflections recorded within the profile. The reflection characteristics and morphology of the PR1 horizon in the 100 MHz GPR data are comparable (e.g. high amplitude and continuous reflection) to that in the 40 MHz profile (Fig. 3) and they only differ in elevation by~2 m. This is likely a product of the higher resolution of the 100 MHz GPR system identifying the boundary of the unit more accurately (e.g. Neal, 2004). The lower frequency 40 MHz GV7 will likely 'bundle' reflections as a result of decreasing resolution with decreasing radar frequency (Neal, 2004). Over the first 1.25 km of the radargram a 10 m thick package of subparallel reflections directly overlies PR1 (Fig. 4 C), with the uppermost reflection of the package corresponding with a moderately continuous sub-parallel reflection identified within the first 2.5 km of SKCP. This package of sub-parallel reflections is most continuous 0.5-1.25 km along the radargram. The two packages of chaotic reflections located directly beneath the surface elevation highs are better resolved in the 100 MHz profile. The radar reflections beneath the first of the two surface elevation highs (0.25-1.25 km along the radargram) (Fig. 4) are characterised by mostly discontinuous planar to wavy reflections, with indications of occasional trough-like structures. At the top of this package there is evidence of reflections sub-parallel to the surface morphology. Reflections beneath the second surface elevation highs (2-3.4 km along the radargram) are generally more undulating and show less consistency with the surface morphology. Reflection geometries are mostly discontinuous to moderately continuous wavy and undulating (Fig. 4 D). Cross-trough reflections can also be identified 2.25-2.5 km along the radargram at an elevation range of 60 -80 m a.s.l. and suggest the presence of a large channel body.

SKC2
A high amplitude, laterally continuous reflection is identifiable in the first 1 km of radargram SKC2-5 at an elevation of 43-45 m (Fig. 5). We are confident that this reflection is PR1 due to its similar stratigraphic location and reflection characteristics. It is subhorizontal and moderately continuous with a high amplitude (Fig. 5). The first 1 km of the radargram is characterised by a package of high amplitude, very continuous sub-horizontal  (Fig. 5). These stacked V-shaped reflections are coincident with a large depression visible in the surface geomorphology (e.g. Fig. 1 B). The main package of V-shaped reflections has a thickness of~30 m and a width of~180 m (Fig. 5 Bi). Low-amplitude reflection free zones are apparent in the vicinity of the large V-shaped reflection package (1.25-1.50 km along Fig. 5), and are coincident with similar reflection free zones in the SKCP radargram (e.g. Fig. 3 D). Reflection characteristics vary considerably between 1.25-2.48 km along the radargram. Reflections here are moderately continuous, vary in amplitude, and are often of a trough-shaped to undulating morphology. They contrast markedly with the horizontal continuous reflections between 0-1 km along the radargram. Hummocky to subhorizontal onlapping reflections are located 1.57-2.25 km along the radargram (Fig. 5 Bii) at an elevation range of 72-87 m a.s.l. (depth of 17-30 m) with a package thickness of~13 m. Radargram SKC2-4 is a south-to-north trending profile acquired south of Route 1 (Fig. 1 A) with the 40 MHz GV7 (Fig. 6). Reflections can be seen at a depth range of up to 50 m, though the highest amplitude reflections are located in the upper 25-30 m of the profile. PR1 is present throughout the entirety of the profile and is moderately continuous with a high amplitude and a sub-horizontal geometry. Material thickness above PR1 is between 5-10 m, thickening in an ice-proximal direction (Fig. 6). Reflections located above PR1 are chaotic and primarily characterised by hyperbolic diffractions, whilst reflections below PR1 are characterised by an undulating and wavy morphology, and are mostly continuous with a high amplitude (Fig. 6).

SKEC
In the eastern most ice-proximal radargram (Fig. 7) a high amplitude, continuous, gently up-sandur dipping reflection is visible and dissects two distinctive radar reflection packages. Given the characteristics of this reflection (e.g. high amplitude and continuous) and its stratigraphic position between two differing radar facies units we determine that this is a continuation of PR1 (Fig. 7). Reflections above the extensive PR1 form a wedge which tapers in a proximal to distal direction from roughly 40 m thick to 5 m thick. The reflections above PR1 are highly chaotic and characterised by a mix of planar and oblique reflections dipping both up-and downsandur. Some evidence of deformed and disrupted reflections can be seen in the most proximal parts of the radargram (0-0.5  km along the radargram) (Fig. 7). Reflections below PR1 are generally less chaotic with a planar to wavy geometry. Some oblique and potentially deformed reflections are also evident within the lower package of reflections 4.5-4.8 km along the radargram at an elevation range of 45-60 m a.s.l. (Fig. 7 Bi).

Major bounding reflection (PR1)
The radargrams presented within this study are all characterised by the major bounding reflection PR1 (Figs 3, 5, 6, 7, 8). The reflection is mostly laterally continuous and of high amplitude. PR1 is visible in all but one radargram, but its extent in each radargram varies, largely due to the increasing thickness of the overlying materials in more ice-proximal regions leading to the attenuation of the transmitted radiowave. In 3D, the unit above PR1 resembles a wedge with the thickest portions (45-50 m) in more ice-proximal locations (Fig. 8). In more ice-distal locations of the survey the material thickness above PR1 reduces to between 5-10 m. In proximal-to-distal radargrams PR1 dips gently up sandur with an average slope of ∘ 0. 2 . The maximum dip angle ( ∘ 0.34 ) is observed in the SKC1-1 radargram and the minimum dip angle ( ∘ 0.08 ) was observed in the combined SKC2 radargrams (Fig. 8 C-F). PR1 is much more variable in the across sandur profile orientation (Fig. 8 B) with pronounced sections of increased material thickness coinciding with localised elevation highs emanating from the LIA push moraines, visible in the surface geomorphology (Figs 1 B, 8 A). In more eastern parts of the SKCP radargram (Fig. 3) PR1 has a much less variable morphology and a more uniform overlying deposit thickness (Fig. 8). Further detailed descriptions of PR1 associated with individual radargrams and the relationship with additional reflections are discussed below.

Radargram interpretations
The thick, high amplitude and continuous nature of PR1 and its presence throughout the majority of the survey lines suggest that it is a major interface with a strong dielectric contrast potentially indicating an unconformity or bounding surface within the subsurface of central to eastern proximal regions of Skeiðarársandur. Extensive and bright radar reflections are often indicative of a major erosional surface or a distinctive change of the sediment properties (e.g. sediment type, compaction, orientation) of the subsurface (Neal, 2004). There are a few potential hypotheses which can explain this surface. These include: (i) PR1 represents a former sandur surface prior to the maximum LIA advance of Skeiðarárökull; (ii) PR1 is a major erosional or depositional surface linked to a jökulhlaup (e.g. Stórahlaup in 1861) that inundated the entire sandur (Þórarinsson, 1974); or (iii) PR1 is a subglacial decollement surface. We prefer the interpretation of PR1 as a former sandur surface following inundation by a sandur-wide jökulhlaup flow. This is based upon its extensive nature across the survey area, its high amplitude, and the thick nature of the radar reflection. These indicate that it is a major interface and change in sediment (& physical) properties (Jol and Bristow, 2003;Neal, 2004;Reynolds, 2011;Robinson et al., 2013). Furthermore, a subglacial decollement or erosional surface seems unlikely due to the lack of a deformational signature, evidencing glaciotectonic processes, below the PR1 surface. We suggest that sedimentation above PR1 occurred as Skeiðarárökull was at a maximum LIA extent, as ice advance is known to influence proglacial aggradation (Maizels, 1979;Marren, 2002b). Based on the interpretation of a sandur-wide jökulhlaup-related unconformity we can provide a minimum age estimate of the surface, the most recent of which is the Stórahlaup in 1861 (Þórarinsson, 1974, Table 2). If PR1 represents a non-jökulhlaup-related former sandur surface, then a maximum age of sedimentation can be inferred and associated with the first LIA advance of Skeiðarárökull (e.g. 1784) (Jóhannesson, 1985;Björnsson et al., 2003).
The gentle up-sandar dip of PR1 in down-sandur radargrams (Fig. 8) is difficult to reconcile with this interpretation, as it would be expected that the sandur surface was flat or dipped gently down-sandur towards the coast. However, there are a number of possible explanations why PR1 has this geometry: (i) compaction of sediments below the surface by increased weight of overlying sediments; (ii) down-draw of the surface as a result of isostatic adjustment due to ice growth during the LIA (e.g. Whitehouse, 2018); (iii) dewatering of jökulhlaup deposits from the region (e.g. Duller et al., 2014); (iv) meltout of any historic widespread buried ice deposits/units (e.g. Harrison et al., 2019); or (v) the former presence of an icemargin parallel drainage network (e.g. Fig. 1 B). Furthermore, slight inaccuracies in the radiowave velocity used for timedepth conversion may result in the 'false' dip of a reflection (Neal, 2004;Reynolds, 2011).

Basin Fill deposits
The package of stacked moderately continuous, planar to wavy, sub-parallel reflections located in the first 2 km of SKCP  (stratigraphically below PR1) indicate boundaries of individual depositional units with thicknesses of 3.5-6 m (Figs 3, 9). The radar systems deployed during this survey cannot resolve bedforms or structures within these units which makes it difficult to determine the processes responsible for their deposition. It is also difficult to assess the continuity of these reflections over the entirety of the survey area as reflections below PR1 in SKCP are limited beyond the first 4 km along the radargram (Fig. 3). The consistency of these reflections and their parallel nature to depth are, however, characteristic of basin fill architecture (Table 6), and potentially an indicator of deposition in a former glaciolacustrine environment (e.g. Roberts et al., 2018). The consistency of the distance between reflections, and their sub horizontal nature, suggest a consistent period of aggradation.

Evidence of surge-related deposits
The continuous sub-parallel, planar reflections identified above PR1 in the first 2 km of the SKCP and SKCP-100-A radargrams are most likely representative of stacked individual sediment units (Figs 3,4,9). The total radar package thickness is~10 m with individual radar unit thicknesses of between 2-5 m. This unit of stacked sub-parallel reflections are topped by a~20 m thick package of sediments that are indicative of outwash fan architecture (e.g. chaotic, planar to wavy and cross-trough reflections) ( Fig. 9 and Table 3). The combination of these reflection packages is suggestive of surge-related outwash deposits on a large-scale (e.g. Russell et al., 2001;van Dijk, 2002;van Dijk and Sigurðsson, 2002). Stacked and sub parallel reflections are an indication of horizontally bedded sands deposited during glacier advance by diurnal meltwater flows (Krüger, 1997;Maizels, 1993a;Russell et al., 2001). However, the unit thicknesses outlined by these reflections (e.g. 2-5 m) are comparable to jökulhlaup unit thickness (e.g. Maizels, 1993bMaizels, , 1997Marren and Schuh, 2009). Therefore, these cannot be ruled out as jökulhlaup deposits, related to sandur-wide jökulhlaup deposition (Þórarinsson, 1974). The upper unit of chaotic, planar to wavy and cross-trough reflections (Table 3) are interpreted as representing a poorly sorted and matrix supported outwash unit deposited during a surge-termination flood event (e.g. Russell et al., 2001). The strong reflection (above PR1) that separates these packages, visible in the first 2.5 km of SKCP and SKCP-100-A (Figs 3, 4) is associated with an erosional surface that separates the lower reflection unit of sands and the upper chaotic reflection unit indicative of outwash fan deposition. The presence of an erosional boundary between diurnal meltwater deposits and flood-related deposits has previously been identified in surge-related fans at Skeiðarársandur (Russell et al., 2001;van Dijk, 2002;van Dijk and Sigurðsson, 2002). Furthermore, outwash fans at the Gígjukvísl gap ( Fig. 1) have previously been hypothesised to be a result of surge related processes, but on a much larger-scale than what was evident from the 1991 surge (Russell et al., 2001). The absence of kettle-holes on the sandur surface of the fans at the Gígjukísl gap ( Fig. 1 B) add further evidence that these are surge-related deposits as opposed to jökulhlaup-related sedimentation (Russell et al., 2001).

Buried ice and jökulhlaup deposits
Sub-horizontal, parallel and moderately continuous reflections (Table 3) are identified within radargram 9). These reflections delineate the boundaries of individual sediment units each with a thickness of roughly 2.5 m. Jökulhlaup events are the most likely interpretation for the deposition of these sediment units. The individual unit thickness of~2.5 m is comparable to known jökulhlaup unit thicknesses reported at Skeiðarársandur (Maizels, 1993b(Maizels, , 1997Marren and Schuh, 2009  during the LIA (Þórarinsson, 1974;Ives, 2007;Blauvelt et al., 2020). The continuity of these sub-horizontal reflections is interrupted by the chaotic package of radar reflections, highlighted by large trough-shaped and highly oblique reflections (Table 5) located 1-2 km along the SKC2-5 radargram (Fig. 9). Similar structures are also present 5-7.5 km along the SKCP radargram, which intersects the SKC2-5 radargram at the location of the large meltout landforms visible in the surface geomorphology (Fig. 1). The combination of highly oblique and trough-shaped reflections are indicative of post-depositional meltout structures formed following the degradation of buried ice masses (Maizels, 1992;Fay, 2002b;Marren and Schuh, 2009;Harrison et al., 2019;Blauvelt et al., 2020). The presence of these next to visible kettle holes and meltout structures in the surface geomorphology reinforces this interpretation (Blauvelt et al., 2020). Furthermore, this adds support to the interpretation of the sub-horizontal reflections in SKC2-5 as boundaries of jökulhlaup-related sediment units, as large ice-blocks are known to be transported and deposited during high-magnitude outburst flood events (Maizels, 1992;Fay, 2002a, b;Russell et al., 2006;Marren and Schuh, 2009;Burke et al., 2010a;Dunning et al., 2013;Harrison et al., 2019;Blauvelt et al., 2020). Ice-block obstacle marks and meltout structures are considered key elements of jökulhlaup activity within the sedimentary record (Fay, 2002a, b;. Our data supports the model presented by Blauvelt et al. (2020) of a jökulhlaup landsystem where a portion of the icemargin broke off and was transported onto the sandur by jökulhlaup flood waters before being buried by the aggrading sandur. Where buried ice is actively present in the subsurface radar reflection are often absent and zones of low-amplitude will be recorded in the radargram (e.g. Fig. 10 and Table 5) as a result of the relative homogeneity of the ice (Moorman et al., 2003;Brandt et al., 2007). Low amplitude zones in the SKCP radargram ( Fig. 10) are indicative of buried ice (6-10 m in thickness with a horizontal extent of up to 100 m). These zones are coincident with the large trough structures and highly oblique reflections (Fig. 10), diagnostic of faults and depressions associated with ice meltout structures (Price, 1971;McDonald and Shilts, 1975;Maizels, 1992;Smellie et al., 2016;Harrison, 2018;Blauvelt et al., 2020). This evidence provides further support that buried ice is present at these locations.

Evidence of glaciotectonism
There is limited evidence of reflections suggestive of glaciotectonic processes. Where these deformed and disrupted reflections (Table 4) are present they are of a limited extent and lack coincident/cross-cutting GPR data. A package of hummocky and oblique radar reflections (Figs 7 Bi, 11), situated stratigraphically below PR1, can be identified in two down-sandur orientated radargrams (SKC1-1 & SKEC3). These reflection characteristics are interpreted as deformed sediment units linked to push moraine formation (e.g. Bennett, 2001;Overgaard and Jakobsen, 2001;Jakobsen and Overgaard, 2002). The location of these reflection packages are a natural extrapolated extension of the buried Sandgígjur moraines located in the western portions of Skeiðarársandur (Harrison et al., 2022) and are coincident with a moraine location reported on historical maps of Skeiðarársandur and Skeiðarárökull (Þórarinsson, 1939). If these features represent a moraine ridge it is possible that it signifies a central-and eastern-sandur representation of the Sandgígjur moraines, and may be linked to a pre-LIA position of Skeiðarárökull (Harrison et al., 2022). However, further geophysical and sedimentological (e.g. borehole constraint) evidence is crucial for further identification and interpretation of these features. Hummocky to sub-horizontal onlapping reflections located in SKC2-5 (Fig. 5 Bii) may represent thrust blocks linked to push moraine formation and could correspond with the LIA maximum (1890-95) or 1904 ice-margin of Skeiðarárökull (Fig. 1). However, it is also possible that these structures are representative of backfilling into conduits within the glacier during jökulhlaup sedimentation. Esker formation as a result of infilling of a conduit during a jökulhlaup was reported from the 1996 November jökulhlaup deposits at Skeiðarársandur (e.g. Burke et al., 2008). Given the lack of coincident radar lines or sediment exposures in the region it is difficult to offer a definitive interpretation of these hummocky and onlapping reflections. However, given their coincidence to the 1904 icemargin limits of Skeiðarárökull our preferred interpretation is push-moraine related deposits (Fig. 1).

Evidence of glaciofluvial braid plain architecture
Reflections below PR1 in most southerly regions of the proximal-distal orientated profiles (e.g. Figs 6, 9) show a distinctive contrast compared to reflections above PR1.
Reflections show similarities to braid plain architecture identified within other GPR investigations of proglacial outwash environments and are likely representative of largescale structures or surfaces related to changing of channel locations and the formation of in-channel bedforms (Beres et al., 1995(Beres et al., , 1999Møller and Vosgerau, 2004;Hansen et al., 2009). Undulating, wavy and oblique reflections likely represent large-scale glaciofluvial braidplain components, such as longitudinal gravel bars and bedforms, scour fills, and downstream accretion macroforms (Miall, 1985(Miall, , 1992. The interpretation of this sediment package as a glaciofluvuial braidplain is predominately supported and constrained by historical accounts. Details from early settlement times in Iceland (ca. 10 th Century) indicate that Vatnajökull was much smaller than its present day extent, with a non-glaciated valley between its eastern and western portions, likely coincident with Skeiðarárjökull and the Grímsvötn depression (Þórarinsson, 1974;Nummedal et al., 1987;Ives, 2007). During this period the Skeiðarárökull ice margin is thought to have been 10 km in land of its current position, with meltwater emanating from a single central proglacial river (Þórarinsson, 1974;Ives, 2007). During later LIA-related expansion of Skeiðarárökull, farmlands on Skeiðarársandur were abandoned as a result of erosion and aggradation of the sandur by jökulhlaups (Ives, 2007).

Reconstruction of the late Holocene proximal proglacial landsystem at Skeiðarársandur
The interpretation of radar data collected within this study for reconstruction of the ice-proximal landsystem at Skeiðarársandur is underpinned by two critical assumptions; (i) the large extant push-moraines ( Fig. 1) are the limit of LIA glacier advance (Molewski and Olszewski, 2000;Björnsson and Pálsson, 2008;Hannesdóttir et al., 2020), and (ii) the PR1 bounding reflection is a boundary between pre-LIA sedimentation and LIA to present-day sedimentation.

Pre-LIA glaciofluvial deposition
Below PR1, specifically in the Háöldukvísl region, sediments are dominated by reflections that are interpreted to be representative of large bedforms and barforms formed by a glaciofluvial braidplain system (e.g. Miall, 1985Miall, , 1992Beres et al., 1995Beres et al., , 1999Møller and Vosgerau, 2004;Hansen et al., 2009). The interpreted palaeo-braidplain comprises at least 40 m of the sandur subsurface below PR1 (Fig. 6) suggesting the presence of a major aggradational meltwater system that drained the central portions of the sandur plain. Historical accounts suggest that a major braidplain was located in the centre of the sandur and was the primary outlet responsible for draining meltwater from Skeiðarárökull prior to its LIA readvance (Þórarinsson, 1974;Nummedal et al., 1987;Ives, 2007). The~40 m thick braidplain we identify in proximal regions of Skeiðarársandur (e.g Fig. 6) may be attributed to a glaciofluvial delta topset that has formed as a response to rapid glacier recession (Lønne, 1995;Tuttle et al., 1997;Bataller et al., 2021). Given the lack of basal contact within any of the radargrams it is difficult to establish the absolute thickness of this radar unit. Seismic investigations in proximal locations of Skeiðarársandur reveal a unit of consolidated/compacted sediments that exist beneath 70-100 m of unconsolidated sediments interpreted to be of Holocene age (Gudmundsson et al., 2002). Therefore, it is plausible that we have only imaged the upper 40 m of a much thicker glaciofluvial unit located in the central proximal regions of the sandur, subsequently capped by more recent (e.g. LIA) glaciofluvial sediments. The formation of the Gardermoen Delta in Norway has been attributed to be a response of persistent maximum ablation-driven meltwater drainage during rapid deglaciation (Tuttle et al., 1997). However, the increased thickness of the potential glaciofluvial braidplain radar unit we identify (e.g. Fig. 6) is likely a product of increased jökulhlaup frequency at Skeiðarársandur (e.g. Fig. 2 & Table 2) as opposed to maximum ablation-related meltwater drainage, driven by deglaciation (Tuttle et al., 1997). Roughly five 1996magnitude jökulhlaups per century are believed to be required, on top of sedimentation from ablation-driven meltwater flow and surge-related outwash events, to account for the sediment accumulation at Skeiðarársandur (Gudmundsson et al., 2002;Snorrason et al., 2002;Russell et al., 2006).

Pre-LIA ice-margin position
Identification of folded reflections that are indicative of morainic and glaciotectonic processes suggest a possible icemargin position down-sandur of the large LIA moraines that may be associated to a pre-LIA position of Skeiðarárjökull (Figs 11, 12). The interpreted moraine-related folded and oblique reflections in the SKC1-1 and SKEC3 radargrams (Fig. 11) are stratigraphically below the PR1 surface limit the down-sandur extent of the PR1 reflection. The identification of these possible moraine features in the radar data may be a lateral continuation of both a moraine limit proposed in the central part of the sandur by Þórarinsson (1939), and the Sandgígjur moraine located in western Skeiðarársandur (Harrison et al., 2022) (Fig. 12). Furthermore, this proposed pre-LIA ice-margin of Skeiðarárjökull may have integrated with outlet glaciers on the east of the sandur, which had Holocene extents that emerged much further onto Skeiðarársandur than what occurred during the LIA (Gudmundsson, 1997). It has previously been suggested that morainic features may be preserved beneath the thick jökulhlaup and glaciofluvial sandur successions (Haraldsson and Palm, 1980;Haraldsson, 1981;Gudmundsson et al., 2002), and large subsurface morainic features has recently been confirmed at depth at Skeiðarársandur (Harrison et al., 2022). However, the absence of similar reflections within the subsurface in central proximal areas of Skeiðarársandur, particularly in the downsandur-orientated SKC2 radargrams (Figs 5,6), suggests that preservation of these moraines located beyond the LIA limit is low within proximal regions of the sandur. This low preservation potential is likely increased where the sandur is prone to a relatively persistent frequency of high magnitude jökulhlaups, which have been prevalent in the central regions of the Skeiðarársandur during the LIA (Þórarinsson, 1974;Ives, 2007;Blauvelt et al., 2020, see Fig. 2). The lack of identification of extensive moraine features in the subsurface of central Skeiðarársandur may therefore be a reflection of the central parts of the sandur being disproportionately impacted by jökulhlaup events, at least relative to its western parts.

LIA to present-day sedimentation and post-depositional modification
The PR1 boundary likely reflects a boundary between pre-LIA and sedimentation in the LIA and post-LIA. The continuity and brightness of the reflector (relative to other reflections in the entire dataset) over large portions of the GPR survey suggests it is a major unconformity. We suggest that formation is a result of the Stórahlaup in 1861 (Þórarinsson, 1974). This would mean that the sediment wedge above PR1, that caps the thick glaciofluvial braidplain, is a product of deposition from 1861 onwards. It is important to note that this assumption is not constrained by physical evidence (e.g. absolute dating) but is consistent with historical records (e.g. Þórarinsson, 1974;Ives, 2007). Sandur aggradation and the formation of steep proximal outwash fans, such as we observe in our radar data, are a common sandar-related landform that develop during advanced ice-margin positions (Maizels, 1979;Marren, 2002b). The presence of a thick sediment wedge with a steep surface gradient above PR1, tapering in thickness from~40 m to~5 m in a proximal to distal direction (Fig. 8), is likely a product of an increasingly limited terrestrial accommodation space leading to a proximal steepening of the sandur (Maizels, 1979;Marren, 2002a, b;Roussel et al., 2018). The presence of this sediment wedge, formed of coalesced outwash fans of alluvial-fan type sedimentation, stratigraphically above glaciofluvial braidplain architecture suggests that outwash fan and alluvial fan type sedimentation is a product of recent glacial advance (e.g. last 250 years) influenced sedimentation. It is not necessarily a major component of the sandur-wide sediment architecture to depth, as previously suggested (e.g. Boothroyd and Nummedal, 1978;Zieliński and van Loon, 2003). If our assumption of PR1 being the product of the 1861 jökulhlaup is correct this would mean the majority of the sediment wedge can be linked to sedimentation associated to LIA advances in 1873 and 1890-95 (Jóhannesson, 1985;Björnsson et al., 2003;Björnsson and Pálsson, 2008), and the period of relative icemargin stability up until the mid-1940s (Klimek, 1973;Gomez et al., 2000). High levels of sandur aggradation are common during periods of ice-margin advance and stability (e.g. Maizels, 1979;Marren, 2002b;Storms et al., 2012;Roussel et al., 2018). The LIA-related proximal outwash-fan sedimentation is spatially heterogeneous: (i) at the Gígjukvísl gap surgetermination related outwash occurred; whilst (ii) jökulhlauprelated deposition occured in the Háöldukvísl and more eastern ice-proximal regions. The west-to-east transition of LIA sedimentation form surge-related outburst floods to jökulhlaup events shows the heterogeneity of sandur sedimentation even in a relatively short timescale (e.g. LIA to present-day). Identification of radar architecture resultant from surge-related outburst floods supports previous interpretations of outwash fans emanating for the large LIA moraines at the Gígjukvísl gap being surge-related as opposed to jökulhlauprelated (Russell et al., 2001). Surge fans often form along the ice-margin in locations distal to the coarse-grained braidplain that usually drains the ice margin (Marren, 2005). The prevalence of these fans only in the regions of the Gígjukvísl gap may therefore be driven by the subsurface palaeobraidplain in the Háöldukvísl region (Fig. 9). Distinguishing ice-marginal fans from surge-related fans can be difficult, however, as there are relatively few variations in their architecture (Krüger, 1997;Zieliński and van Loon, 2000;Russell et al., 2001;van Dijk, 2002;van Dijk and Sigurðsson, 2002;Marren, 2005). However, the known surge history of Skeiðarárökull (e.g. Björnsson et al., 2003, see Table 1) can provide some support for our interpretation. It has also been hypothesised that many Pleistocene terminoglacial outwash fans are actually related to surge processes as opposed to 'normal' glaciological conditions (Þórarinsson, 1969;Drozdowski, 1987;Marren, 2005). Assuming that the PR1 surface is the product of the 1861 jökulhlaup, we hypothesise that the surge-related fans at the Gígjukvísl gap may be a result of the 1873 surge (Björnsson et al., 2003) and possible further aggradation during the LIA maximum in 1890-95 (Björnsson and Pálsson, 2008).
East of the possible surge-related outwash fans (Fig. 9) that emanate from the large LIA moraines, sandur aggradation above PR1 (up to 50 m thick in the most proximal regions, Fig. 8) is attributed to jökulhlaup activity. The presence of jökulhlaup-related outwash deposits is evidenced by the chaotic nature of radar reflections, in combination with the Figure 12. Possible pre-LIA extent of Skeiðarárjökull ('The Sandgígjur stage') that incorporates evidence from the SKC1-1 and SKEC3 radargrams (Fig. 11) and the Sandgígjur moraine (Harrison et al., 2022). Approximate locations of Holocene ice-margins of Skaftafellsjökull (black dash-dot lines) from Gudmundsson (1998). Extent of the 1890/95 ice-margin of Skeiðarárjökull from Hannesdóttir et al. (2020). Background satellite image is a September 2019 cloud-free composite acquired from Planet (PlanetTeam, 2017). [Color figure can be viewed at wileyonlinelibrary.com] presence of ice-block meltout structures, a common diagnostic feature of jökulhlaup activity in the sedimentary and geomorphic record (Russell and Knudsen, 1999;Marren, 2005;Marren and Schuh, 2009;Wells et al., 2022). However, it remains possible that the chaotic unit of reflections above PR1 are a product of deformation linked to glacial advance (e.g. Fig. 7). In competent materials (e.g. gravel-dominated glaciofluvial sediment), pressure from an actively advancing icemargin can be transferred into the proglacial foreland, resulting in chaotic sediment structures linked to deformation (e.g. Boulton et al., 1999;Benediktsson et al., 2010). However, given the highly kettled nature of the topography in the most eastern regions of the survey (Fig. 1) and the well documented jökulhlaup history, we prefer the interpretation of the chaotic unit to be representative of jökulhlaup-related activity. Although, it is plausible that post-depositional modification, influenced by glacier readvances during the LIA, may have contributed to an increase in the chaotic nature of the sediment unit. Which is supported by the presence of highly deformed reflections linked to glaciotectonism (e.g. Table 4) in the most proximal parts of the SKEC3 radargram (Fig. 7).
Ice-block related meltout structures, and the evidence of preserved ice-blocks within the subsurface, are more widespread than what is evident from surface geomorphology. Our data is consistent with the jökulhlaup-landsystem model suggested for the Háöldukvísl region (Blauvelt et al., 2020). The widespread nature of radar architecture consistent with buried ice deposition suggests a large piece of glacier ice was transported and deposited on the sandur plain before being buried by the aggrading sandur during waning stage flood flows and further jökulhlaup events (Blauvelt et al., 2020). The sub-horizontal reflections in 9), down-sandur of ice-block meltout structures, are consistent with burial by jökulhlaups following initial ice emplacement. The emplacement of ice in the Háöldukvísl area can be attributed to the 1903 jökulhlaup (Table 2) where it was reported that a~1 km long and 150 m thick portion of glacier ice was fractured from the ice-margin and emplaced in the central sandur alongside numerous 'house-sized' ice blocks (Þórarinsson, 1974;Blauvelt et al., 2020). The identification of a~5 km 2 area of the subsurface characterised by ice-block meltout structures and reflections indicative of faulting and post-depositional modification supports these accounts (e.g. Figs 9, 10). Our data also suggest that ice-blocks up to 100 m in length and 10 m in thickness (Fig. 10) may be preserved within the subsurface. The size of these preserved buried ice features corresponds to the scale of ice deposition linked to the 1903 jökulhlaup (Þórarinsson, 1974), as well as corresponding to ice-block sizes reported from more recent jökulhlaups (Tómasson, 1996;Russell and Knudsen, 1999;Roberts et al., 2000;Fay, 2002a, b;Roberts, 2005;Russell et al., 2006). Buried and preserved ice-masses are common within the sandar of south-eastern Iceland, and given the thickness of the overlying sediments (approx. 20 m) above the ice-blocks we identify (Fig. 10) it is possible they could survive for up to 100s of years (Everest and Bradwell, 2003).

Terrestrial sediment flux during glacier advance
Given a lack of evidence of large-scale deformation by glaciotectonic processes within the outwash sediments it can be assumed that the majority of sedimentation above PR1 occurred once the glacier was at a relatively advanced LIA position and in a non-surging state. Sediment aggradation is a common occurrence during an advanced glacier position (Maizels, 1979;Marren, 2002a) and sediment supply is often increased during an advanced ice-margin position (e.g. Antoniazza and Lane, 2021). At Skeiðarársandur, this has likely been further enhanced by the reactivation of the main Grímsvötn calderas by the 1783 Laki eruption (Þórarinsson, 1974;Guðmundsson, 1992). Resulting in an increase in the magnitude and frequency of jökulhlaup activity from 1783 onwards (e.g. Fig. 2 and Table 2).
Assuming PR1 represents a former sandur surface, we can provide minimum and maximum sedimentation rates between the first (1784) and last (1890-95) LIA advances of Skeiðarárökull and the present-day (Jóhannesson, 1985;Björnsson et al., 2003;Björnsson and Pálsson, 2008). We can also assume major down-sandur aggradation was limited from 1940 onwards due to ice-margin retreat associated with decoupling of the glacier and sandur surface (Klimek, 1973;Gomez et al., 2000), leading to the rerouting of meltwater into entrenched channel systems (e.g. Gígjukvísl gap and meltwater system). Based on an average sediment thickness above PR1 of 32.3 m (range: 8.9-57.8 m) we can approximate a rate of deposition for the central proximal region of Skeiðarársandur of 0.2-0.65 m − a 1 (i.e. minimum estimate: 1784-1940; maximum estimate: 1890-1940). However, it is unlikely that the sediments above PR1 were deposited at a constant rate. Instead, sedimentation is likely the product of large outburst floods in the region.
Thirteen jökulhlaups between 1861 and 1938 are believed to have impacted the central portions of Skeiðarársandur (Þórarinsson, 1974;Ives, 2007;Blauvelt et al., 2020, Table 2), specifically around the Háöldukvísl gap region (Fig. 1). This would equate to an average jökulhlaup unit thickness of 2.48 m to account for the sediment thickness we observe above PR1. This is within the range of LIA jökulhlaup unit thicknesses predicted by Maizels (1993b). Furthermore, the sub-horizontal reflections observed in radargram SKC2-5 (Fig. 5) are interpreted as~2.5 m thick sediment units, adding further support to the idea that they are jökulhlaup related deposits. The increased sediment thickness of the deposits above PR1 in the more eastern locations of the sandur (e.g. Fig. 8) can be attributed to a greater sediment flux available due to a longer periood of coupling of the glacier and sandur surface (Hannesdóttir et al., 2020;Hauser and Schmitt, 2021). In addition, sediment flux as a result of subglacially-induced jökulhlaups, following the reactivation of Grímsvötn in 1783 (Guðmundsson, 1992), disproportionately affected eastern portions of Skeiðarársandur compared to the central and western sections of the sandur.

Conclusions
• The proximal regions of Skeiðarársandur are composed of a complex, heterogeneous, and often chaotic sediment architecture. The PR1 surface identified throughout the majority of the radargrams is the only major bounding reflection visible in the GPR survey of the central, proximal parts of Skeiðarársandur. This surface splits the GPR data into upper and lower depositional packages (Figs 3,5,9). The lack of other extensive bounding surfaces across the dataset indicates the heterogeneity of sandur environments and suggests sandur-wide depositional and/or erosional impacts can be rare within the sediment record. • We identify a large palaeo-meltwater system in central portions of Skeiðarársandur that likely predates the LIA. Evidence of this beneath a proximal wedge of LIA outwash fan type deposition suggests that the majority of Skeiðarársandur to depth (up to 40 m) is the depositional product of a glaciofluvial braidplain, formed following rapid deglaciation during the Holocene.
• Radar data indicate a maximum approx. 50 m thick proximal outwash fan terrace (Fig. 8), which is comparable to numerical model derived predictions (e.g. 47 m thick jökulhlaup-related proximal terrace) of the Skeiðarársandur sediment succession (e.g. Maizels, 1993b). The thicker and more extensive LIA-related (e.g. deposits located above PR1) ice-proximal outwash sedimentation can be associated to the length of glacier and sandur coupling. Thickness of the deposits above PR1 are greatest in more eastern locations where Skeiðarárjökull has retreated from its LIA maximum at a much slower rate than in more central and western sandur locations (e.g. Hannesdóttir et al., 2020;Hauser and Schmitt, 2021). • We can attribute the majority of aggradation above PR1 to be due to jökulhlaup-related deposition when Skeiðarárjökull was advanced and coupled to the sandur surface. This is most pronounced in the Háöldukvísl (central) and easternsandur regions. The formation of a steep proximal sediment wedge fits with models of sandur aggradation during an advanced ice-margin state. • Radar data across outwash fans close to the Gígjukvisl gap provides supporting evidence that these outwash fans may be a result of a surge-related outwash processes. This suggests that, in some cases, large ice-marginal outwash fans may be a result of surge-related processes as opposed to normal glaciological conditions. • The signature of ice-block deposition and post-depositional modification of jökulhlaup deposits is a major component of the proximal component of sandur sediment architecture. The degradation of ice-blocks emplaced by jökulhlaups has a major impact on the preservation potential of proximal sandur architecture and can impact estimates of sediment flux and aggradation from these high magnitude meltwater flows. Identification of active, in situ buried ice within the subsurface of proximal regions of Skeiðarársandur requires further geophysical research. • We propose an extensive, sandur-wide, pre-LIA ice marginal limit ('The Sandgígjur stage') of Skeiðarárökull (Fig. 12), based on observations reported here (e.g. Fig. 11) and in previous work (e.g. Harrison et al., 2022). The limited observation of morainic features in the central and eastern portions of the sandur is most likely a product of the increased frequency of high-magnitude jökulhlaups in the region. This, which is in contrast to the large buried and preserved Sandgígjur moraines in western Skeiðarársandur, may reflect a variation in magnitude and frequency regime across the ice-margin of Skeiðarárökull.

Data Availability Statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Supporting information
Additional supporting information can be found in the online version of this article.