Elevation Changes of the Fennoscandian Ice Sheet Interior During the Last Deglaciation

The dynamics and paleo‐glaciology of ice sheet interiors during the last deglaciation are poorly constrained, hindering ice sheet model reconstructions. We provide direct evidence of Fennoscandian Ice Sheet (FIS) interior behavior during deglaciation through surface exposure dating. Our results demonstrate early thinning of the FIS, prior to the Younger Dryas (YD, 12.8–11.7 ka). Interior thinning in central Norway was concurrent with retreat along the coastline, exposing ice‐free mountainous tracts, potentially as early as 20–15 ka. The FIS then formed moraines in these ice‐free tracts during the YD. This is contrary to current hypotheses advocating a landscape fully covered by cold, inactive ice during this period. Present empirical and model reconstructions fail to capture rapid interior downwastage, increasing uncertainties in ice sheet volume estimates and sea level contributions.


Introduction
Reconstructions of ice sheet behavior from the last deglaciation (≤21 ka) until the onset of the Holocene at 11.7 ka provide fundamental, long-term insights into ice sheet-climate interactions, internal ice dynamics (Ullman et al., 2015), sea level rise contributions (Ritz et al., 2015), and mechanisms of rapid deglaciation (DeConto & Pollard, 2016). These are vital to constrain and calibrate ice sheet (Lecavalier et al., 2014), glacial isostatic adjustment (Whitehouse, 2018), and earth system models (Vizcaino et al., 2015). By combining onshore and offshore geomorphological evidence with geochronological data, former ice margins and their retreat chronology during deglaciation can be reconstructed Hughes et al., 2016;Lecavalier et al., 2014;Margold et al., 2018;Stroeven et al., 2016). Conversely, geological evidence defining paleo-ice sheet surface elevation is relatively sparse. Although this is being progressively addressed (Ballantyne, 2010;Barth et al., 2019;Bierman et al., 2015;Brook et al., 1996;Koester et al., 2017), reconstructions often rely on indirect indicators of thickness such as glacial isostatic adjustment modeling or far-field relative sea level records (Lambeck et al., 2014;Peltier et al., 2015). As a result, ice thickness and interior dynamics during deglaciation are poorly constrained, limiting estimates of glacial isostatic adjustment and ice sheet sea level contributions in the paleo record and for future projections of sea level contributions.
At its Last Glacial Maximum (LGM) peak, the Fennoscandian Ice Sheet (FIS: 3.9 × 10 6 km 3 )  was similar in volume to the present-day Greenland Ice Sheet (~3.0 × 10 6 km 3 ) (Morlighem et al., 2017) and coalesced with the Barents Sea Ice Sheet (BSIS) and British-Irish Ice Sheet (BIIS) to form the Eurasian Ice Sheet Complex (EISC). The FIS's LGM configuration was complex, with asynchronous, asymmetrical marginal behavior Rinterknecht et al., 2018). By the Younger Dryas (YD: 12.9 ka), and likely much earlier, the FIS disconnected from the BSIS and BIIS and was primarily land based . Terrestrial retreat within the high relief alpine mountains of central and southern Norway (~2,000 m above modern mean sea level [a.m.s.l.]) has allowed FIS retreat patterns to be recorded in well-preserved glacial landforms and sediments, providing an unrivaled and relatively untapped opportunity to reconstruct ice thickness and thinning patterns during deglaciation.
After several decades of FIS reconstructions, ice margin positions, deglacial patterns, and the spatiotemporal distribution of subglacial thermal regime (cold or warm based) have been established ( Figure 1) (Cuzzone et al., 2016;Kleman et al., 2008;Stroeven et al., 2016). However, reconstructions of ice thickness since the LGM remain subject to altitudinal discrepancies of 1,000-1,500 m in central and southern Norway (Linge et al., 2006), precluding any accurate understanding of interior FIS dynamics. Early reconstructions were based on undated ice marginal landforms, including tilted ice marginal meltwater channels (Gjessing, 1965;Mannerfelt, 1945), providing clear evidence of vertical downwastage of the FIS interior. Numerical models (Boulton et al., 2001;Patton et al., 2016Patton et al., , 2017 and empirically derived reconstructions (Hughes et al., 2016;Stroeven et al., 2016) suggest radial recession to independent, cold-based ice domes over southern and northern Norway by 10 ka. This widely accepted ice-dome recession hypothesis has been challenged by geomorphological and paleoecological evidence of ice-free tracts during deglaciation (Dahl & Nesje, 1992;Paus et al., 2011Paus et al., , 2015Sernander, 1905). Moraines in front of high-altitude cirques (1,200-1,600 m a.m.s.l) in central and southern Norway (Sernander, 1905) indicate the existence of past cirque glaciers, detached from the inland ice sheet. Based on paleo-equilibrium line altitude (ELA) estimates, these are thought to have existed during the YD (Dahl et al., 1997). Without a formal chronology, it is unclear if the cirque glaciers and the inland ice sheet were contemporaneous. This concept, based on the aforementioned hypotheses, is referred to here as the Sernander-Dahl Hypothesis. Here we use detailed glacial-geological mapping, surface exposure dating, and comparison to ice sheet model results from the hypothesized cold-based center of the southern dome of the FIS to test this hypothesis and examine how the interior configuration of the ice sheet changed during deglaciation.  (Dahl et al., 1997), and white circle shows Mt. Elgåhogna (Goehring et al., 2008). White dashed line demarcates the boundary of the minimum cold-based FIS extent at the Last Glacial Maximum (LGM), and the black dashed line indicates the boundary of the minimum cold-based FIS extent at the YD .

Study Site
Snøhetta (62.320°N 9.268°E, 2,286 m a.m.s.l.) is the largest mountain in Dovrefjell, a sandstone and mica schist massif (Sigmond et al., 1984) in central Norway (Figure 1). The flanks of Snøhetta support three polythermal cirque glaciers, all of which are fronted by terminal, potentially ice-cored moraines. Beyond these moraines, there is no geomorphological evidence for previous, larger cirque glacier extents, and west east trending glacial striations outboard of the northern cirque glacier moraine could only have been created by FIS ice moving along Stropelsjødalen ( Figure 2). As a result, all other moraines are interpreted to have been deposited by the FIS.
The massif displays abundant evidence of ice sheet glaciation and is heavily dissected by glacial valleys and high-altitude cirques (%3E1,500 m a.m.s.l.). High elevation ridges contain weathered regolith mantles and tors, often found in association with erratic boulders, with occasional exposed, glacially molded bedrock surfaces (supporting information Figure S3). Below this, frost-shattered regolith continues, with frequent subglacially eroded erratics and glacially abraded bedrock. Bedrock-incised lateral meltwater channels formed during FIS downwasting are found at 1,500-1,600 m a.m.s.l. We mapped four valley systems surrounding Snøhetta (see Figures 2 and S1). Each valley is characterized by two distinct FIS terminal moraines, with a series of smaller, tightly spaced (100-500 m) transverse ridges between them ( Figure S1). These are interpreted as ribbed moraines, often formed in subglacial thermal regimes varying between cold and warm based (Dunlop & Clark, 2006;Hättestrand & Kleman, 1999). We targeted two of these FIS moraines, at the heads of Stropelsjødalen (Stropelsjødalen moraine) and Svånådalen (Bandranden moraine). Due to geomorphological similarities, we hypothesize that the ages of these two systems are the same as other nearby valleys, where the same depositional pattern is reproduced.

Surface Exposure Dating
To constrain ice sheet thinning, we sampled boulder, bedrock, and tors in Snøhetta for 10 Be surface exposure dating using hammer and chisel, drill, or rock saw. Boulder samples were taken from stable erratic boulders demonstrating clear evidence of subglacial erosion ( Figure S3). Bedrock samples were taken away from surfaces that may have been previously covered by sediment or seasonal snow cover and from upper horizontal rock surfaces away from edges, joints, and cracks. The location, altitude, site description, surface dip, sample thickness, sample geometry, shielding, and level of surface weathering were recorded in the field.
Sample preparation and 10 Be/ 9 Be ratio measurements were performed at SUERC, as part of routine Be runs. 10 Be concentrations are based on 2.79·10 − 11 10 Be/ 9 Be ratio for NIST SRM4325. Exposure ages were calculated using the iceTEA online interface , based on the CRONUScalc calculation framework , the LSDn scaling model (Lifton et al., 2014), and the global 10 Be production rate . Ages quoted in the text include internal uncertainties, for ages including production rate uncertainty, see Table S1. Use of the Western Norway production rate (Goehring et al., 2012) produces ages 3.3 -4.3% younger and within the reported external uncertainties (9.6%), thus not affecting our interpretations. Long-term erosion rates from this region are unknown, but the hard, crystalline sample lithologies and freshness of bedrock surfaces suggest that corrections for postglacial erosion are unnecessary, in agreement with other studies from Norway (Gump et al., 2017). Due to the high elevation of samples (1,358-1,800 m a.m.s.l.), ages were corrected for elevation change using the iceTEA online interface , using the ICE-6G model, based on glacial isostatic adjustment (Peltier et al., 2015). Seasonal snow cover may affect cosmic ray attenuation, decreasing the reported exposure age compared to actual. No long-term snow cover data are available, but snow records from http://www.senorge.no/ suggest 4 month yr − 1 cover of 80 cm. Past snow cover is unknown, but as a basic calculation using present snow rates, production rates would be reduced by~7%, modifying our ages but not altering our interpretations.

Ice Sheet Model Reconstructions
We compare our geochronological and geomorphological findings to newly extracted data from a 3-D model reconstruction of the last EISC, constrained and validated against geophysical data sets (Auriac et al., 2016;Hubbard, 2006;Patton et al., 2016Patton et al., , 2017. The simulation used here represents a best fit in terms of its ability to reproduce the broad-scale patterns of retreat as interpreted from the geological record (e.g., Stroeven et al., 2016). The sensitivity of the ice sheet growth trajectory toward the LGM in response to a range of internal and external parameters and forcing perturbations is detailed within these previous studies.

Results
Twenty-one cosmogenic surface exposure ages ( 10 Be) are used to reconstruct FIS interior behavior. Twelve samples are from perched erratic boulders, ice-molded bedrock, and tors on a ridge from Skardkollen to Snøhetta ( Figure 2) and cover an altitudinal range of 1,554 -1,800 m a.m.s.l. Ages range from 20.3 ± 0.6 to 11.9 ± 0.4 ka after the removal of one outlier (EXPO6; 7.7 ± 0.2 ka). Nine samples from bedrock (n = 4) and boulders (n = 5), are closely clustered from 12.8 ± 0.4 to 11.2 ± 0.4 ka. The three older ages are from bedrock (EXPO10: 18.1 ± 0.5 ka; EXPO18: 20.3 ± 0.6 ka) and an erratic which sits on a bedrock tor with a 5 ka younger exposure age (EXPO11: 17.8 ± 0.5 ka; tor age EXPO12: 12.6 ± 0.4 ka). Samples EXPO10 and EXPO18 may contain inherited nuclides, and EXPO11 most likely contains inherited nuclides. This is common in areas of unerosive, cold-based ice cover (see Goehring et al., 2008). While we do not eliminate these samples from our study, we acknowledge that they must be treated with some caution.
Boulders were sampled from FIS moraines in Stropelsjødalen and Svånådalen (Figure 2), respectively. Four boulders from Stropelsjødalen moraine (SSD) yield an error-weighted mean age of 11.8 ± 0.2 ka. Three boulders from Bandranden moraine (BR) yield an error-weighted mean age of 12.0 ± 0.2 ka, after removal of one outlier (BR4; 9.2 ± 1.1 ka). 10 Be ages from the SSD and BR moraines overlap within 1-sigma internal uncertainty and indicate deposition of FIS moraines at the end of the YD.
Our 10 Be ages provide direct evidence for emergence of high elevation land (%3E1,570 m a.m.s.l.) from ice sheet cover prior to the YD, potentially as early as 20 -15 ka. FIS margins in this region were glaciologically active, producing moraines during the YD. These ages provide a minimum age control for the onset of ice thinning, as exposure of the region's highest peaks are not captured in our ages.
Each of the four mapped valleys (Figures 2 and S1) contains an up valley terminal moraine (including SSD and BR), ribbed moraine deposited after the YD (constrained by ages on SSD and BR moraines), and a second terminal moraine, marking the last known position of the inland ice sheet in this region (see Figures 2 and S3). Similar innermost moraines are found in Dovrefjell (Sollid & Sørbel, 1994) and Norwegian fjords and coastal valleys (Sollid & Torp, 1984;Stroeven et al., 2016) and are associated with the Preboreal Oscillation (PBO) (~11.47-11.35 ka) (Rasmussen et al., 2014). Given their morphostratigraphic position and similarity to these moraines, we hypothesize that the innermost moraines from this region also date to the PBO ( Figure S2, Kaldvell, Svånålegeret, Kjelsung, Grøna, and Mjogsjobekken). Using moraine height and location, and our chronology (external uncertainties included), we calculate hypothesized ranges of mean annual thinning and retreat rates between the YD and PBO of 0.23-1.68 m yr − 1 and 8.6-43.4 m yr − 1 , respectively.

A New Late Glacial Ice Sheet Configuration in Southern Norway
Together, our data provide the first directly dated evidence for ice-free tracts within the central FIS during deglaciation and a thin YD FIS in the Snøhetta region. Thinning during this period is also recorded at Elgåhogna,~150 km east of Snøhetta (Goehring et al., 2008) (Figures 1 and 3). For consistent comparison with our ages, we have recalculated their 10 Be ages using our input parameters. Their chronology supports our findings, providing further evidence for FIS thinning as early as 20-15 ka and rapid thinning from 13 to 10 ka. The similarity in timing of sample exposure at Elgåhogna and Snøhetta suggests that interior FIS downwastage took place over a %3E150 km 2 region (Figure 3).  Goehring et al. (2008). Ages are grouped as %3E 13.5 ka (n = 11; black outlined circles) and %3C13.5 ka (n = 26; blue outlined circles). Our new hypothesized YD ice surface and the modeled (Patton et al., 2017) YD ice surface are shown in blue and red dashed lines, respectively. Location of the A-A′ transect is shown in Figure 1.
To allow this new hypothesized ice configuration during deglaciation, the ice surface would be tilted eastward between Snøhetta and Elgåhogna. At the YD this ice surface would lie below the mean altitude of cirque floors in Rondane and Sølnkletten (Figure 3), allowing the contemporaneous existence of FIS-independent cirque glaciers and confirming the Sernander-Dahl Hypothesis (Dahl et al., 1997;Sernander, 1905). When viewed alongside evidence of YD cirque glacier activity and YD advance of marine terminating glaciers in western Norway (Bakke et al., 2009;Briner et al., 2014;Lohne et al., 2012), these ages provide evidence for synchronicity between the FIS interior and margin, suggesting an external forcing acting across the entire southern central FIS.
Our new proposed FIS configuration shows startling similarity to the geometry of the Cordilleran Ice Sheet (CIS), where extensive thinning and mass loss during the Bølling-Allerød (B-A) led to the emergence of mountainous regions from the ice sheet and permitted CIS and independent glacier regrowth (Menounos et al., 2017). This suggests that multiple ice sheets experienced systematic interior drawdown at this time, in addition to retreat.

Modeled Dynamics of FIS Thinning
Our reconstruction is compared to thinning data extracted from published numerical ice sheet model reconstructions (Patton et al., , 2017 to further test the plausibility of the Sernander-Dahl Hypothesis and its implications for deglacial patterns beyond our study region. Model results project the LGM ice surface at a maximum elevation of 2,050 m a.m.s.l. and place the southern FIS's central spreading axis 30 km east of Snøhetta, with ice persisting into the Early Holocene, in line with contemporary deglaciation reconstructions . Full glacial conditions stabilized in the region by 32 ka, with a maximum-modeled regional ice surface of 1,850 m a. m.s.l., meaning that Snøhetta and other major peaks within the Dovrefjell range persisted as nunataks throughout the LGM. Initial interior thinning began at~25 ka (Figure 4), despite FIS mass gain until 23 ka. This is driven by late expansion of the FIS along its eastern sectors until 17 ka (E. Larsen et al., 1999;Rinterknecht et al., 2018), facilitating large scale but slow drawdown from the central ice dome. Thinning was then enhanced by breakup of the Norwegian Channel Ice Stream (NCIS) at~20 ka (Svendsen et al., 2015) (see Figure 4a). Large volumes of meltwater discharge of Scandinavian origin are also reported in the Channel River during this period, from~23 -16.5 ka (Toucanne et al., 2015). Despite mass loss through this period, net contribution of the FIS to eustatic sea level rise from 25 to 17.5 ka is only 0.39 m sea level equivalent (SLE).
Initial FIS thinning is synchronous with reconstructions from other ice sheets (Figure 4a). North American Ice Sheet (NAIS) volumes decreased from~25 ka (Patton et al., 2017;Peltier, 2004;Peltier et al., 2015). Data and modeling suggest that BIIS thinning began at 21 ka and accelerated at 19 ka (C. D. Clark et al., 2012). This synchronicity of thinning across multiple ice sheets is likely to have been forced by increasing highnorthern-latitude insolation (P. U. Clark et al., 2009) and climatic warming (Masson-Delmotte et al., 2005), amplified by albedo, changing atmospheric circulation patterns and other positive feedback mechanisms (Bailey et al., 2018).
A second, more pronounced, phase of FIS thinning occurs in central and southern Norway from 15.3 ka, leading to large volumetric losses during the B-A. Modeling results demonstrate widespread surface melt, ice sheet thinning, and exposure of upland terrain (Figure 4b), with a drastic reduction in FIS volume and  Patton et al., 2017;Peltier, 2004;Peltier et al., 2015). NA ice sheet is scaled to a separate y-axis. (b) Ice volume changes from the FIS model Patton et al., 2017) from 30 to 5 ka (black line) and the FIS thickness from the Snøhetta grid cell.  Deschamps et al., 2012), the FIS contributed~2.7 m SLE. This phase of interior downwasting is contemporaneous with ice margin retreat in southwest Sweden (~16 -15.2 ka: Cuzzone et al., 2016) and NCIS retreat (~16 ka: N. K. Larsen et al., 2012), highlighting the potential for FIS-wide drawdown and thinning driven by ice margin retreat. Scandinavian proxy records from this period reconstruct July temperatures of 12-14°C proximal to the FIS margin, indicating a strong thermal summer forcing of rapid ice sheet melt (Schenk et al., 2020).
Rapid thinning during the B-A is also reported for the NAIS, with ice loss at~14.5 ka (Barth et al., 2019;Corbett et al., 2019;Koester et al., 2017;Menounos et al., 2017), and is cited as responsible for North America saddle collapse at 15.7 -14.0 ka (Gregoire et al., 2012) (Figure 4). Greenland and Antarctic ice sheets also responded to B-A warming, with major Greenlandic ice stream retreat from the shelf edge at~15 cal. ka BP (Dowdeswell et al., 2014;Cofaigh et al., 2013) and Antarctic thinning at~15 ka (Bentley et al., 2010;P. U. Clark et al., 2009). Following widespread thinning, YD cooling led to renewed ice sheet thickening (Figure 4b), after which final disintegration of the FIS followed swiftly.
The modeled reconstruction of deglaciation supports our surface exposure ages, demonstrating the early emergence of high-altitude topography in southern central Scandinavia, close to the central ice divide as early as 20 ka. However, model resolution (10 km) is too coarse to resolve the topographic detail across the study region, resulting in oversimplified pattern of margin retreat and interior thinning. Topography has been shown to have a large impact on FIS models (Åkesson et al., 2018), and a finer grid resolution would better resolve arterial pathways within the modeled bed topography, increasing the sensitivity of the modeled ice sheet interior to external drivers such as the B-A. We hypothesize that this would lead to more rapid modeled thinning and an improved correspondence between the land exposure patterns with our findings.

Implications for FIS Thinning
Our calculated mean annual thinning rate range of 0.23-1.68 m yr − 1 from the YD to PBO overlaps with our EISC model, which indicates 0.16-0.5 m yr − 1 of thinning during this period (Figure 4b). Based on this overlap, the modeled thinning rates support our hypothesis that the inner moraines are PBO in age. As no other thinning rates exist for the interior of the FIS, it is possible that our YD-PBO thinning rates were in part driven by the effect of regional mountains on ice sheet meteorology, similar to "ablation islands" in Antarctica (Bintanja, 1999).
Our estimated YD-PBO thinning rate is similar to paleo-thinning rates from other ice sheets. Thinning of ≥0.98 m yr − 1 and up to 4 m yr − 1 were reported during the B-A in the southeast NAIS (Barth et al., 2019) and Cordilleran Ice Sheet (Menounos et al., 2017), respectively. Our thinning rate also falls within the range of present-day thinning rates in West Antarctica and Greenland (Smith et al., 2020), and from Pine Island Glacier (1.12 and 1.67 m yr − 1 ) during the Early Holocene (Johnson et al., 2014), and elsewhere in Antarctica .
In combination with our geomorphological observations, our estimated thinning rates suggest that this region of the FIS was not cold based during deglaciation. It is probable that thinning of this magnitude results from a transition to warm-based conditions, or other interior ice-dynamic changes. However, the dynamic processes that triggered downwastage and thinning across the FIS interior during deglaciation remain understudied.
Investigating the dynamics and basal footprint of paleo-ice sheets offers an alternative route to inform on the efficacy and magnitude of such processes. Given that ice sheet volume reconstructions rely primarily on constraints of areal extent (Hughes et al., 2016;Ullman et al., 2015), these new data on interior FIS thickness have important implications for estimates of ice sheet volume and contribution to both global sea level change and glacio-isostatic loading. These data provide, for the first time, a connection between interior and marginal FIS ice sheet dynamics, which is significant for an improved understanding of the magnitude and rates of drawdown processes of both paleo-ice sheet and contemporary ice sheet. The analysis presented here provides further impetus to examine sensitivity of the FIS to rapid climate change during deglaciation periods.

Conclusion
Our study demonstrates that interior thinning of the central and southern FIS during the last deglaciation occurred earlier than previously thought, in accordance with the Sernander-Dahl Hypothesis. Areas previously considered to be ice covered and cold based until the early Holocene were ice-free prior to the YD and potentially as early as 20 -15 ka. Both thinning and margin retreat appear to have been synchronous at both the B-A and following the YD. With previously published ages and geomorphological evidence from cirque glaciers, this supports the premise of ice-free high elevation massifs in the center of the FIS, while its outer margins persisted on the western coast of Norway.

Data Availability Statement
All data from this study are available in the supporting information and have been uploaded to the Open Science Framework (https://osf.io/jgb9w/).