By far the biggest pain in the neck in the overall field of cosmogenic-nuclide geochemistry is the need to compute nuclide production rates due to muon interactions. Production rates due to muons are quite small relative to total production rates at the surface, but they’re not small enough that you can ignore them completely. In the subsurface, they’re the majority of production, so are actually important for applications such as burial dating and depth profile dating, but they’re complex to calculate and hard to verify due to a general lack of subsurface calibration data. It’s much more difficult to obtain useful calibration data for production rates due to muons than for surface production rates due to spallation, mainly because there aren’t really any geologic processes that completely remove the entire preexisting muon-produced nuclide inventory, which extends to tens of meters depth, all at once, leaving a blank slate. And if there were, production rates are really low, so the blank slate would have to sit around for a really long time  before you could do the measurements accurately, in which case something else bad would probably happen to it. So the only hope for decent calibration data is to collect tens-of-meters-long cores from bedrock surfaces that have been sitting around at extremely low erosion rates for millions of years, and there aren’t a whole lot of such things.

The last couple of years, largely due to work in the CRONUS-Earth and CRONUS-EU projects, has either improved or unimproved this situation by generating some new methods of computing muon production rates and also by collecting several new bedrock cores. This increases the overall complexity of the situation as a whole, but also helps a bit by providing enough data to make a decent start at evaluating which methods of computing muon production rates work best. In the last few months I’ve gotten tired of being perpetually confused on this subject, so I’ve waded through all these data and put together an extraordinarily long, tedious, and pedantic document on the subject of benchmarking muon production rate calculation methods against the available calibration data set. I imagine that I will submit this to an actual journal at some stage, but for now it is here (the PDF file). Warning: it is long, tedious, and pedantic. However, if you want my best shot at a synoptic view of all reasonably accessible knowledge that is relevant to production by muons, it might be potentially useful. Of course, it hasn’t been reviewed in detail by anyone else, so it might also be wrong. The MATLAB code that does all the calculations is also available at the same address (the zip-file).

Presumably, about six people on Earth will be interested in this document. However, I imagine that all six of those folks are also regular readers of this blog. You know who you are.  As noted, this will probably be sent off to a journal at some stage, in which case some of you will also have the truly unpleasant and unfortunate duty of reviewing it. Enjoy!

Oh yeah. The 2%-of-production-98%-of-suffering logo above is also in that directory, in case you want to put it on a T-shirt. For example, from CafePress:

A large amount of the scientific literature on cosmogenic-nuclide geochemistry uses the term “production rate scaling scheme” to refer to one of a number of computational methods used to calculate how cosmogenic-nuclide production rates vary geographically and over time. Of course, the word “scheme” is intended to have its fairly common scientific meaning of a formal mathematical procedure for carrying out some task, as one would refer to a numerical integration scheme or an implicit solution scheme. As most people are aware, however, that’s not really what’s connoted by “scheme” in the rest of the English language. An article in Physics Today from 2011 called attention to this with the following table:

This list is absolutely right on target. In the admittedly unlikely event that a non-cosmogenic-nuclide geochemist were to read the cosmogenic-nuclide geochemistry literature, they might legitimately be mystified as to why you had to choose one of several dastardly plots, whilst rubbing one’s hands together with an evil cackle, to compute nuclide production rates. I’m probably the worst past offender in this regard, starting with this paper, and I’m not proposing to give up terms like “positive feedback” that are specifically meaningful because they are derived from the presence of positive or negative signs in mathematical descriptions of the phenomenon, but the sad fact is that “scheme” doesn’t have a whole lot of redeeming qualities. It really does sound like a devious plot, and it doesn’t have any special meaning that’s not equally well communicated by simpler words like “method” or “procedure.” We don’t need it, so let’s not use it. Henceforth, I suggest just talking about production rate scaling methods, procedures, or possibly algorithms; of these, I think I prefer “methods.” If anyone has any better ideas, let me know, but “scheme” is now out as far as I’m concerned.

This post is mostly about production rate scaling models, but in addition is neat because it’s an actual useful application of the ICE-D Antarctic exposure age database. It’s a pretty simple one, but hey, mighty oaks from little acorns grow. This is another immensely long post, but the short summary is that a look at samples from Antarctica with very high Be-10 and Al-26 concentrations seems to indicate that the relatively recent, and extremely complicated, production rate scaling model of Nat Lifton and colleagues (the “LSD” scaling method) works better at high latitudes than the simpler and very widely used scaling method based on the work of Devendra Lal in the 1980’s. This is interesting because so far it has been quite difficult to determine whether or not there is any substantive difference in performance between these scaling methods.

We start with some recent Be-10 measurements by Gordon Bromley of the University of Maine on sandstone boulders, presumably glacially transported, at the same site in the southern Transantarctic Mountains where we recently observed what appears to be the highest concentration of cosmogenic He-3 ever observed in a terrestrial sample. These Be-10 measurements are interesting because, even after a reasonably thorough effort to make sure that no mistakes were made, they appear to show that these samples contain an impossible amount of Be-10. What does impossible mean? Basically, if a surface is left undisturbed and exposed to the cosmic-ray flux long enough, the concentration of cosmic-ray-produced Be-10 will eventually build up to a level high enough that the rate that Be-10 is lost by radioactive decay is equal to the production rate. At this point, the amount of Be-10 can’t increase any more. So this relationship provides a limit to the maximum possible concentration of Be-10 that can be present in quartz at a particular site. Generally this condition is referred to as either ‘production-decay equilibrium,’ or, more commonly, just ‘saturation.’ In math, accumulation of Be-10 in a non-eroding surface is governed by the differential equation:

$\frac{dN}{dt} = P - N\lambda$

where N is the Be-10 concentration (atoms/g), P is the production rate at the site (atoms/g/yr), and $\lambda$ is the decay constant for Be-10 (4.99e-7 /yr). So saturation occurs when $dN/dt = 0$, or when $N = P/\lambda$. Thus, if we know the production rate at some site, we equivalently know the saturation concentration.

At Gordon’s sample site, the production rate (calculated with the Antarctic atmosphere model and scaling scheme of Stone(2000) and the “primary” production rate calibration data set of Borchers and others (2015), the Be-10 production rate is 38.2 atoms/g/yr, which means the saturation concentration for Be-10 is 7.65 x 10^7 atoms/g. However, the sample actually contains 8.5 x 10^7 atoms/g. Inconceivable!  If we accept this production rate estimate, this is an impossible Be-10 concentration.

OK, that seems weird, but it’s not actually all that  unusual in Antarctica. To show this, look at nearly all known Be-10 and Al-26 measurements from Antarctica as they are represented in the ICE-D:ANTARCTICA database (making sure they are all correctly normalized to the ’07KNSTD’ and ‘KNSTD’ standardizations, respectively). This data set looks like this:

What I have done here is just plot nuclide concentration vs. elevation. Because production rates increase with elevation, the nuclide concentration associated with a particular exposure age also increases with elevation, so the envelope of possible nuclide concentrations that one could potentially observe widens as elevation increases. There are about 1200 Be-10 measurements and 500 Al-26 measurements represented here.

Now, because we are using an atmosphere model and scaling scheme that vary only with elevation (magnetic cutoff rigidity is pretty much effectively zero in Antarctica for all practical purposes), predicted saturation concentrations for Be-10 and Al-26 are also only a function of elevation, so we can plot them on these axes, as follows:

Here the black line shows predicted saturation concentrations. It is evident that there are lots of samples in Antarctica, mostly at high elevation above about 2000 m, whose Be-10 concentrations appear to be above saturation. Gordon’s samples are sitting off to the right of the line at about 2300 m elevation. Some measurements are up to 12% above saturation concentrations calculated in this way.

There are four possible reasons for this. One, the measurements are somehow messed up. This is unlikely, because the measurements represented here are from several different laboratories and/or AMS facilities, and it’s hard to come up with any reason why they would all be similarly spurious. Certainly it is suspicious that the maximum amount that some samples exceed predicted saturation concentrations (12%) is about the same as the difference between the KNSTD and 07KNSTD standardization (11%). However, I spot-checked most of these samples, most of the measurements are well documented, and I am pretty sure that this error has not been made.

The second reason is that we might be overestimating the atmospheric pressure at high-elevation sites, and therefore underestimating the production rate and the saturation concentration. That possibility we can evaluate by using a totally different atmosphere model, one derived from the ERA40 reanalysis and adapted for use in estimating production rates by Nat Lifton (you can read about it here). To do this, we need slightly different axes, because in the ERA40-derived atmosphere model, atmospheric pressure varies with location as well as elevation. So instead of plotting just raw nuclide concentrations, we will plot the ratio of the measured nuclide concentration at a particular site to the saturation concentration predicted for that site. For samples above saturation, this ratio will be greater than one.

The two left-hand panels use the Stone (2000) Antarctic atmosphere, and the two right-hand panels use the ERA40-derived atmosphere. Basically, there is no difference. Both atmosphere models with this scaling scheme indicate that there are a bunch of samples at high elevations that exceed saturation concentrations for both Be-10 and Al-26. The fact that two independent atmosphere models agree on this point would tend to indicate that our problem cannot be explained just by a bad atmosphere approximation.

The third reason is that the sites where we observe greater-than-saturated nuclide concentrations might have changed elevation. If the samples had spent a significant amount of their exposure history at an elevation higher than their present elevation, they could potentially have higher nuclide concentrations than they could ever achieve at their present elevations. Elevation changes could be caused by local effects (for example, all the samples are on the hanging wall of an active normal fault), or by regional effects (for example, changes in dynamic topography during the several-million-year exposure history of these samples; or glacial isostasy, if the ice sheet surrounding sample sites had continuously thickened during the period of exposure). The amount by which the offending samples exceed predicted saturation concentrations requires approximately 300-400 meters of elevation change in the last few million years. We can evaluate some of these possibilities by looking at where these samples are actually located in Antarctica. Here they are:

The white dots are Be-10 measurements exceeding saturation, and the gray ones are for Al-26. There are some more white dots hidden behind some of the gray dots. These samples are widely distributed around the high mountains of East Antarctica, which makes it unlikely that concentrations above saturation are due to local faulting or mass-movement effects at a particular location. The background image in this figure is a model for changes in dynamic topography in the last 3 Ma from a recent paper by Austermann and others; the color scale shows total uplift in the past 3 Ma in meters. Samples exhibiting Be-10 concentrations above saturation are not located in areas that have experienced signficant regional subsidence, at least due to this process, over long time periods. In fact a few of them are located in areas that are predicted to have experienced significant uplift. So we can probably exclude local faulting as a blanket explanation for the many samples in Antarctica that have Be-10 concentrations above saturation, and we can also exclude long-wavelength elevation change due to mantle dynamics. What about glacial isostasy? If the ice sheet in the vicinity of these samples had continuously thickened over the past several million years, then the elevation of the samples would have lowered during that time due to isostatic compensation, potentially leading to above-saturation nuclide concentrations. However, we need this change to be large — order 1000 m of ice to obtain the needed few hundred meters of isostatic response — and this is not consistent with sea level records or ice sheet modeling.

For example, the above figure shows changes in bed elevation from a 5-million-year-long Antarctic ice sheet simulation by David Pollard and others, at two locations that are representative of where we see Be-10 concentrations above saturation: the central Transantarctic Mountains (red), and the East Antarctic marginal highland in the area of the Sor Rondane mountains (blue). At both of these sites, this model run predicts a long-term decrease in land surface elevation over the last ~4 Ma, presumably driven by Antarctic ice sheet expansion due to lower sea levels in the Pleistocene. However, the magnitude of the decrease is only several tens of meters, much too small to account for the above-saturation Be-10 concentrations. Thus, long-term lowering of sample sites due to glacial isostasy does not appear to be a good explanation for above-saturation cosmogenic-nuclide concentrations either.

The fourth possibility, of course, is that our production rate scaling is incorrect. I’ve left it for last, of course, because it’s the most likely. We don’t have any Be-10 or Al-26 production rate calibration data at high elevation near the poles (see here), so we are relying on the scaling method of Stone (2000), which is of course just based on the work of Devendra Lal, to extrapolate from production rate calibration sites at low elevation to sample sites at high elevation. We can test whether inaccuracy in this scaling method explains the apparently-above-saturation samples, by trying a different scaling method. Specifically, let’s try the more recent (and much more complicated) scaling method of Lifton and others (2015), now commonly known as the ‘LSD’ scaling method. At low magnetic cutoff rigidity (i.e., in polar regions), this scaling method predicts a larger altitude dependence than Lal-based scaling methods. Because magnetic field variability is out of the picture at polar latitudes, this elevation dependence is the primary difference between the two scaling methods. Here are the results, calculated in the same way as for the four-panel figure above (compare with the first two panels of that figure):

Using the LSD scaling method completely fixes the problem (for both atmosphere models, although only one is shown here). Here is another view that copies the first figure that I started with above, showing raw concentration vs. elevation, compared with saturation concentrations. The solid lines are saturation concentrations predicted using the Stone/Lal scaling method as shown above, and the dashed lines are saturation concentrations predicted using the LSD scaling method. If we use the LSD scaling method, there are no measured Be-10 or Al-26 concentrations in Antarctica that are impossibly high.

So, overall, nuclide concentrations in many high-elevation Antarctic surfaces are inconsistent with production rates estimated using the Lal/Stone scaling method, but are consistent with production rates estimated using the LSD method. This would appear to imply that LSD scaling performs better for this data set, and, by implication, very likely for high-elevation sites in polar regions generally.

This is interesting for several reasons.

First, it would tend to indicate that one should probably use LSD rather than Lal-based scaling to compute exposure ages from old, high-elevation glacial deposits in  Antarctica. It makes a pretty big difference: almost 20% at 2500 m elevation. For example, you might find a moraine at this elevation that belonged to the mid-Pliocene warm period (ca. 3-3.3 Ma) if you used Lal-based scaling, but was coeval with the onset of Northern Hemisphere glaciation (ca. 2.7 Ma) if you used LSD scaling. Certainly some potential for confusion there.

Second, although the LSD scaling method is very much more complex than the Lal-based methods, and there are some physical arguments that indicate that LSD should be more accurate than Lal, so far it has been very difficult to show that there is any practical difference in performance between the two. This  issue is important mainly because of the complexity of the LSD scheme — it requires a lot more computation, so code to implement it is more complicated and very much slower. If we don’t have to use it, we’d rather not. The calibration exercise of Borchers and others (2015) showed that for available Be-10 and Al-26 calibration data, there was basically no detectable difference in performance, and that’s been taken (certainly by me) to indicate that for most practical purposes, it’s fine to use the simpler method. However, the Antarctic exposure-age data set would tend to indicate the opposite: that the Lal-based methods, when calibrated with existing calibration data sets that lack coverage at high altitudes in polar regions, significantly underestimate production rates at these sites.

Third, if we accept that the LSD scaling method is correct, we’ve gone from the conclusion that there are a lot of Be-10 and Al-26-saturated surfaces in Antarctica, to the conclusion that there aren’t any. No known measurements from Antarctica quite overlap with predicted saturation values. In one sense, this is not really a huge surprise, because for a site to reach true production-decay saturation, the surface erosion rate must be exactly zero for several million years. Really exactly zero erosion for millions of years is a big ask — this doesn’t even occur in outer space — and is probably unlikely to occur in reality, even at -60° C in the upper Transantarctic Mountains. Again, if we accept that the LSD scaling method is correct, we can estimate steady-state erosion rates (that is, maximum limits on the erosion rate) for the samples that are closest to saturation; for the Be-10 measurements this is in the range of 1-2 cm per million years; for the Al-26 measurements it is 3-4 cm per million years. It is not clear why Be-10 measurements are closer to saturation than Al-26 measurements, although because they are all very close to saturation, this calculation is quite sensitive to the calibrated reference production rates for both nuclides, so this could easily be explained by inaccuracies in those values. In any case, these are plausible values for the lowest erosion rates in Antarctica. Not zero, but extremely slow.

So the summary here is that it appears that the LSD scaling method is consistent with the overall data set of Be-10 and Al-26 concentrations in Antarctic surfaces, and the Lal-Stone scaling method isn’t. One final point, however, is that Stone (2000) also considered this issue (with fewer data; see Figure 4 in that paper) and concluded that the Lal scaling method did not, in fact, undershoot saturation concentrations. The difference, I am fairly sure, is that new calibration data imply lower production rates than the calibration data available at the time of that paper; if currently available production rate calibration data had been used in that paper they would have led to the same conclusion I come to here, that Lal-based scaling underestimates saturation concentrations at high elevations in Antarctica. However, there are a lot of moving parts involved in comparing the two calculations (Be-10 standardizations, production rates, decay constants, etc….) and that claim could probably use further investigation.

In a previous post I pointed out that we (that is, the noble gas mass spectrometry facility at BGC) might have recently analysed the sample with the highest concentration of cosmogenic Ne-21 ever observed in a terrestrial sample. By ‘terrestrial sample,’ of course, I mean a rock that has accumulated its entire cosmogenic-nuclide inventory on Earth — because cosmogenic-nuclide production rates are so much higher outside the Earth’s atmosphere, meteorites have orders of magnitude higher concentrations of cosmic-ray-produced radionuclides than possible on Earth. So meteorites, or extraterrestrial-dust-laden ocean-floor sediments from the central Pacific, don’t count. In any case, no one has shown up to tell me about anything with a higher Ne-21 concentration, so the claim may be true.

The point of this post is that it now appears that we may also have analysed the sample with the highest concentration of in-situ-produced cosmogenic helium-3 ever observed in a terrestrial rock. This sample is near 2300 m elevation at Roberts Massif, a nunatak at the head of the Shackleton Glacier in the Transantarctic Mountains at 85.5° S latitude. In fieldwork last year, Gordon Bromley and his colleagues collected this sample from that site; the sample is a loose clast, coarsely faceted and therefore possibly glacially transported, of the local bedrock, a Jurassic mafic intrusive rock known as the Ferrar Dolerite. Here is a picture (photo credit: Gordon Bromley).

Two separate analyses at BGC of pyroxene separated from this sample yielded 1.39 x 10^10 and 1.35 x 10^10 atoms/g He-3. Two adjacent smaller loose rocks — this one and this one — had slightly lower, but still pretty unreasonable, He-3 concentrations near 1.25 x 10^10 and 1.07 x 10^10 atoms/g, respectively. As surface rocks go, that’s a lot of helium-3. Note, however, that even though He-3 is extremely expensive at present, you would still have to mine something on the order of 1000 metric tons of this pyroxene to have \$1 worth of He-3. It’s not ore-grade.

Pyroxene in this lithology is known to contain only about 6 x 10^6 atoms/g of non-cosmogenic He-3, presumably inherited from initial mineral crystallization, so the measured concentration in these samples is essentially all cosmic-ray-produced. 1.39 x 10^10 atoms/g cosmogenic He-3 at this site implies an apparent exposure age of approximately 13 Ma. As discussed in the previous post, this is not the oldest exposure age ever measured; that record still belongs to samples from the Atacama Desert, analysed by Tibor Dunai and colleagues, with apparent exposure ages exceeding 30 Ma. However, those samples are at much lower elevation where production rates are also much lower, so even though they are older, the total cosmic-ray dose they have experienced is smaller.

In any case, a relatively slapdash literature search reveals that the nearest competitor appears to be another sample of pyroxene from the Ferrar at Mt. Fleming, also a high-elevation site in the Transantarctic Mountains. Joerg Schaefer and colleagues found 6.9 x 10^9 atoms/g cosmogenic He-3 in this sample, which is still plenty but is almost a factor of two less than in the Roberts Massif sample. Are there any other contenders?

The point, once again, is that these high-elevation rock surfaces in Antarctica have enjoyed a stunning and expansive view of nearly the entire late Cenozoic evolution of the Antarctic ice sheets. Frankly, a lot of that was probably pretty boring, but if the East Antarctic Ice Sheet really collapsed in past warm-climate periods, they watched.

This post is about the current proliferation of online exposure age calculators. Summary: good thing. Details: complicated.

Here is a short summary of the situation. Actually, it’s a long summary. But I can’t make it any shorter.

2008-v2. The initial online exposure age calculator at hess.ess.washington.edu, also sometimes referred to by the deeply unfortunate and ill-conceived nickname “Balculator” after the first author of the accompanying paper, dates from 2008 and reflected current thinking at that time about how best to compute exposure ages. The initial funding for me to develop it (a Linux server and several months of postdoctoral salary) came from the then-just-beginning, NSF-funded “CRONUS-Earth” project, so it has also been widely known as the “CRONUS calculator;” however, as discussed below, this is no longer quite accurate. To avoid confusion in the rest of this post, I’ll refer to it here as “the 2008-v2 calculator,” because at the time it entered widespread use it was denoted by this version number, which is often reported in papers and other downstream uses.

Aspects of this calculator that have subsequently turned out to be important are twofold. First, it includes five different production rate scaling methods, two based on the work of Devendra Lal (“Lal-based schemes”) and the rest on later interpretation of more recent neutron monitor data (“neutron-monitor-based schemes”). This  probably excessive diversity was simply because at the time no one knew which one worked better. Second, reference production rates for Be-10 and Al-26 were based on a compilation of currently available calibration data (the “2008 global calibration data set”).

Between 2008-2015 there were several changes to this calculator that were the result of discoveries that various parameter values in version 2 were incorrect.

2008-v2.2. This revision (in 2009) made the important change that users were now required to enter information about isotope ratio standardizations used for Be and Al AMS measurements. The details of this are a giant can of worms that you can learn about here. Interestingly, this turned out to be a really effective example of incentive-based behavioral engineering in that users were successfully induced to solve a messy and unpleasant problem (standardization-related confusion) by being presented with an incentive to do so (the opportunity to outsource all the hard work of exposure age calculations to the online calculator). Also, a subsequent 2.2.1 revision changed a few parameter values, most notably the Be-10 decay constant, that had only minor effects.

2008-v2.2 with alternate production rate calibration data sets. Shortly after 2008, it became clear, from new production rate calibration data that were being generated, that the 2008 global calibration data set was simply wrong, predicting Be-10 and Al-26 production rates that were approximately 10% too high. The reasons for this are diverse and not entirely clear in some cases, but the overall result is quite clear. Although this was basically a non-issue for most applications of erosion rate calculations, it was a big issue for exposure-dating (for example, because the ages of the Younger Dryas and the Antarctic Cold Reversal differ by about 10%). To deal with this issue, I added features to the 2008-v2.2 calculator, mainly the capability to enter arbitrary production rate calibration data, use them to generate a best-fitting value of the reference production rate, and then use this production rate value to determine exposure ages at unknown-age sites. In addition, I carried out this exercise for a variety of published alternatives to the 2008 global calibration data set and made them available as distinct data entry pages. Thus, although the erroneous 2008 calibration data set was still presented as the default means of calculating exposure ages, the availability of other calibration data effectively solved this problem for most practical purposes.

This version (2008-v2.2 with various calibration data sets) has been widely used between ca. 2009 and the present.

Results of the CRONUS-Earth project.  Starting in 2014, results of the now-completed “CRONUS-Earth” project began to become available. Of these results, the ones most relevant to the issue of online exposure age calculators are as follows.

1. Calibration data. A lot more calibration data now exist, and they continue to indicate that the 2008 data set gives incorrect results.

2. Scaling schemes. It is now clear that the neutron-monitor-based scaling schemes are inaccurate, for reasons to do with how neutron monitors work. Nat Lifton and colleagues have put together a new class of scaling scheme based on particle transport models (“Sato-based schemes” after the author responsible for the particle transport modeling) that works better. Both Lal-based and Sato-based schemes appear to work indistinguishably well for most practical purposes, although they are not indistinguishable — they make different predictions for some situations.

3. Muons. We have better estimates of muon interaction cross-sections. This is only marginally relevant for surface exposure dating, but has a nonzero effect on erosion rate calculations.

Of these results, only (1) above has been dealt with at all in the 2008 calculator framework. Thus, at this point we have the situation where the 2008-v2.2 calculator includes incorrect default calibration data, scaling schemes that we know to be inaccurate, and muon interaction cross-sections that we know to be wrong. The calibration data issue can be fixed by use of alternate data sets, but not the others.

CRONUS-2016. In 2015-2016 (the accompanying paper is dated 2016), Shasta Marrero, Brian Borchers, and Rob Aumer put together a new online calculator, intended as one of the capstone products of the overall CRONUS project, that is here. Thus, this product is now the “CRONUS-Earth online exposure age calculator,” and, if one can learn from history, will presumably in future also be known as the “Shastalator” or “Marrerolator.”

One implication of this is that it appears that the 2008-v2.2 calculator is no longer the “CRONUS calculator.” Therefore, in homage to Prince, the 2008-v2.2 calculator will henceforth be known as “The online exposure age calculator formerly known as the CRONUS-Earth online exposure age calculator.”

The CRONUS-2016 calculator has the following features.

1. Even more scaling schemes. Sato-based schemes (n = 2) have been added to the existing Lal-based (n = 2) and neutron monitor based (n = 3) schemes from 2008. Thus, n = 7.
2. More complex code. More physical processes are represented in the code.
3. More nuclides. It is possible to compute exposure ages not only for Be-10 and Al-26 measurements, but also C-14, He-3, and Cl-36.
4. No shortcuts with respect to numerical precision. The 2008 code takes many numerical shortcuts to speed up calculations that maintain acceptable accuracy for surface exposure dating purposes, but sacrifice some capabilities, such as, for example, the ability to accurately calculate subsurface production rates. These shortcuts are not taken in the 2016 calculator.
5. Default production rate calibration data derived from recent work.

Items 1-4 above were done for an important reason: the code that the online interface is running was designed also for quantitative testing of all known scaling schemes against all known calibration data (see this paper), and clearly in carrying out such an exercise it is important to take numerical imprecisions and unphysical shortcuts out of the picture. However, the downside of these features is that the code is quite slow; instead of immediate return of results via the web server as in the 2008-v2.2 calculator, the input page of the CRONUS-2016 calculator starts an often-fairly-time-consuming compute job on the server, which emails you the results at a later time when the job is complete.

2008-v2.3. At present, despite the existence of the CRONUS-2016 online calculator, there is still quite a lot of user demand for the 2008-v2.2 calculator, presumably for the following reasons. One, folks are used to it, so it is like a pair of nice fuzzy old socks. Two, it is faster than the CRONUS-2016 code. Three, it facilitates use of non-default calibration data. Possibly, four, it has a web service API. This is all fine, but as noted above, the 2008-v2.2 code perpetuates three major errors: incorrect default calibration data, incorrect scaling schemes, and incorrect muon interaction cross-sections. Thus, I have put together a version 2.3 that mitigates two of these. Specifically, the default production rate calibration data set is now the same (the “CRONUS primary data set” of Borchers and others, 2016) as is used in the CRONUS-2016 calculator, and the muon interaction cross-sections are also updated to reflect calibration data also from that paper. What’s not fixed is that the neutron-monitor-based scaling schemes (now obsolete) are still there, and the Sato-based scaling schemes are not there. However, if you can ignore the neutron-monitor-based scaling schemes, this update removes everything from the 2008-v2.2 code that we actually know to be incorrect. Thus, for most practical purposes the 2008-v2.3 calculators can be used for surface exposure dating with equivalent accuracy to the CRONUS-2016 calculators. Version 2.3 is now the default version at hess.ess.washington.edu. However, it will not be updated further, but, hopefully, instead replaced by:

v3. The v2.3 update doesn’t fix the scaling scheme issue, and in addition it would be nice if the 2008-v2.2 calculators could be used at least for the other primarily spallogenic nuclides (C-14, He-3, and Ne-21), even if not Cl-36 yet. To deal with these issues, there now exists a developmental version 3 of the online exposure age calculators formerly known as the CRONUS calculators, which is here. The design principle of the v3 calculator is to do only exposure-age and erosion-rate calculations for surface samples, but to do them as fast as possible while maintaining acceptable accuracy for these purposes. Thus, it includes a highly streamlined version of  Sato-based scaling that works by interpolation of precalculated gridded data and runs a couple of orders of magnitude faster than the code in the CRONUS-2016 calculator. Muon production systematics are also highly simplified. In addition, the v3 calculator ingests Be-10, Al-26, C-14, Ne-21, and He-3 data for various mineral targets (simultaneously, and also does the multiple-nuclide plots, which is fun). At the date of this writing, I believe it works correctly, but it hasn’t been extensively checked and one should expect that it will be modified fairly often in the near future.

The ICE-D database. One motivation for developing a fast v3 calculator is to facilitate the development of online databases of both exposure-age data and production rate calibration data, such as those here and here. As discussed in other posts, the idea of these databases is that they pretty much make the current approach to online exposure age calculators obsolete, because raw observational data such as nuclide concentrations lives behind the scenes in a database and exposure-age calculations happen dynamically and transparently when viewing the data. This completely deals with the problem of comparing published exposure-age data that were calculated inconsistently. But for it to work, the calculations have to be fast. Hence the need for v3. At present, the developmental v3 code is the back end for the ICE-D database.

Another goal of this database development is to improve how we deal with production rate calibration data. The idea of the ICE-D production rate calibration database is that it represents a compilation of up-to-date production rate calibration data that is generally believed to be complete and accurate, and can be available to whatever online calculator system wants to use it. This overall principle was one of the motivations for…

CREp. This is a new online exposure age calculator put together by Leo Martin, PH Blard, and their colleagues at CRPG in France, that is here. It includes Lal-based and Sato-based scaling schemes and has a large range of options for production rate calibration data sets that are derived from the ICE-D calibration database. This is a good step towards solving the data-assimilation problem for production rate calibration data (which continues to be steadily generated)…theoretically, as we add new data to the ICE-D database, it should be immediately available to the CREp online calculator. CREp also has by far the most attractive user interface (although unfortunately this is not a very high bar).

To summarize. The current situation features at least four different options for online exposure age calculations. One and two, the online calculators formerly known as the CRONUS calculators — the 2008-v2.3  and v3 calculators at hess.ess.washington.edu. Three, the CRONUS-2016 calculator. Four, CREp. Five, the ICE-D database, although as noted above that is just using the v3 code as the back end.

I think this is great.

I should note that there is not universal agreement on this being great. Some members of the CRONUS project envisioned something very different — that an eventual capstone result of the overall project would be a single online calculator that was used by all researchers everywhere as a de facto universal standard for exposure-age calculations. That is the vision described in this paper, which suggests the appointment of an international advisory committee to oversee a single online calculator and issue recommendations for best practices. I think this is a bad idea. In my view, the main value of the CRONUS project is that it generated a really large quantity of valuable information that is relevant to exposure-age calculations and production rate scaling. I think the best possible outcome of generating all this data is that a wide variety of people, whether associated with the original project or not, use it in whatever way they think is best to do whatever research they want to do, and to improve the state of the overall science. Ideally, instead of one committee-approved online calculator, we should have a whole bunch of online calculators that can compete based on their own merits. And that’s what seems to be happening. I think that’s really good.

However, this situation does create a couple of new problems.

One is more general. The question is, how are these things actually going to compete on their own merits? At present, certainly they can compete on the basis of ease of use and general attractiveness. However, really we would like for them to compete on the basis of some sort of quantitative performance metrics that reflect their speed, accuracy, and usefulness for whatever the needed application is. For this to happen, there needs to be some way to easily and transparently quantify the relative performance of different options. This is one important potential function of the ICE-D calibration database: to make the data needed for benchmarking and performance assessment of calculation methods in an easy and straightforward way. Actually making that process straightforward and transparent is more complicated, but at least that is a start.

The second problem is more specific. If I need to calculate some exposure ages right now, what do I do? The answer to this is fairly simple. It doesn’t matter very much which of the four options you use. There are only two things you need to remember. One, don’t use the neutron-monitor-based scaling schemes. Use the Lal- or Sato-based ones. These both work and are effectively equivalent for nearly all practical purposes. Two, which production rate calibration data set you use is much more important than which exposure age calculator you use. If you use the different calculators with the same calibration data, you should get indistinguishable (i.e., differing at less than measurement uncertainty) results in nearly all cases. Much more important, as always, is to completely record all the raw observations needed to compute the exposure ages in your papers, so that readers can recalculate the exposure ages themselves with different methods or calibration data as needed.

In a previous post I described the ICE-D: Antarctica project to collate exposure-age data from Antarctica. This post is about something similar that I am working on with Pierre-Henri Blard  and his colleagues at CRPG. Basically, we’ve now done the same thing with production rate calibration data, which is more useful and more important for several reasons that I will now enumerate in too much detail. If you can’t wait that long, here is the link:

ICE-D: Production rate calibration data

First, some review of what is going on here. To compute exposure ages from cosmogenic-nuclide concentrations, we need to know what the production rate of that nuclide is. We determine this with two steps. We start with some sort of a “scaling scheme,” which is just a model of how production rates vary with location, elevation, and time. Then we need to fit that scaling scheme to a “production rate calibration data set,” which is a set of measurements of cosmogenic-nuclide concentrations in rock surfaces whose exposure ages we already know from some sort of independent evidence. So we measure the average production rate during whatever the exposure time was at some sites whose exposure age we already know, and then we use a scaling scheme to combine data from multiple locations and then apply the results to sites whose exposure age we would like to measure.

The other thing that falls out of this process is the ability to evaluate how well a scaling scheme works. If we can do a good job of fitting the scaling scheme to a set of real calibration data from different locations, elevations, and ages, then we conclude that the scaling scheme is doing a good job. Consult a paper by Brian Borchers and others to see how this works. Basically, this is how we decide which scaling scheme we should be using if we want to compute exposure ages in the most accurate way possible.

Obviously, production-rate calibration data are quite central to this whole process. Another thing that is important is that this whole issue just became a lot more complicated. Expect this to become more blog fodder in the future, but the summary is that the recently completed CRONUS-Earth project, as well as some other efforts, have resulted in a proliferation of scaling methods, online exposure age calculation schemes, and calibration data. In addition, the rate that new production-rate calibration data get generated is reasonably high — a handful of papers describing new calibration data appear each year — which means that the total data set is now rather larger than what was included in the main calibration-and-evaluation project that was done a couple of years ago as part of CRONUS and is described in the paper by Borchers and others. And it’s growing.

So this creates a couple of problems.

One is just a data assimilation problem. What we really want is for the production rate calibration process to be, basically, self-updating — when new calibration data are generated, they are incorporated into our current best estimate of production rates. We, in large part meaning me, have done a terrible job of this in the past few years; there really weren’t any systematic updates of how well scaling models fit the entire set of calibration between the 2008 paper by myself and others and the 2015 Borchers paper. That’s embarrassing, because it was clear very shortly after the 2008 paper was published that the calibration data set in that paper isn’t very good. The data assimilation problem isn’t particularly computationally difficult; we know how to do this. However, there are some obstacles, the first one being the fact that the calibration data needs to be all in one place so that software can easily access the current data set, whatever it may be, in a fast and low-hassle fashion.

The other is the issue of how to decide which of the various scaling and age-calculation options to use, which of course is wound up in the data-assimilation problem because the correct answer will likely change as new data appear. A paper by Fred Phillips and others, that summarizes some of the where-do-we-go-next discussions at the end of the CRONUS project, envisions the existence of an ongoing international committee that will decide what the best way to compute exposure ages is, and prescribe this method to folks who want to do exposure-age calculations. The idea seems to be that this committee would evaluate various proposed scaling schemes and calibration data and determine which were acceptable, which were not, and which one is best.  I’m not a co-author of this paper, and I think this is a bad idea. It creates disincentives for folks who are not on the committee to engage with trying to make the overall field of exposure-dating better, and incentives for people on the committee (who will, of course, be selected because they are responsible for the present state of the art, whatever it is) to maintain the status quo. In my view, this approach would impede progress. As an Economist subscriber who lives near Silicon Valley, of course, I think this is a software problem and not a governance problem: if you create software tools to make it easy to evaluate calibration data and scaling schemes, then people can figure things out for themselves without any help from a committee, and the best-performing calculation methods will float to the top because they are the best-performing. Of course, I don’t have this software yet. But regardless, the first thing that needs to happen here is to make it possible for this software, or anything else, to easily get the entire set of existing calibration data, whatever it is and whenever you want it.

So that’s the initial problem that a calibration data database needs to solve: putting all the calibration data we know about in one place and delivering it to anyone or anything who wants to use it. At present, this problem is not solved. What we have currently is the usual terrible situation in geochronology where various researchers are all maintaining their own mutually inconsistent Excel spreadsheets that each include some fraction of the existing calibration data, are variably up-to-date, and contain different errors and omissions. This situation, of course, maximizes confusion, redundant work, potential for error, and general hassle; and minimizes accuracy, repeatability, and transparency. We can do better. Specifically, the ICE-D calibration data project aims to replace the inconsistent-spreadsheet situation with an online database of all known production rate calibration data that has the following properties: (1) it is generally believed to be reasonably complete and accurate; (2) it only exists in one place so that there are not multiple inconsistent versions; and (3) it can easily be ingested into any software that wants to use these data to do calculations, whether elaborate online exposure age calculator frameworks or odd bits of MATLAB code running on one’s local machine.

So here are the details:

What. The ICE-D: Production Rate Calibration database project.

Where. It’s here: http://calibration.ice-d.org.

Who. I and Pierre-Henri Blard of CRPG/CNRS are at present collaborating on putting it together. If you are involved with collecting production rate calibration data, or you think some data are missing or wrong, you should help too. That does require some knowledge of relational databases and MySQL. Give one of us a call.

Disclaimer. This project is not part of the CRONUS-Earth project.

What, in more detail. As with the ICE-D:Antarctica database, data live in a MySQL database hosted on the Google Cloud SQL service. There is a  front end running in Python on Google App Engine. At present, the front end provides various browser interfaces (by location, publication, etc.) to look at data associated with individual samples or sites. It looks very similar to the Antarctica one. Some interesting features are as follows:

Nonprescriptiveness. The database is organized such that samples can be grouped into “calibration data sets.” For example, the database contains many beryllium-10 measurements from calibration sites that were not included in the CRONUS calibration exercise described in the Borchers and others paper. However, it’s possible to access just the samples that were used in that study as a distinct “calibration data set.” The idea here is to make it possible to replicate previous calibration studies using the data that were used in that study, even as the overall data set grows. It’s also to make the project non-prescriptive: the database should contain all data that can plausibly be described as decent calibration data, but you should be able to decide which data you want to use.

No more spreadsheets. What we want here is for any software to be able to ingest the most up-to-date calibration data from anywhere. This is the fun part. For example, say I am running MATLAB on my laptop and I want to ingest some calibration data to do some kind of calculations. I can get the current up-to-date data for the “Primary Be-10 calibration data set” described in the Borchers et al. paper noted above by using the ‘urlread’ function of MATLAB as follows:
 urls = ['http://calibration.ice-d.org/cds/4']; s = urlread(urls); l1 = '<!-- begin v3 --><pre>'; l2 = '</pre><!-- end v3 -->'; ss = s((strfind(s,l1)+length(l1):strfind(s,l2)-1)); 
What you just did was to read the entire calibration data set (102 measurements) in online exposure age calculator v3 input format into a string variable. The following then parses it into a useful data structure:
in = validate_v3_input(ss); 
That uses the ‘validate_v3_input’ function from the version 3 online exposure age calculator, which isn’t in really great shape yet so it’s not posted anywhere, but if you want a copy let me know. The point is, you don’t need the spreadsheet any more. You just need an internet connection and five lines of code.

Where is it going?  A couple of improvements and applications are really going to happen. First, Pierre-Henri and  his colleagues at CRPG are working on a MATLAB-based online exposure age calculator that will use this database to update calibration data as needed. This is in testing and it is happening. Second, the database is reasonably complete right now for Be-10, Al-26, and He-3 calibration data, but nothing else. In-situ-produced carbon-14 in quartz will happen — that is just a matter of data entry. Including chlorine-36 data is part of the plan, but is more difficult simply because a lot more data need to be recorded. That’s more speculative, but there’s been a bit of progress. Other items on the wish list include making it look less like a static database and more like Wikipedia, so more people can contribute. That’s certainly feasible but a lot more programming work; we’ll need help for that. Finally, we’ll need means to use these data in whatever computational environments people are using for exposure-age calculations. The above example shows how it works for MATLAB, but that could be a lot smoother and, obviously, it would be useful to do the same thing for stuff like Python, R, and, yes, Excel if you really must.

Summary. You don’t need to keep the incomprehensible spreadsheet updated any more. Done with spreadsheets.

This post is about elevation/atmospheric pressure models (or, “atmosphere models”), and which one should be used for exposure-dating purposes. These are important because the reason cosmogenic-nuclide production rates change with elevation is that the air pressure changes with elevation; at higher elevation there is less atmosphere between you and the extraterrestrial cosmic-ray flux. Thus, the parameter that is actually controlling the production rate is not the elevation, but the atmospheric pressure. The difficulty this presents is that while it is relatively easy to determine the elevation of an exposure-dating sample site by direct measurement, it is not at all obvious how one might directly measure the mean atmospheric pressure at that site over the entire exposure history of the sample. Mostly we deal with this by (i) measuring the elevation, and then (ii) converting that to a mean atmospheric pressure by reference to some sort of an atmosphere model that relates elevation to atmospheric pressure and is based on the characteristics of the modern atmosphere.

Early on in the development of production rate scaling schemes (like in the ’90’s), we generally accomplished step (ii) using the ICAO standard atmosphere model. This consists of a single formula relating elevation above sea level to atmospheric pressure that is intended to be a decent approximation for the part of the Earth’s atmosphere lying in temperate latitudes commonly transited by aircraft. Unfortunately for exposure-dating use, although the standard atmosphere is reasonably accurate in an average sense, it fails to capture fairly large spatial variations in the global atmospheric pressure distribution. The following plots show this by comparing the standard atmosphere model to measured mean barometric pressures at various places on Earth. The data set for comparison here is this, which is a set of long-term means of various climate data measured at weather stations globally.

This first figure shows observed mean barometric pressures as a function of elevation from the data set of station measurements (red dots) compared with the standard atmosphere approximation (black line). As noted, the standard atmosphere does pretty well on average, but there is quite a lot of scatter. In particular, as noted some time ago in this paper among others, it does a terrible job in the deep southern hemisphere and Antarctica. To see how terrible, here is the above plot showing only station data south of 60° S.

The elevation-pressure relationship is significantly and systematically different in Antarctica, and the standard atmosphere approximation doesn’t work at all.

The important question, of course, is how important is the scatter from the perspective of production rate estimation? To attempt to answer this, try the following. For all the weather stations in the global data set, calculate the production rate scaling factor for the actual measured mean atmospheric pressure at each station. Here and henceforth in this post I will do this with the ‘St’ scaling scheme in the online exposure age calculators. Then, use the standard atmosphere approximation to estimate the atmospheric pressure at the site of each station, and calculate the production rate scaling factor based on this estimated pressure. Then compare the production rate computed using the actual observed atmospheric pressure to that calculated using the estimated atmospheric pressure. For the standard atmosphere and just considering stations north of 60° S, this yields the following:

What is being plotted here is the production rate calculated for the pressure estimated from site elevations and the standard atmosphere approximation, divided by the production rate calculated for the actual measured pressure. A value above 1 means the standard atmosphere overestimates the production rate relative to measured mean pressure; a value below 1 means an underestimate. Some of the few rather extreme values in this plot (like 40% errors) most likely have to do with errors in the weather station data set that I did not remove by a quick screening, so we won’t get too worried about them. However, if we disregard the worst 2% of the data as outliers most likely due to errors, exclude data south of 60° S, and take the mean and standard deviation of the remaining results we find that using the standard atmosphere approximation even outside the high southern latitudes systematically overestimates production rates across the board by about 1 % (with a standard deviation of 3.2%). In addition, as evident in the above plot there is a significant systematic drift toward larger overestimates at higher elevations. So the summary here is that if we just apply the standard atmosphere approximation to estimate production rates at arbitrary locations globally, we will incur a systematic error at high elevations and an uncertainty in the production rate estimate of at least 3%, depending on how you interpret the distribution of the residuals.

In the paper describing the 2008 online exposure age calculators, we showed that you could do better than this by using a spatially variable atmosphere model. We used the same formula relating pressure to elevation as used in the standard atmosphere approximation, but made two input parameters to the formula — sea level pressure and temperature — spatially variable. The MATLAB function that implements this method (NCEPatm_2.m — downloadable from here) looks up long-term average sea level temperature and pressure values computed by the NCEP reanalysis  (specifically, the 2008 version thereof) at one’s sample location, and uses them as input parameters to the standard atmosphere equation. The equivalent plot as above for this method is as follows:

Note that the y-axis scale has changed. The red circles are production rate residuals using the standard atmosphere as plotted above in the previous figure, and the blue dots use the 2-spatially-variable-parameter model based on the NCEP reanalysis. The residuals are smaller across the elevation range and the overall systematic offset as well as the systematic drift at high elevations is corrected. Note that this figure is the same as Figure 1 in the paper describing the online exposure age calculators linked above (although the y-axis is inverted here). Much of the apparent improvement visible here, however, is just from bringing in the serious outliers — the standard deviation of the residuals for the entire data set, calculated in the same way as described above, is 2.1% in contrast to 3.2%, which is a bit better but not all that much different. So if we apply this atmosphere model to arbitrary exposure-dating sites globally, we get rid of the systematic overestimate and the uncertainty in production rates attributable to this source drops to about 2%. This is pretty good,  but note that this is still a significant contributor to production rate uncertainties. Scatter in globally distributed production rate calibration data sets based on geologic calibration data is typically in the 6-8% range; uncertainty in air pressure estimation could account for a quarter to a third of that total scatter; more if you consider the issue of past changes in the atmospheric pressure distribution.

Now to get to the actual point of this post. A recent paper by Nat Lifton, as well as the recent CRONUS-Earth production rate calibration project described in a paper by Brian Borchers, used a different atmosphere model developed by Nat Lifton. It is similar to the 2-spatially-variable-parameter model based on the NCEP reanalysis used in the 2008 online exposure age calculators and described above, but uses data from the more recent ERA40 reanalysis and adds latitude-dependent variation in a third parameter related to the lapse rate. The MATLAB code is available as part of the supplement to the Lifton paper; I’ll try to add a link here in future. The question here is how this compares to the 2008 NCEP-based scheme. To look at this, replicate the above figure such that instead of comparing the standard atmosphere to the spatially variable NCEP scheme, we are comparing the 2008 NCEP scheme to the 2014 Lifton/ERA-40 scheme.

Here the blue points are the same as the blue points in the previous figure (using 2008 NCEP model) and the green points use the 2014 ERA40 model. There are some differences, but the two models agree very well in total performance, both having insignificant systematic offset and 2.1% scatter. Thus, although to my knowledge it is generally believed that the ERA40 reanalysis is more accurate in a climatological sense than the NCEP reanalysis, this does not appear to make a significant difference from the perspective of global uncertainty in production rate estimates contributed by air pressure estimation, at least given the data available to test this.

Note that there are some regional differences between the two atmosphere models. The following figure shows sites at which the production rate based on the ERA40 model is less than that based on the NCEP model in blue and sites where the opposite is true (P(NCEP) < P(ERA40)) in red.

These differences are actually quite small (the standard deviation of P(NCEP)/P(ERA40) is 0.4% and the largest symbols in the plot above correspond to a difference of about 2%), but it is possible that although there appears to be no significant difference in the performance of the two atmosphere models when evaluated based on the entire data set, there may be some percent-level differences in performance in certain regions. I haven’t investigated this in much detail, but based on the map above it might be worth looking into for some mountainous regions, specifically the central Andes, Himalayas, and central Asia.

One more thing for completeness. Here is the relative performance of various atmosphere models for the deep southern hemisphere and Antarctica (all stations south of 60° S).

The red data are for the standard atmosphere, which we already know does a terrible job here. Blue data, spatially variable NCEP-based scheme; green data, spatially variable ERA40-based scheme. The black data are for the single-formula Antarctic atmosphere approximation described in the paper by John Stone referenced above (the antatm‘ formula in the online exposure age calculators, which is actually based on this paper). Of the methods that don’t do a terrible job (antatm, ERA40, NCEP), none display a significant systematic offset, and the scatter in the residuals (considering all data; no effort to remove outliers) is 3.5% (antatm), 4.6% (NCEP), and 4.0% (ERA40). So it appears that the best air pressure approximation for Antarctica continues to be the simple single-formula antatm’ approximation, although both spatially variable schemes perform nearly as well.  Note that although this scatter for Antarctic sites appears to be significantly larger than for the data north of 60° S, there are not very many station measurements for Antarctica, so this could be the result of only a couple of errors in the data set. Don’t take those scatter estimates too seriously without investigating the input data in more detail.

Summary:

Even without considering changes in the atmospheric pressure distribution over the entire exposure history of exposure-dating samples, uncertainties in approximating mean atmospheric pressure from sample elevations probably contribute at least 2% uncertainty to production rate estimates, which is a significant fraction of the scatter evident in fits to geological calibration data.

All atmosphere models that are more sophisticated than the standard atmosphere (NCEP-based scheme from Balco et al., 2008; ERA40-based scheme from Lifton et al., 2014; simple Antarctic model from Stone, 2000) display similar performance when tested against a global data set of measured mean station pressures (obviously, I am only testing the Antarctic model against Antarctic stations as described above). There appear to be some significant differences between the ERA40 and NCEP-based models for certain regions, but I haven’t looked into that in enough detail to determine which is better.

One more thing:

If the NCEP-based scheme of Balco et al. (2008) and the ERA40-based scheme of Lifton et al. (2014) perform basically the same, which one to use? As discussed above, there is a possibility that ERA40 is better in some mountain areas. This would be worth looking into. The other important thing is, which one is faster? In MATLAB R2014b, the ERA40 model (ERA40atm.m written by Nat Lifton) is faster than NCEPatm_2.m by approximately a factor of 2, the reason for which is unclear but seems to have something to do with differences in interpolation schemes between versions of MATLAB. So it might not be faster in Octave; that would be worth looking into.  In any case, there is no reason not to use ERA40atm.m henceforth.