| 21 Jul 2026
iceHIST v1.0: a cosmogenic-nuclide model suite for testing glacier and ice-sheet histories
Abstract. Cosmogenic-nuclide inventories in glaciated landscapes record changes in glaciers or ice sheets, and other surface processes. Scenario-based models can extract glacier and ice-sheet history information from these inventories, but many approaches remain specialised to a single site, data type, scenario structure, or inference method. Here I present iceHIST v1.0, a model suite for testing hypotheses about glacier and ice-sheet history using in situ 10Be, 14C and/or 26Al. iceHIST supports depth profiles (e.g. bedrock cores), surface samples as transects, and combined core-plus-transect datasets, including modern subaerial and subglacial samples. Four scenario families define the tested histories: OneStage, DynamicPeriods, TimeSeries and IceSurface. These scenarios encompass simple cover histories, repeated exposure and burial intervals, thresholded time-series histories, and explicit or imported ice-surface histories. For each scenario, a common forward model integrates nuclide production over depth and time, accounts for radioactive decay, ice cover, and sample thickness, and includes options for snow shielding, subaerial and subglacial erosion, and inherited inventories. Model predictions are compared with data using Monte Carlo screening, Markov chain Monte Carlo (MCMC) posterior exploration, or likelihood weighting of fixed imported histories. Synthetic demonstrations illustrate surface and depth-profile recovery, time-series threshold recovery, thinner-than-present lowstand testing, a Monte Carlo-to-MCMC workflow, and ranking of imported ensemble histories. Applied in this manner, iceHIST can determine how long glaciers were larger or smaller than today, whether ice sheets thinned below current levels and subsequently re-thickened, and identify which glacier or ice-sheet simulations best match cosmogenic-nuclide observations.
General comments:
iceHIST is a suite of models designed to test the compatibility of various ice cover histories against measured cosmogenic nuclide inventories. I believe the author has correctly identified a gap in the current published works. While other works have endeavored to predict and model cosmogenic nuclides in order to reveal possible ice-cover histories, this is the most unified and flexible tool yet with a highly customizable setup. A code/model like this is very valuable to the community, both for non-modelling quaternary geologists to test their interpretations of glacial histories, but also for other cosmo-modelers that have been keeping personal, poorly documented code repositories to solve these kinds of problems. Hence, iceHIST presents itself as a novel tool that is ready to be used by the scientific community, with a clear, transparent framework that allow for reproducibility of results. This integration of multiple models makes for a welcomed advance in cosmogenic nuclide simulation/analysis and hence it is well worth publishing a few revisions.
Specific comments:
1) Examples in section 7 are missing explanations and documentation. All examples are missing important input parameters. For example, input bounds on the parameters are not reported, and neither are fit or acceptance criteria. In this very manuscript, the author himself strongly emphasizes that such information should always be reported. Lines 666-667: “MC results should be reported with the scenario family, data mode, sampled bounds, … fit or acceptance criterion…”.
It is also not mentioned in the examples section what uncertainty is prescribed the synthetically generated nuclide concentrations that are used as data.
These shortcomings together make it hard to properly judge whether the nice fits in the examples could be because of too tightly chosen bounds, or low assumed nuclide uncertainty.
A minor comment regarding the examples that confuses me a little is that in Fig. 6a the ‘true’ scenario is 121-125 ka, and the ‘best’ fit is also 121-125 ka, yet these two lines look offset?
2) In section 6.2 I’m missing an explanation of how the new parameter vectors are proposed/ perturbed from the previous accepted parameters. And are the bounds imposed in the MC simulations also used in the MCMC simulations, or are they completely free?
3) I do not think it is clear exactly what scenarios employ partial shielding versus complete shielding during ice cover. The IceSurface is well explained how the ice thickness affects the production throughout. At first read, I thought maybe the DynamicPeriods used 100% shielding during cover, but then I looked at Figure 2, and saw the term ‘Relative cover’ on the y-axis. How should this be understood? Is this a factor that is multiplied by the production rate? I’m guessing not since the cover thickness and production rate do not scale linearly. Anyways, I would like to see this better explained in the text.
4) Regarding the inheritance. In only very rare cases is it possible to give a good prediction of the inherited nuclides in any sampling situation. Even guessing the inheritance bounds could be quite tricky. For example, in the OneStage scenario inheritance of nuclides could completely compensate for the exposure duration. Hence, for making assumptions, I think there are two end-member scenarios that would be to most natural to assume, a) that the inheritance is 0 (this is implemented), b) that the samples are fully saturated at the start of a simulation at their given depth/position. I think implementing this directly would be very valuable for testing scenarios. And all the codes are already there to make the calculation I’m sure.
5) I’ve tested some of the codes from the repositories, which I downloaded from the zenodo link.
I think there is a lack of discussion of performance of the code, since many people will be hoping to run the codes on laptops I presume. When running these codes, there are no time estimates. There is not even any progress bars indicating the running time, particularly the MCMC perturbations can take a very long time to run. My suggestions: A progress bar is a must for these types of codes, for general accessibility. Secondly, I would like the run-times reported for the examples presented in the manuscript, alongside the computer setup used.
On this note, I was also wondering about the lack of parallel computing. As far as I can see it is never applied. It’s generally an easy implementation in matlab, unless the code architecture prevents it. MC methods are highly parallelizable and using 8 cores instead of 1, well, it is 8 times as fast. Even for the MCMC methods, the different chains could run on separate cores. Perhaps, the author could leave a comment on this choice.
Another final comment on the examples provided in the repository. Some of them were not running correctly directly on run, which I think they should. I encourage the author to test all again before publication:
I did not run all codes, but some example errors:
Model_icesurface_core_transect.m gave this error, which I was unable to fix:
” Error using regexprep
The 'STRING' input must be either a char row vector, a cell array of char row vectors, or a string array.
Error in save_icesurface_results>normalise_data_type_label (line 194)
data_type = regexprep(data_type,'[^A-Za-z0-9_]+','_');
Model_icesurface_imported_ensemble_weights.m gave an error on run with prod_max_depth = 100. I was able to fix the error after reading the code comments and changing it to 1000.
It has to be said that I am using the 2025b version currently, so I guess there is a small chance that could be the problem. If everything runs smoothly on the authors end.
Technical comments:
Lines 163 and 165 are almost identical.. “This hierarchy is…” Line 238 again repeats this message.
Section 3.1 also sounds quite repetitive of section 2.3.
Line 348: “For bedrock-core samples, zs,i is the core-top or bedrock-surface elevation.” Which is it, I assumed it is bedrock-surface, since the core could be missing parts of the top.
Lines 412-414. You use the word “thickness” for two different quantities in the same sentence.
Line 422: “Missing measurements should be entered as NaN..” I noticed in the supplied example .xlsx sheet that 0 values were entered, so maybe this is not correct?
Lines 464: I think units for t, d, and Z could be specified here.
Line 481: Rock density unit is missing, unit of t is still also not specified.
Equation (4): Is semi-colon intended between t and S?
Line 646 and 725: the symbol q is defined as two different things in section 6, I suggest changing the symbol of one.
/J. Nørgaard