Episodic Neoglacial expansion and rapid 20th century retreat of a small ice cap on Bafﬁn Island, Arctic Canada, and modeled temperature change

. Records of Neoglacial glacier activity in the Arctic constructed from moraines are often incomplete due to a preservation bias toward the most extensive advance, often the Little Ice Age. Recent warming in the Arctic has caused extensive retreat of glaciers over the past several decades, exposing preserved landscapes complete with in situ tundra plants previously entombed by ice. The radiocarbon ages of these plants deﬁne the timing of snowline depression and glacier advance across the site, in response to local summer cooling. Erosion rapidly removes most dead plants that have been recently exposed by ice retreat, but where erosive processes are unusually weak, dead plants may remain preserved on the landscape for decades. In such settings, a transect of plant radiocarbon ages can be used to construct a near-continuous chronology of past ice margin advance. Here we present radiocarbon dates from the ﬁrst such transect on Bafﬁn Island, which directly dates the advance of a small ice cap over the past two millennia. The nature of ice expansion between 20 BCE and ∼ 1000 CE is still uncertain, but episodic advances at ∼ 1000 CE, ∼ 1200, and ∼ 1500 led to the maximum Neoglacial dimensions ∼ 1900 CE. We employ a two-dimensional numerical glacier model calibrated using the plant radiocarbon ages ice margin chronology to assess the sensitivity of the ice cap to temperature change. Model experiments show that at least ∼ 0.44 ◦ C of cooling over the past 2 kyr is required for the ice cap to reach its 1900 CE margin, and that the period from ∼ 1000 to 1900 CE must have been at least 0.25 ◦ C cooler than the previous millennium, results that agree with regional temperature reconstructions and climate model simulations. However, signiﬁcant warming since 1900 CE is required to explain retreat to its present position, and, at the same rate of warming, the ice cap will disappear before 2100 CE.


Introduction
Although summer insolation in the Northern Hemisphere has declined steadily through the Holocene, favoring cryosphere expansion, glaciers worldwide are currently losing mass (Stocker et al., 2013). Globally, summer temperatures have been increasing over the last 60 years, a trend that is more pronounced at high latitudes due to strong positive feedbacks (Serreze and Barry, 2011). Since summer temperatures are the dominant control on glacier mass balance in the Canadian Arctic (Koerner, 2005), the reversal of late Holocene cooling (Kaufman et al., 2009) has caused recent retreat of ice caps and glaciers in the region. This continued shrinkage of the cryosphere reinforces the need for records of past glacier and climate change to provide context for contemporary warming. most evidence of previous advances. Other archives, such as lake sediment records, provide continuous records of climate, but typically require glacial activity to be inferred indirectly. "Threshold" lakes record a glacier entering and exiting a well-defined catchment but provide the timing of a glacier margin at a single point (Briner et al., 2010). Lichenometric surface exposure dating, often used on glacial landforms, provides relative age information, but conversion to absolute exposure age has large uncertainties (e.g., Osborn et al., 2015;Rosenwinkel et al., 2015). Furthermore, even though significant advances have been made improving analytical measurements for cosmogenic radionuclide (CRN) exposure ages of moraine boulders, the technique can be compromised by issues of nuclide inheritance and post-depositional moraine destabilization (Crump et al., 2017;Pendleton et al., 2017). Consequently, moraines can be difficult to date precisely, and moraine records by nature are discontinuous.
Recent ice-margin retreat of cold-based ice caps on Baffin Island is exposing preserved landscapes. The radiocarbon ages of in situ plants on these preserved lands record the time when snowline lowering and/or ice advance covered the site, killing and entombing the plants beneath the glacier. Commonly, the re-exposure of dead plants leaves them highly susceptible to erosion, most often by meltwater flowing along the margins of cold-based ice, but also by blowing snow in winter. The ephemeral nature of these re-exposed dead plants means that their ages are representative of the most recent snowline advance at that location; multiple exposures and burials are unlikely (Miller et al., 2013).
However, in certain topographic settings, where plant removal by meltwater has been inefficient, dead plants may remain for decades. A transect of dead plant radiocarbon ages perpendicular to a receding ice margin can provide a more reliable and more continuous record of ice advance than is possible from moraine records . The occurrence of long-preserved ice-killed tundra plants under optimal circumstances hundreds of meters beyond retreating ice margins provides an opportunity to derive near-continuous records of ice-margin advance through the Neoglacial period. Here we present a chronology of ice-margin advance from Divide Ice Cap (informal name), a small mountain ice cap in southeastern Baffin Island, derived from radiocarbondated dead plants, and utilize numerical modeling to estimate the changes in summer temperature required to reproduce the observed record of ice margin advance.

Neoglaciation in the eastern Canadian Arctic
There is substantial evidence of repeated glacier advances on Baffin Island during Neoglaciation (∼ 5 ka through the Little Ice Age; Miller et al., 2005) as summer temperatures declined, but absolute chronologies remain sparse and imprecise. Foundational mapping and initial chronologies derived from lichenometry from nested moraines show that some glaciers were approaching their Neoglacial maximum dimensions as early as 3.5 ka (Davis, 1985;Miller, 1973). Radiocarbon ages of tundra plants killed by snowline depression, and only now re-emerging from beneath cold-based ice, show that glaciers began to expand as early as ∼ 5 ka, followed by episodic intensifications culminating in the Little Ice age (Anderson et al., 2008;Miller et al., 2012Miller et al., , 2013Margreth et al., 2014). Although most studies conclude that glaciers reached their maximum Neoglacial dimensions during the LIA, Young et al. (2015) produced a CRN chronology of nested moraines in northeastern Baffin Island that suggests that Naqsaq glacier advanced to a position similar to that of its LIA extent at ∼1050 CE. However, moraine chronologies hinting at the nature of pre-LIA temperature change are sparse, and the chronology from an individual site can be influenced by non-climatic factors; thus, other data sets that provide direct evidence of Neoglacial ice advance are needed.

Field setting: Divide Ice Cap
Divide Ice Cap (DIC) is situated ∼ 60 km southwest of the settlement of Qikiqtarjuaq on eastern Baffin Island ( Fig. 1; Glacier Atlas of Canada 46204L11/12, 46204M247/248). The ice cap currently mantles the northern slope of a local summit between ∼ 1550 and 1180 m a.s.l., spilling down into a local saddle, where it separates into two outlet glaciers flowing down either side of the saddle, orthogonal to the main ice cap flow. A vegetation trimline high on the southfacing slope on the opposite side of the saddle defines the maximum LIA limit (Fig. 1). The uniform nature and clarity of this trimline as well as the southerly aspect of the slope all suggest that the trimline is indeed the result of an ice margin and not persistent snow. Aerial imagery from 1960 shows the ice margin ∼ 15 m below this trimline (Fig. 1;National Air Photo Library, 1960).
The 2015 ice margin lies 60 m below the LIA trimline. Our sample transect runs from the modern ice margin upslope toward the LIA trimline and is aligned with the topographic drainage divide where ice growth upslope is driven by internal deformation (Fig. 2). Between the ice margin and the trimline, the land surface is carpeted by till boulders amidst finer-grained morainal debris. The presence of preserved plants still in growth position exposed by the retreating ice margin is evidence for a high degree of subglacial landscape preservation, indicating a cold-based, non-erosive basal regime for DIC at our transect.
Until now, the majority of studies utilizing dates on previously ice-entombed plants have relied on plants collected close to the ice margin at multiple ice caps to build a composite data set that highlights periods of regional cooling (e.g., Miller et al., 2013;Margreth et al., 2014). However, because meltwater erosion at the location of our transect was inefficient, in situ dead plants remain preserved from the ice margin to the trimline, although with decreasing abundance at higher elevations as the trimline is approached (Fig. 1). Miller et al. (2017) recently employed a similar method on Svalbard by collecting a transect of six dead, in situ mosses up to 250 m away from the 2013 CE ice margin. The radiocarbon ages of these plants provide evidence of an episodically advancing ice margin over the last ∼ 1200 years.

Vegetation collection
In August 2015, when seasonal melt was at its maximum, 11 in situ dead mosses were collected along a ∼ 200 m transect between the 2015 CE ice margin (1185 m a.s.l.) and the LIA trimline (1238 m a.s.l.; Fig. 1). Within ∼ 5 m of the ice margin, all exposed vegetation was dead. Regrowth and recolonization of plants began ∼ 30 m from the ice margin, with dead plants increasingly rare and living plants increasingly abundant up the transect. A single ice margin plant was also collected on the western margin of the ice cap in August of 2014. During field collection, plants are closely inspected for spontaneous regrowth, which is discernable in the field based on plant color (La Farge et al., 2013). Woody plants were avoided because of their higher survival potential and because their stems have an average radiocarbon age that is older than the actual kill date. In the previous year, we collected in situ dead vegetation at the western margin of the ice cap, closer to the summit (sample #12; Fig. 1). Between 1 and 3 filaments from each sampled plant were washed with deionized water and freeze dried before being graphitized at the Laboratory for AMS Radiocarbon Preparation and Research (NSRL). Graphite targets were measured at the W. M. Keck Carbon Cycle Accelerator Mass Spectrometry Laboratory at the University of California, Irvine. Final radiocarbon dates were calibrated using OxCal 4.2.4 (Bronk- Reimer et al., 2013) and reported here as the weighted mean, 1σ , and 2σ uncertainties.

Glacier model
We employ a two-dimensional finite difference numerical model in conjunction with the in situ vegetation to further investigate the evolution of the Divide Ice Cap over the last ∼ 2000 years. The model, modified from Kessler et al. (2006), calculates an ice thickness on a two-dimensional terrain model governed by explicit equations for ice flux and mass balance using shallow-ice physics (e.g., Blatter et al., 2011), an approach similar to Åkesson et al. (2017).

Terrain production
Prior to model implementation, a terrain model of the bedrock surface is required. This involves removing (to the best approximation) the modern DIC from an existing digital elevation product. ASTER digital elevation data from 2011 CE was used as a base for the two-dimensional terrain model, resampled to a 60 m pixel size. The ASTER data product was retrieved from the online Data Pool, courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, https: //lpdaac.usgs.gov/data_access/data_pool. The resulting surface was smoothed using a 7 × 7 mean filter to remove artifacts in the raw data that would lead to instabilities in the model. Modern ice, including DIC and the ice on the surrounding summits, had to be removed to create an ice-free terrain to model upon. Using the best approximation of basal shear stress (τ b ) to be ∼ 100 kPa (Haeberli, 2016), current ice thicknesses (H ) were calculated following Cuffey and Paterson (2010): where ρ is the density of ice (917 kg m −3 ), g is gravitational acceleration (9.81 m s −2 ), and θ is the surface slope of the modern ice surface. Calculated thicknesses were then subtracted from the modern terrain surface to produce an icefree surface. Eventual model runs show that calculated and modeled modern ice thickness are the same within 10 % of each other. Kessler et al. (2006) drove ice formation with a climate dictated by an equilibrium line altitude (ELA), a mass balance gradient with elevation, and a maximum positive balance (maximum accumulation). However, the overall low accumulation rates of our high-latitude site and the variable aspect of the DIC, which increases the influence of solar radiation during the melt season, necessitates a different approach (Benn and Evans, 2010). Ice core records, observational studies in the eastern Canadian Arctic, and previous modeling work suggest a maximum accumulation of 0.3 m w.e. yr −1 (meters water equivalent per year) throughout the Holocene (Hooke, 1976;Serreze et al., 1995;Anklin et al., 1998;Mair et al., 2005), which, given the limited elevation range of the DIC location, is applied as the annual accumulation across the whole model surface.

Mass balance
Although precipitation records are sparse near the study site, ice core records from the summit of the Greenland Ice Cap show that, regionally, precipitation varied by only ∼ 6 % over the last 1200 years (Alley, 2004). Given the relatively low accumulation rates, summer temperature is likely the dominant driver of glacier advance and retreat at DIC (Koerner, 2005).
Temperature index melt models generally capture the majority of summer melt due to the strong relationship between air temperature and components of the surface energy balance (Braithwaite and Olesen, 1989;Lang and Braun, 1990), and have been successfully applied on larger scales in the eastern Canadian Arctic ). However, the relatively low accumulation rates, small elevation range, and variable aspect and slope at DIC along with the asymmetric trimlines surrounding summits (Fig. 1) suggest that the glacier system is sensitive to small temperature changes and that incoming solar radiation is an important factor in the overall glacier mass balance. Additionally, previous studies of the orientation of cirque glaciers on Baffin Island (Williams, 1975) and summertime snow patch distribution in the eastern Canadian Arctic (Lauriol et al., 1986) emphasize the role of solar energy in enhancing melt on southerly aspects. Wind redistribution of snow could be another factor in the more positive mass balance of northerly aspect of the Divide Ice Cap. However, available meteorological wind direction data from Clyde River to the north and Iqaluit to the southwest both show the prevailing winds to be north-north westerly (Gearheard et al., 2010;Nawri and Stewart, 2010). While these observations do not rule out the possibility that local winds at DIC differ, and move snow from south-facing to north-facing slopes, they strongly suggest that the dominant wind patterns cannot explain the asymmetric mass balance. In areas of overall low accumulation and relatively low relief where the snowline is only just below mountaintops, solar radiation modulated by slope and aspect can exert a strong control on the annual pattern of accumulation (Benn and Evans, 2010).
Given the above concerns regarding melt driven by temperature and solar radiation, we calculate the summer melt rate (M, mm day −1 ) using a radiation modified positive degree day melt model for air temperatures (T ) above 0 • C (Hock, 2005;Jonsell et al., 2012;Kane et al., 1997): The melt contribution from air temperature is calculated using the product of a degree day melt factor (m T , 6.3 mm day −1 • C −1 ; Braithwaite, 1981) with positive degree days (PDD) over the terrain surface. Using a prescribed sea-level mean annual temperature (MAT), a MAT at elevation is calculated using a near surface lapse rate of −4.9 • C km −1 (derived from Canadian Arctic glaciers; Gardner et al., 2009). Annual temperature cycles are then calculated at all elevations using an amplitude of 20 • C (based upon daily temperature records from Dewar Lakes meteorological data from 1959 to 2015) around the MAT. Integration of the portion of the curve where T > 0 • C provides the PDD for each location on the terrain surface. Calculation of the melt contribution from radiation employs the product of a radiation melt factor (m R , 0.036 mm day −1 (W m −2 ) −1 ; obtained via model calibration; see Sect. 5.3.1. below) with that portion of the incident solar radiation (R; W m −2 ) that is not reflected from the surface (albedo = a = 0.5; Benn and Evans, 2010). Solar radiation (R) for the DIC latitude and elevation is calculated following Kustas et al. (1994) and Kumar et al. (1997). In this case, radiation is also modulated based on slope and aspect (Cuffey and Paterson, 2010;Hock, 2005;Kumar et al., 1997). Melt from solar radiation only takes place when T > 0 • C, which is calculated from the annual temperature cycle from the PDD component. Net mass balance is then calculated as the winter mass balance (accumulation) minus the summer mass balance (melt).

Modeling ice surfaces
Driven by the above mass balance, the model, following Kessler et al. (2006), calculates an ice surface on the supplied 2-D terrain surface using explicit equations for ice flux and mass conservation (Eqs. 1 and 8, respectively, in Kessler et al., 2006). Using a shallow ice approximation and the recommended coefficient for Glen's flow law for polar ice of 3.5 × 10 −25 Pa −1 s −1 (Cuffey and Paterson, 2010), ice discharge is driven by the shear stress associated with ice thickness and surface slope. Field observations of the highlypreserved land surface, including vegetation still in growth position, justify a no-slip basal boundary condition. Sliding was therefore disallowed, and ice moves only by internal deformation.

Plant ages
The single in situ sample collected on the western margin of DIC produced a calibrated age of 26 CE + 29/−21 (Fig. 1, Table 1). Calibrated radiocarbon ages of plant material from the transect become younger away from DIC, ranging from 942 CE + 41/−36 and 1029 + 33/−3 at the 2015 ice margin to 1780 CE + 111/−165 near the trimline. The ages of four samples closest to the 2015 ice margin, although collected up to ∼ 20 m away from the margin, overlap at ±2σ (Table 1). A single plant collected near the LIA trimline returned a postbomb age; this likely reflects random death of a plant that recolonized after deglaciation (Fig. 1, Table 1).

Neoglacial Divide Ice Cap expansion and 20th century recession
The kill date at the western margin of DIC indicates ice expansion across that location ∼ 2 ka and continuous ice cover until 2014 CE. The ages near the 2015 CE ice margin (henceforth the 1000 CE margin) define DIC expansion at ∼1000 CE. Since the ages define the time when ice advanced through that location, the distribution of ages requires that the rate of ice margin expansion increases around ∼1000, 1200,   (Hua et al., 2013). n/a stands for "not applicable".
www.clim-past.net/13/1527/2017/ and 1500 CE (Fig. 3). These periods of accelerated ice expansion with periods of enhanced regional cooling and ice expansion shown in the compilation of dated ice-entombed mosses from Baffin Island from Miller et al. (2013) and Margreth et al. (2014;Fig. 3). This correlation between the transect ages and periods of regional expansion shows that DIC responded similarly to regional cooling as the rest of the cryosphere on Baffin Island. Given that each transect plant age represents the time when the ice margin advanced through that position (not necessarily a standstill such as a moraine deposit would represent), the distribution of plant ages, at minimum, requires a change in the rate of ice margin expansion between samples sites (dashed line; Fig. 3). The resolution of our transect does not allow for the definitive detection of any periods of ice margin standstill, but the change in rate of ice margin advance between transect ages indicates that DIC episodically advanced over the past ∼ 1000 years. Possible explanations for the change in rate of ice margin advance include continued expansion, ice margin stasis, or ice recession. However, the spatial and temporal distribution of the transect data and their correspondence with regional ice-expansion records supports our interpretation that DIC expanded episodically between ∼1000 and 1900 CE. The proximity of the highest elevation radiocarbon sample to the trimline suggests that DIC approached its Neoglacial maximum position sometime between ∼1600 CE and the end of the 19th century. Aerial photographs from 1960 show that the ice margin had already retreated ∼ 15 m vertically from its maximum extent (trimline), indicating abandonment of its Neoglacial maximum position prior to 1960 (Fig. 1).

Ice cap modeling
Using the above spatial and temporal constraints on ice margin movement over the last ∼ 2000 years, model simulations can be used to explore the climate sensitivity of the Divide Ice Cap system, and to estimate the temperature changes required to reproduce the observed advance and retreat cycle.

Model calibration
Prior to running full simulations of the last ∼ 2000 years, the model must be calibrated for an appropriate solar radiation melt factor (m R ). Traditionally this value is calibrated using in situ solar flux and melt rate data. However, since that information is lacking here, a different approach must be taken.
Given the transect chronology observed, the simplest (and most likely) history involves continuous ice advance (though likely at varying rates) from ∼ 26 BCE to 1900 CE followed by modern retreat from ∼ 1900 to present. Given this scenario, the model was calibrated to the observed transect chronology and run to maximum Holocene extent conditions (or at ∼ 1900, culmination of the LIA) using a range of accepted m R (Jonsell et al., 2012) values. In other words, we iteratively tested m R values to most closely match the max- imum LIA dimensions. With this approach, we find that an m R value of 0.036 mm day −1 (W m −2 ) −1 provides the best reproduction of Holocene maximum ice extent (Fig. 4). Values of m R above or below 0.036 mm day −1 (W m −2 ) −1 produced far too much or too little ice, respectively (Fig. S1 in the Supplement), making this value our closest approximation.

Full simulation
Using the calibrated model from above, it is possible to reproduce the full history of advance and retreat at Divide Ice Cap over the last ∼ 2000 years. Assuming a stable ice cap with a margin just inside sample #12, the model requires an average minimum summer cooling of 0.19 • C persisting for 1000 years to grow ice across the sample #12 location and reach the 1000 CE margin in the allotted ∼ 1000 years. A subsequent additional cooling of ∼ 0.25 • C is required for ice to advance across the 1000 CE margin to the LIA trimline by 1900 CE. It is of course likely that some decades were much colder than others; our modeling suggests that a minimum average cumulative summer cooling of ∼ 0.44 • C over the ∼ 1900 years is required to advance the observed DIC ice margin through the constraints provided by the kill dates. Building upon the ice expansion simulations, the model was used to investigate recent warming and ice retreat. DIC likely began receding around ∼1900 CE, retreating from its maximum extent to its modern position over the past ∼ 100 years. Under a simple linear warming scenario of 0.028 • C yr −1 , the modeled DIC retreats from the LIA trimline to its 2015 position in 105 years (Fig. 5). This rate of warming is slightly higher than the 0.014 • C yr −1 (since 1959 CE) documented from the interior Dewar Lakes weather station ∼ 320 km to the northwest, but significantly lower than the recent warming rate from Qikiqtarjuaq (0.87 • C yr −1 1995-2009 CE). The modeled rate of warming amounts to ∼ 2.7 • C of cumulative warming over the last century, which is similar to values recorded elsewhere in the eastern Canadian Arctic (Bekryaev et al., 2010;Stocker et al., 2013). Compared to the simulations of glacier advance up through the LIA, it would appear that the last century of elevated warming has reversed the last ∼ 1000 years of ad-vance from slow cooling. Continuing under the same rate of warming, DIC will disappear by ∼ 2100, and sooner if local warming accelerates.
The temperature changes reported here are minimum values due to assumptions made in the model. The best available data was used for both model initial conditions and mass balance forcings, but limitations in these data introduce uncertainty to the model (see Supplement). However, the general agreement of ice cap simulations to observational data in the region support the first-order model results presented here.

Regional comparisons
The timing of episodes of ice expansion found here agree well with those found by Miller et al. (2017) at ∼1100 and 1500 CE, using a similar transect method in Svalbard. The episodic nature of ice advance in Svalbard and in this study, in the face of steadily decreasing Northern Hemisphere insolation, indicates additional mechanisms and feedbacks are likely significantly influencing local climate. Also, similar to the Svalbard study, plant radiocarbon ages and aerial imagery suggests that DIC reached its maximum Neoglacial extent during the LIA and that warming since the early 1900s has reduced glacier dimensions to a smaller size than any time since 1000 CE. Across the North Atlantic, episodic cryosphere expansion despite steadily declining insolation is not unique. Balascio et al. (2015) found that glaciers in southeast Greenland periodically advanced throughout the late Holocene at similar times as DIC. The authors pointed to large ice rafting events and changes in the Atlantic Meridional Overturning Circulation as the likely driver of regional coolings. Glacier records from the North Atlantic also show episodic advances over the past ∼ 1000 years, suggesting strong climate connections across the region (Solomina et al., 2015). Additionally, the magnitude of temperature change from the glacier model agrees well with a reconstruction of Arctic temperatures over the past 2000 years from Kaufman et al. (2009;Fig. 6).
The magnitude of Neoglacial glacier advances and derived climate interpretations presented here differ somewhat from those based on the Naqsaq valley moraine suite farther north on Baffin Island. Based on sets of tightly clustered 10 Be ages on moraines boulders in a nested moraine complex,  suggests that the Naqsaq glacier reached a position similar to that of its LIA extent at ∼1050 CE. This led the authors to suggest that the climate was similar, with little or no additional cooling from ∼1050 CE through the LIA. This interpretation differs from the observed significant expansion of DIC from ∼1000 CE to its LIA maximum and from the model simulations requiring that the same time period must have been, on average, ∼ 0.25 • C cooler than the preceding millennium.
This discrepancy between the two records may be explained by non-climatic factors affecting the deposition and preservation of moraines. Glaciers with a majority of their area at high elevations relative to the glacier tongue may not respond uniformly to equilibrium line altitude (ELA) fluctuations (e.g., Pedersen and Egholm, 2013). At the site studied by Young et al. (2015), the majority of the glacier area resides on a high plateau, with a narrow outlet glacier occupying a deeply incised valley, terminating in an open, wide valley floor. Although the moraine chronology developed by Young et al. (2015) itself is convincing, it is possible that the terminal position is relatively insensitive to small late Holocene ELA changes. This is because once the ELA drops below the plateau and is within the narrow outlet glacier, each incremental ELA lowering increases the accumulation area only slightly. In contrast, DIC has a more symmetric hypsometry (i.e., the glacierized area is more evenly distributed over glacier elevation). Consequently, we expect that the correlation between ELA change and glacier dimension response should be more linear than for Naqsaq Glacier. Both the plant radiocarbon dates and glacier model simulations require increased summer cooling between ∼1000 CE and the LIA.  claus et al., 2016). The past2k CESM simulation suggests that the period from 0 to 1000 CE was, on average, ∼ 0.24 • C cooler than 1 CE control conditions (Fig. 6), and that the period from 1000 to 1900 CE was on average ∼ 0.30 • C cooler than the preceding millennium, meaning that the period from 1000 to 1900 CE was, on average, ∼ 0.54 • C cooler than 1 CE conditions, and that the second millennium CE temperature decline was dominated by cooling through the LIA (Fig. 6). The period encapsulating the LIA, from 1250 to 1850 CE, was on average 0.40 • C cooler than 1 to 1000 CE. In general, these temperature trends agree well with the glacier model outputs in this study (Table 2). Both the CESM simulation and Kaufman et al. (2009) Arctic temperature reconstructions show less overall warming following the LIA up to present than the cumulative warming from the glacier model. It is possible that the imposed linear warming in our glacier model does not accurately capture the nonlinear warming that is likely taking place in the region, or that the CESM simulation and temperature reconstruction may not capture local factors influencing warming at DIC. However, the overall agreement between the climate and glacier models and the similarity of both to regional temperature reconstructions (Fig. 6, Table 2) support the contention that the record of DIC expansion derived from the death ages of entombed plants faithfully records the average cooling that occurred in the region over the past 2000 years.

Conclusions
We have used radiocarbon-dated in situ tundra plants exposed by retreating ice margins to construct a spatially and temporally constrained record of ice cap expansion over the past 2 kyr. DIC grew between ∼ 26 BCE and 1000 CE, then advanced episodically at ∼1000, 1200, and 1500 CE, reaching its Neoglacial maximum during the LIA. The LIA was terminated by the warming of the 20th century.
Glacier model simulations show that a minimum average cooling of ∼ 0.44 • C is required to match the radiocarbon constrained pattern of ice expansion over the last 2000 years, with the period from 1000 to 1850 CE being on average ∼ 0.25 • C colder than the preceding ∼ 1000 years. A climate model simulation for the past 2 kyr driven by natural and anthropogenic forcings, as well as reconstructed Arctic temperatures show similar summer temperature decreases, reinforcing the glacier modeling conclusions. Both the radiocarbon record and climate model simulations indicate that the coldest interval of the past ∼ 2000 years was during the LIA (1250-1850 CE). Glacier model simulations matching observed ice-cap retreat since 1900 CE suggest that a cumulative warming of 2.7 • C over the last ∼ 100 years has reversed ∼ 1000 years of ice cap expansion under only ∼ 0.25 • C cooling, suggesting modern warming is unprecedented over the past 2 kyr.
At the present rate of warming, DIC will likely disappear before ∼ 2100. Future collection of in situ plants exposed as DIC continues to retreat will both extend the record of ice cap advance and provide more constraints for modeling of ice cap activity and climate perturbations during Neoglaciation. Data availability. All data, model code, and primary model input and output are available in a PANGEA data repository (https://doi.org/10.1594/PANGAEA.881156).