Impact of glacial isostatic adjustment on cosmogenic surface-exposure dating

Calculating cosmogenic-nuclide surface-exposure ages is critically dependent on a knowledge of the altitude of the sample site. Changes in altitude have occurred through time as a result of glacial isostatic adjustment (GIA), potentially altering local nuclide production rates and, therefore, surface-exposure ages. Here we assess the impact of GIA on surface-exposure dating by calculating global time-dependent production rates since the Last Glacial Maximum using surface elevations that were cor- rected and uncorrected for GIA. We ﬁ nd that the magnitude of the GIA effect is spatially and temporally variable. Nuclide production could be reduced by up to 50% in the interior of large ice masses (in North America, Scandinavia and West Antarctica) at times of maximum glacial isostatic depression. Although smaller, the effect is still signi ﬁ cant at ice sheet margins, where nuclide production is reduced by > 5% and potentially > 10%, making exposure ages older in those areas. Away from the ice sheet margins, land surfaces can be isostatically elevated, which can increase nuclide production by > 5% and, therefore, make exposure ages younger. Areas that were more recently exposed or that are distal to large ice masses will generally be less affected. Importantly, we ﬁ nd that the effect at the primary 10 Be production calibration sites is < 1%. Applying a GIA correction to surface-exposure data may help resolve mismatches between some chronologies, but not necessarily in all regions, implying that additional factors may need to be considered. Past atmospheric changes can amplify or reduce the impact of GIA on nuclide production, and the combined effects should be fully accounted for in the future. These time-dependent in ﬂ uences on surface-exposure dating have potentially important implications for interpreting chronologies and for using the data to constrain ice sheet models.


Introduction
Cosmogenic-nuclide surface-exposure dating is a widely used approach for constraining the timing of past geomorphic events and, in particular, for reconstructing ice margin history during the Quaternary. The technique is heavily dependent on a knowledge of the nuclide production rate and how this rate has varied through time (Gosse and Phillips, 2001;Lifton et al., 2014). The altitude of a rock surface exerts a fundamental control on nuclide production and, if not appropriately estimated, can lead to an inaccurate surface-exposure age. Over glacial-interglacial cycles, the expansion and recession of ice sheets caused the altitude of the land surface to change through a process called glacial isostatic adjustment (GIA) (e.g. Peltier et al., 2015;Whitehouse, 2018). The aim of this paper is to provide a first global assessment of the impact of GIA on nuclide production and the implications for surfaceexposure dating.

Cosmogenic nuclide production with altitude
Altitude determines the atmospheric pressure at a location, which controls the site-specific nuclide production rate and, therefore, the surface-exposure age. Nuclide production reduces approximately exponentially with atmospheric depth, which is dependent on atmospheric pressure and acceleration due to gravity (Desilets and Zreda, 2003;Dunai, 2000;Stone, 2000). Primary cosmic rays collide with oxygen and nitrogen nuclei near the top of the atmosphere, where a cascade of protons, neutrons and other secondary particles is produced. Energy is then lost as the cascade passes through increasingly denser air, causing progressively fewer interactions by secondary cosmic rays to occur towards the Earth surface (Gosse and Phillips, 2001). Consequently, nuclide production rates increase with altitude (Lal, 1991;Lal and Peters, 1967;Sato et al., 2008).
Scaling models are used to calculate the site-specific production rate as a function of altitude (i.e. atmospheric pressure), as well as latitude (i.e. geomagnetic cutoff rigidity for incoming cosmic rays). By convention, scaling factors are referenced to conditions at sea level (1013.25 hPa or 1033.2 g cm 2 ) and high geomagnetic latitude (0 GV at 90 N/S), often referred to as 'SLHL'. While scaling models account for time-dependent changes in the geomagnetic field and solar modulation (Balco et al., 2008;Desilets and Zreda, 2003;Dunai, 2000;Lifton et al., 2005Lifton et al., , 2014, they assume that 1) the atmospheric pressure at sea level and 2) the position of a sample site relative to SLHL have remained constant through time. The total mass of the atmosphere above sea level, and therefore the atmospheric pressure at sea level, should have been effectively constant over a glacial-interglacial cycle; changes in global mean atmospheric pressure at the altitude of present-day sea level have been estimated to be < 1 hPa (M eli eres et al., 1991), which corresponds to a <1% change in the nuclide production rate (Stone, 2000). Thus, the first assumption likely holds true. However, the position of a sample site, notably the elevation of a site relative to sea level, may have varied significantly in the past. The largest changes in surface elevation over glacial-interglacial timescales occur from GIA.

Glacial isostatic adjustment (GIA)
GIA describes the viscoelastic response of the solid Earth to time-dependent changes in ice and ocean loading over the course of a glacial-interglacial cycle. Two factors control the spatial and temporal characteristics of surface elevation change due to GIA: the history of global ice sheet change and the rheology of the solid Earth, in particular the mantle viscosity. Geological, geomorphological, geodetic and seismic evidence has been used to constrain these two components in combination with physically-based modelling (e.g. Lambeck et al., 1998;Peltier et al., 2015;van der Wal et al., 2015;Whitehouse et al., 2012a;Whitehouse et al., 2012b), but uncertainties remain due to data gaps and nonuniqueness. For example, there exists a trade-off between the details of ice sheet change and the mantle viscosity that can be invoked to explain ongoing, GIA-related solid Earth deformation.
Numerical models can be used to predict the magnitude and extent of surface elevation change during a glacial-interglacial cycle. In addition to accounting for the relaxation time of the mantle, such models also account for the damping effect of the elastic lithosphere. As a result, subsidence (>10 2 m) is predicted both beneath and adjacent to the major ice sheets at the Last Glacial Maximum (LGM), while uplift (<10 2 m) is predicted in broad 'peripheral bulge' regions surrounding the ice sheets ( Fig. 1). Due to the viscoelastic nature of the mantle, it takes many thousand years for the solid Earth to return to a state of equilibrium following deglaciation.
Some previous cosmogenic nuclide studies have tried to account for changes in surface elevation resulting from GIA, either using local relative sea level records (e.g. Goehring et al., 2012;Rinterknecht et al., 2006;Young et al., 2013) or model-derived estimates of GIA (e.g. Cuzzone et al., 2016;Suganuma et al., 2014;Ullman et al., 2016). However, there has not yet been a systematic global assessment of the pattern and magnitude of potential GIA effects.

Methodology
Here we evaluate the potential global effects of GIA on cosmogenic nuclide production and corresponding surface-exposure ages since the LGM. Time-varying surface-elevation change due to GIA was determined using the global ICE-6G ice model (Peltier et al., 2015) and a three-layer approximation of the VM2 Earth model (see Supplementary Material). Nuclide production was calculated using the 'LSD' time-dependent scaling model (Lifton et al., 2014), which was modified to allow for time-varying changes in atmospheric pressure (Supplementary Material). We focused on the production of beryllium-10 ( 10 Be), as it is the most widely-used nuclide, but the effects that we discuss will be similar for other nuclides.
Our objectives were three-fold. Firstly, we assessed the impact of global surface elevation change on production rates during the last deglaciation. This was achieved by calculating nuclide production at 20 ka, 15 ka, 10 ka and 5 ka BP for the modern surface topography and a surface topography corrected using GIA model output (Supplementary Material). Our calculations were done for grid cells inside as well as outside of former ice sheet margins in order to highlight the effects at potential ice-free areas in the ice sheet interior (i.e. nunataks). We prescribed zero shielding for terrain and surface cover, zero erosion, a rock density of 2.6 g cm 3 , and a rock surface thickness of 1 cm. For each interval, the effects from GIA were determined by considering the difference between the global nuclide production grids that were corrected and uncorrected for past elevation change.
Secondly, we examined the effect of GIA at production rate calibration sites. Applying a GIA correction to calculate timedependent nuclide production will be inherently inaccurate if the calibration sites used to determine the production rate have been affected by GIA. We looked at the effects at the six primary 10 Be calibration sites of Borchers et al. (2016) by calculating nuclide production for each sample within the dataset since the independently dated 'true' ages of the sites, and compared the results to production calculated without the GIA correction. For these calculations, the published location (latitude, longitude and modern elevation), density, thickness and shielding of each sample were used.
Thirdly, we evaluated how GIA-corrected nuclide production impacts the corresponding surface-exposure ages. This was carried out using published 10 Be surface-exposure data that was previously compiled by Jakob Heyman (http://expage.github.io/). To highlight the potential effects of GIA, we included all samples with (uncorrected) exposure ages between 5 ka and 25 ka, and did not exclude samples that were considered outliers in the original publications. We have assumed that all published sample elevations were reported relative to the geoid, which is necessary for accurate atmospheric pressure estimation and thus nuclide production calculation, rather than to the ellipsoid, which is typically measured by GPS. All 10 Be concentrations were normalised to 07KNSTD (Nishiizumi et al., 2007) prior to calculating surface-exposure ages, which was then performed using the CRONUScalc framework ) and LSD scaling model for both the modern elevation and GIA-corrected time-dependent elevation for each sample. We have developed a tool, available as part of iceTEA, to allow users to assess the potential effects of GIA on their own surface-exposure data (http://ice-tea.org; Jones et al., 2019).

Results
The impact of GIA on nuclide production reflects the magnitude and spatial pattern of modelled surface elevation change (Fig. 2). As might be expected, the largest effects are in the vicinities of major ice volume change (Laurentide, Cordilleran, Fennoscandian, Greenland and West Antarctic ice sheets), but the magnitude of the effects elsewhere are potentially significant. Nuclide production reduces by as much as 50% relative to present at the centre of ice masses near to the time of maximum ice expansion, at approximately 20 ka. At the same time, production at the peripheral bulge of large ice masses increases by up to 8% when land surface elevations are corrected for GIA (Fig. 2). Regions with smaller ice masses have less significant effects. In Iceland and southern Patagonia, areas may have had nuclide production reduced by >10% early during deglaciation, but in the vicinity of the British-Irish Ice Sheet and in New Zealand, the effects of GIA on production were <5% (Fig. 2). The ice masses in the Himalayas, European Alps, Caucasus and Andes were not included here (Supplementary Material), but based on this assessment, their effects on nuclide production should not be more than 10% and were probably less than 5%. The spatial distribution of the production rate effects evolves through time as the ice masses retreat to their modern configurations (Fig. 2), highlighting the need to use timedependent scaling frameworks to accurately account for these changes.
Surface elevation changes at the primary 10 Be production rate calibration sites of Borchers et al. (2016) were relatively small as the sites were all sufficiently distal to the large ice masses to be negligibly impacted by GIA (Fig. 3). For all sites other than PPT, the difference between nuclide production that was corrected and uncorrected for GIA effects is <0.3%. The samples at PPT were located at the edge of the peripheral bulge of the Cordilleran and Laurentide ice sheets, but this only increases the time-averaged production rate by 2.6%. The true impact at this site may be even smaller as local isostatic effects associated with the presence of Lake Bonneville could have partially offset the regional GIA signal (Adams and Bills, 2016). Overall, this primary 10 Be production rate calibration dataset is influenced by the effects of GIA by <1%, which is well within the uncertainty of the measured sample concentrations (Fig. 3).
Existing surface-exposure ages, recalculated using the GIAcorrected surface elevations and published 10 Be production calibration dataset, become both younger and older than ages uncorrected for changes in surface elevation change ( Fig. 4;   Fig. 2. The difference in 10 Be nuclide production between a modern surface topography and a surface corrected for GIA at 20 ka, 15 ka, 10 ka and 5 ka. Only the land surface above present-day sea level is shown. Red areas highlight greater nuclide production due to higher surface elevations relative to today, while blue areas reveal lower production values due to isostatically depressed surface elevations. Dotted and solid lines link areas of ± 5% and 10% production rate difference, respectively. The largest effects are in the vicinities of the Laurentide, Cordilleran, Fennoscandian, Greenland and West Antarctic ice sheets. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Supplementary Data). The largest offsets between corrected and uncorrected ages are in Scandinavia and eastern Canada, where the correction of existing samples makes exposure ages up to 6,300 years and 37% older. For Scandinavia, if we exclude samples that are considered unreliable estimates of deglaciation (Hughes et al., 2016), then the maximum exposure age correction is 2,600 years and 15%. The effect is relatively smaller away from the isostatically depressed centres of the former Fennoscandian and Laurentide ice sheets, as well as for smaller ice masses. In Greenland, Patagonia and Antarctica, exposure ages become up to~7%,~7.5% and~13% older, respectively, while in the vicinity of the British-Irish Ice Sheet and in New Zealand, the correction produces age offsets of <2.5% and <1% respectively (Fig. 4). Beyond the former ice margins, in the region of a peripheral bulge, the GIA correction has an opposite effect, making exposure ages up to~5% younger than the uncorrected ages (Fig. 4).
Whether corrected exposure ages become older or younger, and the magnitude of the age correction, depends on the time that the sample became exposed relative to the GIA response. Even in areas of large isostatic depression, samples that become exposed after the main period of elevation change will have a minor or no age correction (small circles overlying large circles in Fig. 4). In some cases, an area that was once isostatically depressed but then became part of a peripheral bulge during ice sheet retreat can produce exposure ages that are both older and younger relative to ages that are uncorrected for GIA (e.g. west Greenland and northeastern USA; Fig. 4).

Discussion
Our assessment has highlighted the potential magnitude and spatial pattern of the impact of GIA on surface-exposure dating, and how these effects have varied through time. The greatest effects occur near to the centre of former large ice masses, however, most of these areas likely became exposed towards the end of deglaciation, after most isostatic rebound had taken place. Ice-free areas in the interior of ice sheets (i.e. nunataks) that became exposed earlier during deglaciation, such as in Antarctica, would have experienced a longer duration of isostatic rebound and, therefore, potentially larger exposure age corrections. The effects of GIA can still be significant beyond the ice sheet margins, even many hundreds of kilometres away from the major ice masses. Crucially, exposure ages in these areas can become older, younger or be unchanged when corrected for GIA, with the form of the effect dependent on the timing of exposure and the local GIA response. Based on our assessment, exposure ages that do not account for GIA are likely incorrect by >5% in some regions, and for some sites this error would exceed the total sample uncertainty (~8.7 ± 4.1%, based on the global dataset used in this study).
Testing whether the GIA-corrected rather than uncorrected exposure ages are better estimates of the 'true' age requires comparison to independent chronologies. We looked at two example areas in Scandinavia and North America, which are regions where the GIA correction is greatest. The oldest onshore postglacial sediments in Norway are from northernmost Andøya, where basal radiocarbon ages from lake sequences have produced minimumlimiting ages for deglaciation of 26.0e22.0 cal ka BP (Hughes et al., 2016;Vorren et al., 2013). Cosmogenic 10 Be surfaceexposure ages from locations proximal to the lake sites (23.4 ± 2.5 to 18.7 ± 2.0 ka) are generally younger than the radiocarbon ages (Nesje et al., 2007;Stroeven et al., 2016) and, therefore, seem irreconcilable. Correcting these exposure ages for GIA increases the apparent exposure ages by~600e1100 years, however, for most of these samples, this correction is still not enough to be consistent with the radiocarbon chronology. An offset between radiocarbon and 10 Be surface-exposure ages is also well demonstrated in northwestern Ontario, Canada. Here radiocarbon ages from basal lake sediments range from~11.5 to 12.5 cal ka BP (Bajc et al., 2000;Teller et al., 2005), but the uncorrected exposure ages date deglaciation tõ 13.6 ± 0.1 ka (calculated here for samples FR-1-15, FR-6-15 and FR-11-15;Leydet et al., 2018). Our GIA correction increases this offset by another~1000 years to give a surface-exposure date of~14.6 ka, however, the offset could also be partly attributed to the minimumlimiting nature of deglacial radiocarbon ages. Such techniquespecific and site-specific geochronological uncertainties will often preclude robust data comparisons and, instead, a regional comparison may provide a more suitable evaluation of the GIA correction.
In Antarctica there is a longstanding puzzle about the timing of major deglaciation. Far-field estimates of sea-level change and some marine geological data indicate that a significant volume of Meltwater Pulse 1 A (MWP-1A) at~14.5 cal ka BP may have been sourced from the Antarctic ice sheets (Clark et al., 2002;Liu et al., 2016). In contrast, surface-exposure chronologies that record interior and coastal ice sheet thinning imply that the majority of ice volume change in those regions occurred in the Early-mid Holocene (e.g. Bentley et al., 2010;Hein et al., 2016;Jones et al., 2015;Small et al., 2019;Spector et al., 2017). When GIA corrections are applied (Fig. 4), the age difference between the chronologies of ice sheet thinning and the timing of MWP-1A can be reduced in some places by~1e2 ka, approximately half of the total discrepancy.
These comparisons highlight the fact that additional factors may need to be considered to fully resolve differences between chronologies. In particular, the accuracy of the exposure age correction depends on our ability to quantify Earth rheology and past ice sheet change. Using Antarctica as an example, we compare GIA corrections that have been computed using different Earth and ice models (Argus et al., 2014;Peltier et al., 2015;Whitehouse et al., 2012b) and find that the choice of Earth rheology makes a relatively minor impact (<1%), but the choice of ice history model could potentially account for up to~15% difference in the nuclide production rate (see Supplementary Material).
The impact of GIA on surface-exposure ages may also be   4. Existing surface-exposure ages corrected for GIA in the regions of formerly large ice masses. Red (blue) symbols highlight ages that become older (younger) when accounting for time-dependent GIA effects. The size of the symbol represents the magnitude of the offset (years) between the corrected and uncorrected mean age, while a darker (lighter) shade represents a greater (smaller) percentage difference. The central-north parts of Scandinavia were not deglaciated until the Holocene, but we include older, and likely anomalous, ages in this area to highlight the potential age correction if these sites were exposed (e.g. as nunataks) early in the deglaciation. The data shown here can be found in the Supplementary Data. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) amplified or reduced with the inclusion of time-dependent atmospheric changes. Atmospheric compression due to cooling and a decrease in atmospheric pressure from katabatic winds would have increased nuclide production by up to 10% at the LGM (Staiger et al., 2007). This would make exposure ages at some ice-marginal locations relatively younger, potentially reversing the GIA effect. However, changes in the pressure distribution were highly spatially variable, with the displacement of atmospheric mass in some regions causing production rates to decrease by~5e10% at the LGM, such as for sections of the Fennoscandian, Laurentide and Cordilleran ice sheet margins (Staiger et al., 2007). This reduction in nuclide production would make exposure ages older, potentially adding to the effects from an isostatically depressed land surface at that time. Although Staiger et al. (2007) do not provide an evaluation of the atmospheric effects through deglaciation, it is apparent that changes both in the atmosphere and from GIA significantly influenced past nuclide production. Future efforts should focus on developing time-dependent production scaling schemes that account for these spatially and temporally varying effects in order to provide more accurate surface-exposure age estimates. Finally, surface-exposure ages provide crucial constraints for ice sheet models (e.g. Briggs and Tarasov, 2013;Tarasov et al., 2012;Whitehouse et al., 2012a). However, there is some circularity in correcting exposure ages using ice models that are partially constrained by these ages. The circularity feeds into GIA corrections, which require an ice-load history, as well as atmospheric corrections, which need to use climate simulations that incorporate timedependent changes in ice sheet volume. Resolving this issue will be a challenge for the surface-exposure dating and ice-sheet modelling communities.