Bering Glacier Research Articles

Evolution of a Surge Cycle of the Bering‐Bagley Glacier System From Observations and Numerical Modeling

AbstractThe Bering‐Bagley Glacier System (BBGS), Alaska, Earth's largest temperate surging glacier, surged in 2008–2013. We use numerical modeling and satellite observations to investigate how surging in a large and complex glacier system differs from surging in smaller glaciers for which our current understanding of the surge phenomenon is based. With numerical simulations of a long quiescent phase and a short surge phase in the BBGS, we show that surging is more spatiotemporally complex in larger glaciers with multiple reservoir areas forming during quiescence which interact in a cascading manner when ice accelerates during the surge phase. For each phase, we analyze the simulated elevation‐change and ice‐velocity pattern, infer information on the evolving basal drainage system through hydropotential analysis, and supplement these findings with observational data such as CryoSat‐2 digital elevation maps. During the quiescent simulation, water drainage paths become increasingly lateral and hydropotential wells form indicating an expanding storage capacity of subglacial water. These results are attributed to local bedrock topography characterized by large subglacial ridges that dam the down‐glacier flow of ice and water. In the surge simulation, we model surge evolution through Bering Glacier's trunk by imposing a basal friction representation that mimics a propagating surge wave. As the surge progresses, drainage efficiency further degrades in the active surging‐zone from its already inefficient, end‐of‐quiescence state. Results from this study improve our knowledge of surging in large and complex systems which generalizes to glacial accelerations observed in outlet glaciers of Greenland, thus reducing uncertainty in modeling sea‐level rise.

Application of a structured decision‐making process in cryospheric hazard planning: Case study of Bering Glacier surges on local state planning in Alaska

AbstractSurging glaciers are glaciers that experience rapidly accelerated glacier flow over a comparatively short period of time. Though relatively rare worldwide, Alaska is home to the largest number of surge‐type glaciers globally. However, their impact on the broader socioecological system in the state is both poorly understood and under‐researched, which poses a challenge in developing appropriate sustainability decisions in Alaska. We investigated how the surge patterns of the Bering Glacier in Alaska have potentially devastating effects on the local ecological biodiversity of its watershed via a structured decision‐making analysis of the different possible consequences. Specifically, this analysis was conducted to explore the various outcomes of a Bering Glacier surge particularly if humans have an increased presence near the glacier due to the area potentially becoming a state park. This work explored the benefits of applying a risk and decision analytical framework in a cryosphere context, to better understand the socioeconomic impact of glacier surges. This is a novel approach in which a decision analysis tool was used to better understand an environmental sustainability challenge, offering an innovative method to support the achievement of the United Nations Sustainability Development Goals in Alaska. We therefore emphasise the need for integrated biophysical and socioeconomic analyses when it comes to understanding glacier hazards. Our research highlights the importance of understanding and researching biophysical changes as well as using a structured decision‐making process for complicated hazard planning scenarios, exemplified via glaciated regions in Alaska, in order to create adaptation strategies that are sustainable and encompass the range of possible outcomes.

Interrogating a Surging Glacier With Seismic Interferometry

AbstractSeismic interferometry can be spectacularly sensitive to temporal changes within the Earth. Zhan (2019, https://doi.org/10.1029/2019GL082411) analyzed 12 years of data from seismic stations situated near Bering Glacier, Alaska, that spanned a 2008–2011 surge. Interferometrically isolated glacier‐crossing surface waves revealed an anisotropic seismic wave speed decrease during the surge. 1‐D and 3‐D modeling suggests pressurized englacial water flow being stored and guided by a network of basal crevasses, thereby sustaining high water pressure at the bed and supporting the high surge velocities. This result is demonstrative of the potential for seismic interferometry and other cryoseismological methods to continuously monitor and reveal changes within dynamic glacial‐hydrological systems.

Late Pleistocene marine resources from the Bering Glacier Foreland and human coastal migration in the northern Gulf of Alaska region

AbstractRecent research on the Bering Glacier forelands in the northern Gulf of Alaska provides new insights into late Pleistocene/early Holocene shorelines, providing a favorable route for human migration as early as ~16,000 cal yr BP. This route included an irregular coastline with embayments and islands offering protection from the open ocean; edible marine invertebrates dating from 15,000 to 5,500 cal yr BP; and marine vertebrates dating as early as 16,000 cal yr BP. The latter included walrus (Odobenus rosmarus), bearded seal (Erignathus barbatus), and ringed seal (Phoca cf. hispida), all associated with pack ice conditions unlike those present today. While this ecosystem could have supported humans migrating along the coastline, and coastal refugia may have existed elsewhere in the region, coastal archaeological sites in the northern Gulf of Alaska and southwest Alaska are no older than ~9,500 cal yr BP. This suggests that the earliest sites have been eroded or destroyed, that the earliest migrants ignored available marine resources, and/or that these migrants did not use a coastal route. In contrast, the earliest archaeological sites in southeast Alaska date to ~12,500 cal yr BP, suggesting migration from interior Alaska to the coast somewhere east of the Copper River delta.

Late Quaternary glacial dynamics and sedimentation variability in the Bering Trough, Gulf of Alaska

Abstract Ice dynamics, tectonic setting, and sediment supply are the key parameters controlling the architecture of high-latitude margins and the formation of trough mouth fans (TMFs). Current understanding of these archives of paleo–ice streams is based on studies of ice sheets adjacent to stable, passive margins, while the behavior of active, convergent glacier-influenced margins remains relatively unconstrained. We integrate high-resolution seismic data and chronology from Integrated Ocean Drilling Program Expedition 341 cores in southeast Alaska across the actively converging Yakutat terrane margin to examine the late Quaternary evolution of the Bering Glacier, the largest outlet glacier of the poorly understood Cordilleran Ice Sheet (CIS). We interpret at least eight glacial advances to the shelf break since the end of the mid-Pleistocene transition, showing a more dynamic CIS than hitherto realized. During the past ∼130 k.y., the temperate, meltwater-charged Bering Glacier delivered ∼925 km3 of sediment to the shelf and slope, providing one of the highest rates of sustained sediment accumulation (5–10 m/k.y.) ever reported globally. Rapid formation of a TMF, reaching ∼600 m thick in ∼130 k.y., emphasizes the extreme sediment flux that can be produced by wet-based glacial systems, and its critical role in the development of high-latitude margin stratigraphy. TMF formation despite initially steep, tectonically controlled slopes in this active setting reflects an autogenic shift in the evolution of the Bering Trough, suggesting that major transitions between sedimentary regimes need not reflect some externally driven change in climate variability.

Melt water input from the Bering Glacier watershed into the Gulf of Alaska

The annual runoff from the melting of large glaciers and snow fields along the northern perimeter of the Gulf of Alaska is a critical component of marine physical and biological systems; yet, most of this freshwater is not measured. Here we show estimates of melt for the watershed that contains the largest and longest glacier in North America, the Bering Glacier. The procedure combines in situ observations of snow and ice melt acquired by a long-term monitoring program, multispectral satellite observations, and nearby temperature measurements. The estimated melt is 40 km3 per melt season, ± 3.0 km3, observed over the decadal period, 2002–2012. As a result of climate change, these estimates could increase to 60 km3/yr by 2050. This technique and the derived melt coefficients can be applied to estimate melt from Alaska to Washington glaciers.

Low‐frequency radar sounding of temperate ice masses in Southern Alaska

We present the Warm Ice Sounding Explorer (WISE), a low‐frequency (2.5 MHz) radar for the sounding of temperate ice. WISE deployment in southern Alaska in 2008 and 2012 provides comprehensive measurements of glacier thickness, reveals deep valleys beneath glaciers and the full extent of zones grounded below sea level. The east branch of Columbia Glacier is deeper that its main branch and remains below sea level 20 km farther inland. Ice is 1000 m deep on Tazlina Glacier. On Bering glacier, two sills separate three deep bed depressions (>1200 m) that coincide with the dynamic balance lines during surges. The piedmont lobe of Malaspina Glacier and the lower reaches of Hubbard Glacier are entirely grounded below sea level 40 and 10 km, respectively, from their termini. Knowledge of ice thickness in these regions helps better understand their glacier dynamics, mass balance, and impact on sea level.

Active tectonics of the St. Elias orogen, Alaska, observed with GPS measurements

We use data from campaign and continuous GPS sites in southeast and south central Alaska to constrain a regional tectonic block model for the St. Elias orogen. Active tectonic deformation in the orogen is dominated by the effects of the collision of the Yakutat block with southern Alaska. Our results indicate that ~37 mm/yr of convergence is accommodated along a relatively narrow belt of N‐NW dipping thrust faults in the eastern half of the orogen, with the present‐day deformation front running through Icy Bay and beneath the Malaspina Glacier. Near the Bering Glacier, the collisional thrust fault regime transitions into a broad, northwest dipping décollement as the Yakutat block basement begins to subduct beneath the counterclockwise rotating Elias block. The location of this transition aligns with the Gulf of Alaska shear zone, implying that the Pacific plate is fragmenting in response to the Yakutat collision. Our model indicates that the Bering Glacier region is undergoing internal deformation and could correspond to the final stage of offscraping and accretion of sediments from the Yakutat block prior to subduction. Predicted block motions at the western edge of the orogen suggest that the crust is laterally escaping along the Aleutian fore arc.

Mass balance, runoff and surges of Bering Glacier, Alaska

Abstract. The historical net, ablation and accumulation daily balances, as well as runoff of Bering Glacier, Alaska are determined for the 1951–2011 period with the PTAA (precipitation-temperature-area-altitude) model, using daily precipitation and temperature observations collected at the Cordova and Yakutat weather stations, together with the area-altitude distribution of the glacier. The model mean annual balance for this 61 yr period is −0.6 m w.e., the accumulation balance is +1.4 and the ablation balance is −2.0 m w.e. Average annual runoff is 2.5 m w.e. Periodic surges of this glacier transport large volumes of ice to lower elevations where the ablation rate is higher, producing more negative balances and increasing runoff. Runoff from Bering Glacier (derived from simulated ablation and precipitation as rain) is highly correlated with four of the glacier surges that have occurred since 1951. Ice volume loss for the 1972–2003 period measured with the PTAA model is 2.7 km3 w.e. a−1 and closely agrees with losses for the same period measured with the geodetic method. It is proposed that the timing and magnitude of daily snow accumulation and runoff, both of which are controlled by the glacier's area-altitude distribution and are calculated with the PTAA model, can be used to determine the probability that a glacier will surge.

A continental shelf sedimentary record of Little Ice Age to modern glacial dynamics: Bering Glacier, Alaska

The Bering Glacier System is the world's largest surging temperate glacier with seven events occurring over the past century under a range of north Pacific climatic conditions. Onshore records reveal changes in glacial termini positions and evidence of late Holocene glacial advances, but the Little Ice Age (LIA) record of potential glacial surging and associated flooding has not been examined. A 13.6m-long jumbo core collected on the adjacent continental shelf reveals a 600-yr-long record of sedimentation associated with changing glacifluvial discharge. The chronology is based on 210Pb geochronology and five radiocarbon dates, and the core can be separated into three distinct lithologic units based on the examination of X-radiographs and physical properties: (1) an uppermost unit dating from ∼125calyrBP to the present characterized by bioturbated mud interbedded with laminated, thick (5–20cm) low-bulk density clay-rich beds; (2) a middle unit dating from ∼120–400calyrBP that includes numerous interlaminated-to-interbedded low- and high-bulk density beds with infrequent evidence of bioturbation; thick laminated clay-rich beds are rare; (3) a lowermost unit that predates ∼400calyrBP and is composed of rare laminated beds grading down into mottled to massive mud. In each of these units, the laminated lithofacies from this mid-shelf location indicates both flood deposition and likely sediment transport in the wave-current bottom-boundary layer. The thick low-density, clay-rich beds in the uppermost unit correlate with historic outburst floods associated with known surge events. Based on previous terrestrial studies, the terminus was at its Holocene Neoglacial maximum extent close to the modern coastline at some point in the middle to late stages of the LIA in southern Alaska (100–350calyrBP). During the LIA, preservation of bioturbated intervals is rare while laminated intervals are common. This style of interbedding indicates frequent (<10yr recurrence interval) event-scale mud deposition, suggesting that frequent summer flooding and redistribution by winter storms were more prevalent during the LIA rather than the outburst flooding typical of the past century. Rare event-scale bedding indicative of outburst flooding and possible surge events is found within the middle unit, and may correspond to periods with similar climatic trends as in the 20th century. The infrequent deposition of event layers in the lowermost unit could be attributed to the less frequent flooding and/or enhanced diversion of glacial drainage to the eastern terminus instead of present day Seal River. The thickness and depositional frequency of event-scale bedding can be related to Gulf of Alaska tree-ring proxy temperature reconstructions, where more numerous event bed formation occurs when there are more frequent, higher-amplitude temperature excursions. These frequent fluctuations may have prevented the decadal-long periods of positive mass balance required to enable numerous surge events during this period.

Elevation changes of Bering Glacier System, Alaska, from 1992 to 2010, observed by satellite radar altimetry

Satellite radar altimetry has been extensively used to estimate ice-sheet elevation change rates over Antarctica and Greenland from repeated measurements of surface elevation changes through either crossover or collinear (repeat-track) analysis. However, there has been no attempt to use satellite radar altimetry over mountain glaciers because the instrument often loses lock or does not provide usable measurements over the rough terrains and steep slopes due to its large footprint size. In this study, we demonstrate the feasibility of satellite radar altimetry to estimate elevation change rates over the Bering Glacier System in Alaska for the period of 1992–2010 from TOPEX/Poseidon and Envisat radar altimeter measurements. We are able to observe the 1993–1995, and 2008–2011 surge events from the T/P and Envisat time series, respectively. We also observe the accelerated elevation decreases in 2002–2007, after slightly negative or near nil elevation changes in 1996–2001, which are related to the temperature and snow depth variations. Our method can be applied to other wide (>7km) glaciers worldwide, and provide new insights into the behavior of glaciers responding to climate change.

Examination of the interplay between glacial processes and exhumation in the Saint Elias Mountains, Alaska

The combination of large, temperate glaciers and rapid crustal convergence in the Saint Elias Mountains (southeastern Alaska, USA, and Yukon Territory and British Colombia, Canada) provides an exceptional opportunity to study the interactions between the tectonic and surface processes that have shaped most active orogens on Earth during much of the Quaternary. This research first provides a review of thermochronometric data sets recording exhumation under two major glacier systems of the Saint Elias Mountains, the Bagley-Bering and the Seward-Malaspina systems. These data sets are integrated over the single glacier systems and used in conjunction with glaciological data to investigate the interactions of glacial erosion and tectonics. Despite their proximity, the glaciological processes and geological settings of these two glacial systems differ significantly. On the east side of the orogen, sediments eroded from bedrock underneath the Malaspina Glacier reflect regions of rapid erosion under the slowly moving Seward Ice Field. Because the Seward Ice Field overlies a localized zone of major faulting and rapid exhumation, the strained and fractured bedrock is primed for erosion. On the west side, the Bering Glacier is the primary outlet for the Bagley Ice Field, which covers half of the crest of the orogen; however, few if any of the sediments at its terminus originate from under the Bagley Ice Field. Sediment transport is likely hindered by subglacial freeze-on processes that reduce the sediment-carrying capacity of subglacial rivers, though glacial surges are likely exceptions that deposit sediment far beyond the active margin of the glacier. Our study concludes that the widely invoked concepts of glacial erosion should be used with caution, as oversimplification can fail to account for important site-specific differences in geologic and glacial conditions.

The propagation of a surge front on Bering Glacier, Alaska, 2001–2011

Bering Glacier, Alaska, USA, has a ∼20 year surge cycle, with its most recent surge reaching the terminus in 2011. To study this most recent activity a time series of ice velocity maps was produced by applying optical feature-tracking methods to Landsat-7 ETM+ imagery spanning 2001-11. The velocity maps show a yearly increase in ice surface velocity associated with the down-glacier movement of a surge front. In 2008/09 the maximum ice surface velocity was 1.5 ±0.017 km a-1 in the mid-ablation zone, which decreased to 1.2 ±0.015 km a-1 in 2009/10 in the lower ablation zone, and then increased to nearly 4.4 ± 0.03 km a-1 in summer 2011 when the surge front reached the glacier terminus. The surge front propagated down-glacier as a kinematic wave at an average rate of 4.4 ±2.0 km a-1 between September 2002 and April 2009, then accelerated to 13.9 ± 2.0 km a-1 as it entered the piedmont lobe between April 2009 and September 2010. The wave seems to have initiated near the confluence of Bering Glacier and Bagley Ice Valley as early as 2001, and the surge was triggered in 2008 further down-glacier in the mid-ablation zone after the wave passed an ice reservoir area.

The propagation of a surge front on Bering Glacier, Alaska, 2001–2011

AbstractBering Glacier, Alaska, USA, has a ∼20 year surge cycle, with its most recent surge reaching the terminus in 2011. To study this most recent activity a time series of ice velocity maps was produced by applying optical feature-tracking methods to Landsat-7 ETM+ imagery spanning 2001-11. The velocity maps show a yearly increase in ice surface velocity associated with the down-glacier movement of a surge front. In 2008/09 the maximum ice surface velocity was 1.5 ±0.017 km a-1 in the mid-ablation zone, which decreased to 1.2 ±0.015 km a-1 in 2009/10 in the lower ablation zone, and then increased to nearly 4.4 ± 0.03 km a-1 in summer 2011 when the surge front reached the glacier terminus. The surge front propagated down-glacier as a kinematic wave at an average rate of 4.4 ±2.0 km a-1 between September 2002 and April 2009, then accelerated to 13.9 ± 2.0 km a-1 as it entered the piedmont lobe between April 2009 and September 2010. The wave seems to have initiated near the confluence of Bering Glacier and Bagley Ice Valley as early as 2001, and the surge was triggered in 2008 further down-glacier in the mid-ablation zone after the wave passed an ice reservoir area.

Bering Glacier surge 2011: analysis of laser altimeter data

AbstractThe Bering Glacier–Bagley Icefield system in Alaska is currently surging (2011). Large-scale elevation changes and small-scale elevation-change characteristics are investigated to understand surge progression, especially mass transport from the pre-surge reservoir area to the receiving area and propagation of the kinematic surge wave as manifested in heavy crevassing characteristic of rapid, brittle deformation. This analysis is based on airborne laser altimeter data collected over Bering Glacier in September 2011. Results include the following: (1) Maximal crevasse depth is 60 m, reached in a rift that separates two deformation domains, indicative of two different flow regimes. Otherwise surge crevasse depth reaches 20–30 m. (2) Characteristic parameters of structural provinces are derived by application of geostatistical classification. Parameters include significance and spacing of crevasses, surface roughness and crevasse-edge curvature (indicative of crevasse age). A classification based on these parameters serves to objectively discriminate structural provinces, indicative of surge progression down-glacier and up-glacier. (3) Elevation changes from 2011 and 2010 altimetry show 40–70 m surface lowering in the reservoir area in lower central Bering Glacier and 20–40m thickening near the front in Tashalich arm. Combining elevation changes with results of crevasse profilometry and pattern analysis, the rapid progression of the surge can be mathematically–physically reconstructed.

Bering Glacier and Bagley Ice Valley surge 2011: crevasse classification as an approach to map deformation stages and surge progression

AbstractThe dynamics of a surge is manifested in the crevasse patterns: literally, deformation state frozen in ice. This basic observation is utilized as the concept of an automated approach to map and analyze deformation stages and progression of surge kinematics. The classification method allows imagery to be used as geophysical data and is applied to aerial observations (photographic and video imagery, GPS data) collected in September 2011 during the surge of the Bering Glacier–Bagley Ice Valley system, Alaska, USA. As the third dimension that complements two-dimensional imagery, ice-surface elevation is observed using aerial laser altimetry. The classification method builds on concepts from signal processing, geostatistical data analysis and neural networks. Steps include calculation of generalized directional vario functions from image data and composition into feature vectors. The vario function operates as an information filter that retains spatial characteristics at an intermediate scale that captures crevasse spacing, anisotropy and other generalized roughness properties. Association of feature vectors to crevasse classes and hence deformation types employs a connectionist algorithm. In general, the connectionist–geostatistical classification allows the mapping of kinematic changes in crevassed glaciers.Dedicated to the memory of Austin Post

Surge dynamics on Bering Glacier, Alaska, in 2008–2011

Abstract. A surge cycle of the Bering Glacier system, Alaska, is examined using observations of surface velocity obtained using synthetic aperture radar (SAR) offset tracking, and elevation data obtained from the University of Alaska Fairbanks LiDAR altimetry program. After 13 yr of quiescence, the Bering Glacier system began to surge in May 2008 and had two stages of accelerated flow. During the first stage, flow accelerated progressively for at least 10 months and reached peak observed velocities of ~ 7 m d−1. The second stage likely began in 2010. By 2011 velocities exceeded 9 m d−1 or ~ 18 times quiescent velocities. Fast flow continued into July 2011. Surface morphology indicated slowing by fall 2011; however, it is not entirely clear if the surge is yet over. The quiescent phase was characterized by small-scale acceleration events that increased driving stresses up to 70%. When the surge initiated, synchronous acceleration occurred throughout much of the glacier length. Results suggest that downstream propagation of the surge is closely linked to the evolution of the driving stress during the surge, because driving stress appears to be tied to the amount of resistive stress provided by the bed. In contrast, upstream acceleration and upstream surge propagation is not dependent on driving stress evolution.

Icequake locations and discrimination of source and path effects with small aperture arrays, Bering Glacier terminus, AK

The complex source processes associated with both glacier calving and the breakup of icebergs, combined with commonly heterogeneous periglacial seismic velocity structure, can result in complicated seismic records. Key features of the waveforms, which are typically characterized by low‐amplitude or emergent first‐arrivals and long‐duration, narrow‐band codas, have been attributed to either source processes or propagation path effects. This uncertainty must be addressed in order for seismic data to be effective for studying the calving process as it relates to terminus dynamics. In this study, we use sets of 3‐element arrays of 3‐component geophones and infrasound sensors to locate calving and iceberg breakup events and isolate path effects in the seismic records obtained near the Bering Glacier terminus in the summer of 2010. Using waveform correlation, we treat each array as an antenna and determine the direction to the source and apparent velocity of the wavefield across the array. The initial few (∼3) cycles ofPwaves recorded from an array, beam formed to identify coherent arrivals for each event, are useful for deriving a propagation azimuth and apparent velocity, allowing for location of events using a small number of similar arrays. We locate 125 calving and iceberg breakup events near the terminus with this method. We also demonstrate that the longer‐lived narrow‐band coda is not coherent across individual arrays, suggesting that the narrow‐band coda observed at the Bering Glacier is attributable to a path effect rather than to the source process. The large number of iceberg breakup events that we located has important implications for other calving glaciers where icebergs are present, and calving rates may be erroneously overestimated from the seismic data if their contribution is not taken into account.

The Suckling Hills Fault, Kayak Island Zone, and accretion of the Yakutat microplate, Alaska

[1] The Suckling Hills and Kayak Island are isolated mountain blocks located along strike from each other within the foreland of the St. Elias orogen in southern Alaska. These blocks preserve an erosional surface that was deformed by slip on northwest-dipping reverse faults in the Pleistocene. We suggest that the Suckling Hills Fault and Kayak Island Zone form a segmented fault network that links with the Bering Glacier structure to the north. This fault network separates the central Yakataga fold and thrust belt from complex, multiply deformed structures in the western syntaxis. Ongoing accretion of the Yakutat microplate to North America results in translation of structures of the fold and thrust belt into the western syntaxis. The composite Suckling Hills Fault, Kayak Island Zone, and Bering Glacier structure may have formed because the older structures of the fold and thrust belt were unfavorably oriented within the western syntaxis region. This pattern of deformation provides a template for understanding the complex deformation within the core of the western syntaxis and predicts refolding and straightening of the western syntaxis margin with continued accretion. This study provides an analog for structural overprinting and changing deformation patterns through time in orogenic corners.

The sedimentary architecture of outburst flood eskers: A comparison of ground-penetrating radar data from Bering Glacier, Alaska and Skeithararjokull, Iceland

We present ground-penetrating radar (GPR) profiles that reveal the sedimentary architecture of an esker deposited during a surge-associated outburst flood at the Bering Glacier, Alaska. The wide, up-flow end of the esker contains a transition from large backset beds to large foreset beds interpreted to reflect composite macroform development in an enlarged part of the conduit. By contrast, the narrow, down-flow portion of the esker is dominated by plane beds interpreted to have been deposited where the conduit was constricted and the flow was faster. A previously studied outburst esker at Skeidararjokull, Iceland, has a similar morphology and stratigraphic architecture. This suggests that outburst floods generate distinct depositional signatures in eskers, both in terms of morphology and sedimentary architecture. Identification of these distinct signatures in ancient eskers will help assess the paleohydraulic conditions under which ancient eskers formed and, by extension, the nature of meltwater drainage systems beneath the Laurentide and Eurasian ice sheets.