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
Be and
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
O proxy glacial record. In this model, a given
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
Be and
Al concentrations considering fixed values of glacial erosion and subaerial weathering. Knudsen et al. [
10] described a method to solve not only the
O cutoff value but also the glacial and interglacial erosion rates from a set of multiple cosmonuclide concentrations that can include
Be,
Al,
C, and/or
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
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
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 (Be, Al, Ne, He, Cl, and/or 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
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 He, Be, C, Ne, Al, and 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]:
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
for stable isotopes, Equation (
1) results in
.
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 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
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 () and elevation () 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 Be. All other productions are scaled accordingly.
The contribution of muons to the total surface
Be production (
) is calculated as
This approximation is based on the
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
Be and
Al muon production rates, respectively [
17], the uncertainty of the calculated muon contributions based on Equation (
2) should be at least a 7% for
Be and 14% for
Al.
The calculation of the muon contributions for other nuclides are based on the following ratios:
. See
Figure 2.
, according to Heisinger and Nolte [
19].
, according to Balco and Shuster [
20].
, consistent with Blard et al. [
21].
, 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
Be and
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
(
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 () considered for each isotope is:
Be:
. See
Figure 2.
Al:
. See
Figure 2.
Cl:
, according to Heisinger and Nolte [
19].
Ne:
, according to Balco and Shuster [
20].
He:
, consistent with Blard et al. [
21].
C:
, according to Heisinger and Nolte [
19].
The uncertainties of these approximations are within the uncertainties described in
Section 2.3 for
Be and
Al.
All these data can be changed in the file constants.m.
2.5. Densities
A density of
g
is considered for ice [
23], and a density of
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 (
) based on the input conditions (weathering
w, glacial erosion rate, and maximum and current ice levels) for each time range (
) defined by the climate curve and for each sample. The model concentration is calculated as
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 is calculated by adding all the for all production types and time ranges.
The effect of the glacial erosion rate is ignored inside each time range (), as usually is much more sensitive to 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:
where
and
are the sample concentrations and their uncertainties derived from the apparent surface exposure ages (
Section 2.1),
are the model concentrations corresponding to sample
i (
Section 2.6), and
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
, and models fitting the data within a 2
confidence level are defined by the ones with
. 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 (
) for poor fittings and high
. The formula
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 , this confidence range is increased to 4, 8, 16, etc. 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
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
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.
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 (
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
Be and
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
Al/
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 (
) for any synthetic concentration ratio under the ice-sheet, which is probably an overestimation of the
Al/
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.