Next Article in Journal
Age of Initial Submarine Volcanism in the Paleo-Tsushima Basin and Implications for Submarine Volcanism in the Opening Stage of the Japan Sea in Northern Kyushu
Next Article in Special Issue
Cosmogenic Exposure Dating (36Cl) of Landforms on Jan Mayen, North Atlantic, and the Effects of Bedrock Formation Age Assumptions on 36Cl Ages
Previous Article in Journal
High Preservation Potential Volcaniclastic Sedimentation in the Serravallian Sequence of the Amantea Basin (Coastal Chain, North-Western Calabria)
Previous Article in Special Issue
Glacial Erosion Rates Determined at Vorab Glacier: Implications for the Evolution of Limestone Plateaus
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The NUNAtak Ice Thinning (NUNAIT) Calculator for Cosmonuclide Elevation Profiles

Scottish Universities Environmental Research Centre, University of Glasgow, East Kilbride G75 0QF, UK
Geosciences 2021, 11(9), 362; https://doi.org/10.3390/geosciences11090362
Submission received: 15 July 2021 / Revised: 20 August 2021 / Accepted: 23 August 2021 / Published: 26 August 2021
(This article belongs to the Special Issue Cutting Edge Earth Sciences: Three Decades of Cosmogenic Nuclides)

Abstract

:
Cosmogenic nuclides are widely used to constrain the landscape history of glaciated areas. At nunataks in continental polar regions with extremely arid conditions, cosmogenic nuclides are often the only method available to date the ice thinning history of the glacier. However, the amount of cosmogenic isotopes accumulated at the surface of nunataks depends not only on the length of time that rock has been exposed since the last deglaciation but also on the full history of the surface, including muon production under ice, exposure during previous interglacials, subaerial weathering rate, glacial erosion rate, and uplift rate of the nunatak. The NUNAtak Ice Thinning model (NUNAIT) simulates the cosmonuclide accumulation on vertical profiles, fitting the aforementioned parameters to a set of multi-isotope apparent ages from samples taken at different elevations over the ice-sheet surface. The NUNAIT calculator is an easy-to-use tool that constrains parameters that describe the geological history of a nunatak from a set of surface exposure ages.

1. Introduction

Quantifying the changes in the thickness of the Greenland and Antarctic ice sheets is key to understanding future sea-level rise [1]. Cosmogenic nuclides are widely used for the quantification of glacial chronologies. However, the climatic interpretation of the existing cosmonuclide data sets requires accounting for geologic processes that cause apparent exposure ages on glacial landforms to differ from the age of deglaciation [2].
Nunataks, the mountains emerging from polar ice sheets, have been used as vertical dipsticks that record past changes in the thickness of the polar ice sheets (e.g., [3]). Cosmogenic signatures at the surface of nunataks are the result of the intermittent exposure of the surfaces to cosmic radiation through the glacial cycles (e.g., [4]), glacial erosion (e.g., [5,6]), and the subaerial weathering of these surfaces (e.g., [7]). Therefore, the abundance of one or more cosmonuclides in one of these surfaces can be explained by the combination of multiple possible scenarios [8].
Stroeven et al. [4] modelled the accumulation of 10 Be and 26 Al in tors. The model they used is based on complex exposure-burial histories forced along the ice-free/ice-covered conditions provided by a marine oxygen isotope δ 18 O proxy glacial record. In this model, a given δ 18 O cutoff value defines when the surface of the tor is exposed or shielded from cosmic radiation. Li et al. [9] developed a method to solve the cutoff value in the marine oxygen isotope record that satisfies a set of 10 Be and 26 Al concentrations considering fixed values of glacial erosion and subaerial weathering. Knudsen et al. [10] described a method to solve not only the δ 18 O cutoff value but also the glacial and interglacial erosion rates from a set of multiple cosmonuclide concentrations that can include 10 Be, 26 Al, 14 C, and/or 21 Ne data.
The models described by Stroeven et al. [4], Li et al. [9], and Knudsen et al. [10] are designed to be applied on a single site, and therefore one cutoff δ 18 O value can be solved at a time. To solve the elevation of the ice surface during the glaciations, several samples should be used to obtain an elevation profile of δ 18 O cutoff values, which would allow the reconstruction of the ice sheet thickness with time. If the elevation of the ice surface is known, a second iteration of modelling would allow calculating how cosmonuclides accumulate during glacial times by reconstructing the muonic production cross-section under the ice sheet for any time.
In summary, the interpretation of cosmonuclide concentrations from nunataks often requires accounting for the effects of surface and subglacial erosion, glacial dynamics, and tectonic activity. The models described in the literature are focused on solving one or two of the parameters that emulate these processes, usually using data from a single sample. Here, I describe an easy-to-use method to solve up to five parameters that emulate the glacial history, surface erosion, and tectonic uplift using a set of surface exposure ages.
The NUNAtak Ice Thinning (NUNAIT) calculator presented here solves (1) the elevation history of the ice surface, (2) the glacial erosion rate, (3) the subaerial weathering rate, and (4) the nunatak uplift rate from a multi-sample (elevation profile) and multi-isotope ( 10 Be, 26 Al, 21 Ne, 3 He, 36 Cl, and/or 14 C) data set. The calculator does not require the input of production rates, as the default inputs are not cosmogenic concentrations but apparent surface exposure ages, and approximate muon cross sections are calculated using the latitude and elevation of the sampling sites.

2. Method Details

Here, I present a set of MATLAB®/GNU Octave© scripts that form the NUNAIT calculator and their mathematical descriptions. All scripts needed to run the NUNAIT calculator (Supplementary Materials) are freely accessible at https://github.com/angelrodes/NUNAIT (accessed on 19 August 2021).
When running the script START.m, the user is asked to run the calculator or select previous data to display the text and output.
If the first option is selected (Run simulation), two types of files can be selected:
  • A .csv file containing basic input data;
  • A .mat file containing full input data, including apparent concentrations and apparent production rates. A *_sampledata.mat is generated every time a .csv is processed. This allows, for example, changing the distribution of the production rates before running the simulations by editing the *_sampledata.mat file.
If the second option is selected (Display results), a .mat file containing previously calculated data is required. This type of file is generated at the end of each fitting session with the same name as the input file and _model.mat.

2.1. Input Data

Site data have to be inputted in individual comma separated files (.csv) for each measurement. Some examples of input files are included in the folder “Examples”. The input file contains the following headers (first line) that we recommend are not changed:
  • name: Sample name without spaces or symbols.
  • lat: Latitude used to calculate the muon contributions (decimal degrees).
  • site_elv: Elevation of the sample above sea level (m).
  • isotope: Mass of the cosmogenic isotope. Currently accepting 3, 10, 14, 21, 26, and 36 for 3 He, 10 Be, 14 C, 21 Ne, 26 Al, and  36 Cl, respectively.
  • base_level: Current elevation of the glacier surface above sea level at the sampling site (m). This is used to calculate the ice position through time.
  • apparent_years: Apparent surface exposure age calculated with any cosmogenic calculator, any scaling scheme, and any production rate reference.
  • dapparent_years: External uncertainty of the previous age.
Apparent concentrations (C) are calculated from apparent surface exposure (T) ages following Lal [11]:
C = 1 λ · 1 e λ T
where λ is the decay constant of the isotope considered. The values of λ are stored in constants.m.
Note that the concentrations described in Equation (1) are scaled to site production rates. Therefore, they should be expressed in time units (years).
To reduce computing time, conditional statements are avoided in the code by considering all cosmonuclides radioactive. To do this, stable isotopes are assigned values of λ corresponding to 100 times the age of the Earth. As  T < < 1 / λ for stable isotopes, Equation (1) results in C T .
The calculated concentrations, together with the muon relative contributions described in Section 2.3 and Section 2.4, are stored in a .mat file with the same name as the original .csv file and the suffix _sampledata. If the user needs to change the calculated concentrations or relative production rates, this file can be modified and used as an input file.

2.2. Climate Curves

The scripts make_climatecurves.m and make_climatecurves_ant.m generate a time series of δ 18 O values that will be used to calculate the vertical position of the glacial surface over the samples.
The curves from Lisiecki and Raymo [12] and Zachos [13] are combined and scaled with NGRIP data [14] or Five-core data [15] in Antarctica.
All records are arbitrarily scaled to the LR04 stack data [12]. As the δ 18 O values generated will finally be transformed into elevations by the model described in Section 2.6, the choice of one data set as reference is irrelevant.
To reduce the number of calculations and the computing time while representing the ice changes relevant to the cosmogenic accumulation, the data are interpolated for ages every 10 years for the last century, every 100 years until 20 ka, every 200 years until 50 ka, every 500 years until 100 ka, and every 1% increase for ages older than 100 ka. The resulting simplified curve is shown in Figure 1.

2.3. Muon Contributions

The function muon_contribution.m generates the muon contribution and its uncertainty based on latitude ( l a t ) and elevation ( e l v ) for a given nuclide. If either latitude or elevation is not a number, a global average is given. A single value of latitude and elevation is used to calculate the contribution of muons to the total surface production of 10 Be. All other productions are scaled accordingly.
The contribution of muons to the total surface 10 Be production ( R μ ( 10 Be ) ) is calculated as
P μ ( 10 Be ) P t o t a l ( 10 Be ) = 1 100 · 1.29 + l a t 900 + 1.056 · e l a t + 1 30.31 2 · 0.1 + 0.9 · e e l v 2000
This approximation is based on the 10 Be production at 1678 sites equally distributed on land areas according to ETOPO1_Bed_g_geotiff.tif [16] and calculated using P_mu_total_alpha1.m and stone2000.m from Balco [17] and Balco et al. [18], respectively. The fitting of this approximation is shown in Figure 2. This formula fits the original data within a 5% standard deviation.
Considering that P_mu_total_alpha1.m fits the empirical data available within a ∼5% and a ∼13% for the 10 Be and 26 Al muon production rates, respectively [17], the uncertainty of the calculated muon contributions based on Equation (2) should be at least a 7% for 10 Be and 14% for 26 Al.
The calculation of the muon contributions for other nuclides are based on the following ratios:
  • R μ ( 26 Al ) / R μ ( 10 Be ) = 1.4587 . See Figure 2.
  • R μ ( 36 Cl ) / R μ ( 10 Be ) = 3.2720 , according to Heisinger and Nolte [19].
  • R μ ( 21 Ne ) / R μ ( 10 Be ) = 4.086 , according to Balco and Shuster [20].
  • R μ ( 3 He ) / R μ ( 10 Be ) = 1 , consistent with Blard et al. [21].
  • R μ ( 14 C ) / R μ ( 10 Be ) = 8.2767 , according to Heisinger and Nolte [19].
As the uncertainties of these ratios are unknown, this script assigns a conservative 20% uncertainty for muon contributions calculated using Equation (2) to cover both the uncertainties at the surface and the subsurface extrapolations described in Section 2.4.
All these data can be changed in the files constants.m and muon_contribution.m.

2.4. Muon Cross Sections

To simulate production under ice and rock surfaces, muon production was approximated as three exponential functions of depth [22]. A total of 1678 10 Be and 26 Al muon production rates generated using P_mu_total_alpha1.m[17] were analysed to fit three exponential decays with attenuation lengths of 850, 5000, and 500 g cm 2 (Figure 3). These attenuation lengths correspond to 75% of the fast muon, 25% of the fast muon, and the negative muon productions at the surface, respectively.
The share of surface fast muon production with respect to the total muon production ( P μ f a s t / P μ t o t a l ) considered for each isotope is:
  • 10 Be: P μ f a s t / P μ t o t a l = 0.32069 . See Figure 2.
  • 26 Al: P μ f a s t / P μ t o t a l = 0.22282 . See Figure 2.
  • 36 Cl: P μ f a s t / P μ t o t a l = 0.0620 , according to Heisinger and Nolte [19].
  • 21 Ne: P μ f a s t / P μ t o t a l = 1 , according to Balco and Shuster [20].
  • 3 He: P μ f a s t / P μ t o t a l = 0.32069 , consistent with Blard et al. [21].
  • 14 C: P μ f a s t / P μ t o t a l = 0.0672 , according to Heisinger and Nolte [19].
The uncertainties of these approximations are within the uncertainties described in Section 2.3 for 10 Be and 26 Al.
All these data can be changed in the file constants.m.

2.5. Densities

A density of ρ i c e = 0.917 g cm 3 is considered for ice [23], and a density of ρ = 2.65 cm 3 for bedrock.

2.6. Nunatak Accumulation Model

The nunatak accumulation model (nuna_model.m) considers the depth of the sample under the bedrock surface (z) and the thickness of the ice on top of the surface ( z i c e ) based on the input conditions (weathering w, glacial erosion rate, and maximum and current ice levels) for each time range ( Δ t ) defined by the climate curve and for each sample. The model concentration is calculated as
C i = P λ + w · ρ / Λ · e ( z · ρ + z i c e · ρ i c e ) / Λ · 1 e Δ t · ( λ + w · ρ / Λ ) · e λ · t
where P is the production rate considered (spallation and each of the muon types), Λ is the attenuation length for the production rate considered, λ is the decay constant of the nuclide, and t is the age corresponding to the end of the time range defined by the climate curves.
The final concentration for each sample C m o d e l is calculated by adding all the C i for all production types and time ranges.
The effect of the glacial erosion rate is ignored inside each time range ( Δ t ), as usually C i is much more sensitive to z i c e than to the change in position of the sample under the bedrock surface due to glacial erosion during ice-covered periods.

2.7. Model Fitting

The fitting of the model described in Section 2.6 is performed by the script fit_nuna_model.m.
The script asks the user to set maximum and minimum values for the parameters to be fitted: ice-free weathering rate, glacial erosion rate, ice-thinning since maximum glacier extension, deviation of the current ice surface, and uplift rate. As weathering and erosion rates are simulated in logarithmic space, minimum values of 0.1 mm/Ma are assumed.
It also allows changing the fit type. With a value of 0, the script will try to fit the model to the data normally. With a value of 1, models with concentrations below the sample concentrations will be ignored. With a value of 2, models with concentrations above the sample concentrations will be ignored. A value of 3 is used to represent the models within the stated parameter limits ignoring the sample concentrations. If fit type 3 is used, the script assumes that all generated models fit the data.
The script selects the climate reference based on the average latitude of the samples. The Antarctic curves described in Section 2.2 are used for latitudes south of 55 S.
The degrees of freedom ( ν ) are calculated by subtracting the number of parameters with an initial range greater than 0 from the number of data in the input file (Section 2.1). A minimum ν of 1 is always considered.
The script calculates concentrations corresponding to the sample positions for random parameter values between the parameter limits. Randomisation of the weathering and erosion rate values is performed logarithmically. A combination of random parameter values is computed in each iteration. The goodness of fit is defined by the chi-squared function:
χ 2 = i = 1 n C m o d e l C i σ C m o d e l 2 + σ C i 2 2
where C i and σ C i are the sample concentrations and their uncertainties derived from the apparent surface exposure ages (Section 2.1), C m o d e l are the model concentrations corresponding to sample i (Section 2.6), and σ C m o d e l is the model uncertainty corresponding to the uncertainty of the muon produced concentration (Section 2.3) plus the minimum age spacing of the climate data (10 years, as described in Section 2.2).
Models fitting the data within a 1 σ confidence level are defined by the ones with χ 2 χ m i n . 2 + ν , and models fitting the data within a 2 σ confidence level are defined by the ones with χ 2 χ m i n . 2 + 2 · ν . Note that these formulas do not fully represent the chi-squared distribution described in Rodés et al. [24] (Section 2.2.1). The method described in Rodés et al. [24] often yields infinite values when computing maximum fitting values ( χ m a x . 2 ) for poor fittings and high ν . The formula χ m a x . 2 = χ m i n . 2 + n · ν is an approximation to the method described by Avni [25] for high degrees of freedom.
After a learning cycle of 3000 iterations (consts.minmodelstoconverge in constants.m), the limits of the randomised parameters start converging. Initially, the new limits converge to the models that fit the data within a 2 σ confidence level, and within a 1 σ confidence level for the last 1/3 of the total iterations. If the number of models fitting the data within this confidence level is lower than n c o n v , this confidence range is increased to 4 σ , 8 σ , 16 σ , etc. n c o n v is initially equal to the desired fitting models and decreases exponentially with time. The new parameter limits are calculated every 100 iterations (consts.convergencestep in constants.m) from the models fitting the desired confidence level and expanded by 10% of the range to avoid missing fitting values at the limits of the 1 σ range.
The script runs iterations until one of the following conditions are met:
  • The simulations reach the maximum number of models to calculate:
    consts.maxnmodels = 50,000 in constants.m.
  • There are more than the desired fitting models that fit the data within the 1 σ confidence level:
    consts.targetnmodelsonesigma = 300 in constants.m.
To represent the results as probability density distributions of the parameters, the relative probability corresponding to each model is calculated as
P ( χ 2 ) ν χ 2 · e χ 2 / ( 2 · ν )
which has a similar shape as the cumulative distribution function of the chi-squared distribution but can be computed avoiding zeros for high values of χ 2 and ν .
Finally, a set of fake samples covering a wide range of altitudes and all the fitted nuclides is generated. The parameter values of the models fitting the original data within a 1 σ confidence level are used to generate altitudinal concentration profiles within the fake sample’s data. Maximum and minimum concentration profiles are generated and used to plot the scatter of the fitting models.

2.8. Data Representation

A summary of the results is outputted in the command window by the script display_ results.m.
Three figures are generated by plot_results.m as graphical output:
  • The probability distribution of the models for each of the parameters.
  • A representation of the ice surface evolution and the altitudinal trajectories of the samples (if no uplift rate is considered, these will be horizontal lines). Uncertainties corresponding to all 1 σ models are also represented.
  • Altitudinal profiles of the apparent exposure ages for all 1 σ models and all nuclides, and the actual apparent exposure ages of the samples (model vs. data).
An example of the full graphical output generated by the NUNAIT calculator is shown in Section 3.2.

3. Examples

Two natural examples of input files are included in the folder “Examples”. Inputs can be generated from new or published data. Data from published data can be easily generated from the ICE-D: ANTARCTICA database [26] (http://antarctica.ice-d.org, accessed on 19 August 2021), that compiles a large number of cosmogenic data sets, including updated exposure ages, organized by sites. The only datum required by the NUNAIT calculator that is missing in the ICE-D database is the current elevation of the ice surface. This value could be easily guessed by checking the lowest sampling site in the set, which often coincides with the ice surface.

3.1. Marble Hills

Marrero et al. [27] reported the first cosmogenic nuclide-derived erosion rates for carbonate rocks in Antarctica. Erosion rates were derived from carbonate bedrock samples at the Marble Hills field site (Ellsworth Mountains). I generated the exposure ages required by the NUNAIT calculator using the CRONUS online calculators [28] with the data included in Marrero et al. [27] (Appendix A).
In the original paper, Marrero et al. [27] calculated an apparent 36 Cl erosion rate of 0.22 ± 0.02 mm/ka from the samples above the elevation of ∼550 m above the present ice-surface, and apparent 36 Cl erosion rates >4 mm/ka at lower elevations. Marrero et al. [27] interpreted that samples below 550 m over the present ice surface require complex exposure-burial histories to explain their composition.
As shown in Figure 4, the NUNAIT calculator predicts a subaerial weathering rate between 0.52 and 0.84 , 2 to 4 times higher than the one calculated by Marrero et al. [27]; a glacial erosion rate below 37 mm/ka, concordant with the cold-based glacial processes expected in Marble Hills; and a maximum ice extension of 176–232 m above the present ice surface. These results suggest that the data from samples above 232 m could be compatible with simple exposure conditions, and the scatter of these data could be explained by inhomogeneous weathering ratios of the continuously exposed surfaces.
The fit type 1 was used to constrain the minimum weathering rate that is compatible with these data (Figure 5), implicitly assuming that faster apparent weathering rates due to different lithologies, slopes, etc., can produce shorter apparent exposure ages in some surfaces. A (minimum) subaerial weathering rate between 0.19 and 0.21 m/Ma was obtained using this setup. These values agree with the value of 0.22 ± 0.02 mm/ka obtained by Marrero et al. [27].

3.2. Mount Hope

The longest nunatak elevation profile showing a wide range of cosmogenic isotope exposure ages in the ICE-D ANTARCTICA database is probably the site HOPE (Mt. Hope, Beardmore Glacier, Southern Ross Sea). Samples from the HOPE site appear in Spector et al. [29] and Spector [30]. The data set contains 10 Be, 26 Al, 21 Ne, 3 He, and  14 C exposure ages.
According to Spector et al. [29], Mt. Hope (836 m) remained ice-covered until 14.4 ± 0.5 ka. Several kilometres upstream from this position, two lateral moraines at 1050 and 1200 m mark the maximum elevation of ice during the Last Glacial Maximum.
Using only the cosmogenic exposure ages from the bedrock at Mt. Hope, the NUNAIT calculator yields a maximum elevation of ice during the Last Glacial Maximum of ∼1065 m above sea level (Figure 6), which seems to be in good agreement with the position of the lateral moraines described by Spector et al. [29].
Figure 6 shows that the optimum fit of the NUNAIT model mimics the distribution of the data for each isotope but does not fit the ratios between isotope data well.
As for the Marble Hills’ data, the NUNAIT model fits this data set for very low weathering and glacial erosion rates and predicts a maximum uplift rate of 15 m/Ma.

4. Discussion and Conclusions

The NUNAIT and previous models described by Stroeven et al. [4], Li et al. [9], and Knudsen et al. [10] are based on the same principle: using a climate proxy ( δ 18 O record) to solve complex exposure-burial histories that fit the surface cosmogenic nuclide data. Previous models focused on solving the problem for sets of multiple isotope data from single samples. The method presented here focuses on solving the same problem but considering data from all the sampled sites on the nunatak, and yielding results that are consistent with the whole data set. Thus, the results obtained using the NUNAIT calculator are expected to be less precise than the results based on single samples but more robust, as the model can consider more possible scenarios by randomizing erosion rates, the position of the current ice surface, or uplift.
The NUNAIT calculator requires an input of cosmogenic data as apparent exposure ages with no surface erosion rate considered. The models fitting the data are also expressed as apparent exposure ages in the output. The use of input/output in this intuitive format has some advantages and disadvantages.
The user needs to calculate apparent exposure ages using a local or online exposure age calculator (e.g., [18,28,31]), allowing the user to consider any calculator, production rate reference, and scaling factor. This simplifies the use of the NUNAIT calculator, as the information about the production of cosmogenic isotopes (production rate, shielding, self-shielding factors, radiogenic produced concentrations, etc.) is implicitly included in the input data.
Although the user does not need to deal with production rates, the NUNAIT model works internally with scaled concentrations and constant production rates. Equation (1) assumes that the average production rate for the apparent (minimum) age equates to the constant production rate. This introduces differences with the time-dependent production models typically considered for the calculation of surface production rates. According to Balco et al. [18] (Figure 3 and 4), these differences should not exceed 10% of apparent exposure ages for most altitudes in polar regions. However, this uncertainty should be represented in the input data by the external uncertainty of the apparent exposure ages.
As the model is fitted using external uncertainty of the apparent exposure ages (Equation (4)), the fittings provided by the NUNAIT calculator are more sensitive to the spatial distribution of cosmogenic concentrations than to the ratios between different isotopes (e.g., Section 3.2), in contrast with the methods described by Stroeven et al. [4], Li et al. [9], and Knudsen et al. [10]. This effect is intentional and seeks to reflect the uncertainties of the cosmogenic surface production rate ratios realistically (e.g., [32,33,34]).
As the default input does not include any information on the muon contributions, these values need to be estimated as shown in Section 2.3. This approximation introduces an uncertainty of a similar magnitude as the one derived from the scaling scheme. The simplification of the muon cross-sections described in Section 2.4 introduces an additional uncertainty of 5% in the muon production rates under the ice sheet. The uncertainty of the muon produced cosmonuclides should also include the scatter of the global data available for the calibration of muon production under the surface, which is ∼5% and ∼14% for 10 Be and 26 Al, respectively, according to Balco [17] (Table 1). The NUNAIT calculator incorporates these uncertainties by considering a 20% uncertainty for all muon-produced concentrations.
According to the data summarized in Balco [17], the best predictions of the muon-produced 26 Al/ 10 Be ratios fit the empirical data within ∼20% uncertainty (Figure 7). For other isotope pairs, the existing empirical data about their production rate ratios at depth are more scarce (e.g., [35]). Therefore, we should assign an uncertainty greater than 20% to our modelled concentration ratios at great depths. As the NUNAIT model considers a 20% uncertainty for all muon-produced concentrations, it assumes a 28% uncertainty ( 20 2 + 20 2 ) for any synthetic concentration ratio under the ice-sheet, which is probably an overestimation of the 26 Al/ 10 Be uncertainties. However, similar to the surface predictions, this overestimation of the uncertainty makes the model less sensitive to the ratios between isotopes and therefore relatively more sensitive to the spatial distribution of the data.
The NUNAIT model considers a constant ice density for the ice column covering the samples during glaciations. This value can be adjusted by the user, and the effect of its uncertainty is not expected to exceed the 20% uncertainty considered for muon-produced concentrations.
When uplift is considered in the model, it is assumed to be a constant rate. Isostatic rebound is not emulated by the NUNAIT model. This should not greatly affect the distribution of glaciated and ice-free elevations through time, as the isostatic rebounds are expected to be coupled with the changes in the elevations of the ice surface. However, a constant fast uplift could result in surfaces accumulating cosmogenic nuclides at slower rates in the past due to the reduced production at lower elevations. This effect is not yet considered by the NUNAIT model. Therefore, this model could overestimate the concentration of stable isotopes in highly uplifted areas that have not been glaciated in the past.
During ice-free periods, the NUNAIT model considers a homogeneous weathering rate along with the elevation profile. This might not be very realistic for areas with intense periglacial processes that produce increased erosion rates in local areas (e.g., rock falls). When fitting the model to data from areas with evident periglacial processes, the minimum fitting type should be selected (fit type 1 described in Section 2.7).
Marrero et al. [27] described a systematic difference between bedrock and boulder samples, with boulder samples yielding systematically lower erosion rates. By default, the NUNAIT model considers a homogeneous erosion rate under the ice. Therefore, when fitting data from erratic boulders that could have been preserved during glacial periods, and hence maintaining a higher surface cosmonuclide concentration than the bedrock, the maximum fitting type should be selected (fit type 2 described in Section 2.7).
The examples in Section 3 show that the NUNAIT calculator yields results that go beyond the typical observations deduced from surface exposure dating, such as glacial erosion rates and uplift rates. Therefore, the NUNAIT calculator is presented as an easy-to-use tool that will help glaciologists to interpret cosmogenic data from nunataks, where exposure histories are usually complex. Moreover, the methods described in Section 2 can be used to develop new cosmogenic-based tools with intuitive and simplified inputs and outputs.

Supplementary Materials

All scripts discussed in Section 2 and the data discussed in Section 3 are freely accessible at https://github.com/angelrodes/NUNAIT (accessed on 19 August 2021, subject to the terms of the GNU General Public License, version 3, as published by the FreeSoftware Foundation).

Funding

This research received no external funding.

Data Availability Statement

All scripts and example data are freely accessible at https://github.com/angelrodes/NUNAIT (accessed on 19 August 2021).

Acknowledgments

I would like to thank Laura Medina-Ruiz for proofreading the manuscript and the three anonymous reviewers for their constructive insights during the review of the article.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Alley, R.B.; Clark, P.U.; Huybrechts, P.; Joughin, I. Ice-Sheet and Sea-Level Changes. Science 2005, 310, 456–460. [Google Scholar] [CrossRef] [Green Version]
  2. Balco, G. Contributions and unrealized potential contributions of cosmogenic-nuclide exposure dating to glacier chronology, 1990–2010. Quat. Sci. Rev. 2011, 30, 3–27. [Google Scholar] [CrossRef]
  3. Stone, J.O.; Balco, G.A.; Sugden, D.E.; Caffee, M.W.; Sass, L.C.; Cowdery, S.G.; Siddoway, C. Holocene Deglaciation of Marie Byrd Land, West Antarctica. Science 2003, 299, 99–102. [Google Scholar] [CrossRef] [Green Version]
  4. Stroeven, A.P.; Fabel, D.; Hättestrand, C.; Harbor, J. A relict landscape in the centre of Fennoscandian glaciation: Cosmogenic radionuclide evidence of tors preserved through multiple glacial cycles. Geomorphology 2002, 44, 145–154. [Google Scholar] [CrossRef]
  5. Davis, P.T.; Bierman, P.R.; Marsella, K.A.; Caffee, M.W.; Southon, J.R. Cosmogenic analysis of glacial terrains in the eastern Canadian Arctic: A test for inherited nuclides and the effectiveness of glacial erosion. Ann. Glaciol. 1999, 28, 181–188. [Google Scholar] [CrossRef] [Green Version]
  6. Stroeven, A.P.; Fabel, D.; Harbor, J.; Hättestrand, C.; Kleman, J. Quantifying the erosional impact of the fennoscandian ice sheet in the torneträsk–narvik corridor, northern sweden, based on cosmogenic radionuclide data. Geogr. Ann. Ser. A Phys. Geogr. 2002, 84, 275–287. [Google Scholar] [CrossRef]
  7. Nishiizumi, K.; Kohl, C.; Arnold, J.; Klein, J.; Fink, D.; Middleton, R. Cosmic ray produced 10Be and 26Al in Antarctic rocks: Exposure and erosion history. Earth Planet. Sci. Lett. 1991, 104, 440–454. [Google Scholar] [CrossRef]
  8. Bierman, P.R.; Marsella, K.A.; Patterson, C.; Davis, P.; Caffee, M. Mid-Pleistocene cosmogenic minimum-age limits for pre-Wisconsinan glacial surfaces in southwestern Minnesota and southern Baffin Island: A multiple nuclide approach. Geomorphology 1999, 27, 25–39. [Google Scholar] [CrossRef]
  9. Li, Y.; Fabel, D.; Stroeven, A.P.; Harbor, J. Unraveling complex exposure-burial histories of bedrock surfaces under ice sheets by integrating cosmogenic nuclide concentrations with climate proxy records. Geomorphology 2008, 99, 139–149. [Google Scholar] [CrossRef]
  10. Knudsen, M.F.; Egholm, D.L.; Jacobsen, B.H.; Larsen, N.K.; Jansen, J.D.; Andersen, J.L.; Linge, H.C. A multi-nuclide approach to constrain landscape evolution and past erosion rates in previously glaciated terrains. Quat. Geochronol. 2015, 30, 100–113. [Google Scholar] [CrossRef] [Green Version]
  11. Lal, D. Cosmic ray labeling of erosion surfaces: In situ nuclide production rates and erosion models. Earth Planet. Sci. Lett. 1991, 104, 424–439. [Google Scholar] [CrossRef]
  12. Lisiecki, L.E.; Raymo, M.E. A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records. Paleoceanography 2005, 20, PA1003. [Google Scholar] [CrossRef] [Green Version]
  13. Zachos, J. Trends, Rhythms, and Aberrations in Global Climate 65 Ma to Present. Science 2001, 292, 686–693. [Google Scholar] [CrossRef]
  14. North Greenland Ice Core Project Members. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. Nature 2004, 431, 147–151. [Google Scholar] [CrossRef]
  15. Buizert, C.; Sigl, M.; Severi, M.; Markle, B.R.; Wettstein, J.J.; McConnell, J.R.; Pedro, J.B.; Sodemann, H.; Goto-Azuma, K.; Kawamura, K. Abrupt ice-age shifts in southern westerly winds and Antarctic climate forced from the north. Nature 2018, 563, 681–685. [Google Scholar] [CrossRef]
  16. Amante, C. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis; National Geophysical Data Center: Washington, DC, USA, 2009. [Google Scholar] [CrossRef]
  17. Balco, G. Production rate calculations for cosmic-ray-muon-produced 10Be and 26Al benchmarked against geological calibration data. Quat. Geochronol. 2017, 39, 150–173. [Google Scholar] [CrossRef]
  18. Balco, G.; Stone, J.O.; Lifton, N.A.; Dunai, T.J. A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements. Quat. Geochronol. 2008, 3, 174–195. [Google Scholar] [CrossRef]
  19. Heisinger, B.; Nolte, E. Cosmogenic in situ production of radionuclides: Exposure ages and erosion rates. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. At. 2000, 172, 790–795. [Google Scholar] [CrossRef]
  20. Balco, G.; Shuster, D.L. Production rate of cosmogenic 21Ne in quartz estimated from 10Be, 26Al, and 21Ne concentrations in slowly eroding Antarctic bedrock surfaces. Earth Planet. Sci. Lett. 2009, 281, 48–58. [Google Scholar] [CrossRef]
  21. Blard, P.H.; Braucher, R.; Lavé, J.; Bourlès, D. Cosmogenic 10Be production rate calibrated against 3He in the high Tropical Andes (3800–4900 m, 20–22° S). Earth Planet. Sci. Lett. 2013, 382, 140–149. [Google Scholar] [CrossRef]
  22. Granger, D.; Smith, A. Dating buried sediments using radioactive decay and muogenic production of 26Al and 10Be. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. At. 2000, 172, 822–826. [Google Scholar] [CrossRef]
  23. Harvey, A. Properties of Ice and Supercooled Water. In CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, USA, 2019. [Google Scholar]
  24. Rodés, A.; Pallàs, R.; Braucher, R.; Moreno, X.; Masana, E.; Bourlés, D.L. Effect of density uncertainties in cosmogenic 10Be depth-profiles: Dating a cemented Pleistocene alluvial fan (Carboneras Fault, SE Iberia). Quat. Geochronol. 2011, 6, 186–194. [Google Scholar] [CrossRef]
  25. Avni, Y. Energy spectra of X-ray clusters of galaxies. Astrophys. J. 1976, 210, 642–646. [Google Scholar] [CrossRef]
  26. Balco, G. Technical note: A prototype transparent-middle-layer data management and analysis infrastructure for cosmogenic-nuclide exposure dating. Geochronology 2020, 2, 169–175. [Google Scholar] [CrossRef]
  27. Marrero, S.M.; Hein, A.S.; Naylor, M.; Attal, M.; Shanks, R.; Winter, K.; Woodward, J.; Dunning, S.; Westoby, M.; Sugden, D. Controls on subaerial erosion rates in Antarctica. Earth Planet. Sci. Lett. 2018, 501, 56–66. [Google Scholar] [CrossRef]
  28. Marrero, S.M.; Phillips, F.M.; Borchers, B.; Lifton, N.; Aumer, R.; Balco, G. Cosmogenic nuclide systematics and the CRONUScalc program. Quat. Geochronol. 2016, 31, 160–187. [Google Scholar] [CrossRef] [Green Version]
  29. Spector, P.; Stone, J.; Cowdery, S.G.; Hall, B.; Conway, H.; Bromley, G. Rapid early-Holocene deglaciation in the Ross Sea, Antarctica. Geophys. Res. Lett. 2017, 44, 7817–7825. [Google Scholar] [CrossRef]
  30. Spector, P. Antarctic Glacial History Inferred from Cosmogenic-Nuclide Measurements in Rocks. Unpublished. Ph.D. Dissertation, University of Washington, Washington, DC, USA, 2018. [Google Scholar]
  31. Martin, L.; Blard, P.H.; Balco, G.; Lavé, J.; Delunel, R.; Lifton, N.; Laurent, V. The CREp program and the ICE-D production rate calibration database: A fully parameterizable and updated online tool to compute cosmic-ray exposure ages. Quat. Geochronol. 2017, 38, 25–49. [Google Scholar] [CrossRef]
  32. Goethals, M.; Hetzel, R.; Niedermann, S.; Wittmann, H.; Fenton, C.; Kubik, P.; Christl, M.; von Blanckenburg, F. An improved experimental determination of cosmogenic 10Be/21Ne and 26Al/21Ne production ratios in quartz. Earth Planet. Sci. Lett. 2009, 284, 187–198. [Google Scholar] [CrossRef] [Green Version]
  33. Balco, G.; Shuster, D.L. 26Al-10Be-21Ne burial dating. Earth Planet. Sci. Lett. 2009, 286, 570–575. [Google Scholar] [CrossRef]
  34. Phillips, F.M.; Argento, D.C.; Balco, G.; Caffee, M.W.; Clem, J.; Dunai, T.J.; Finkel, R.; Goehring, B.; Gosse, J.C.; Hudson, A.M.; et al. The CRONUS-Earth Project: A synthesis. Quat. Geochronol. 2016, 31, 119–154. [Google Scholar] [CrossRef] [Green Version]
  35. Fernandez-Mosquera, D.; Hahm, D.; Marti, K. Calculated rates of cosmic ray muon-produced Ne in subsurface quartz. Geophys. Res. Lett. 2010, 37, L15403. [Google Scholar] [CrossRef]
Figure 1. δ 18 O glacial proxies. Combination of scaled δ 18 O curves from Lisiecki and Raymo [12], Zachos [13], NGRIP [14], and Buizert et al. [15], depicted with colours. Black lines show the simplified curves used by NUNAIT for latitudes north (A) and south (B) of latitude 55 S.
Figure 1. δ 18 O glacial proxies. Combination of scaled δ 18 O curves from Lisiecki and Raymo [12], Zachos [13], NGRIP [14], and Buizert et al. [15], depicted with colours. Black lines show the simplified curves used by NUNAIT for latitudes north (A) and south (B) of latitude 55 S.
Geosciences 11 00362 g001
Figure 2. 10 Be and 26 Al surface muon contributions. (A) Percentage of 10 Be muon production rates with respect to the total muon production rate generated using P_mu_total_alpha1.m [17] for 1678 land sites, and the approximation calculated using Equation (2) for the same sites. (B,C) Share of 10 Be and 26 Al fast muon production with respect to the total muon production at the surface. (D) Ratio between the 26 Al and the 10 Be muon shares.
Figure 2. 10 Be and 26 Al surface muon contributions. (A) Percentage of 10 Be muon production rates with respect to the total muon production rate generated using P_mu_total_alpha1.m [17] for 1678 land sites, and the approximation calculated using Equation (2) for the same sites. (B,C) Share of 10 Be and 26 Al fast muon production with respect to the total muon production at the surface. (D) Ratio between the 26 Al and the 10 Be muon shares.
Geosciences 11 00362 g002
Figure 3. 10 Be muon cross sections. Fast and negative muon production rates scaled to surface values calculated using P_mu_total_alpha1.m [17], for 1678 land sites from ETOPO1_Bed_g_geotiff.tif [16], and random depths between 0 and 100 m below the surface (blue dots). Red lines represent the exponential decay approximations used in this work.
Figure 3. 10 Be muon cross sections. Fast and negative muon production rates scaled to surface values calculated using P_mu_total_alpha1.m [17], for 1678 land sites from ETOPO1_Bed_g_geotiff.tif [16], and random depths between 0 and 100 m below the surface (blue dots). Red lines represent the exponential decay approximations used in this work.
Geosciences 11 00362 g003
Figure 4. Marble Hills NUNAIT results. Results of fitting the NUNAIT model to Marble Hills’ data. Right graph: data are depicted in red, and the models fitting the data within 1 σ confidence level are plotted in blue. The best-fitting model, with a reduced chi-squared value of 7.2, is shown as a black line. Left graphs: probability distribution of the tested parameters.
Figure 4. Marble Hills NUNAIT results. Results of fitting the NUNAIT model to Marble Hills’ data. Right graph: data are depicted in red, and the models fitting the data within 1 σ confidence level are plotted in blue. The best-fitting model, with a reduced chi-squared value of 7.2, is shown as a black line. Left graphs: probability distribution of the tested parameters.
Geosciences 11 00362 g004
Figure 5. Marble Hills NUNAIT results using fit type 1. A reduced chi-squared value of 23 was obtained for the best fitting model. Note that when using fit type 1, models producing apparent exposure ages shorter than sample data are discarded. This results in no data being shown above a certain threshold in the first two graphs on the left column (ice-thinning and weathering plots).
Figure 5. Marble Hills NUNAIT results using fit type 1. A reduced chi-squared value of 23 was obtained for the best fitting model. Note that when using fit type 1, models producing apparent exposure ages shorter than sample data are discarded. This results in no data being shown above a certain threshold in the first two graphs on the left column (ice-thinning and weathering plots).
Geosciences 11 00362 g005
Figure 6. Full graphical output of the NUNAIT calculator for Mt. Hope data set. Left graphs show the probability distribution of the fitting parameters. Right graphs show the best model, the models fitting the data within 1 σ , and the sample exposure ages in black, blue, and red, respectively. The best fit yielded a reduced chi-squared value of 78.7. The bottom graph shows the evolution of the ice surface and the position of the samples according to the model best fit and one sigma results.
Figure 6. Full graphical output of the NUNAIT calculator for Mt. Hope data set. Left graphs show the probability distribution of the fitting parameters. Right graphs show the best model, the models fitting the data within 1 σ , and the sample exposure ages in black, blue, and red, respectively. The best fit yielded a reduced chi-squared value of 78.7. The bottom graph shows the evolution of the ice surface and the position of the samples according to the model best fit and one sigma results.
Geosciences 11 00362 g006
Figure 7. Subsurface 26 Al/ 10 Be ratios. (A) Plot of 26 Al/ 10 Be concentration ratios shown in Balco [17] (coloured circles with error bars) and cross sections predicted by P_mu_total_alpha1.m (black lines), also from Balco [17]. All concentration ratios (measured) scaled to their predicted ratios (B) are plotted as a camel-plot with 68% of its area at ∼ 100 ± 20 % (C).
Figure 7. Subsurface 26 Al/ 10 Be ratios. (A) Plot of 26 Al/ 10 Be concentration ratios shown in Balco [17] (coloured circles with error bars) and cross sections predicted by P_mu_total_alpha1.m (black lines), also from Balco [17]. All concentration ratios (measured) scaled to their predicted ratios (B) are plotted as a camel-plot with 68% of its area at ∼ 100 ± 20 % (C).
Geosciences 11 00362 g007
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Rodés, Á. The NUNAtak Ice Thinning (NUNAIT) Calculator for Cosmonuclide Elevation Profiles. Geosciences 2021, 11, 362. https://doi.org/10.3390/geosciences11090362

AMA Style

Rodés Á. The NUNAtak Ice Thinning (NUNAIT) Calculator for Cosmonuclide Elevation Profiles. Geosciences. 2021; 11(9):362. https://doi.org/10.3390/geosciences11090362

Chicago/Turabian Style

Rodés, Ángel. 2021. "The NUNAtak Ice Thinning (NUNAIT) Calculator for Cosmonuclide Elevation Profiles" Geosciences 11, no. 9: 362. https://doi.org/10.3390/geosciences11090362

APA Style

Rodés, Á. (2021). The NUNAtak Ice Thinning (NUNAIT) Calculator for Cosmonuclide Elevation Profiles. Geosciences, 11(9), 362. https://doi.org/10.3390/geosciences11090362

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop