Late Holocene canyon-carving floods in northern Iceland were smaller than previously reported

Catastrophic floods have formed deep bedrock canyons on Earth, but the relationship between peak discharge and bedrock erosion is not clearly understood. This hinders efforts to use geological evidence of these cataclysmic events to constrain their magnitude – a prerequisite for impact assessments. Here, we combine proxy evidence from slackwater sediments with topographic models and hydraulic simulations to constrain the Late Holocene flood history of the Jökulsá á Fjöllum river in northern Iceland. We date floods to 3.5, 1.5 and 1.35 thousand years ago and confirm that flow peaks during these events were at most a third of previous estimates. Nevertheless, exposure ages suggests that nearby knickpoints retreated by more than 2 km during these floods. These findings support a growing consensus that the extent of bedrock erosion is not necessarily controlled by discharge and that canyon-carving floods may be smaller than typically assumed. Flood events which carved canyons along the Jökulsá á Fjöllum river in northern Iceland occurred 3.5, 1.5 and 1.35 thousand years ago and were at most a third of the size of previous estimates, according to stratigraphic analyses and flood simulations.

C atastrophic floods can reshape entire landscapes in a matter of hours to days. This raises fundamental questions about the genesis of fluvial landforms. For example, were these features formed gradually or by abrupt events? And how much water was involved? Geological evidence suggests that such outburst floods stem from the rapid release of glacial meltwater. During deglaciation, radiative warming triggered the build-up and drainage of glacier-dammed lakes along the perimeter of wasting ice sheets 1 . In geologically active glaciated regions, volcanic heat episodically melts large volumes of ground ice 2 .
These catastrophic events carved deep canyons that have been extensively investigated to reconstruct past flow rates 3 . Yet, after decades of research, process-form relationships remain contested. For example, recent modeling work suggests that erosion rates are also determined by other factors than discharge 4 . If true, paleoflood discharge estimates may be revised down considerably. This has major implications for the use of geological evidence of floods to better understand pressing unknowns in planetary science such as the sensitivity of Earth's climate to meltwater pulses or the extent of hydrological activity on Mars 5,6 .
The analysis of sediments from depositional zones perched above flood-carved canyons can help narrow these uncertainties. When sufficiently sheltered from channel flow velocities, rising water may pond in these sites and allow deposition of suspended sediments. If preserved during subsequent discharge peaks, such so-called slackwater deposits can record the timing and frequency of floods over thousands of years 7 . When then combined with hydraulic simulations, these flood level (paleostage) indicators can be used to constrain the magnitude of flooding. While successfully used for decades 8 , recent methodological developments have advanced the potential of slackwater deposits to reconstruct past floods and predict future flood hazards. This progress includes the use of new statistical and scanning approaches to fingerprint the signature of slackwater sedimentation with unmatched fidelity 9 . In addition, geochronological advances now allow us to precisely date flood sediments using a range of independent methods 10 . Finally, a new generation of observation-calibrated hydraulic models can be parameterized with centimeter-scale Digital Elevation Models (DEMs) to produce highly accurate peak discharge estimates 11 .
This study harnesses the potential of these advances to deepen our understanding of the links between flood magnitude and bedrock erosion by targeting slackwater sediments deposited by the Jökulsá á Fjöllum river in northern Iceland (Fig. 1). Since the last deglaciation, numerous glaciovolcanic outburst floods or jökulhlaups have extensively modified the local landscape as evidenced by large-scale erosional and depositional features 12 . Because of the regular occurrence of these events and the similarities of associated landforms with those observed on Mars 13,14 , the Jökulsá á Fjöllum catchment has been extensively studied in recent decades 2 . But despite a wealth of geomorphological, sedimentological and geochronological evidence, the timing and magnitude of past floods and their role in landscape evolution remain uncertain and much-debated [15][16][17] . This ambiguity is bestillustrated by peak discharge estimates from the upper reaches of the watershed that differ by two orders of magnitude 18,19 and references therein.
We constrain the Late Holocene (past 4.5 ka) glacial outburst flood history of the Jökulsá á Fjöllum by investigating sediments from Ástjörn-a lake that is uniquely suited to trap slackwater sediments during extreme events because of its sheltered position above the river channel 7 . In our analyzed core, we identified flood successions using a combination of physical properties (density and organic content), granulometry, Computed Tomography (CT) and elemental geochemistry (XRF). To estimate flood magnitude, we utilized grain size End-Member Modeling Analysis (EMMA), Structure-from-Motion (SfM) photogrammetry and available hydraulic simulations [20][21][22] . To date each event, we combined tephrochronological analysis of visible basaltic ash horizons with radiocarbon dates that bracket flood deposits.
We refine the Late Holocene flood history of the Jökulsá á Fjöllum by confirming that Ástjörn was flooded three times, at 3.5, 1.5, and 1.35 cal. ka BP. We also show that the peak discharge of these events was at least three times smaller than previously reported 19 . Our findings support the emerging view that the magnitude of canyon-carving floods has been overestimated.
Setting. The 200 km long Jökulsá á Fjöllum river drains the northern part of Vatnajökull, Europe's largest ice cap, before entering the Nordic Seas ( Fig. 1). Since deglaciation around 10 ka BP 23 , its watershed has witnessed multiple jökulhlaupsoutbursts of glacial meltwater that were triggered by sub-glacial volcanic eruptions 17 . These cataclysmic floods incised various channels, cataracts and canyons along the river's course 19,24 and references therein. These include the max. 100 m deep Jökulsárgljúfur gorge, which is situated in the lower reaches of the river and features a few perched lakes that can help constrain the timing and magnitude of past floods by trapping slackwater sediments 25 . We target one of these basins, Ástjörn ( Fig. 2; 37 m a.s.l.), which is a particularly promising flood stage indicator for two reasons. First, its setting: situated close to ultimate base (sea) level in channels eroded into lava flows beneath and in a catchment uniformly affected by post-glacial uplift and volcanic rifting 16,26 , the distance between Ástjörn and the river has remained stable since both features were formed during the Early Holocene 17,27 . This is a prerequisite for robust millennial-scale sediment-based flood magnitude reconstructions 28 . Also, as the lake is separated from the Jökulsá á Fjöllum channel by a min. 60 m high (108 m a.s.l.) escarpment (Fig. 3 degrees junction angle with the river channel reduces flow speed to warrant effective deposition and preservation of suspended slackwater sediments 7 . Secondly, the availability of published flood simulations that allow us to determine the discharge required to inundate the lake [20][21][22] . While approaching these levels during the largest recorded flood of 1725-1726 CE ( Fig. 2A), the lake has not been inundated during historical times 14 . The Holocene flood history of the Jökulsá á Fjöllum watershed remains loosely constrained. Exposure dates ( 3 He) of fluvially sculpted surfaces and tephra markers or radiocarbon ( 14 C) ages from overlaying soil profiles identify three broad periods with evidence of floods around 9-7, 4-3, and 2-1 ka BP 15,19,29 . In between such catastrophic events, sediment transport in the watershed is dominated by glacigenic suspended load during late summer discharge peaks of~500 m 3 s −1 30,31 . As shown in Fig. 2A, the river barely overtops its banks during these seasonal floods. In addition, prevalent katabatic southwesterly winds easily entrain silty sediment from the Dyngjusandur dust plume area to the south (Fig. 1) 32 . The lake catchment is, however, sheltered from eolian processes by a well-developed vegetation cover that is dominated by birch and willow woodland ( Figs. 2A and 3A). Ástjörn has no outlet and presumably drains through the subsurface, while the lake's only inlet is a small brook that occupies an over-dimensioned paleoflood channel that carries no sediment and enters the basin across a 10 m high headwall at its southern end (Fig. 3B) 12 .

Results and discussion
The Late Holocene evolution of Ástjörn. We selected AST-P2-18 as master core for this study owing to its greater (450 cm) length and recovery of surface sediments (see "Coring" under "Methods"). Investigation of a processed GPR profile across our coring site reveals a sharp facies change from continuous reflectors to reflection-free at 5 m sediment depth (Fig. 2D). This transition marks the sediment-bedrock boundary 33 , and suggests that we retrieved the entire infill of Ástjörn. Field observations support this notion as bouncing of the hammer weights during coring suggested that an impenetrable surface was reached 34 . Based on this evidence, we claim that the lake sits in an overdeepened bedrock basinlikely an extension of the paleoflood channel to the South ( Fig. 3B; see "Setting"). Visual logging and multi-proxy stratigraphic analysis of core AST-P2-18 reveal 9 (numbered from the top) units that comprise 4 facies (Figs. 4-6): peat beds, organic lacustrine sediment, minerogenic slackwater deposits and intercalated overbank sediments.
First, the peat deposits of units 9 (base-440 cm) and 3 (212-178 cm) whose organic character is reflected by exceptionally  high (45-100%) Loss on Ignition (LOI) percentages (Fig. 6). As evidenced by the widespread occurrence of roots and stems in both deposits, a boggy woodland occupied Ástjörn at the time of deposition. Unit 3 is bracketed by 14 C ages that suggest rapid accumulation of up to 0.5 cm yr −1 (Fig. 5): similar growth rates have been reported for other Late Holocene sub-Arctic peatlands 35,36 . Because of the sharp contacts with adjoining facies (Fig. 6), we hypothesize that abrupt reorganizations in sub-surface drainage of the basin allowed these conditions. As demonstrated by 37 , single floods can extensively modify the surface profile of Icelandic floodplains like the Axarfjörður sandur that borders Ástjörn (Fig. 3A). We note that material retrieved in our core catcher shows that unit 9 contains an additional~10 cm of sediments, which minimizes the likelihood that this peat horizon was redeposited 7 , and strengthens our confidence in the accuracy of the reported basal age of 4418-4589 cal. yr BP (see "Core chronology" and Supplementary Table 1). Second, the dark-colored and coarse-textured clastic sediments of units 8 (440-342 cm) and 4 (248-212 cm). Closer examination of CT imagery reveals intercalated, laterally discontinuous organic horizons that consist of lumps of peat and fragments of roots or stems (Fig. 6). The thickest of these max. 2 cm lenses are captured by LOI peaks and Titanium (Ti) minima. The observed alternations have also been reported in similar channel-marginal basins up-stream 38 , and are attributed to overspill flood deposition in line with evidence from other bedrock river canyons 7 . Under such circumstances, clastic sediments settle from suspension in overbank flows during discharge peaks, while light organics settle on top as the water recedes. Often, this dateable material is eroded from older deposits 39 , which may explain the outlying age LuS 15023 in unit 8 (Supplementary Table 1 and Fig. 5). Based on the parallel orientation of separate sediment beds and the fine sand-dominated size distribution of particles (Figs. 6 and 7C), we argue that both units were deposited in the lower flow regime of seasonal floods. This interpretation is supported by the near-identical grain size distribution of catchment samples from a seasonally-flooded channel on the adjacent Axarfjörður sandur (Figs. 2A and 7C). The absence of buried soils and high accumulation rates also hint at frequent inundation.
Third, distinctly colored units 7 (342-304.5 cm), 5 (257.5-248 cm), and 2 (178-89.5 cm) that range from dark brown at their base to beige toward the top (Fig. 6). The uniformly dense (DBD; max. 1.2 g cm −3 ) and minerogenic character of these sediments, reflected by high (~1.4) Total Scatter Normalized (TSN) Ti ratios and low LOI (~2%), notably set them apart from the seasonal flood deposits of AST-P2-18. In addition, mean grain size data reveal a distinct normal grading from basal fine sands to upward-fining coarse siltdominated caps in each unit (Fig. 6). In similar settings, such sequences characterize geologically instantaneous slackwater deposition during flood events 40,41 : sand first settles from inundating currents as flow velocity drops while the finer silts only fall out of suspension when waters pond. As discussed before, light organic detritus settles last and 14 C ages from this flood-eroded material may be older than the time of deposition. This could explain the anomalous ages of samples LuS 14881 (unit 7), LuS 14877 (unit 2), and LuS 15020 (unit 2), justifying our decision to identify them as outliers (Supplementary Table 1 and "Core chronology"). Finally, the light brown-colored sediments of units 6 (304.5-257.5 cm) and 1 (89.5-0 cm). Elevated (10-15%) LOI values and high scattering ratios demonstrate a high organic content, while low (~1) TSN Ti ratios suggest minimal minerogenic input. The mean grain size (~30 µm) of clastic sedimentation falls in the medium silt fraction, which is considerably finer than other core facies (Fig. 6). With the notable exception of the dark visible ash layers targeted for tephra analysis (see "Methods"), all measured parameters and core photos show that both units are structureless (massive) and homogeneous (also see Supplementary Fig. 1). Based on this sediment signature and its overlap with modern deposition, we argue that units 1 and 6 were laid down when the Ástjörn basin was occupied by a lake that did not receive clastic material from the river or sandur. By deriving low accumulation rates, the presented chronology also indicates slow background sedimentation during these intervals (Fig. 5). Minerogenic input that entered the lake during these quiescent phases was likely wind-blown: their medium silt-dominated size distribution matches that of sediment sourced from nearby Dyngjusandur -Iceland's largest source of dust 42 (Fig. 1). The katabatic southwesterly winds that prevail during the snow-free summer season frequently blow large plumes across Ástjörn 32 .
Following from the above, we argue that the Late Holocene evolution of Ástjörn was characterized by multiple sharp transitions between terrestrial, lacustrine and fluvial sedimentation (Figs. 5 and 6). Shortly after the onset of (peat) accumulation in the basin prior to 4.5 cal. ka BP (unit 9), seasonal overbank floods deposited a sequence of organic-minerogenic couplets (unit 8). The graded sandy-to-silty minerogenic sequence of unit 7 marks the first phase of slackwater deposition in Ástjörn. Following this event, organic lacustrine conditions similar to today prevailed (unit 6), until the basin was inundated again (unit 5). The subsequent two centuries were marked by rapid (~0.5 cm/ yr) overbank accretion (unit 4) and peat accumulation (unit 3). A third and final flood deposited the massive slackwater deposit of unit 2, before lacustrine background sedimentation resumed until the present (unit 1).
Core chronology. All sampled radiocarbon (n = 10; Supplementary Table 1) and tephra (n = 4) age ties were included in our linearly interpolated Clam-generated chronology (Fig. 5) 43 . 14 C ages were calibrated with the Intcal13 curve 44 . We eschewed a Bayesian approach as the abrupt shifts between stratigraphic units in core AST-P2-18 indicate highly variable sedimentation rates; this restricts the ability of such models to robustly parameterize accumulation rate priors 45 . Based on visual correspondence between piston core AST-P2-18 and gravity core AST-G2-18 (see "Coring" under "Methods"), we argue that no sediments were lost during coring; we thus assigned a zero age (2018 CE) to the core top. Basal age LuS 14882 shows that the sediment infill of Ástjörn covers the Late Holocene (past 4.5 ka). A number of inverted 14 C ages suggest contamination with reworked material: to avoid mixing carbon sources, we exclude all 14 C ages derived from bulk organic material from our chronology (Supplementary Table 1). Correlation of the analyzed basaltic tephra horizons AST 1-4 to known eruptions is therefore key to a correct identification of outliers. To do so, we characterized the major elements data (expressed as oxides) using the approaches described for Icelandic tephra by 46 . Based on minimal tailing, a homogenous layer thickness and an angular shard morphology ( Supplementary  Fig. 1), we argue all four horizons derive from primary air fall. All data, except for a sub-population (n = 7) of shards in chemically bi-modal layer AST4 ( Supplementary Fig. 3), reveal a tholeiitic basaltic composition (Supplementary Data 1): this restricts their provenance to Iceland's North (NVZ), West (WVZ) or East (EVZ) volcanic zones (Fig. 1). To correlate our horizons to active volcanoes in those zones, we relied on key discriminatory bi-plots established by 47 and 48 . These notably include TiO 2 vs. MgO to distinguish between VAK (Veiðivötn-Bárðarbunga, Askja, Krafla) and TGK (Thordarhyrna, Grímsvötn, Kverkfjöll) sources, and K 2 O vs. FeO to (better) separate a Veiðivötn-Bárðarbunga provenance from other TGK edifices (Fig. 4). Based on this assessment and with the help of Reference Compositional Fields (RCFs) for tephra from the foregoing systems 47,48 , each tephra horizon could be attributed to particular volcanoes. Even more so, with the additional support from our calibrated radiocarbon ages, we correlate each marker to a specific eruption (Fig. 4).
Horizon AST1 consists of two populations: the largest and most homogeneous (1-1; n = 15) has a Veiðivötn-Bárðarbunga affinity (Fig. 4a). While a Reykjanes Volcanic Belt (RVB) provenance cannot be excluded based on geochemical evidence, this option is ruled out by our radiocarbon-based chronology. No explosive eruptions of this system have been recorded during the last millennium, when AST1 was deposited 49 . Linear interpolation between the core top (2018 CE) and the rangefinder 14 C age at the base of unit 1 (LuS 15020; Supplementary Table 1) yields a 180 cal. yr BP age for horizon AST1 (Fig. 5). This estimate is consistent with a well-dated regional marker from Bárðarbunga; the 233 cal. yr BP V-1717 tephra 50 , which was dispersed to the North across our field area 51 . This correlation is further supported by the nearidentical geochemistry of reference material (Fig. 4b). We should note that the geochemistry of a 1477 CE eruption from the same system is near-identical; however, this age falls far outside the constraints provided by our 2018 CE core top zero age and subjacent 14 C ages (Fig. 5). While minor (n = 9) sub-set AST 1-2 is mixed, most shards show a geochemical affinity with the TGK Grímsvötn and Kverkfjöll volcanoes ( Supplementary  Fig. 2b). We favor a correlation with the former as the latter has not erupted during historical times 52 . Following from the above, we assign the reported 233 cal. yr BP age of V-1717 to AST1 50 .
With the exception of three higher silica (SiO 2 ≥ 52 wt. %) outliers, the analyzed shards from horizon AST2 (n = 27) reveal a homogeneous geochemistry that matches the composition of the Kverkfjöll volcano (Fig. 4C). This edifice has the lowest eruption frequency of all TGK volcanoes, 1 per millennium during the investigated Late Holocene 52 ; this greatly aids source identification. To achieve this, we rely on tephra data from the afore-mentioned Kárahnjúkar soil profile 52 . This deposit contains a 1325 cal. yr BP old Kverkfjöll tephra (Kári1-113) that may be deposited by the last known eruption of this volcanothe coincident Lindahraun event e.g. 53 . This age is also in broad agreement with the~1500 cal. yr BP estimate for AST2 derived from linear interpolation between the core top and radiocarbon dating sample LuS 15020 at the base of unit 1 (Fig. 5). We confirm this correlation with two lines of geochemical evidence based on AST2 and Kári1-113 (n = 6) major element glass data: (1) the highly similar values of key Kverkfjöll discriminators TiO 2 and MgO (see Fig. 4C), and (2) Similarity Coefficients (SCs) ≥ 0.95 (Supplementary Data 1), calculated on oxides with >1 wt. % (n = 7) after 54,55 . As the chronology of the Kárahnjúkar profile is well-constrained by 21 known regional tephra markers, we assign its 1325 cal. yr BP age estimate to AST2 while also applying the 250 yr uncertainty margin recommended for this record 56 .
AST3 is also geochemically homogeneous and analyzed shards (n = 42) have a composition very similar to that of the TGK Grímsvötn volcano (Fig. 4D). Assuming instantaneous deposition of flood deposit unit 5 (see "The Late Holocene evolution of Ástjörn"), linear interpolation between included 14 C ages LuS 14879 and LuS 15022 ( Fig. 5 and Supplementary Table 1) suggests an age of~1900 cal. yr BP. This places AST3 in a period characterized by a low eruption frequency of the highly active Grímsvötn system 52,57 , narrowing its likely source down to two candidates: the 1698 cal. yr BP G-Svart tephra 57,58 , or the 2436 cal. yr BP Layer 578-579 ash 57 . Comparison with oxide data from both these eruptions reveals that the geochemistry of AST3 is indistinguishable from G-Svart (Fig. 4D).
AST4 contains three glass populations. Most shards (n = 31) display a compositional affinity with either the VAK Askja (AST4-1, n = 17) and Krafla (AST4-2, n = 12) volcanoes ( Fig. 4E and Supplementary Fig. 3). The~3100 cal. yr BP age that we derive for AST4 through linear interpolation between AST3 (G-Svart) and 14 C sample LuS 15022 agrees with known eruptions of both volcanoes that have been dated but lack geochemical fingerprints. Askja, which produces tholeiitic magma although no basaltic tephras have been attributed to this system 59,60 , experienced fissure eruptions between 2900 and 3500 cal. yr BP 61 . Perhaps more significantly, the only known Holocene explosive basaltic eruption of Krafla occurred around 2850 ± 250 cal. yr BP: this so-called Hverfjall event dispersed ash in the direction of Ástjörn 62 . We consider the third small (AST4-3, n = 7) sub-set of AST4 shards, which likely derive from Katla ( Supplementary Fig. 3), as outliers. In light of the above, we cannot confidently link AST4 to one specific eruption, but plot the concurring (and overlapping) ages of the afore-mentioned Krafla and Askja events in our chronology for reference (see Fig. 5).
Timing and magnitude of flood events. By precisely dating the three slackwater deposits of units 7, 5, and 2 in the Ástjörn basin, this study refines the Late Holocene outburst flood chronology of the Jökulsá á Fjöllum catchment. Previous efforts primarily relied on tephra horizons that solely provide minimum or maximum age estimates for floods because of their irregular stratigraphic distribution 63 . Also, a dearth of reliable provenance indicators for some of these ash markers raises the possibility of miscorrelation in an environment where volcaniclastics are omnipresent and easily redistributed by katabatic winds 12,17,32 . Here, we combine robust geochemical tephra fingerprints with 14 C ages that bracket flood deposits to overcome these limitations and capitalize on the strengths of both methods. This approach identifies Late Holocene floods around (1) 3500 ± 500 cal. yr BP -based on the 95% confidence range of our Clam-derived age-depth model at the dated upper contact of unit 7, (2) 1500 ± 100 cal. yr BP-based on extrapolating the linear fit between plant macrofossil-derived 14 C age LuS 15022 and the G-Svart (AST3) tephra to the basal contact of unit 5, and (3) 1350 ± 50 cal. yr BP -based on the calibrated 2σ range of peat macrofossil-derived 14 C age LuS 15021 taken at the base of unit 2 (Supplementary Table 1 and Fig. 5). The use of available basal ages, which are less likely to be affected by reworked flood-eroded organic material 7 , is justified by the absence of erosive contacts: our 63.5 µm resolution CT imagery reveals that unit transitions are sharp but conformable (Fig. 7D). Good agreement between these age constraints and estimates from the top of slackwater deposits further strengthens confidence in the presented flood history. Compared to previous reconstructions e.g. 15,17,63 , our results show that the 1-2 ka BP flood identified by many workers actually comprises 2 closelyspaced events. This discovery helps resolve recent cosmogenic evidence of knickpoint retreat and terrace abandonment after 1.5 ka BP at the up-stream Dettifoss waterfall in greater detail 15 (Fig. 1). Our findings also confirm previous evidence of extensive flooding around 4 ka BP from exposure ages and flood deposits capped by Hekla 4 ash in sediment sections along the lower Jökulsá á Fjöllum 15,17,29,38 . Finally, as GPR and field evidence suggest that master core AST-P2-18 covers the entire sediment infill of the lake (Fig. 2D), the presented 4.5 cal. ka BP basal age provides a minimum age estimate for the last flood that was sufficiently powerful to remove all sediment from Ástjörn 17 .
Available flood simulations with the GeoClaw flow model by 64 allow us to constrain the magnitude of past floods in the catchment. Using a Manning's roughness coefficient of 0.05 following the recommendations of 65 for the Jökulsá á Fjöllum watershed, [20][21][22] show that waxing floods first enter the lake from the North when flow exceeds 20,000 m 3 s −1 (Fig. 2A). We should note that this simulation prescribes a 37 m a.s.l peak stage (present-day lake level) while our SfM-generated DEM shows Ástjörn is separated from the adjacent sandur plain by a 44 m a.s.l levee (Fig. 2C). However, by raising questions about the stability of this unconsolidated landform, the stratigraphy of master core AST-P2-18 justifies the use of such a conservative discharge estimate. Notably, accumulation of overbank deposits (unit 8) and peat (unit 3) prior to flooding suggests a more effective exchange of water between river and lake. Relying on output from the same model setup and cross-sections from our SfM-generated DEM (Fig. 2), we calculate that overtopping of the 108 m a.s.l spill-over at the catchment's southeastern perimeter requires a discharge in excess of 130,000 m 3 s −1 21, 22 . The presence of flow- aligned boulder lags, flood-carved bedrock channels and the 10 m high headwall at the lake's southern edge reveal that such events are characterized by catastrophic high-energy flow regimes (Fig. 3) 17,18 . However, the discussed fine sand-dominated signature of analyzed flood deposits in AST-P2-18 are indicative of low-energy backflooding (Fig. 6) 41 . Also, the CT imagery of Fig. 7 shows no obvious evidence of changes in flow direction and regime like erosive contacts. Following from the above, we argue that all three Late Holocene floods inundated Ástjörn from the adjacent sandur plain to the North-constraining discharge peaks to 20,000-130,000 m 3 s −1 .
To assess the relative magnitude of each flood, we analyze the Grain Size Distributions (GSDs) of slackwater units 7, 5, and 2 (Fig. 7). Because of the coupling between flow speed and sediment competence, paleohydrologists often rely on the abundance of coarse sediments to do so 9 . However, this approach is not suitable for the Jökulsá á Fjöllum watershed as observational data reveal that there is no clear-cut relationship between flood discharge and grain size 30 . Summer discharge peaks are dominated by silt-dominated glacigenic suspended load from the foreland of Vatnajökull, whereas coarser sandsized material is only available in weathered upland areas and mostly mobilized during spring snow melt 30,42 . End-Member Modeling Analysis (EMMA; see "Methods") permits us to unmix the contribution of these different sediment sources and derive a robust predictor of flood magnitude after 66 . As can be seen in Fig. 7B, our analysis identifies 3 significant End Members (EMs) that together explain 97% of the data variance. EM 2 dominates all slackwater flood deposits and has a modeled GSD that is near-identical to samples from these units ( Fig. 7A and C). The mean of these distributions also overlaps with the observed range of silty sediments that dominate modern seasonal glacigenic floods 30 . Based on this evidence, as well as the notion that past flooding mobilized the same sediment sources because of their glaciovolcanic origin 17 , we argue that EM 2 best captures flood intensity. To assess the relative magnitude of each event, we developed a so-called Flood Magnitude Index (FMI) by (1) normalizing EM 2 abundances per unit to account for the previously discussed changes in levee height (flood threshold) after 28 , before (2) calculating definitive integrals (area under the curve) from these z-scores as slackwater sediment thickness also reflects flood discharge and duration 25 . This analysis suggests that the most recent 1.35 cal. ka BP flood was far greater than any other Late Holocene event (FMI = 109), with a twice smaller flood magnitude around 3.5 cal. ka BP (FMI = 45), and about one-tenth of this strength at 1.5 cal. ka BP (FMI = 10).
The upper 130,000 m 3 s −1 limit of our model-derived discharge range is at least three times lower than reported lower-end 400,000 m 3 s −1 peak estimates for Late Holocene Jökulsá á Fjöllum floods 19,63 . We attribute this mismatch to dating uncertainties: i.e., past workers fitted modeled water surfaces to undated wash limits e.g. 17 , while our results provide chronological constraints on contemporaneous flooding (Fig. 5). However, surface exposure dates of flood-carved bedrock surfaces reveal that knickpoints at the nearby Dettifoss waterfall retreated more than 2.5 km during these events (Fig. 1) 15 . Taken together, this evidence underscores the highly erosive nature of comparatively modest floods. In doing so, our work implicitly supports simulations and observations that highlight that erosion rates in bedrock canyons can be controlled by other factors than discharge 4,67-69 . Jointed volcanic flows like those found in the Jökulsá á Fjöllum watershed allow shear and drag from flood waters to topple basalt columns with relative ease: calculations suggest that this happens during exceptional seasonal discharge maxima that exceed 500 m 3 s −1 16,70 . Our work thus strengthens these and other studies that suggest that the magnitude of canyon-carving floods may have been over-estimated 69,71 . This could have implications for the use of geological evidence left by these events to assess their impact. Locally, a downward revision of the magnitude of millennial floods along the Jökulsá á Fjöllum reduces the probability of infrastructure damage or the loss of life. Globally, the lowering of flood-derived runoff volumes may challenge assumptions about the sensitivity of ocean circulation to freshwater fluxes. Finally, if comparatively modest floods can be highly erosive, formation of bedrock canyons requires less running water than thought.

Methods
Geomatics. To assess lake depth and sediment thickness prior to coring, we performed a Ground Penetrating Radar (GPR) survey. For this purpose, we used a Malå ProEx system connected to a 50 MHz shielded antenna dragged behind a zodiac in coring tubes during surveying. Following acquisition, data was processed in version 1.4 of the RadExplorer software to improve the contrast of our GPR profiles using prescribed bandpass filtering, DC removal and time-zero adjustment routines after 33 . To map the altitude of flooding thresholds in the Ástjörn catchment with high accuracy, we created a~10 cm resolution Digital Elevation Model (DEM) using Structure from Motion (SfM) photogrammetry on imagery captured with Unmanned Aerial Vehicles (UAVs). For this purpose, we flew a DJI Phantom 4 Pro drone that covered regular double grids at 100 m altitude and took images with 70% overlap using the Pix4Dcapture software. Collected imagery was subsequently processed following the four steps of the conventional workflow in version 1.4.5 of the Agisoft Photoscan suite. First, all photos (n = 831) were aligned at high accuracy using adaptive model fitting, limiting key-and tie-points to 40,000 and 10,000, respectively. Secondly, we generated a high-quality dense point cloud in mild depth filtering mode. To further improve accuracy, camera alignment was automatically optimized and obvious outliers were manually removed. Next, we used this point cloud to build a 2.5-D height field mesh with 10,000,000 faces and referenced it with handheld GPS-derived (Garmin 62 S) Ground Control Points (GPCs). Finally, we built a DEM projected to WGS 84/UTM zone 28 N and exported it in GeoTIFF format for further analysis using the 3D Analyst toolbox in v. 10.4 of Esri ArcGIS to draft the final figures.
Coring. Sediment cores were recovered from Lake Ástjörn in September 2018. After surveying bathymetry and sediment thickness with Ground Penetrating Radar (GPR; see "geomatics"), we selected a max. 5 m deep flat section in the central part of the basin to avoid disturbance (Fig. 2B). By coring in the center of this area, we stay clear of the potentially erosive plunge pool-like feature fronting the headwall at the lake's southern end (Fig. 3B). At the same time, this site is sufficiently close to the levee that is first over-topped during floods to maximize the contrast between background and slackwater sediments 9 . For this purpose, we relied on two systems: a hammer driven piston corer to extract long sequences and a Uwitec gravity corer to warrant preservation of an intact sediment/water interface (zero year; 2018 CE). Two long cores were extracted: piston cores AST-P1-18 (~300 cm) and AST-P2-18 (450 cm). Owing to its greater length, we focus this study on the latter archive and a gravity core (AST-G2-18; 140 cm) taken from the same location (Fig. 2B). Following fieldwork, cores were split, logged and analyzed at the EARTHLAB facility in Bergen.
Stratigraphy. Following visual inspection, we performed a series of nondestructive scanning techniques on master core AST-P2-18. To determine minerogenic input, we mapped sediment geochemistry with an X-Ray Fluorescence (XRF) core scanner. To measure heavier elements with high sensitivity, we fitted the instrument with a Molybdenum (Mo) tube set to 27 kV and 27 mA. In total, 8717 elemental profiles were acquired at 500 μm resolution. We excluded elements with a Signal-to-Noise ratio (SNR; μ/σ) lower than 2 after 72 . Moreover, all XRF data are expressed as Total Scatter Normalized (TSN) ratios after 73 to account for changes in organic and water content. To visualize sediment structure in 3-D, we employed a ProCon-X-Ray CT-ALPHA Computed Tomography (CT) scanner that was operated at 125 kV and 600 μA with a 267 ms exposure time to generate 63.5 μm resolution 16 bit scans. Reconstructed scans were subsequently processed using the ThermoFisher Amira-Avizo software to generate 704*6001 pixel 2-D (XY; ortho) slices and rendering a 20 mm 3 3-D (voxel) volume. Following logging and scanning, we carried out destructive analyses to determine down-core variations in organic content, density and grain size distribution. First, 0.5 cm 3 cubes of sediment were extracted every 10th centimeter throughout core AST-P2-18, increasing to 2 cm resolution to cover lithological transitions in detail (n = 89). Sediments were then dried overnight at 105°C and heated for 4 h at 550°C to measure Dry Bulk Density (DBD) and Loss On Ignition (LOI; a measure of organic content) after 74 . Next, we assessed the grain size distribution of these samples on a Malvern Mastersizer 3000 linked to a Hydro SV dispersion unit. Each sample was measured in triplicate to monitor analytical precision; we subsequently analyzed averages that were sufficiently reproducible (within ISO 13320 limits).
Dating. We sent 10 samples from core AST-P2-18 for Accelerator Mass Spectrometer (AMS) radiocarbon dating at the Radiocarbon Dating Laboratory in Lund (LuS Supplementary Table 1). Because of a scarcity of dateable material, we had to submit bulk samples in certain sediment sections. Non-bulk (aquatic plants and peat) samples were extracted by wet-sieving using a 250 µm mesh before submitting stems or leaves (macrofossils). To determine the timing of major shifts in depositional regime, we focused these efforts on stratigraphic transitions. To strengthen this chronology, we also extracted and identified four basaltic tephra horizons (AST1-4). This method has significant geochronological potential in our study area owing to its proximity to several highly active volcanoes (Fig. 1) that have produced numerous well-characterized tephra horizons e.g. 52,75,76 . We confined our efforts to discrete visible (centimeter-scale) basaltic horizons in organicrich sediment units (Fig. 6). As a result, we do not report any rhyolitic regional markers in this study. The sharp contrast in color (darker) and composition between basaltic ash and surrounding sediments aids visual identification 77 . In contrast, minerogenic units post challenges owing to the basaltic composition of most Icelandic bedrock. We subsequently complemented visual identification of horizons with XRF and CT data: elevated Titanium (Ti) and Manganese (Mn) counts highlight basaltic ash 78 , while the CT density values of tephra are much higher than those of organic material 79 (Supplementary Fig. 1). Next, we extracted 1 g of wet sediment in 0.5 cm wide slices from four discrete horizons (Supplementary Table 2). Examination of these samples under a light microscope (×200) revealed that the material solely comprised basaltic tephra, eliminating the need for density separation 80 . We did, however, float off any organics by centrifuging our samples in distilled water for 5 min at 2500 rpm. We only retained the >125 μm fraction for further analysis following subsequent sieving. Next, tephra shards were mounted on frosted microscope slides and embedded in epoxy resin. Samples were then ground with 1 mm silicon carbide paper and polished with 1μm diamond suspension to expose shards for geochemical characterization, which was performed using wavelength electron microprobe (WDS-EMP) analysis. For this purpose, a Cameca SX100 was used at the Tephra Analysis Unit of the University of Edinburgh. The instrument was operated at an accelerating voltage of 15 kV, with a variable beam (8.8 μm diameter) current of 2 nA (Na, K, Mg, Al, Si, Ca) or 80 kA (P, Ti, Mn, Ti). A secondary glass standard (BCR2g) was analyzed in between and within runs to monitor analytical precision (Supplementary Data 1). Measurements with total oxide wt. % <94.5 or >102.5 were rejected.
Statistics. To generate a chronology for core AST-P2-18 using our age tie-points (see "dating"), we used version 2.3.2 of the Clam package in R 43 . We also employed End-Member Modeling Analysis (EMMA) to decompose sample Grain Size Distributions (GSDs) after e.g. 66 ; this approach has shown significant promise to derive information about the magnitude of past floods 81 . To this end, we relied on the AnalySize tool by 82 in version 9.3 of MATLAB. We specifically used the nonparametric HALS-NMF algorithm, which has accurately reproduced a broad range of artificial datasets 83 . Additional statistical analyses include the calculation of metric Folk and Ward measures for raw grain size data with the GRADISTAT software by 84 , as well as correlation and regression approaches in version 15 of StataSE.

Data availability
All stratigraphic data from investigated sediment core AST-P2-18 that is presented in our main figures has been shared through the DataverseNO repository. Here, the files can be accessed using the following DOI -https://doi.org/10.18710/F8XGVF. Received: 25 September 2020; Accepted: 24 March 2021;