Elevation/atmospheric pressure models
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.
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.