Controls on jökulhlaup-transported buried ice melt-out at Skeiðarársandur, Iceland: Implications for the evolution of ice-marginal environments

The full-text may be used and/or reproduced, and given to third parties in any format or medium, without prior permission or charge, for personal research or study, educational, or not-for-pro t purposes provided that: • a full bibliographic reference is made to the original source • a link is made to the metadata record in DRO • the full-text is not changed in any way The full-text must not be sold in any format or medium without the formal permission of the copyright holders. Please consult the full DRO policy for further details.


Introduction
Glaciers are frequently used as indicators of climate change as they respond dynamically to changes in the climate driven components of their mass balance. Knowledge of the former extent of glaciers can be used to reconstruct palaeo-climate and to define the former position of contemporary glaciers. Distinctive assemblages of landforms and deposits at modern glacier-margins have stimulated the development of models which can be used to reconstruct ice-marginal processes. Models such as the ice-marginal landsystem (e.g. Krüger, 1994;Evans and Twigg, 2002;Evans et al., 2019) were developed from the detailed investigation of contemporary glacier margins and have been used to reconstruct palaeo-glacier margins in the Quaternary record (e.g. Evans et al., 1999). By definition, many studies of contemporary icemarginal processes provide only a snapshot of the geomorphological and sedimentological evolution of an ice-marginal zone, as they do not take account of post-depositional landscape modification processes and are not always able to constrain the influences of previous landsystems on landscape evolution. Palimpsest landscapes comprising a number of superimposed landsystems allow landform overprinting and the potential for landsystem legacy to persist for periods of 10 1 -and Stroeven, 1997; Schomacker and Kjaer, 2007;Korsgaard et al., 2015).
Ice-marginal and proglacial geomorphology can be modified by a number of post-depositional processes such as aeolian deflation and deposition, fluvial erosion and deposition, periglacial and paraglacial slope processes (e.g. Ballantyne, 2002;Mountney andRussell, 2006, 2009). Melt-out of buried glacier ice has been well-documented from modern ice-margins (Price, 1969;Schomacker and Kjaer, 2007;Tonkin et al., 2016) as well as being interpreted from the Quaternary record (Eyles et al., 1999;Fard, 2003). Buried ice melt-out has been invoked to account for a number of distinctive landforms such as 'kame and kettle' topography and 'hummocky moraine', both of which result from topographic inversion as buried ice melts, causing slow collapse of overlying sediment (e.g. Everest and Bradwell, 2003;Lukas et al., 2005;Bennett and Glasser, 2009).
Despite widespread acknowledgement of the importance of buried ice within former glacier margins, relatively little attention has been paid to the process of ice emplacement and how this may determine buried ice distribution and melt-out styles. Studies of buried ice melt out also tend to have focused on the immediate ice proximal areas of proglacial outwash plains. Similarly, there are scant data on the rates of ice melting beneath thick debris mantles; exceptions include McKenzie (1969) and Schomacker (2008). Buried ice is known to have survived for decades to centuries (e.g. French and Harry, 1990;Evans and England, 1992;Everest and Bradwell, 2003;Schomacker, 2008) however there have been few detailed studies of melt-out rates over these timescales.
Processes occurring in the ice-marginal zones of glaciers and ice sheets are complex. Subaerial processes re-work glacially-deposited debris and the melt-out of buried ice leads to collapse structures and topographic inversion (Price, 1969;Bennett and Glasser, 2009). Even a thin (N0.01 m) layer of debris covering glacier ice can provide sufficient insulation to retard ablation (e.g. Lister, 1953;Østrem, 1959;Young, 1981, 1982;Nicholson and Benn, 2006) and ablation can be very slow under thick debris mantles. Very slow melt rates can therefore permit the survival of buried glacier ice for long periods of time. The sustained collapse of overlying sediment due to buried ice meltout is a significant post-depositional modification process in deglaciated landscapes with ice-cored topography (Ballantyne, 2002). The correlation between climatic parameters and melt rates of buried ice bodies is weak however, suggesting that both burial processes and topography play a key role in the rates of ice melting (Nicholson and Benn, 2006;Schomacker, 2008).
Glacier ice can also be buried 'in situ' by supraglacial sediment deposition on top of an active or stagnant glacier margin (e.g. Russell and Knudsen, 2002;Schomacker and Kjaer, 2007;Schomacker et al., 2006). During jökulhlaups, ice blocks become detached by englacial hydrofracturing, meltwater conduit collapse and ice cliff collapse (Roberts et al., 2000;Roberts, 2005). Ice-blocks up to 10 2 m in diameter are known to have been washed from glacier margins by jökulhlaups on to outwash plains (sandar) and subsequently either partially or completely buried by sandur aggradation (Tómasson, 1996;Russell and Knudsen, 1999;Roberts et al., 2000;Fay, 2001Fay, , 2002aFay, , 2002bRoberts, 2005;Russell et al., 2006). Melt out of ice blocks transported by the 2010 Eyjafjallajökull eruption generated jökulhlaups is thought to have played a major role in post depositional landscape evolution (Harrison et al., 2019).
During the November 1996 jökulhlaup,~8.3 × 10 6 m 3 of ice was removed from Skeiðarárjökull (Fay, 2002a(Fay, , 2002bRussell et al., 2001aRussell et al., , 2005. Sections of the snout of Skeiðarárjökull were fractured in situ into blocks up to 200 m × 400 m in size (Roberts et al., 2000). Ice blocks as large as 45 m in diameter were transported from the glacier margin (Fay, 2002a(Fay, , 2002bRussell et al., 2005). Some of these ice blocks were deposited as a linear jökulhlaup flow-parallel cluster, resulting in a single coalesced kettle hole approximately 130 m wide and 40 m long (Fay, 2002a;Knudsen, 1999, 2002). The largest accumulation of ice blocks was over 1 km in length with a width of up to 300 m (Fay, 2002a;Russell and Knudsen, 2002). Ice blocks transported and buried by such single high-magnitude jökulhlaups are known to persist for 10 1 -10 2 yr −1 (e.g. Fay, 2002aFay, , 2002bEverest and Bradwell, 2003;Russell et al., 2005).
The aims of this paper are to: (1) determine the origin of a series of actively developing depressions within the proglacial area of  (Thórarinsson, 1974). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 3. Proximal to distal profiles of the glacier and the sandur, demonstrating the retreat of the margin and the base level lowering of drainage within the proglacial depression, and assumed lowering of groundwater table. The 1997The , 2003The , 2007 profiles are derived DEMs from imagery acquired by Landmaelingar Íslands, Loftmyndir ehf. and NERC ARSF, respectively. The 2018 profile is derived from ArcticDEM. Skeiðarársandur, southeast Iceland; (2) evaluate the mode and significance of buried ice emplacement and subsequent melt for depression development; and (3) explain their significance for sandur evolution. To fulfil these aims we quantify the decadal evolution of the depressions and characterise their sub-surface sedimentary architecture and relate to the wider record of jökulhlaups on Skeiðarársandur.

Study area
Harðaskriða is located 3.6 km from the current margin of Skeiðarárjökull within the central zone of Skeiðarársandur, a 1300 km 2 outwash plain fed by Skeiðarárjökull in southeast Iceland (Fig. 1). Skeiðarárjökull is a temperate, surge-type, outlet glacier of Vatnajökull ice cap with a 23 km wide piedmont snout (Björnsson, 1998). Skeiðarársandur has a strong maritime climate, where the maximum depth of winter freezing is only of the order of centimetres (Douglas and Harrison, 1996;Thórhallsdóttir, 1996) making the presence of permafrost impossible. Skeiðarárjökull and Skeiðarársandur have been subject to repeated high-magnitude jökulhlaups generated both by subglacial volcanic eruptions and the drainage of subglacial and ice-marginal lakes (Thórarinsson, 1974;Björnsson, 1992Björnsson, , 1997. The Harðaskriða area of Skeiðarársandur last experienced meltwater flow during the 1922 jökulhlaup, after which glacier margin recession created incised proglacial channels at the Háöldukvísl and Gígjukvísl; subsequently a proglacial trench developed, allowing all subsequent meltwater to be routed westward (Galon, 1973) (Figs. 2, 3 and 4; Table 1). The Harðaskriða area is now part of an elevated sandur surface which was unaffected by the large November 1996 jökulhlaup (Snorrason et al., 2002). The Harðaskriða area was however inundated by eleven large jökulhlaups between 1861 and 1938 when Skeiðarárjökull was at its Little Ice Age maximum extent (Thórarinsson, 1974;Björnsson, 1997;Glaciorisk, 2005) (Table 1). These high frequency high-magnitude jökulhlaups resulted in significant aggradation on the eastern and central proximal areas of Skeiðarársandur with accumulated elevations of up to~125 m above sea level (a.s.l.) compared with elevations of below 90 m a.s.l. for equivalent-aged sandur surfaces to the west (Blauvelt, 2013). Historic accounts indicate a number of large jökulhlaups associated with glacier margin disruption and ice block release in the Harðaskriða area between 1861 and 1938 (Thórarinsson, 1974;Glaciorisk, 2005) (Fig. 2; Table 1). The 1897 and 1903 jökulhlaups can also be specifically linked to the Harðaskriða area (Thórarinsson, 1974). A 1 km long, 150 m high piece of the glacier margin was washed out on to the sandur during the 1903 jökulhlaup possibly associated with detachment along a large, ice flow transverse, hydrofracture generating a large 'embayment' (Roberts et al., 2000;Roberts, 2005).
The Harðaskriða comprises outwash surfaces characterised by a well-developed channel and bar pattern supporting numerous ice block obstacle marks up to 10 m in diameter (Fig. 4). These outwash surfaces are however disrupted by four large depressions the largest of which has maximum dimensions of 604 × 108 m and a depth of 13 m (Figs. 5, 6 and 7). Although there are historic accounts of high magnitude jökulhlaup processes and immediate impacts at Harðaskriða there has been no detailed examination of the development of the resulting landforms and deposits.

Topographic survey and DEM generation
DEMs for 1968 and 2007 were produced using stereophotogrammetry on digital aerial imagery. Vertical aerial photographs collected for DEM creation were obtained from the National Land Survey of Iceland, Landmaelingar Íslands. Images were scanned by Landmaelingar Íslands with an Eversmart Jazz+ Scitex scanner, at a resolution of 2000 dpi and delivered as tagged image format (tif) files. Camera calibration documentation and flight lines drawn on 1:100,000 topographic maps were provided for all photographs taken after 1954. Photographs from 1945, taken by the U.S. Air Force, were purchased from Landmaelingar Íslands, who provided known flight elevations and camera focal lengths. Colour digital images acquired in 2007 by NERC ARSF (IPY07/13) were also used in this study. Digital elevation models (20 m resolution) and photo mosaics were purchased from Loftmyndir HF for 1997 and 2003 datasets.
A differential GPS (dGPS) survey was conducted in 2007 between the glacier margin and Iceland's ring road (Fig. 1). Transects of four of the Harðaskriða depressions were surveyed (Fig. 7). Additionally, large-scale, persistent features across Skeiðarársandur including kettle holes, boulders and ridges that were visible on all historical aerial photographs were utilised for ground control points (GCPs). Survey points were collected using a Thales ProMark III unit, corrected to Icelandic Roads Authority survey sites. BAE's SocetSet 5.5 (Ngate) software was used to generate the DEMs. Triangulation (interior and exterior orientation) was accomplished for all photosets. Once an internal coordinate system was established within the photographs, the control points measured in the field could be used to relate the image to the ground (absolute orientation). The ISN93 coordinates of the GCPs were used to identify points on the images. Once the x, y and z values of GCPs were identified on both (or more) images, SocetSet then performed point measurement automatically, using digital image matching (Baily et al., 2003).
Systematic errors and random errors were evaluated by comparing apparent elevation differences between the DEMs and ground control points measured with the dGPS. The location of the check points and ground control points are summarised in Table 2. Systematic errors are given as root-mean-square error (RMS) measures and the 95th percentile limit is given for random errors a technique commonly used in DEM quality analysis (Schiefer and Gilbert, 2007). All units are in metres above sea level (m a.s.l.). Following manual clean up ('post pushing') of the study area, all check points fell under 1 m. Due to the comparatively small scale and dynamic terrain of the study area, no additional registration was applied. Table 1 Chronology of glacier margin fluctuations and jökulhlaups  from the drainage of Grímsvötn subglacial lake draining from Skeiðarárjökull. Information in this table is sourced from Thórarinsson (1974), Björnsson (1997), Glacier risks database (2005)  Little Ice Age Maximum One of the largest jökulhlaups in river Skeiðará with a duration of four days. Peak discharge reached after two days accompanied by tremendous noise. The next day ice blocks were covering the sandur all the way to the ocean. Ice blocks were up to 20 m in diameter and very tightly packed. South of the main flood outlet, icebergs, up to 20 m high, covered a 7 km wide region. Mud was spread over the flooded area and was unusually thick. Melting icebergs and quicksand made it difficult to cross the Skeiðarársandur for several months after the flood.

Yes
Little Ice Age Maximum This jökulhlaup was smaller than the one in 1892 and had a duration of 10 days with a 6 day rising stage. It burst from Skeiðarárjökull near Harðaskriða, a steep moraine south of central Skeiðarárjökull. Ice blocks up to 20 m in diameter were spread over a 6 km wide area between Harðaskriða and the Skeiðará.
(Potential source of jökulhlaup transported ice blocks to the Harðaskriða area).

Yes
Large jökulhlaup with a duration of 4 days reaching discharge peak very quickly and covering more of the western outwash plain than usual. Ice blocks were carried all the way to the ocean. A large piece of the glacier margin detached during the 1903 jökulhlaup approximately 1 km in length and up to 150 m in height; it was also documented that a fracture of similar size and length developed up glacier located where the floodwaters burst from the glacier margin.
(Potential source of jökulhlaup transported ice blocks to the Harðaskriða area).

No
Recession Large jökulhlaup of 12 days duration. The flood was focussed on the eastern part of Skeiðarársandur. Large amounts of ice detached from snout of Skeiðarárjökull with ice blocks described as being the size of houses. 1922 Yes Recession This large jökulhlaup had a duration of 14 days and had its main outlet on the eastern side of Skeiðarársandur. Rising stage discharge increased slowly for 6 days before any ice blocks were observed being transported downstream. Followed by 8 days of recession. The jökulhlaup waned within one day.
(Potential for aggradation in the Harðaskriða area).

Yes
Recession A large jökulhlaup with a volume of 4.5 km 3 with a 16 day duration and a peak discharge of 25,000 to 30,000 m 3 s −1 . The eastern most outlet was 2.5 km wide and ice blocks carried to ocean.
(Potential for aggradation in the Harðaskriða area).

Yes
Recession A large jökulhlaup with a volume of 4.7 km 3 , a peak discharge of 25,000 to 30,000 m 3 s −1 and a 16 day duration. Almost all Skeiðarársandur was flooded with ice blocks covering the sandur, although they were smaller than those generated by the 1934 jökulhlaup.
(Potential for aggradation in the Harðaskriða area).

No
Recession A small jökulhlaup with a duration of 17 days and a volume of 1.4 km 3 . Relatively small ice blocks released from the glacier snout.
(No impact on the Harðaskriða area due to the formation of a proglacial trench diverting flow in a westward direction).

No
Recession A jökulhlaup with a duration of 12 days, a volume of 2.6 km 3 and a peak discharge of 10,000 m 3 s −1 . Jökulhlaup flowed in the Skeiðará, Sandgígjukvísl and Núpsá river and produced relatively small ice blocks.
(No impact on the Harðaskriða area due to the formation of a proglacial trench diverting flow in a westward direction).

No
Recession A jökulhlaup with a duration of 17 days, a volume of 2.2 km 3 and a peak discharge of 5000 m 3 s −1 .
(No impact on the Harðaskriða area due to the formation of a proglacial trench diverting flow in a westward direction).

No
Recession A jökulhlaup with a duration of 14 days, a volume of 3.2 km 3 and a peak discharge of 10,000 m 3 s −1 . Jökulhlaup flowed in the Skeiðará and Sandgígjukvísl with peak discharges of 6000 m 3 s −1 4000 m 3 s −1 respectively.
(No impact on the Harðaskriða area due to the formation of a proglacial trench diverting flow in a westward direction).
Elevation differences were used to provide an approximate estimate of the volumetric loss for four of the major Harðaskriða depressions (Fig. 8) and the elevation loss of the adjoining glacier (Fig. 3). Whilst the poor quality of the 1945 images precluded the production of a 1945 DEM surface to quantify subsidence over the last sixty-two years, an attempt was made to provide as close an approximation as possible. Comparing surfaces constructed from two subsequent time periods (before and after) is often utilised as a cost-effective method to quickly quantify large-scale volumetric changes due to melt out, subsidence, flooding, human interference or other causes (Schiefer and Gilbert, 2007). By removing elevation points that lay within the depressions on the 2007 imagery and generating a triangular irregular network, or TIN, across the missing data points, a 1945 DEM surface could be simulated.
Although the 2018 ArcticDEM was used to provide a general comparison of glacier recession and proglacial fluvial system incision (Porter and 28 others, 2018) (Fig. 3), offsets between the 2018 ArcticDEM and photogrammetrically-derived DEMs, as well as the difficulty in removing bias in this type of terrain, precluded its use for quantification of rates of lowering of the Harðaskriða depressions.

Ground Penetrating Radar survey
Ground-Penetrating Radar (GPR) profiles were collected from depression 4 (Figs. 5 and 6) using the MALÅ ProEx system with a 13 m long (distance Tx -Rx = 6 m), low-frequency (30 MHz) Rough Terrain Antenna (RTA). GPR-lines (all corrected for topography), both outside and across the depression, were collected in 2013 (i.e. 6 years after the GPS surveys). The objectives were to gain insight in the subsurface sediment architecture, deformation or collapse structures, and to investigate the possibility of the presence of remnants of buried ice.
The basic principles of GPR surveying are that electromagnetic waves travel at different velocities dependent on the electrical and magnetic properties of the earth materials, and that incident waves are refracted or reflected on interfaces between materials with contrasting dielectric permittivities. The nature of the signal that is returned to the surface (i.e. its intensity, polarity and propagation velocity) can then be analysed using processing software (ReflexW; cf. Sandmeier, 2012) which allows the reconstruction and modelling of the architecture of the subsurface. Fig. 6 shows the survey plan with a single E-W profile capturing the length of the depression, and three shorter N-S profiles across the depression. Data were collected as a continuous array with additional transects surveyed to connect the long and cross profiles (the radar profiles outside the depression were used for reference only). Assuming an average propagation velocity of c. 0.07 m ns −1 (cf. Cassidy et al., 2003), depth penetration in the sandur sediments was c. 30 m. This is a value typical for velocities in wet, sandy to gravelly materials and may be an underestimation in case of significant quantities of buried ice present in the subsurface.

Morphology of Harðaskriða depressions
All depression measurements are based on their 2007 dimensions. Depression 1 is approximately oval in shape and ranges in width from 164 m (north-south) to 108 m (east-west) (Figs. 7 and 8). The northern and southern rims are characterised by outwardly dipping arcuate and concentric normal faults. Along the southern rim, normal faulting has resulted in the rotation of two large blocks (up to 60 m in length and 15 m wide). The base of this depression is characterised by sagging, uneven terrain, and is divided into two portions of unequal depth. The northernmost part of the depression measured 10 ± 1.64 m in depth, while the southern part of the depression measured 12 ± 1.64 m in depth. Recent satellite imagery (DigitalGlobe, 2016) indicates further depression widening attributed to continued melt-out.
Depression 2 is approximately circular in shape and ranges in width from 89 m (north-south) to 104 m (east-west) and 13 ± 1.64 m deep (Fig. 7). The northern portion contains concentric normal faults and two extensional faults that trend north-south (40 m and 50 m in length). Along the southernmost rim of this depression, normal faulting has resulted in the rotation of two blocks, the largest 60 m long and 13 m wide.
Depression 3 possesses an irregular, elongate morphology that trends east-west. The depression ranges in width from 342 m (eastwest) and 116 m (north-south) and is 12 ± 1.64 m deep (Fig. 7). The margin, while not circular in shape, contains numerous normal faults and recesses that surround the depression. The southern margin is marked by several rotated blocks and steep walls. The margin appears to slump in rotational blocks towards the centre, resulting in 'steps' that dip outwards from the depression. A dirt road observed on the 1945 aerial photographs remains visible on the 2007 aerial photographs. Its original surface, although now undulating, remains discernible as it traverses depression 3, suggesting that subsidence has been gradual in nature.
Depression 4, the widest of the depressions, is similar in shape to depression 3, possessing an irregular shape and trending east-west (Fig. 7). The depression ranges in width from 604 m (east-west) to 148 m (north-south). Similar to the other depressions, the walls are steepest along the southern margin, and the margin is characterised by horst and graben blocks and concentric extensional fractures (Figs. 5a, b and 6). Numerous recesses have developed along the northern, eastern and western margin and possess a relatively gentle, stepped slope, compared to those along the steeper southern margin.
Directly 800 m north of these depressions and visible on the 1965 photographs are three drumlinised, elongate ridges that lead to the elevated sandur (Fig. 9a). These ridges trend north-south and, from east to west, are 176 m, 79 m and 151 m in length and 30, 37, and 34 ± 1.64 m in height respectively. Elevation profiles were extracted of the area immediately adjacent to these ridges (Profile 1) and the prominent ridges (Profiles 2-4) that rose towards the elevated sandur surface (Fig. 9b). While varying in height, the ridges all span an average of 33 m from the base of the depression to the sandur. Later photo series  indicate that these ridges were largely removed by the fluvial erosion of shifting proglacial

Sub-surface structure of Harðaskriða depressions
Using a relatively low radar frequency, and assuming that the subsurface sediment mostly comprises sand and small to medium gravel (with no outsized boulders to cause 'disruptive' hyperbolae) it may be expected that the signal-to-noise ratio is adequate for resolving metre-scale sandur sedimentology down to a depth of c. 30 m.
There are two main sub-horizontal reflectors in the GPR crossprofiles: one at an estimated depth of 6 m and one at c. 20 m below the surface (white solid lines in Fig. 11). A third discontinuous reflector is visible just above the noise which starts at 30 m. As it is right at the detection limit, interpreting this reflector will not be attempted below. The 6 m reflector tends to mirror the surface topography and the 20 m reflector is generally less undulating and continuous across the depression. All three profiles also show shorter sub horizontal reflectors that are thought to represent prominent bedding surfaces. Northward dipping reflectors, a single one in the central cross-profile and two parallel features in the east cross-profile (white dashed lines in Fig. 11), appear to extend down from the 6 m reflector, cross and deflect the 20 m reflector and then connect in a stepped fashion with the reflector at 30 m depth. (See Fig. 12.)  In all three cross-profiles, high-angle linear or curvilinear structures (indicated in red in Fig. 11) intersect, or terminate onto the aforementioned reflectors. They are interpreted as joints or normal faults. Particularly near the margins of the depression, they can be seen to disrupt or offset reflectors. Most structures are outward dipping, but there are also less common, apparently younger, inward-dipping faults. All such fractures seem to have developed to accommodate the flexure in the depression and are attributed to progressive subsidence due to gradual melt-out of buried ice. At the surface around the periphery of the depression, the structures present as stepped ring-structures (Figs. 5a, b and 6), which are very similar to the 'concentric' ring-fractures described for collapsing calderas and other ice-melt phenomena by Branney (1995) and Branney and Gilbert (1995).
Whilst confident about the interpretation of the joint and fault structures, the characterisation of the subsurface materials from the radargrams is more challenging, particularly without the possibility of direct ground truthing. There are good exposures near the Gígjukvísl, 5 km kilometres to the west (see Russell et al., 2001aRussell et al., , 2001b, but they are developed into proglacial surfaces which lack the 'pitted' surfaces diagnostic of jökulhlaup deposits. Fortuitously, the 1996 jökulhlaup cut a 30 m high section into sandur sediments 1.3 km north of the Harðaskriða depressions (Fig. 10) so there is at least some information available on textural heterogeneity of the sandur sediments and the distribution and dimensions of buried ice.
Apart from the aforementioned reflectors, the most conspicuous zones in the GPR profiles are those that seem to be devoid of energy returns. Such zones (indicated in blue in Fig. 11), can be observed mostly away from the centre of the depression and below the 6 m reflector, but there are also a few at greater depths. Assuming that these features are not processing artefacts, they must represent homogeneous materials with a low relative dielectric permittivity ε r . Where a reflector -mostly that at 6 m depth, forms the upper surface of such zones -the polarity is opposite to that of the air wave (phase change of 180°) which would suggest that the overlying sediment has a relatively high ε r . Since ice has an ε r of 3-4 (Brandt et al., 2007) and overlying materials, which can logically assumed to be relatively (wet) jökulhlaup sediments may be expected to have an ε r in the order of 10-30 -and drawing analogies with nearby exposures -the zones are tentatively interpreted as remnants of buried ice.
The observation that the interpreted blocks of buried ice are ubiquitous at the north and south sides of the cross-profiles, but less common in the central parts is compatible with the idea that subsidence has been greatest in the centre of the depression. The inward dipping reflector on the south side in the central and east cross-profiles (Fig. 11, middle and lower panels) may delineate the upper surface of relatively intact buried ice. The deeper parts of this surface may have served as a slip-plane extending into one of the deeper identified normal faults.
Although its strength is variable, it is clear that the reflector at 20 m is more continuous than the 6 m reflector. Interestingly, its polarity is the same as the air wave which suggests that it  represents a contact between a lower permittivity (above) to a higher permittivity material below. The fact that it does not show significant offsets where intersected by high angle faults is taken as evidence that the reflector is not a sedimentary surface. Instead it is proposed that it represents the local groundwater table (wet/saturated sand: ε r = 10-30; Brandt et al., 2007), although it is noted that other studies on the sandur (cf. Cassidy et al., 2003;Burke et al., 2010) have found the groundwater table to be significantly shallower.

Discussion
The average rate of 19.4 ± 2.6 cm of lowering per year between 1945 and 2007 determined from this study is an order of magnitude lower than the 1.88 ma −1 reported for immediate post jökulhlaup ice-melt out within the Gígjökull basin between 2010 and 2016 (Harrison et al., 2019). Higher buried ice melt rate at Gígjökull can be attributed to the simultaneous deposition of smaller ice fragments with jökulhlaup deposits rather than the melt of large isolated blocks.
Combined DEMs, dGPS measurements and GPR surveys reveal that the Harðaskriða depressions experienced the greatest vertical loss within their centres, due to slump in rotational blocks towards the centres, characteristic of 'horst and graben' structures. This process has resulted in 'steps' that have developed along the side of each of the features into the centre. The concentric rings of normal faulting, horst and graben and normal and extensional faults described at the Harðaskriða depressions 1-4 are consistent with observations made at other field sites involving the melt-out of smaller bodies of ice that have been transported by lahars and jökulhlaups (Maizels, 1992;Branney, 1995;Branney and Gilbert, 1995;Olszewski and Weckwerth, 1999). Such features have also been used as indirect evidence of buried bodies of ice at other locations (e.g. Boulton, 1972;Hambrey, 1984;Krüger and Kjaer, 2000;Kjaer and Krüger, 2001;Dickson and Head, 2006).
As the ice bodies buried at Harðaskriða began to melt, the loss of volume and drainage of subsurface water may have resulted in the subsidence of the overlying sediment (McDonald and Shilts, 1975;Maizels, 1992). This sort of subsidence can produce outwardly-dipping arcuate hairline fractures that can elongate into a ring, causing the subsidence of a coherent block of sediment as seen in depressions 1 and 2 (Branney, 1995;Branney and Gilbert, 1995). These overhanging scarps become unstable and collapse along new arcuate faults, resulting in the development of extensional crevasses that may continue to expand along small vertical and normal faults causing some walls to collapse, resulting in keystone graben (Sanford, 1959;McDonald and Shilts, 1975). Continued collapse leads to intersection of arcuate fractures resulting in blocks that tilt and subside into the pit, while mass movements and slumping may accelerate the melting rate of a buried ice body (Johnson, 1992).
At larger collapse pits, such as depressions 3 and 4, irregular topographic margins with embayments also developed (Branney, 1995;Branney and Gilbert, 1995). These features, combined with the steep walls of the depressions and undisturbed nature of the surrounding outwash plain are consistent with bodies of ice that have been surrounded by sediment (Maizels, 1991). The gentle slopes of the northern walls and the steeper slopes of the southern walls are consistent with the development of a 'normal' kettle hole (Maizels, 1992;Olszewski and Weckwerth, 1999), as proglacial outwash would have resulted in the development of gravitational flow on the northern side, while block displacement and subsidence developed on the southern side following melt out.
The existence of large bodies of buried ice N30 m in thickness on Skeiðarársandur have been identified and documented using resistivity studies (Everest and Bradwell, 2003) and confirmed at exposures (Klimek, 1972;Bogacki, 1973;Churski, 1973;Jewtuchowicz, 1973;Russell and Knudsen, 1999;Molewski, 2000). Ridges and detached slabs of dead ice in the eastern and western parts of Skeiðarársandur have also been identified and described and are attributed to deposition by the retreating ice margin (Galon, 1973;Jewtuchowicz, 1973;Wojcik, 1973). Unlike ice-cored ridges, plains or moraines elsewhere on the sandur, the geometry, orientation and size of the bodies of ice that resulted in the Harðaskriða depressions are consistent with other descriptions of isolated blocks of ice emplaced during high-magnitude jökulhlaups (Maizels, 1992;Maizels and Russell, 1992;Branney, 1995;Branney and Gilbert, 1995;Harrison et al., 2019).
A topographic map published in 1904 (Danish Staff Map) depicts several elongated, east-west trending ridges extending across central Skeiðarársandur that appear to be continuations of the 19th century moraines that persist today in the western region of the sandur (Fig. 2). By 1945, aerial photographs reveal that these moraines are no longer visible on the central sandur, and reportedly buried or removed by jökulhlaups (Galon, 1973;Jewtuchowicz, 1973;Wojcik, 1973;Wisniewski et al., 1997;Knudsen et al., 2001). While some of the depressions and landforms correspond to the approximate positions of the 19th century moraines, the largest Harðaskriða depressions are developed approximately 400 m south of this limit, suggesting that they are not related to buried ice bodies contained within the pre-existing 19th century moraines (Fig. 4).
According to Thórarinsson (1974), a large piece of the glacier margin detached during the 1903 jökulhlaup approximately 1 km in length and up to 150 m in height; it was also documented that a fracture of similar size and length developed up glacier located where the floodwaters burst from the glacier margin. During this same flood, house-size ice blocks were emplaced on the sandur and the flood waters "dug into the sand a deep, 'many persons high', steep-sided channel" (Thórarinsson, 1974). Ice blocks, regardless of their original shapes, result in circular depressions, such as Depression 2, however dumbbellshaped pits may form where circular collapse pits from two closely adjacent buried blocks of ice overlap, such as Depression 1 (Branney and Gilbert, 1995). The geometry and orientation of the largest elongated depressions (Depressions 3 and 4) may therefore correspond to the 1 km wide portion of the margin that was detached during the 1903 jökulhlaup described by Thórarinsson (1974). In the absence of evidence of a disrupted glacier snout or ice blocks on the topographic map published in 1904, it is presumed that the field survey that formed the basis for this map pre-dated the 1903 jökulhlaup. Thórarinsson (1974) stated that jökulhlaups in 1913 and 1922 inundated the central sandur with floodwaters and sediment. During later periods of glacier stillstand, meltwater runoff was concentrated in the central part of the sandur, resulting in the formation of wide outwash channels (Galon, 1973). In common with glacier termini elsewhere in Iceland, the margin of Skeiðarárjökull experienced climate-forced recession from their Little Ice Age maximum extents (Thórarinsson, 1943;Sigurðsson, 2005). Recession of Skeiðarárjökull resulted in meltwater drainage from the glacier margin at progressively lower elevations leading to sandur incision (Galon, 1973). As such, subsequent jökulhlaups in 1934 and 1938 did not affect Harðaskriða, as the floodwaters were routed through other channels such as the Háöldukvísl, 1.5 km to the east (Fig. 1). Aerial photographs taken in 1945 show the formation of a proglacial trench and meltwater flow in a westerly direction towards the Gígjukvísl (Figs. 1 and 3).
While the jökulhlaup-transported ice bodies may have been emplaced as early as 1897 and as late as 1922, any melting that occurred during that time is not captured due to a lack of available imagery. The rate of melt of a buried ice body may be affected by a variety of factors, including the amount of sediment within the ice, depth of burial and geothermal heat flux (Nakawo and Young, 1981;Nicholson and Benn, 2006), making it difficult to estimate the initial size of the buried ice body. Ice blocks emplaced and completely buried by the 1903 jökulhlaup would have been further insulated by additional sediment aggradation during the 1913 and 1922 jökulhlaups (Thórarinsson, 1974). That glacier ice buried by November 1996 jökulhlaup deposits has survived for 23 years illustrates the feasibility of buried ice preservation between the 1903 and 1913 jökulhlaups.
It is noticeable that the Harðaskriða depressions are not visible on the 1945 photographs, suggesting that the buried ice has not exhibited high melting rates. It is not until the 1965 photographs, following the retreat of the central lobe of the glacier margin and the subsequent formation of the proglacial trench post-1945, that subsidence is visible. This observation and the sequence of events presented in this study suggests that the melt rate of the buried ice bodies may have been accelerated as a result of the retreat and decoupling of the glacier margin and the associated rise in ambient temperatures and lowering of local groundwater table (Robinson et al., 2008;Levy et al., 2015). This demonstrates the control that glacier margin stability has on post-depositional modification processes, as buried ice bodies may be capable of persisting for much longer periods at a stable or advancing margin, characterised by proglacial aggradation, rather than at a retreating or stagnating margin characterised by proglacial incision.
According to Björnsson et al. (1999) profiles of the surface of Skeiðarárjökull in 1904 were~100 m higher than in 1945, which would have resulted in a steeper ice surface gradient and therefore increased hydraulic gradient during high-magnitude jökulhlaups (Roberts et al., 2000(Roberts et al., , 2001Roberts, 2005). This would have increased the capacity of jökulhlaups to excavate and transport sediment. The elongate, drumlinised ridges observed on the 1965 images on the down-glacier side of the proglacial trench generated by the retreat of the glacier margin are interpreted as conduit-fill eskers created by sediment deposition as meltwater ascended by at least 30 m over a distance of~200 m from the proglacial depression to innundate Harðaskriða (Fig. 9).
The landform and sediment assemblage at Harðaskriða reflect the role of multiple jökulhlaups just after the Little Ice Age maximum extent of Skeiðarárjökull. Initial glacier position before the 1903 jökulhlaup is associated with unconfined proglacial drainage Knudsen, 1999, 2002;Russell et al., 2005Russell et al., , 2006 (Fig. 12a and a(i). Erosion of a 1 km wide ice-walled re-entrant into the snout of Skeiðarárjökull by the 1903 jökulhlaup liberated large ice blocks which were transported by the jökulhlaup onto the sandur for distances of up to 0.5-0.8 km (Fig. 12b). The largest 1903 jökulhlaup-transported ice blocks were probably partially buried as was the case with the largest ice blocks during the 1996 jökulhlaup (Russell and Knudsen, 1999;Fay, 2001, 2002a. Sediment aggradation during the 1913 and 1922 jökulhlaups buried the ice blocks emplaced in 1903 ( Fig. 12c(i). It is likely that the ice blocks had reduced in size by ablation between 1903 and 1913. Continued glacier recession resulted in the abandonment of the Harðaskriða sandur surface between 1933 and 1945 (Fig. 12d). Melt of buried ice results in depressions which have deepened and expanded in surface area between 1968 and 2007 ( Fig. 12d (i). The GPR survey undertaken in 2013 of the largest depression indicates the presence of buried glacier ice which together with the recent satellite observations of depression widening, suggests that the melt out processes are ongoing.

Conclusions
Continued melting of the Harðaskriða ice bodies nearly a century following their emplacement and burial demonstrates that jökulhlaups may continue to be an important control on sandur evolution over decadal to centennial timescales (10 1 -10 2 years). Buried ice meltout associated with the development of the Harðaskriða depressions was enhanced by the lowering of the groundwater table following abandonment of the sandur brought about by glacier margin recession during the second half of the twentieth century. The occurrence of three high magnitude jökulhlaups within an 18-year period following the Little Ice Age glacier maximum extent resulted in significant sandur aggradation and ice block burial, assisting the long term preservation of ice. By contrast, a similar succession of jökulhlaups during a period of glacier margin recession will reduce the potential for jökulhlaup-transported ice blocks to be buried as repeated 'decoupling' of the glacier margin from the its sandur reduces the potential for stacking of jökulhlaup deposits.
Our model of the jökulhlaup landsystem at Harðaskriða and the ability to identify them at other warm-based sediment-rich glaciers that may be subject to some or all the large-scale processes including margin fluctuations, jökulhlaup dynamics and secondary modification may provide a useful analogue for interpreting landforms and strata emplaced by margin fluctuations, jökulhlaups and melt out generated by the retreating continental Pleistocene ice sheets.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.