The NUNAtak Ice Thinning (NUNAIT) calculator for cosmonuclide elevation proﬁles.

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


14
Quantifying the changes in the thickness of the Greenland and Antarctic ice sheets is    54 Herewith I present a set of MATLAB ® / GNU Octave © scripts that form the 55 NUNAIT calculator and their mathematical descriptions. All scripts needed to run the 56 NUNAIT calculator are freely accessible at https://github.com/angelrodes/NUNAIT.

57
When running the script START.m, the user is asked to run the calculator or select 58 previous data to display the text and output.

59
If the first option is selected (Run simulation), two types of files cn be selected: To reduce the number of calculations and the computing time while representing 106 the ice changes relevant to the cosmogenic accumulation, the data is interpolated for 107 ages every 10 years for the last century, every 100 years until 20 ka, every 200 years until 108 50 ka, every 500 years until 100 ka, and every 1% increase for ages older than 100 ka.

109
The resulting simplified curve is shown in Fig. 1. or elevation is not a number, a global average is given. A single value of latitude and 114 elevation is used to calculate the contribution of muons to the total surface production 115 of 10 Be. All other productions are scaled accordingly.

116
The contribution of muons to the total surface 10 Be production (R µ ( 10 Be)) is calcu-117 lated as: 118 P µ ( 10 Be) P total ( 10 Be) = 1 100 · 1.29 + lat 900 This approximation is based on the 10 Be production at 1678 sites equally distributed 119 on land areas according to ETOPO1_Bed_g_geotiff.tif (Amante 2009) and calculated  A. Percentage of 10 Be muon production rates with respect to the total muon production rate generated using P_mu_total_alpha1.m (Balco 2017) for 1678 land sites, and the approximation calculated using equation 2 for the same sites. B and 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.

134
As the uncertainties of these ratios are unknown, this script assigns a conservative 135 20% uncertainty for muon contributions calculated using equation 2.

136
All these data can be changed in the file constants.m.

Muon cross sections 138
To simulate production under ice and rock surfaces, the muon production was 139 approximated as three exponential functions of depth (Granger and Smith 2000). 1678 140 10 Be and 26 Al muon production rates generated using P_mu_total_alpha1.m (Balco 141 2017) were analysed to fit three exponential decays with attenuation lengths of 850, 5000, 142 and 500 g cm −2 (Fig. 3). These attenuation lengths correspond to 75% of the fast muon,

143
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 145 (P µ f ast /P µtotal ) considered for each isotope is:

152
The uncertainties of these approximations are within the uncertainties described in 153 section 2.3 for 10 Be and 26 Al.

154
All these data can be changed in the file constants.m. 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 ice ) 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: where P is the production rate considered (spallation and each of the muon types),

160
Λ is the attenuation length for the production rate considered, λ is the decay constant of type 3 is used, the script assumes that all generated models fit the data.

182
The script selects the climate reference based on the average latitude of the samples.
where C i and σ C i are the sample concentrations and their uncertainties derived from To represent the results as probability density distributions of the parameters, the 221 relative probability corresponding to each model is calculated as: which has a similar shape as the cumulative distribution function of the chi-square    be explained by inhomogeneous weathering ratios of the continuously exposed surfaces.
268 Figure 5. Marble Hills NUNAIT results using fit type 1. A reduced chi-square 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).
The fit type 1 was used to constraint the minimum weathering rate that is compat-269 ible with this data (Fig. 5), implicitly assuming that faster apparent weathering rates  Fig. 6 shows that the optimum fit of the NUNAIT model mimics the distribution of 288 the data for each isotope, but does not fit well the ratios between isotope data.

289
As for the Marble Hills' data, the NUNAIT model fits this data set for very low  1 assumes that the average production rate for the apparent (minimum) age equates 315 to the constant production rate. This introduces differences with the time-dependent 316 production models typically considered for the calculation of surface production rates. available for the calibration of muon production under the surface, which is ∼5% and 335 ∼14% for 10 Be and 26 Al respectively according to Balco (2017, Table 1). The NUNAIT 336 calculator incorporates these uncertainties by considering a 20% uncertainty for all 337 muon-produced concentrations.

338
According to the data summarized in Balco (2017), the best predictions on the muon 339 produced 26 Al/ 10 Be ratios fit the empirical data within ∼20% uncertainty (Fig. 7). For   During ice-free periods, the NUNAIT model considers a homogeneous weathering 363 rate along with the elevation profile. This could not be very realistic for areas with 364 intense periglacial processes that produce increased erosion rates in local areas (e.g. rock 365 falls). When fitting the model to data from areas with evident periglacial processes, the 366 minimum fitting type should be selected (fit type 1 described in section 2.7). go beyond the typical observations deduced from surface exposure dating, like glacial 376 erosion rates and uplift rates. Therefore, the NUNAIT calculator is presented as an easy-377 to-use tool that will help glaciologists to interpret cosmogenic data from nunataks, where 378 exposure histories are usually complex. Moreover, the methods described in section 2