Ground‐penetrating radar (GPR) investigations of a large‐scale buried ice‐marginal landsystem, Skeiðarársandur, SE Iceland

The sedimentary record of Icelandic ice‐contact environments provides critical insights into past glacier margin dynamics and position, relative sea level, and the geomorphic processes that drive the evolution of proglacial environments. This important archive has been little exploited, however, with most glacier and sea‐level reconstructions based on limited sedimentary exposures, coring and surface geomorphic evidence. We report an extensive (42 km of data within a 24‐km2 study area) and deep (reflections recorded at depths up to 100 m) low‐frequency (40 and 100 MHz) ground‐penetrating radar (GPR) survey of the Sandgígur moraines, SE Iceland. GPR profiles reveal a much larger (67 m high) and extensive (1.25 km wide) buried moraine ridge than that suggested by surface topography (typically 125 m wide and 7 m high). These data reveal that the Sandgígur moraines was deposited during a major Holocene re‐advance of Skeiðarárjökull. The moraine ridge is buried by sediments dominated by glacifluvial deposits with an estimated sediment volume of 1.04 km3. We combine GPR‐derived subsurface architecture and the surface morphology to develop a conceptual model detailing the geomorphic evolution of the moraine and surrounding region. These results provide new insights into the Holocene evolution of Skeiðarársandur, identifying the presence of a former major ice‐margin position, as well as a past relative sea‐level limit. Furthermore, we establish that sediment supply and available terrestrial accommodation space are dominant drivers in the formation and evolution of vast sandar environments.

Sandur plains (sandar) are depositional systems beyond the margin of a glacier or ice sheet, typically located in the present-day Arctic and common within the Quaternary record of low latitudes. The development and evolution of these depositional systems are primarily controlled by fluctuations in the mass balance of their feeding ice mass (Krigstr€ om 1962;Þ orarinsson 1972;Maizels 1983). However, the topographic setting of sandar, tectonic processes and changes in global sea level also play an important role in their evolution (Maizels 1993a). Despite considerable research (Krigstr€ om 1962;Church 1972;Maizels 1993a;Gomez et al. 2000;Marren 2002a, b;Zieli nski & van Loon 2002, the primary controls on the formation and evolution of sandur environments during glaciation and rapid deglaciation remain to be established. The controls of sandur sedimentation and geomorphological reworking can be broken into active or passive factors (Weckwerth et al. 2021). Active factors include: (i) changing climate conditions and the resultant glacier margin fluctuations; and (ii) variations in the hydrology of proglacial fluvial systems, especially related to the magnitude and frequency regime of the meltwater system (e.g. flood events vs. non-flood events) (Marren 2005;Carrivick & Heckmann 2017;Strzelecki et al. 2018;Kociuba et al. 2019;Weckwerth et al. 2021). Passive factors influencing sandur sedimentation are generally regarded as: (i) longterm geomorphic processes (e.g. location and distribution of existing proglacial fluvial systems); (ii) calibre of the prevailing sediment (e.g. gravel, sand or clay); and (iii) changes in the geology of the catchment (Weckwerth et al. 2021).
Some of the best examples of sandar are those in Iceland. It has been hypothesized that Icelandic sandar comprise thick alluvial successions containing detailed records of the events that contributed to their formation (e.g. Marren 2005). However, surprisingly limited work has been undertaken to analyse the large-scale sedimentary architecture of these systems, and to reconstruct how they have evolved during the Holocene. It can be argued that more is known about the large-scale architecture of palaeo-sandur and glacigenic depositional environments from ancient glaciations (e.g. Ghienne & Deynoux 1998;Ghienne et al. 2010;Girard et al. 2012) and Quaternary glaciations (Brennand 1994;Brennand & Shaw 1994;Tuttle et al. 1997;Winsemann et al. 2009Winsemann et al. , 2018Lang et al. 2021), than is known about contemporary outwash plains. However, inferences from the sedimentary record lack temporal constraint on the processes that have influenced sandur development.
Glacier margin fluctuations are known to generate proglacial fluvial system responses in the form of sandur morphological change. Changes in glacier mass balance influence meltwater hydrology and its ability to modify sandur environments. This, in turn, is dependent on the coupling of the glacial source to the proglacial depositional system (Marren 2005;Weckwerth et al. 2011Weckwerth et al. , 2019Weckwerth et al. , 2021Knight & Harrison 2014Carrivick & Heckmann 2017;Kociuba et al. 2019). Aggradation and a steepening of sandar are seen to occur during periods of glacier margin stability or advance where the glacier is coupled to the sandur (Maizels 1979;Marren 2002b). In contrast, during periods of glacier retreat, where the glacier is decoupled from the sandur plain, vast braided networks transition into singular entrenched channels with increasing terrace formation, as incision dominates (Galon 1973a, b;Klimek 1973;Marren 2002b;Marren & Toomath 2013). The rates of incision and sandur erosion are controlled by extreme precipitation-driven ablation events which impact short-term changes in the meltwater system (Marren 2005;Beylich & Laute 2015;Lane et al. 2017;Guillon et al. 2018;Kociuba & Janicki 2018;Weckwerth et al. 2021). Retreating ice margins can also lead to the development of lateral drainage networks in topographic lows formed following ice retreat from the previously aggraded sandur surface (Price 1969;Bogacki 1973;Galon 1973a, b;Klimek 1973;Gomez et al. 2000).
The relationship between ice-margin fluctuations and the related aggradation or incision of the sandur, coupledwith the sediment supply of the glacial meltwater systems and occurrence of high magnitude flood events leads to a complex sediment stratigraphy (Maizels 1979(Maizels , 1983Marren 2002bMarren , 2005. Early sandur models (e.g. Boothroyd & Nummedal 1978) do not capture this complexity and predominately represent a snapshot of the sandur surface. Over longer time scales this may be repeated vertically. However, fluctuations in ice-margin position, meltwater outlets and sediment supply may lead to variations in the thickness and lateral extent of the overlying deposits (Marren 2004). The inherent relationship between glacier margin fluctuations and sandur formation means that sandur plains have considerable potential to provide insight into past ice mass history, meltwater dynamics, and sedimentation, over recent geological time scales (e.g. the Holocene).
The sandar of Iceland represent the currently unrealized potential to ascertain the controls on large-scale proglacial sedimentation. The terrestrial history of Icelandic glaciations and their impact on the evolution of the surrounding landscape has largely been based on interpretation of recent events from sedimentary exposures and geomorphological features within the near surface (<20 m). Relatively little is known about the longer-term Holocene glacial history of southeast Iceland, particularly the Late Holocene to Little Ice Age (LIA) (AD 1300-1900) glacial history, which remains poorly constrained, and hypothesized neoglacial activity (1.5-4.2 ka) remains to be tested (Gudmundsson 1997;Geirsd ottir et al. 2019). The evolution of the ice-marginal environments of Iceland since the LIA has been extensively studied. This Iceland-focused work has driven the development of the active temperate glacial land system Evans & Twigg 2002;Chandler et al. 2020), which has been used to inform interpretation of ice-marginal dynamics and processes for Quaternary glaciation (e.g. . Icelandic-derived glacial land system models, however, include little information about the large-scale sedimentary architecture of the proglacial environments. This limits their application to the understanding of the impact of changes in base-level, accommodation space and sediment flux for sandur evolution. It is widely assumed that Icelandic sandar are the product of multiple high magnitude j€ okulhlaups over the last~10 ka following rapid deglaciation from the Younger Dryas maximum (Maizels 1991(Maizels , 1993a(Maizels , 1997Guðmundsson et al. 2002). Sandur systems prone to large magnitude episodic j€ okulhlaups (e.g. Iceland) produce a distinctive facies succession that is dependent upon the nature of the j€ okulhlaup (e.g. exponentiallyrising or linearly-rising stage hydrograph shape) and the available sediment supply (Maizels 1989(Maizels , 1993b(Maizels , 1997Rushmer et al. 2002;Roberts 2005;Rushmer 2006Rushmer , 2007. In contrast, sandar dominated by relatively low magnitude seasonally-produced meltwater (i.e. ablationcontrolled) tend to produce a suite of sedimentary facies driven by channel migration and the formation of large gravel bars and in-channel large-scale dunes (Allen 1983;Miall 1985;Bristow 1996;Marren 2005).
The large-scale sedimentary architecture of sandar environments can be characterized, and its potential unlocked, by borehole or geophysical observations (e.g. Kurja nski et al. 2021). However, such observations constraining the subsurface of Icelandic sandar are sparse and are typically based on shallow near-surface (up to 10 m depth) sediment exposures. A lack of widespread geophysical observations is surprising, given that Icelandic sandar are highly suitable for groundpenetrating radar (GPR) surveys due to limited vegetation coverage, the lack of soil development and the 'GPR-friendly' physical properties of outwash sediments (i.e. sand and gravel ;Neal 2004;Reynolds 2011). These factors are conducive to extensive and deep geophysical characterization of the subsurface. Todate, the majority of GPR surveys conducted in Iceland have focused on specific landforms (e.g. Kjaer et al. 2004;Carrivick et al. 2007;Burke et al. 2008Burke et al. , 2009Burke et al. , 2010aBlauvelt et al. 2020), and surveys targeting the sedimentary architecture of recently emplaced proglacial j€ okulhlaup deposits (e.g. Russell et al. 2001bRussell et al. , 2006Cassidy et al. 2003;Harrison et al. 2019). The influence of glacier margin fluctuations, accommodation space and base-level changes on sandur-wide sedimentation cannot be extracted from the spatially restricted nature of existing studies.
Poor knowledge of the large-scale subsurface architecture of sandur environments limits our ability to constrain the driving mechanisms (e.g. glacier growth and decay, accommodation space, sediment supply, relative sea level) in how these vast depositional systems are formed and evolve during glacial cycles. Therefore, detailed analysis of the sedimentary architecture of modern large-scale unconfined outwash systems and how these reflect important drivers such as glacier recession, sea-level change and episodic sediment flux from enormous high-magnitude glacier outburst floods (j€ okulhlaups) is paramount. Characterization of sandur systems in varying topographic and glaciological settings is crucial to aid the development of more sophisticated models of sandar sediment architecture (Marren 2004).
This paper aims to: (i) evaluate the impact of glacial fluctuations and glacial sedimentation on sandur formation; (ii) establish the primary controls on the development of sandar sedimentary architecture; and (iii) develop a conceptual model for the formation and evolution of sandar environments. To fulfil these aims we characterize the large-scale sedimentary architecture of a large-scale active sandur in southeast Iceland using extensive low-frequency GPR surveys.
Seismic refraction investigations of Skeiðar arsandur suggested a three-layer subsurface sequence of bedrock, overlain by over-consolidated sediments and unconsolidated sediments of a glacifluvial origin (Guðmundsson et al. 2002). The volume of sediment above bedrock across the sandar was estimated at 100-200 km 3 (Guðmundsson et al. 2002). The uppermost unconsolidated sediments, thought to be Holocene in age, increase in thickness from 80 m in proximal regions to 150 m in distal regions of the sandur, and have an estimated volume of up to 100 km 3 , equating to a volumetric sedimentation rate of~1 km 3 per century (or 0.01 km 3 a À1 ) (Guðmundsson et al. 2002). Such estimations of sediment volume and sedimentation rates at Skeiðar arsandur should be viewed with caution, however, due to the limited spatial extent of the seismic data collected and the lack of any borehole constraints on the nature of the sediments (Guðmundsson et al. 2002). J€ okulhlaup sediment yield during a large outburst flood at Skeiðar arsandur in November 1996 was 1.8910 8 m 3 (≥3.1910 11 kg) equating to a sediment flux of 1.8910 6 kg s À1 during the 47-hour-long flood Roberts 2005;Russell et al. 2006). The Holocene sedimentation rate proposed by Guðmundsson et al. (2002) roughly equates to five 1996-equivalent j€ okulhlaups per century (Guðmundsson et al. 2002;Russell et al. 2006).
The topographic relief of the western part of Skeiðar arsandur is dominated by two moraine complexes: (i) the large push-moraine complex best defined where it is dissected by the G ıgjukv ısl river herein referred to as the 'G ıgjukv ısl' moraines; and (ii) the Sandg ıgur moraines located 2 km further down-sandur, which comprise a set of linear, ice-margin parallel, isolated ridge fragments (Fig. 1).
Skeiðar arj€ okull is believed to have reached its Little Ice Age (LIA) maximum position in 1784, 1857, 1871 and 1890-95 (J ohannesson 1985). The first three of these LIA maximum positions coincide with glacier surges in 1787, 1857 and 1873 (Bj€ ornsson et al. 2003) with the 1890-95 position corresponding with the maximum LIA extent across the southern and southeastern margins of Vat-naj€ okull (Bj€ ornsson & P alsson 2008). Þ orarinsson (1939) suggested that the Sandg ıgur moraines represent the maximum LIA extent of Skeiðar arj€ okull in~1750 (Fig. 1C). Þ orarinsson (1939), however, does not provide evidence to support this interpretation, which is at odds with historic accounts suggesting the LIA maximum extent of Skeiðar arj€ okull in the vicinity of the G ıgjukv ısl moraines was reached in 1784, 1857, 1871and 1890-95 (J ohannesson 1985Grove 2004). Sedimentary evidence suggests that at least three glacier advances were required to form the G ıgjukv ısl moraine complex (Molewski & Olszewski 2000) providing further support to the notion that G ıgjukv ısl moraines represent the maximum LIA extent of Skeiðar arj€ okull. Given that the G ıgjukv ısl moraines represent the LIA maximum, the Sandg ıgur moraines have to be pre-LIA in age. That the St oralda moraine at nearby Sv ınafellsj€ okull is ascribed to a neoglacial (4.5-1.5 ka) ice margin (Gudmundsson 1997) supports the interpretation of the Sandg ıgur moraines as representing a pre-LIA neoglacial frontal position of Skeiðar arj€ okull.
The Sandg ıgur moraines are broken into six segments which have a relatively subtle surface geomorphic expression (typically 125 m wide and 7 m high; Fig. 1C). The total extent of the moraine fragments is 3.5 km with each segment ranging between 150-650 m long. Information regarding the formation of the moraines is limited, with no detailed geomorphic or sedimentological investigations present in the literature.

Material and methods
A 40-MHz UTSI Groundvue 7 (GV7) low-frequency ground-penetrating radar (GPR) system (Francke & Utsi 2009;Harrison et al. 2019) was used to collect radar profiles (~42 km of along-track measurements) of the large-scale subsurface architecture (Fig. 1B). The GV7 has a centre frequency of 40 MHz and an emitted bandwidth of 10-80 MHz. Further technical specifications of the GV7 system are reported in Ross et al. (2018). Coincident 100-MHz GPR profiles (using a separate GV7 system) were also acquired collected at specific locations ( Fig. 1B) to provide higher resolution subsurface information. Positional data were acquired with a handheld GPS streamed into the GPR recording software. WARR (wide angle reflection and refraction) surveys suggest radio-wave velocities between 0.13-0.19 m ns À1 for the sediments at the Sandg ıgur moraines. We therefore apply a conservative velocity of 0.13 m ns À1 for time depth conversion, which corresponds with the maximum velocity for regions comprised of unsaturated sands and gravel (Neal 2004;Reynolds 2011). This is supported by the observation of the spring line intersecting the sandur surface in more medial-to-distal reaches of the sandur plain at an elevation of~37 m a.s.l., well below the elevation of our survey (50-80 m a.s.l.). Based upon our radio-wave velocity and the centre frequency of the 40 MHz GPR we achieve a theoretical vertical resolution of 0.8125-1.625 m. This limits our ability to resolve sediment structures and/or units less than 1.625 m in thickness.
Radargrams were processed using standard processing steps for low-frequency GPR data ( Fig. 2) and interpreted in Reflexw geophysical software (http://www.sandmeier-geo.de/reflexw.html). Topographic information used for elevation correction during the processing of GPR profiles was extracted from Arc-ticDEM tiles (Porter et al. 2018) to radar profiles trace headers, after having been converted from WGS 84 ellipsoid to metres above sea level (m a.s.l.) (Fig. 3). Processed data were then imported into OpendTect seismic interpretation software for 3D visualization (https://www.opendtect.org/). Radar stratigraphical descriptions of radar reflections are taken from Neal (2004). Picks of large-scale bounding surfaces were imported into ArcGIS where the Topo-to-Raster tool was used for digital elevation model (DEM) generation. Topo-to-Raster uses the ANUDEM algorithm, which is considered the most appropriate when aiming to preserve stream networks and ridge lines during interpolation (Hutchinson 1989;Arun 2013). DEM differencing between surface and subsurface elevation models was then conducted and volumetric calculations were performed.

Results
The GPR data (Figs 3, 4) reveal that the Sandg ıgur moraines have a complex and chaotic stratigraphy, with the large-scale architecture characterized by strong continuous radar reflections (Fig. 4). The radargrams indicate numerous different radar facies that are linked to different depositional styles. Radar reflections were identified at depths of up to 100 m. Radar profiles presented here are representative of the large-scale architectural features and radar facies within the subsurface across the Sandg ıgur moraines survey area ( Fig. 1). Three major contacts and/or boundaries that delineate changes in depositional setting were identified across the GPR survey and are referred to as reflections RS0, RS1 and RS2.

Radar surfaces (RS)
The RS0 reflection is the lowermost major interface that we identify in the GPR data. This reflection could not be traced in all radargrams and is most prominent where it is located at shallower depths in more distal areas of the survey (Fig. 4). In radargrams where RS0 is identifiable it is continuous and coherent. The RS0 reflection reaches a maximum elevation of À30 m a.s.l. (depth of 80 m), and a minimum elevation of À50 m a.s.l. (depth of 90 m) (Fig. 4E). The reflection is of a low amplitude and is subhorizontal in lines collected in a down-sandur (NE-SW) orientation. Lines across the sandur (NW-SE) show that RS0 has a variable morphology. For example, in the most distal radargram (Fig. 4E), RS0 is characterized by a subhorizontal reflection that transitions into a broad (~750 m) trough or depression 1.25-2 km along the radargram. A~120m-wide and~12-m-high ridge is located within the depression 1.4 km along the survey line. The 3D form of both the trough and the ridge are poorly constrained due to a lack of intersecting radargrams.

830
Devin Harrison et al.
The most extensive bounding reflection within the survey area is RS1. This reflection was identified in all profiles. The extent and continuity of RS1 vary between profiles, however; it is strongest and most coherent in radargrams near the Sandg ıgur moraines (e.g. Fig. 4B). This bounding surface reaches a maximum elevation of 52.4 m a.s.l. (depth of~10 m) and a minimum elevation of À15.4 m a.s.l. (depth of~59 m). The morphology of RS1 is characterized by a prominent ridge that coincides with the surface of the Sandg ıgur morainic ridges (Fig. 4B). RS1 is often truncated in radar data acquired between the moraine ridges and is capped by chaotic reflections (Fig. 4C). The peak elevation of RS1 occurs beneath the Sandg ıgur moraines, lowering to both the north and south of the moraine. Analysis of the radargrams reveals that RS1 is an integral part of the Sandg ıgur moraines, with the reflection continuing beneath the surface expression of the moraine and extending up and down-sandur beyond it (Fig. 4B). Further down-sandur RS1 varies from a mostly subhorizontal reflection to a more hummocky topography, particularly in the more eastern sections of the radar survey (Fig. 4D).
Reflection RS2 is a coherent, high amplitude, continuous reflector that is traceable in all profiles immediately surrounding the moraine ridges ( Fig. 4A-C). The reflector onlaps onto parts of the lower bounding reflection RS1 (Fig. 4B) and in some places extends contiguously down-sandur through gaps in the moraine where the RS1 reflector cannot be traced to the surface (Fig. 4C, D). The reflection intersects the sandur surface~1 km downsandur of the Sandg ıgur moraines (Fig. 4D). RS2 cannot be interpreted as the water table as it truncates lower reflections and where it pinches out to the sandur surface no surface water is present.

Sediment volume and thickness
DEMs of the extensive bounding surfaces RS1 and RS2 provide estimates of sediment thickness and sediment volume for locations surrounding the radar survey. Average sediment thickness above RS1 is 43 m, with a maximum thickness of 68.5 m and minimum of 7.9 m (Fig. 5C). Sediment thickness is similar both up-sandur and down-sandur of the pronounced ridge (e.g. location of the Sandg ıgur moraines; Fig. 5C). However, upsandur of the Sandg ıgur moraines RS1 is located at an elevation of 10-30 m a.s.l., whereas RS1 is mostly below 0 m a.s.l. in the down-sandur direction (Fig. 5). The relief of the subsurface ridge associated with RS1 is 40-50 m and extends in width for~1.25 km (Fig. 5). The volume of sediment between the sandur surface and the RS1 boundary is approximately 1.04 km 3 over an area of 24.09 km 2 (Fig. 5). This equates to 0.043 km 3 of sediment per 1 km 2 .
The average sediment thickness above RS2 is 14.9 m, with a maximum of 50.4 m and a minimum of 1.7 m (Fig. 6). The morphology of RS2 is fan-shaped and generally dips in a NW-SE direction, parallel with the Sandg ıgur moraines (Figs 4A, 6), with its apex located at or close to the SW2 radargram. The slope of the long profile of the RS2 bounding surface is relatively low (0.15°). The thickness of sediment bounded by RS2 and the sandur surface increases along the long profile of RS2, ranging from 10 m at the fan apex to 17 m at the end of the profile (Fig. 6D). In profiles perpendicular to the Sandg ıgur moraines (Fig. 6E, F) the thickness of this unit decreases from~30 m at the LIA limit to~10 m in the proximity of the Sandg ıgur moraines. This trend of decreasing sediment thickness in a down-sandur direction is evident from the elevation difference DEM produced for RS2 (Fig. 6C). The volume of sediment above RS2 is calculated at 0.13 km 3 over an area of 9.45 km 2 where the RS2 surface is visible in the radar data (Fig. 6). This equates to 0.014 km 3 per 1 km 2 .
The volume estimates we generate are highly dependent on the radar velocity of 0.13 m ns -1 being accurate for the sediments around the Sandg ıgur moraines. Using a radar velocity of 0.08 m ns -1 , similar to GPR surveys conducted elsewhere in Skeiðar arsandur (Burke et al. 2008(Burke et al. , 2010a, would alter the volume estimates by roughly 38%. For example, at a radar velocity of 0.08 m ns -1 sediment volume above RS1 would be 0.64 km 3 . Therefore, our estimates of sediment volume above RS1 and RS2 should be taken as maximum estimates.

Radar facies (RF)
Radar facies 1 (RF1). -RF1 (RS1 to RS2) is the most dominant radar facies identified. The entirety of this facies is located above RS1 (Fig. 7). Reflections within RF1 down-sandur of the Sandg ıgur moraines are generally more chaotic and less continuous (Table 1). Upsandur of the Sandg ıgur moraines, reflections are more continuous with fewer diffractions. In radargrams orientated NE-SW (e.g. SW3 and SW4) numerous reflections can be seen onlapping onto the RS1 surface (Figs 4B, 7). Smaller structures (<75 m long and <5 m thick) are present between the more continuous reflections and have a more wavy/undulating geometry (Table 1). However, the low-frequency nature of the 40-MHz GPR system means the true geometry of smallscale structures cannot be fully resolved in these data. Trough-like structures are evident within RF1 and are orientated in both NE-SW and NW-SE directions (Fig. 8). The largest of these trough structures is located down-sandur of the Sandg ıgur moraines in radargram SW7 (Figs 7, 8). This trough, roughly 700 m wide and 25 m deep (Fig. 8), is orientated transverse to the active fluvial systems of the S ula and G ıgjukv ısl. Possible trough structures are also present in radargrams SW1 and SW2 and are spatially coincident with gaps in the Sandg ıgur moraines.
Radar facies 2 (RF2). -RF2 is a localized facies that is found beneath the RS1 surface. It is only found directly beneath the Sandg ıgur moraines (Fig. 7). In most radargrams, RF2 is in contact with the RS1 surface and sandur surface; however, where there are gaps in the Sandg ıgur moraines the RF2 unit is overlain by a unit of chaotic and planar reflections that can be up to 20 m thick (e.g. Fig. 4B, C). RF2 is characterized by large hummocks and concave reflections and is mostly discontinuous, with offsets apparent between reflections (Table 1). Pockets of low amplitude reflections are located throughout RF2 and in some instances have enhanced levels of noise.
Radar facies 3 (RF3). -Facies RF3 is another localized facies. It is positioned below the RS1 reflection and is only recorded in the SW1 and SW3 radargrams (e.g. Fig. 7). RF3 ranges in elevation from 0 to 35 m a.s.l. (depth of 15-50 m). In down-sandur (NE-SW) orientated radargrams (e.g. Fig. 4B) RF3 is dominated by strong, coherent, oblique reflections (Table 1). These reflections dip at an angle of 1.29-2.27° (Fig. 4B) and are truncated by the RS1 bounding surface. When observing RF3 and its related reflections at the intersection of SW1 and SW3 it can be suggested that the radar facies has a conical 3D form.
Radar facies 4 (RF4). -Radar reflections associated with facies RF4 are mostly confined to the SW8 radargram. This facies exists below the RS1 boundary (Fig. 7), which has a hummocky morphology. Reflec- tions within this sediment package generally dip upsandar in a NE-SW orientation (Table 1). Only one line provides a cross-section of this radar facies. In perpendicular orientation RF4 is dominated by discontinuous reflections that have a convex to planar geometry and onlap onto one another.
Radar facies 5 (RF5). -Reflections associated with RF5 are located below the RS1 surface (Fig. 7), with an elevation ranging between À20-20 m a.s.l. The thickness of RF5 is not easy to establish due to the lack of a lower bounding surface in most radargrams. However, in radargram SW7 the RF5 facies is uniform and~10 m thick bounded by RS1 and RS0. RF5 is dominated by a low-amplitude signal with sparse reflections (Table 1). Where reflections are visible, they are moderately continuous but incoherent.
Radar facies 6 (RF6). -Very low amplitude, incoherent, oblique clinoforms (Table 1) can be observed in radargrams SW3 and SW8 at elevations of À30-À60 m a.s.l. (Fig. 7). These structures dip at a moderate angle (2.19-6.87°) down-sandur (NE-SW orientation). The basal contact of this facies is not visible. This is likely due to attenuation of returned radar power as a result of the increased depth of the reflections. It is important to note that our observed dip angles may differ from the true dip of reflections due to changes in water content within the subsurface, especially close to or below present-day sea level.

Radar facies interpretations
RF1 is interpreted as glacifluvial sediments (Table 1) related to outburst flood sedimentation or proximal fan delta deposition (e.g. Jakobsen & Overgaard 2002;Hansen et al. 2009). Diffractions within RF1 are caused by the presence of boulders buried in the sandur. The large trough-like structures present within this unit are interpreted as palaeo-channels. RF2 is an indicator of glaciotectonic deformation that is linked to push moraine formation (e.g. Overgaard & Jakobsen 2001;Jakobsen & Overgaard 2002). Large hummocky and concave reflections indicate anticlinal folds and offset reflections that may represent faulting of the sediments ( Table 1). The presence of this facies directly beneath the Sandg ıgur moraines and within the pronounced ridge of RS1 (Figs. 4B, 7) supports this interpretation. Low amplitude lenses present within RF2 represent a homogenous material. This may be buried ice or finegrained massive deposits. The uniform stacking and oblique nature of the reflections associated with RF3 suggests these are stratified deposits that represent an end-moraine/ice-contact fan ( Table 1). The dip angle of the reflections and possible conical 3D-form is comparable to similar end-moraine fans reported elsewhere in Iceland (e.g. Russell et al. 2001a;Kjaer et al. 2004). The interpretation of RF4 is difficult due to the lack of radargrams coincident with the unit, and the depth at which it is recorded. However, RF4 does show radar characteristics similar to the internal structure of drumlins elsewhere in Iceland (e.g. Woodard et al. 2020). The low amplitude nature of RF5 makes interpretation difficult. Based upon its reflection characteristics, and location beneath RS1 and the morainic facies (i.e. RF2), however, it could represent either subglacial material or an overridden fluvial deposit. RF6 is interpreted as prograding foresets (Table 1) associated with Gilberttype delta formation and sub-aqueous deposition (e.g. Hansen et al. 2009). It is important to stress, however, that the very low amplitude nature of the RF6 reflections and the depth at which they are recorded mean that we  Table 1. have lower confidence in this interpretation. Furthermore, the likely increase in groundwater content at depth may have caused reflection distortion within the RF6 radar facies.

Discussion
The landscape evolution of the Sandg ıgur morainic system can be divided into three specific phases: (i) the pre-moraine landsystem, (ii) formation of the Sandg ıgur moraines, and (iii) aggradation of the sandur plain and burial of the moraine ridge. Based on the radarderived sediment architecture and radar facies of the study area we present a schematic model of landscape evolution (Fig. 9) and assess its significance and implications.

Pre-moraine formation
Radar reflections and radar facies located at depth and distal from the present-day ice margin are difficult to confidently interpret due to a combination of their No radar reflection (NRR).

No visible reflections.
Signal attenuated due to increased depth.

BOREAS
stratigraphical position and the attenuation of the radar signal with depth. Reflections located at depths of 70-100 m should be taken with caution as this is typically viewed as the limit of GPR penetration depths within sands and gravels (Neal 2004;Reynolds 2011). As such interpretation of these deeper reflections makes reconstruction of events prior to the formation of the Sandg ıgur moraines less certain. However, in the right conditions broadband GPRs with high pulse repetition frequencies, such as the UTSI GV7, have achieved penetration depths of up to 120 m (e.g. Francke & Utsi 2009). Furthermore, if buried ice is present within the subsurface, penetration depths would be improved due to the increased radio-wave propagation through low loss frozen materials (Neal 2004;Reynolds 2011). The RS0 surface represents the base of the major reflections and radar units recorded within the GPR survey. RS0 is likely to reflect a major radar-properties (i.e. dielectric) boundary that likely marks a change in depositional setting. The change in reflection characteristics and brightness at depth could indicate that the RS0 surface is a sediment-bedrock boundary. However, similar reflection characteristics can occur between sediment boundaries due to changes in sediment composition, orientation, compaction and the resultant changes in porosity of the substrate (Neal 2004). The Fig. 9. Schematic model of the evolution of the Sandg ıgur morainic system from formation (A), to burial (B), and eventual erosional breaches (C). BOREAS elevation of RS0 is similar to that of a transition between an unconsolidated sediment package and a unit of overconsolidated sediment reported from seismic refraction surveys (Guðmundsson et al. 2002). The presence of reflections below RS0 suggests that it is not bedrock and that it is more likely the over-consolidated sediment unit reported by Guðmundsson et al. (2002).
The progradational foresets within the RF6 radar facies inferred from the down-sandur dipping reflections ( Table 1) provide evidence of deposition into a body of standing water. Deposition related to the formation of radar facies RF6 (Table 1) in a marine environment is favoured over a proglacial lacustrine environment. Proglacial lake formation would require a topographic high distal of the ice margin (e.g. morainic debris, landslide debris or less commonly glacifluvial sediments) to act as a dam (Carrivick & Tweed 2013), and there is a lack of evidence within the radargrams of a topographic high down-sandur of the reflections within RF6. If the oblique reflections characteristic of RF6 are foresets deposited in a marine environment it would suggest a former relative sea level (RSL) limit of approximately À30 m (excluding any effects of GIA) and a former shoreline located beneath the sandur in the immediate vicinity of the Sandg ıgur moraines (Fig. 4B). The potential RSL limit indicated by the reflections within RF6 is comparable to the minimum RSL limit during the Holocene regression (8.9 ka BP) in western and southwest Iceland (Norðdahl & P etursson 2005;Norðdahl et al. 2008). This would support the widely held belief (e.g. Þ orarinsson 1974; Russell & Marren 1999;Guðmundsson et al. 2002;Marren 2002a, b;Russell et al. 2006) that the majority of Skeiðar arsandur was formed during the Holocene, with development beginning at 8.9 ka BP and persisting until present-day. However, this inference should be taken with caution as the depth of the reflections makes confident interpretation difficult.
RF4 is located down-sandur of the Sandg ıgur moraines (Fig. 7) but it is difficult to determine its stratigraphical relationship to the moraine ridge. The up-sandur dipping nature of the reflections in the RF4 package and the hummocky nature of the surface above them are similar to those identified within recently formed drumlins elsewhere in Iceland (e.g. Woodard et al. 2020). If these are drumlins, they must have formed prior to the formation of the Sandg ıgur moraines and would signify a former subglacial environment. Alternatively, these features may be representative of large sand waves deposited in a coastal environment (e.g. Allen 1980) or antidune trains associated with high magnitude outwash events with supercritical flow regimes (e.g. Alexander et al. 2001;Cassidy et al. 2003). Large gravel antidune and/or dune trains preserved in the sedimentary record would require formation and subsequent burial during the same j€ okulhlaup event to facilitate preservation (Marren & Schuh 2009). Similar features have already been reported in other locations on Skeiðar arsandur (Marren & Schuh 2009).

Moraine formation
The RS1 surface is interpreted as a major boundary between: (i) an ice-marginal/subglacial dominated depositional system; and (ii) a glacifluvial proglacial depositional system. The morphology of this surface has a prominent subsurface ridge spatially coincident with the surface expression of the Sandg ıgur moraines (Fig. 4B). The radar data therefore demonstrate that the Sandg ıgur moraines continue into the subsurface and are a much larger landform than indicated by the present surface topography; with the majority of the moraine buried beneath the surface of Skeiðar arsandur. Morainic features have been identified at depth beneath sandur plains elsewhere in Iceland (e.g. Haraldsson & Palm 1980) and in Greenland (e.g. Storms et al. 2012). The maximum height of the Sandg ıgur moraines, including surface and subsurface expressions and based on the radio-wave velocity of 0.13 m ns -1 , is 67 m (Fig. 5). The GPR observations show that the moraine has a maximum width of 1.25 km (e.g. Fig. 4B). These dimensions are comparable to the surface expression of the LIA G ıgjukv ısl moraines (Fig. 1), suggesting that the Sandg ıgur moraines represent a glacial terminus position of Skeiðar arj€ okull of comparable scale and significance. The nature of moraine formation and related processes can be inferred and hypothesized (Fig. 9) from the RF2 and RF3 radar facies. RF3, located on the distal side of the moraine, is interpreted as glacifluvial sediments deposited within an ice-contact outwash fan (Table 1). This suggests that the Sandg ıgur moraines were formed during a temporarily stationary, or slowly oscillating, ice-marginal phase, as an ice-contact outwash fan of this scale would require a long period (years to decades) of time to form (Kr€ uger 1997;Krzyszkowski & Zieli nski 2002;Kjaer et al. 2004;Lukas 2005;Reinardy & Lukas 2009). The full extent of the ice-contact fan inferred from RF3 is unknown. This is because the reflections associated with the fan are truncated by the high amplitude radar reflection of RS1 (e.g. Fig. 4B). This radar stratigraphical relationship is consistent with the RS1 surface representing anunconformity, and hence evidencing erosion of the most distal sections of the (RF3) fan.
The interpretation of RF3 as an ice-contact fan suggests that deformation structures identified within RF2 could be linked to the collapse of the proximal part of the fan during ice-margin retreat, leaving a proximal rectilinear ice-contact face that creates a marked morphological asymmetry (Lukas 2005). However, the relative symmetry of the moraine ridge recorded by the RS1 surface does not support this. It is more likely that RF2 is evidence of glaciotectonic processes and deformation of the sediment from small oscillations of the glacier margin (Fig. 9); such processes and structures are commonly associated with terrestrial ice-contact fans and dump moraines (Lukas 2005). The large hummocky nature of the reflections and offsets between the reflections within RF2 are comparable in form and scale to folded and faulted sediment units generated by glaciotectonism (e.g. Jakobsen & Overgaard 2002;Benediktsson et al. 2010). The absence of deformation within the end-moraine fan unit suggests that there was no major readvance or reoccupation of the moraine following glacial retreat from its position at the Sandg ıgur moraines (Lukas 2005) and that fan development was a result of syntectonic deposition. The presence of lenses of lowamplitude reflectivity within the RF2 unit (e.g. Fig 4B) could indicate the presence of buried glacier ice. Radar unit lenses with few and low-amplitude reflections have been interpreted as pockets of sediment-free massive ice (e.g. buried glacier ice; Moorman et al. 2003). Buried glacier ice is preserved at depth in sandur in southern Iceland (Everest & Bradwell 2003), and has been identified in more-recent morainic sections at Skeiðar arsandur (e.g. Molewski & Olszewski 2000;Russell et al. 2001a). The presence of buried ice at depth in the Sandg ıgur moraines cannot therefore be ruled out.
Moraine development during a glacial surge is another possible glaciological scenario for the formation of the Sandg ıgur moraines. Skeiðar arj€ okull is a known surgetype glacier that is documented to have previously surged seven times, the earliest surge recorded occurred in 1787 and the most recent in 1991 (J ohannesson 1985;Bj€ ornsson et al. 2003). The surges of Skeiðar arj€ okull were predominately associated with the western and central portions of the glacier (Bj€ ornsson et al. 2003) with the limit of the most recent surge (c. 1991) marked by a push moraine and an associated fine-grained, stratified outwash fan (Russell et al. 2001a; van Dijk & Sigurðsson 2002;van Dijk 2002). Push moraines formed during glacier surges are dominated by heavily deformed and glaciotectonized sediments (Russell et al. 2001a;Benediktsson et al. 2010), which is consistent with the presence of RF2 within the Sandg ıgur moraines. The outwash fans associated with the surges of Skeiðar arj€ okull have similarities to other ice-contact/endmoraine fans in the literature (Russell et al. 2001a; van Dijk 2002; van Dijk & Sigurðsson 2002). As a result, RF3 is also consistent with a surge model  for moraine formation (Fig. 4B). However, aside from the potential evidence from the moraine ridge itself, no other evidence suggesting moraine formation by glacial surge is apparent in the GPR data.

Burial of the Sandg ıgur moraines and sandur aggradation
The large proportion of the Sandg ıgur moraines that currently lies beneath the sandur surface indicates that there has been significant sandur aggradation post moraine formation. Sediments above the morainebounding RS1 surface were deposited after the formation of the Sandg ıgur moraines (Fig. 9). The estimated volume of the sediment package above the Sandg ıgur moraines (1.04 km 3 ), is comparable to the sediment yield from ten 1996-equivalent j€ okulhlaups (Guðmundsson et al. 2002;Snorrason et al. 2002). The volume of sediment above the Sandg ıgur moraines represents~1% (1.04 km 3 ) of the Holocene sedimentation (100 km 3 ) approximated by Guðmundsson et al. (2002), accumulating in 1.9% (24.09 km 2 ) of the total sandur area (1300 km 2 ) assuming a radio-wave velocity of 0.13 m ns À1 .
The majority of the sediment above RS1 is interpreted as glacifluvial in origin (Fig. 7). Glacifluvial deposition is the dominant process of sedimentation at Skeiðar arsandur (e.g. Þ orarinsson 1974; Boothroyd & Nummedal 1978;Guðmundsson et al. 2002) and GPR reflection characteristics within RF1 (Table 1) are consistent with sedimentary architectures (e.g. chaotic and planar to wavy moderately continuous to discontinuous reflection geometries) formed byglacifluvial deposition reported in other geophysical investigations (e.g. Beres et al. 1995Beres et al. , 1999Jakobsen & Overgaard 2002;Hansen et al. 2009). The stratigraphy down-sandur of the moraine ridge (i.e. RF1) is dominated by chaotic reflections with numerous diffractions. Chaotic reflection geometries such as these suggest the presence of reworked and poorly sorted, structureless and matrix supported sediments indicative of deposition from sediment-rich flows during the rapid rising stage of a j€ okulhlaup (Brennand 1994;Russell & Knudsen 1999;Marren & Schuh 2009). During rapid rising-stage flood flows in an unconfined proglacial setting there is insufficient time for sediment sorting and development of well-defined bedding to occur (Church & Gilbert 1975;Smith 1986;Miall 1992;Hiscott 1994;Todd 1996;Russell & Marren 1999;Russell et al. 2002). This results in high rates of aggradation and deposition, and the presence of a poorly sorted structureless deposit (Rushmer 2006(Rushmer , 2007. Rapid emplacement of sediment during a j€ okulhlaup is the most likely depositional mechanism of this facies due to the burial and preservation of the interpreted icecontact/end moraine fan (RF3). Some erosion of the fan is evident (e.g. high amplitude nature of the RS1 reflection and truncation of reflections within the fan body). However, the majority of the fan is preserved at depth beneath the interpreted j€ okulhlaup deposits. The lack of any clear boundaries in the sediments that characterize RF1 makes it difficult to determine whether this unit is the product of a single or multiple j€ okulhlaup.
Glacifluvial sediments on the up-sandur side of the moraine ridge below the RS2 surface (i.e. RF1) are typically characterized by less chaotic and more moderately-continuous, planar-to-wavy radar characteristics, than that seen in down-sandur radargrams. Small channel fills and reflections linked to migrating bedforms characterize the unit below RS2 and are predominately orientated parallel (NW-SE) to the Sandg ıgur moraines (Fig. 9B). This suggests that sediment was deposited by meltwater systems that flowed parallel to the palaeo-ice-margin in an iceproximal depression (Fig. 9B), which is similar to the present-day drainage system of Skeiðar arj€ okull. The dip angle (0.15°) and orientation of RS2 (e.g. Fig. 6D) indicate that sediment was deposited from flows running parallel to the Sandg ıgur moraines. The presence of a lateral drainage network is common in landsystem models for retreating ice margins (Price 1969;Bogacki 1973;Galon 1973a, b;Klimek 1973;Gomez et al. 2000). Meltwater drainage parallel to the ice margin requires the lowering of the glacier bed below the elevation of the iceproximal sandur, thus creating a topographic discontinuity that prevents meltwater from draining onto the sandur surface and diverts meltwater to flow parallel to the ice margin . This indicates that the burial of the up-sandur portion of the moraine ridge occurred after a period of down-sandur aggradation once the glacier had retreated from the Sandg ıgur moraines (Fig. 9B).
Sedimentation above RS2 is thought to have occurred when Skeiðar arj€ okull reached its LIA maximum (Figs 1C,9C). The surface geomorphology suggests that fans and fluvial channels emanated from the LIA push moraines and continued to just beyond the Sandg ıgur moraines. In addition, the tapering thickness of sediment above RS2 is consistent with outwash fan formation (Fig. 6E, F). Where reflections are visible within the sediments above RS2 there is evidence of some down-sandur dipping reflections proximally, and more planar continuous reflections distally. This is consistent with outwash fan architecture (Krzyszkowski & Zieli nski 2002). The push moraines associated with the LIA margin ( Fig. 1) are believed to be formed as a result of multiple glacier surges (Molewski & Olszewski 2000;Bj€ ornsson et al. 2003). Large sediment-rich outburst floods have been known to occur at the end of surge events as a result of the reorganization of the subglacial drainage network (Kamb et al. 1985;Bj€ ornsson 1998;Russell et al. 2001a;van Dijk 2002; van Dijk & Sigurðsson 2002;Eisen et al. 2005). Furthermore, the post-surge quiescent phase can allow for the development of supraglacially-fed outwash fans (Russell et al. 2001a;van Dijk 2002; van Dijk & Sigurðsson 2002). The sediments above RS2 could therefore represent surge related outwash fans, which are contrasting to j€ okulhlaup-related outwash fans due to the absence of both heavily kettled topography and large-scale bar and channel patterns (Russell et al. 2001a). The presence of surge related outwash fans in front of the LIA moraines has already been reported from a large, exposed sediment section at the G ıgjukv ısl gap (Russell et al. 2001a). The sediments above RS2 and the potential surge outwash fans at the G ıgjukv ısl gap section are separated by the incised G ıgjukv ısl meltwater system and it is therefore plausible that they may have previously been connected. Outwash fans related to surge events would likely be regularly spaced along the ice margin and form laterally extensive and coalesced belts of outwash fans (Russell et al. 2001a). If the fans above RS2 are surge-related it would add further support to the suggestion that glacial surges are a significant control in the development and formation of glacifluvial landforms and deposits (Russell et al. 2001a; van Dijk 2002; van Dijk & Sigurðsson 2002).
The discontinuous surface expression of the Sandg ıgur moraines and coincident subsurface features suggest that the moraines have been breached and eroded by later outwash floods. The erosion of the Sandg ıgur moraines likely occurred during the period of sedimentation above RS2 and during occupation of the LIA moraines. However, some sections of the moraine may have been eroded prior to the formation of the RS2 surface, as the moderately continuous reflection in the SW5 radargram below RS2 can be seen to truncate a buried portion of the Sandg ıgur moraines (RS1) at depth (Fig. 4C). Therefore, timing of the moraine breaches cannot be conclusively determined. Historic accounts suggest that j€ okulhlaups in the 18th century eroded portions of the Sandg ıgur moraines (Grove 2004). It is thought that these outburst floods were caused by the drainage of the N upsl on ice-dammed lake (Þ orarinsson 1939;Grove 2004). This scenario is plausible in leading to sedimentation behind the Sandg ıgur moraines and the erosional breaches. However, historical accounts show that there have been numerous high-magnitude j€ okulhlaups that have inundated the entire sandur (Þ orarinsson 1974) that would have been capable of eroding portions of the Sandg ıgur moraines.

Implications for sandar evolution
Sandar sedimentary architecture is largely controlled by variations in the glacier margin position (e.g. Marren 2004). Fluctuations of the ice margin influence the regions at which sandar aggradation or incision can take place (Fig. 10). During periods of glacier advance and stability the ice margin is coupled to the extramarginal sandur environment and there is an abundant sediment supply (Fig. 10), promoting continued aggradation beyond the morainic limits. Where climate-driven lowmagnitude high-frequency meltwater regimes dominate stratified ice-contact fans are formed (e.g. Fig. 9A). However, a switch to a high-magnitude low-frequency dominated meltwater regime, for example as a result of increase in j€ okulhlaup prevalence, can result in major down-sandur aggradation of thick j€ okulhlaup units (e.g. Fig. 9B; Marren 2005). As the glacier retreats accommodation space is created, leading to the activation of intramarginal geomorphic activity and aggradation during periods of increased sediment supply (e.g. j€ okulhlaup events). Sediment aggradation within a confined sandur region by ice-margin parallel glacifluvial activity in an ice-proximal depression, following retreat from the Sandg ıgur moraines (Fig. 9B), is coincident with landsystem models for retreating ice margins (e.g. Price 1969;Bogacki 1973;Galon 1973a, b;Klimek 1973;Gomez et al. 2000).
Sediment supply to the proglacial environment peaks during the initial stages of ice-margin retreat (Antoniazza & Lane 2021) promoting aggradation within the confined intramarginal zone (Fig. 10). Furthermore, during periods of rapid glacier recession sediment flux (Fig. 10) to the proglacial outwash environment increases (Marren & Toomath 2013). This has important implications for how sandar evolve during periods with varying rates of glacier margin recession. As a result of higher sediment supply during rapid glacier retreat (e.g. Marren & Toomath 2013), accommodation space between the newly established ice margin and the coupled end-moraine and sandur plain is filled. This leads to the sandur backfilling into available accommodation space as the ice margin continues to retreat (Fig. 9B). However, it is important to note that thiswill onlyoccur where there is a high sediment flux available to the proglacial environment and/or sediment-laden j€ okulhlaups are prevalent. This process of backfilling and aggradation into accommodation space left by the ice mass following retreat can lead to ice-margin parallel draining networks that create complexity within the stratigraphy compared to what may be evidenced in the surface or shallow subsurface of the sandar (e.g. Marren 2004). However, when the rate of retreat reduces, particularly when combined with proglacial lake development and a continued reduction in accommodation space, incision will dominate within the extramarginal (unconfined) sandur environment (Marren & Toomath 2013;Fig. 10). This is predominately the result of the lowering of the Fig. 10. Conceptual diagram showcasing the dominant controls on the evolution of sandar environments and how they vary in response to glacier margin advance and retreat. BOREAS proglacial long profile, due to the decrease of meltwater outlet elevation (e.g. Marren & Toomath 2013), and the reduction of sediment supply to the extramarginal sandur due to ice -marginal proglacial lake systems acting as an efficient sediment trap (e.g. Carrivick & Tweed 2013).
The burial of the Sandg ıgur moraines by thick and extensive glacifluvial deposits demonstrates that high rates of sedimentation, likely the result of j€ okulhlaup events, are a dominant driver in the formation and aggradation of Skeiðar arsandur. Sea level in southern Iceland is known to have been relatively stable since 6 ka BP (Norðdahl & P etursson 2005). Given the relative stability of the sea level we can infer that in this study base level is not a dominant control. Instead, system-scale sandur evolution is driven by ice-proximal changes. High sedimentation rates and aggradation depths above the Sandg ıgur moraines (i.e. average sediment thickness of 43 m) suggest that sediment supply and the available terrestrial accommodation space are more dominant drivers of sandur evolution (Fig. 10).
The identification of sediment supply as a primary control on sandur formation is consistent with numerical modelling of the Holocene evolution of Skeiðar arsandur (e.g. Maizels 1993a). It is likely that large increases of sediment supply are a result of an increased magnitude and frequency of j€ okulhlaup drainage over the sandur. The presence of an unstable proglacial fluvial system as a result of these extreme flood events promotes geomorphic activity and sedimentation in both extra-and intramarginal portions of the sandur environment (Carrivick & Heckmann 2017;Weckwerth et al. 2021), dependent on the state of the glacier margin at the time of sedimentation (Maizels 1979;Marren 2002aMarren , b, 2005Russell et al. 2006) (Fig. 10). Drainage from ice-dammed lakes, in the form of j€ okulhlaups, is known to increase in frequency during periods of ice-margin retreat (Þ orarinsson 1939;Evans & Clague 1994). Therefore, glacier margin retreat from its position at the Sandg ıgur moraines and thinning may have contributed towards the formation and drainage of ice-dammed lakes in the region (e.g. Þ orarinsson 1939; Bj€ ornsson 1992; Grove 2004; Roberts et al. 2005). Furthermore, it is believed that the main calderas of Gr ımsv€ otn became active after the 1783 Laki eruption thus increasing the frequency of volcanogenic j€ okulhlaups in the region that could deliver sediment to the entire sandur (Þ orarinsson 1974;Guðmundsson 1992).

Conclusions
• We identify a large, partially buried, moraine system attributed to the Holocene or neoglacial, terminus position of Skeiðar arj€ okull outlet glacier. The Sandg ıgur moraines are an isolated remnant of a much more extensive moraine system that would have continued eastwards across Skeiðar arsandur. • J€ okulhlaup-related glacifluvial sediments dominate the sediment architecture around and overlying the Sandg ıgur moraines. Erosional breaches of the Sandg ıgur moraines were the result of either j€ okulhlaup flows emanating from the large LIA push moraines, drainage of ice-dammed lakes in the region, or large volcanically-induced j€ okulhlaups that inundated the entire sandur (e.g. Þ orarinsson 1939, 1974Grove 2004;Roberts et al. 2005). • Our GPR data indicate that in proglacial areas prone to high levels of sediment supply, and rapid sedimentation, there is potential for ice-marginal landforms to be preserved in the subsurface. This has important implications for the preservation potential of sandur systems in the geological record. However, preservation potential of ice-marginal landforms remains low where sedimentation has not sufficiently buried the landform. • The stratigraphical architecture of Skeiðar arsandur beneath the Sandg ıgur moraines system includes reflections indicative of progradational foresets (RF6). This suggests deposition into a body of standing water, such as the sea, at an elevation of approximately À30 m a.s.l., prior to subaerial moraine formation and glacifluvial deposition.
Acknowledgements. -This work was supported by a PhD studentship awarded to Devin Harrison through the IAPETUS Natural Environment Research Council Doctoral Training Partnership (NE/L002590/ 1). This research was carried out under RANNIS research declaration number 11/2019. We thank Vatnaj€ okull National Park for permitting the collection of GPR data at Skeiðar arsandur. We acknowledge the ArcticDEM project (Porter et al. 2018) and DEMs provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691 and 1542736. We are grateful to Will Smith for his assistance in the field during data collection in June 2019. We are thankful to Piotr Weckwerth and David Sharpe for providing helpful, supportive and constructive reviews that considerably improved the paper.
Author contributions. -DH and NR collected the GPR data used in this study. DH processed and analysed the data. DH wrote the paper with input from all co-authors.