CosmoCalc: an Excel AddIn for cosmogenic nuclide calculations
Abstract As dating methods using Terrestrial Cosmogenic Nuclides (TCN) become more popular, the need arises for a generalpurpose and easytouse data reduction software. The CosmoCalc Excel addin calculates TCN production rate scaling factors (using Lal, Stone, Dunai and Desilets methods); topographic, snow and selfshielding factors; exposure ages, erosion rates and burial ages; and visualizes the results on bananastyle plots. It uses an internally consistent TCN production equation that is based on the quadruple exponential approach of [Granger and Smith, 2000, NIMB (172) pp. 824828]. CosmoCalc was designed to be as userfriendly as possible. Although the userinterface is extremely simple, the program is also very flexible, and nearly all default parameter values can be changed. To facilitate the comparison of different scaling factors, a set of converter tools is provided, allowing the user to easily convert cutoff rigidities to magnetic inclinations, elevations to atmospheric depths and so forth. Because it is important to use a consistent set of scaling factors for the sample measurements and the production rate calibration sites, CosmoCalc defines the production rates implicitly, as a function of the original TCN concentrations of the calibration site. The program is best suited for ^{10}Be, ^{26}Al, ^{3}He and ^{21}Ne calculations, although basic functionality for ^{36}Cl and ^{14}C is also provided. CosmoCalc can be downloaded along with a set of test data from http://cosmocalc.googlepages.com 1 IntroductionThe first applications of Terrestrial Cosmogenic Nuclide (TCN) geochronology
appeared about 20 years ago (Kurz, 1986; Nishiizumi et al., 1986; Phillips et
al., 1986). The method has rapidly developed since those early days, truely
revolutionizing geomorphology and related fields in the process. TCN dating is no
longer a specialized tool used by a small group of experienced users, but has found an
ever growing base of users who are not necessarily familiar with all the details of the
method. Today we are facing a paradoxal situation. On the one hand, a better
understanding of cosmogenic nuclide production systematics has improved the
accuracy of TCN dating. But on the other hand, many users of the method may be
less familiar with its intricacies than was the case in the pioneering days. An
important example of this situation is that of the production rate scaling
factors. In a landmark paper, Lal (1991) presented a method to calculate
cosmogenic nuclide production rates as a function of latitude and elevation. Lal’s
scaling factors are elegant and easy to use, but overestimate the importance of
muons and are only valid for standard atmosphere. Later authors introduced
several improvements, incorporating atmospheric effects and improved muon
production systematics. The scaling factors of Stone (2000), Dunai (2000),
and Desilets et al. (2003, 2006) more accurately represent the scaling of
cosmogenic nuclide production rates with latitude and elevation, but the
increased sophistication of these methods is an obstacle to their widespread use.
CosmoCalc is an addin to MSExcel developed with the intention to alleviate this
problem. The CosmoCalc interface was designed to be as user friendly and
easytouse as possible. Default parameters are set so that beginning users
only have to make a minimal number of decisions. At the same time, all the
default parameters can be changed so that CosmoCalc is highly customizable
and also experienced users should find it useful. The program as well as a
spreadsheet with test data can be downloaded from the CosmoCalc website
(http://cosmocalc.googlepages.com), which also provides detailed installation
instructions. The addin requires MSExcel 2000 or higher. Because of small
differences between the MSWindows and Apple OSX implementations of Excel, two
versions of CosmoCalc are provided. The functionality of both programs is the
same, but the Macintosh version is significantly slower than the PCversion.
After installing the CosmoCalc addin, a toolbar menu appears (Figure 1) that guides the user through the data reduction and closely follows the outline of this paper:
The following sections will provide more details about these calculations. Thus, the present paper serves as an abridged review of TCN calculations, with an emphasis on the numerical methods that are needed to solve the equations. More details about the physics of TCN production are given in the review article of Gosse and Philips (2001). CosmoCalc is not the first computational tool for TCN calculations. Useful alternatives are CHLOE, an Excel spreadsheet for cosmogenic ^{36}Cl calculations (Phillips and Plummer, 1996) and the CRONUSEarth webcalculator (Balco and Stone, 2007; http://hess.ess.washington.edu/math). CosmoCalc was developed independently from these other tools, except for its topographic shielding correction function, which was translated into VBA from the Matlab code of Balco and Stone (2007). The reader is strongly encouraged to try these other programs. CosmoCalc is optimized for ^{26}Al, ^{10}Be, ^{21}Ne and ^{3}He dating. Because geomagnetic field fluctuations and thermal neutron reactions are ignored, results for ^{36}Cl and ^{14}C may be inaccurate. CosmoCalc can be used as an exploratory tool for these nuclides, but for more accurate results, CHLOE or the spreadsheet of Lifton et al. (2005) are recommended.
2 Production rate scaling factorsIn the simplest case (no shielding or burial), only three pieces of information are
needed to calculate an exposure age or erosion rate: TCN concentration, halflife and
production rate. Production rates are the “achilles heel” of the TCN method.
There exist only a few calibration sites where TCN production rates are
accurately known thanks to the availability of independent age constraints
(e.g., Nishiizumi et al., 1989; Niedermann et al., 1994; Kubik et al., 1998).
These production rates are only valid for the specific conditions (latitude,
elevation, age) of each particular calibration site. To apply the TCN method
to other field settings, the production rates must be scaled to a common
reference at sea level and high latitude (SLHL). Up to 20% uncertainty is
associated with this scaling, constituting the bulk of TCN age uncertainty.
Although several efforts have been made to directly measure TCN production
rate scaling with latitude and elevation using artificial H_{2}O and SiO_{2} targets
(Nishiizumi et al., 1996; Brown et al., 2000; Graham et al., 2007), all currently used
scaling models are based on neutron monitor surveys. The oldest and still most
widely used scaling model is that of Lal (1991). This model is a simple set
of polynomial equations giving the (spallogenic + muogenic) production
rate relative to SLHL as a function of geographic latitude and elevation. In
CosmoCalc, Lal’s scaling factors can be calculated by simply selecting two columns
of latitude (in degrees) and elevation (in meters) data and clicking “OK”.
Lal’s scaling factors use elevation as a proxy for atmospheric depth, assuming a
standard atmosphere approximation. Stone (2000) noted that this approximation is
not valid in certain areas, such as Antarctica and Iceland. To avoid the systematic
errors caused by the standard atmosphere model, Stone (2000) recast the polynomial
equations of Lal (1991) in terms of air pressure instead of elevation. A second
improvement of the Stone (2000) model is the independent scaling of TCN
production by slow (negative) muons (Heisinger et al., 2002). In spite of this
added complexity, the CosmoCalc interface for Stone (2000) scaling factors is
identical to that for Lal (1991) scaling: the user simply needs to provide two
columns of data, one with latitude and one with air pressure (in mbar).
The scaling factors of Stone (2000) can be different for different nuclides,
because the importance of muons depends on the nuclide of interest. Because
most TCN production rate calibration sites are not located at SLHL, it is
crucially important to scale the production rates using the same method as the
unknown sample. This is exactly what CosmoCalc does when the user selects a
nuclide from the scrolldown menu of the scalingform. Thus, the program
“forces” the user to be consistent. The program comes with a set of default
calibration sites, but these can be changed. Also the relative importance of the
production pathways (neutrons, slow and fast muons) can be changed (Section 7).
Behind the scaling models of both Lal (1991) and Stone (2000) lies an extensive
database of neutron monitor measurements, ordered according to geomagnetic
latitude. This ordering implies that Earth’s magnetic field can be accurately
approximated by a simple dipole. To avoid this approximation, Dunai (2000) ordered
the neutron monitor data according to geomagnetic inclination, which also represents
the nondipole field. Just like Stone (2000), also Dunai (2000) incorporates
separate muon scaling and atmospheric effects. However, atmospheric depth
(g/cm^{2}) is used instead of air pressure. Using CosmoCalc it is very easy to
convert air pressure to elevation or atmospheric depth and back (Section 6).
Ultimately, both geomagnetic latitude and inclination are merely proxies for a
more fundamental physical quantity: the geomagnetic cutoff rigidity (R_{c}), which is
the minimum momentum per unit charge (in GV), required for a primary cosmic ray
to reach the atmosphere. Ordering the neutron monitor data according to this
parameter results in yet another set of scaling factors. Unfortunately, at least three
different methods for calculating R_{c} exist in the literature. Dunai (2001) used a
database of horizontal magnetic field intensities and inclinations to estimate the
cutoffrigidity of an equivalent axial dipole field, for which an approximate analytical
solution exists. Desilets and Zreda (2003) used a model based on trajectory tracing of
an axial dipole field, which is done by numerically testing the feasibility of vertically
incident antiprotons to travel from the top of the atmosphere back into space.
Finally, Lifton et al. (2005) fit a cosine function to a database of geomagnetic
latitudes versus trajectory traced cutoff rigidities for the 1955 magnetic
field. These authors consider the scatter around their fit to be a realistic
estimator of the natural variability of R_{c} at any given geomagnetic latitude.
In order to avoid confusion, CosmoCalc currently only implements one of these
three methods, namely that of Desilets and Zreda (2003). The scaling model of
Desilets and Zreda (2003) makes a distinction between slow and fast muons, each of
which scales differently. In an update of their model, Desilets et al. (2006) did a
neutron monitor survey at low latitudes, which were undersampled by previous
surveys. This resulted in a slightly different set of attenuation length polynomials for
spallogenic reactions. The method of Desilets et al. (2003, 2006) is probably the most
accurate of all the scaling models implemented in CosmoCalc. It is also the most
complex model, but this did not change the user interface. Using the scaling
models of Desilets et al. (2003, 2006) is just as easy as that of Lal (1991) in
CosmoCalc. Default values for the relative SLHL production rate contributions of
neutrons, slow and fast muons can be changed in the Settings menu (Section 7).
The scaling models of Dunai (2001), Desilets et al. (2003, 2006), Pigati and Lifton (2004), and Lifton et al. (2005) are a sensitive function of magnetic field intensity and solar activity, both of which are poorly constrained over geologic time. On the one hand, this temporal variability is the biggest downside to R_{c}based models for long exposures (> 20 ka; Dunai, 2001). On the other hand, such models also offer the possibility to correct for secular variation of TCN production rates for short exposures. Instructions for doing so are given by Dunai (2001) and Desilets et al. (2003), provided a local record of paleomagnetic intensity is available. Compiling such a record is something for advanced users and falls beyond the scope of CosmoCalc. The scaling models of Pigati and Lifton (2004) and Lifton (2005) are accompanied by global datasets of magnetic field intensity, polar wander and solar activity and in principle, it would be possible to incorporate these datasets into CosmoCalc. However, because they are very large (4 and 7Mb, respectively) in comparison with CosmoCalc (~ 500kb), the cost of including these scaling models was considered too high. Therefore, researchers working with ^{14}C, where secular variation of the magnetic field is really crucial, should use CosmoCalc only as an exploratory tool, and use the spreadsheets of Pigati and Lifton (2004) and Lifton et al. (2005) for final calculations.
3 Shielding correctionsThe scaling factors discussed in the previous section allow the calculation of TCN production rates at any location on the Earth’s surface, assuming that the sample is a slab of zero thickness taken from a horizontal planar surface. If these assumption are not fulfilled, the SLHL production rates must be multiplied by a second set of correction factors, quantifying the extent to which the cosmic rays were blocked. CosmoCalc implements three such correction factors: topographic shielding, self shielding and snow cover.
3.1 Topographic shieldingTwo kinds of topographic shielding corrections can be distinguished for (a) samples taken from a tilted rather than horizontal surface, and (b) samples that are located in the vicinity of topographic irregularities. CosmoCalc follows the approach of Balco and Stone (2007) (their Matlab function skyline.m) and treats these two effects together using the following equation:
With h(θ) the “horizon” in the azimuthal direction θ, i.e. either the elevation (in radians) of the topography or the sloping sample surface, whichever is greatest. Sometimes, an exponent of 2.3 is used instead of 3.5 in Equation 1 (Staudacher and Allègre, 1993). CosmoCalc treats this exponent as a variable, which can be changed in the Settings form (Section 7). In practice, the integral of Equation 1 is solved by linear interpolation between a finite number of azimuth/elevation measurements. The input needed by CosmoCalc is two mandatory columns of strike and dip (in degrees, where the strike is 90 degrees less than the direction of the dip), followed by an optional series of topographic azimuth/elevation measurements (in degrees). There is no restriction on the total number of measurements, provided they come in multiples of two.
3.2 SelfshieldingCosmic rays are rapidly attenuated as they travel through matter, causing TCN production rates to vary greatly with depth below the rock/air contact. They must be integrated over the actual sample thickness and scaled to the surface production rates before an exposure age can be calculated. Different reaction mechanisms are associated with different attenuation lengths. Gosse and Philips (2001) consider four kinds of thickness corrections, for spallogenic, thermal and epithermal neutrons, and muons. Because selfshielding corrections are generally small, CosmoCalc considers only the spallogenic neutron reactions:
with Λ_{0} the spallogenic neutron attenuation length (default value 160 g/cm^{2}), ρ the rock density (default value 2.65 g/cm^{3}) and z the sample thickness (in cm). Neglecting the remaining pathways makes little difference, with the possible exception of ^{36}Cl, because the latter can be strongly affected by thermal neutron fluxes, which are currently ignored by CosmoCalc.
3.3 Snow coverPerhaps the most popular and powerful application of TCN techniques is the dating of glacial moraines (e.g., Gosse et al., 1995; Schäfer et al., 1999). These features are generally located at high latitudes or elevations, where snow cover poses a potential problem. The snow correction is similar to the selfshielding correction with the important difference that the former is highly variable with time. Given n (e.g., 12 for monthly or 4 for seasonal) measurements of average snow thickness z and density ρ, CosmoCalc computes the snow correction factor S_{c} as follows:
4 Banana plotsBefore calculating an exposure age or erosion rate, it is a good idea to check if the TCN measurements are consistent with a simple or complex exposure history. This can be done with two nuclides (including at least one radionuclide) using a “banana plot” (Lal, 1991). CosmoCalc accomodates two types of banana plot: ^{26}Al^{10}Be and ^{21}Ne^{10}Be. Depending on whether or not a sample plots above, below or inside the socalled “steadystate erosion island” (Lal, 1991), one can decide whether or not to pursue the calculation of an exposure age, erosion rate or burial age. For the construction of the banana plots and the age/erosion calculations of Section 5, CosmoCalc implements a modified version of the ingrowth equation of Granger and Muzikar (2001):
With N the nuclide concentration (atoms/g), P the total surface production rate
(in atoms/g/yr) at SLHL, τ the burial age, ϵ the erosion rate, t the exposure age and
λ the radioactive halflife of the nuclide. Equation 4 models TCN production by
neutrons, slow and fast muons by a series of exponential approximations. The first
term of the summation models TCN production by spallogenic neutron reactions, the
second and third terms model slow muons and the last term approximates TCN
production by fast muons. Thus, F_{0},...,F_{3} are dimensionless numbers between zero
and one, and Λ_{0},...,Λ_{3} are attenuation lengths (g/cm^{2}). The approach of Granger et
al. (2000, 2001) was chosen because of its flexibility. For instance, neglecting muon
production can be easily implemented by setting F_{1},F_{2} and F_{3} equal to zero in
Equation 4. CosmoCalc uses Granger et al.’s (2000, 2001) recommended
values of F_{0},...,F_{3} for ^{10}Be and ^{26}Al, but also offers an alternative choice
of preset values approximating either the alternative parameterization of
Schaller et al. (2001), neglecting the contribution of muons, or only using three
exponentials (for more details, see Section 7). Banana plots with nonzero
muon contributions feature a characteristic crossover of the steadystate and
zero erosion lines which is absent when muons are neglected (Figure 2).
CosmoCalc’s banana plots are normalized to SLHL, meaning that the TCN concentrations of each sample are divided by the cumulative effect of all their correction factors, represented by the “effective scaling factor” S_{e}:
with
where S_{p} is one of the production rate scaling factors of Section 2 and S_{t}, S_{s} and
S_{c} are defined in Section 3. If muon production is neglected, then S_{e} = S_{t} × S_{p} × S_{s}
× S_{c}. However, in the presence of muons, the effective scaling factor S_{e} may deviate
from this value because the relative importance of the different production
mechanisms changes as a function of age, erosion rate, elevation, latitude, sample
thickness and snow cover. The exact form of the function f(S) will be defined in
Section 5.2.4. Note that the topographic shielding correction S_{t} does not
“fractionate” (i.e., change the fractions F_{0},...,F_{3} of) the different production
mechanisms and is placed outside the scaling function f(S). This means that, strictly
speaking, the TCN concentrations should be multiplied by S_{t} prior to generating a
banana plot. The input required by CosmoCalc’s “Banana” function are (1)
the composite scaling factor S for the first nuclide (^{26}Al or ^{21}Ne), (2) the
concentration and 1σ measurement uncertainty of the first nuclide (^{26}Al or
^{21}Ne), both multiplied by S_{t}, (3) S for the second nuclide (^{10}Be) and (4)
the concentration and 1σ measurement uncertainty of the second nuclide
(^{10}Be), also multiplied by S_{t}. Because topographic shielding corrections are
generally small, the systematic error caused by lumping S_{t} together with
the other correction factors is very small. Therefore, if S_{t} > ~ 0.95, say, it
is safe to approximate Equation 5 by S_{e} = f(S_{t} × S_{p} × S_{s} × S_{c}). In this
case, the nuclide concentrations do not need to be premultiplied by S_{t}.
The graphical output of CosmoCalc can easily be copied and pasted for editing in vector graphics software such as Adobe Illustrator or CorelDraw. The yaxis of the ^{26}Al^{10}Be plot is logarithmic by default whereas the yaxis of the ^{21}Ne^{10}Be plot is linear. These defaults can be changed in the “Banana Options” userform. Note that MSExcel (versions 2000 and 2003) only allows logarithmic tickmarks to have values in multiples of ten. To get around this limitation, CosmoCalc uses a “pseudo yaxis”, which cannot be edited by the usual right mousebuttonclick. Hopefully, this limitation will not be necessary in later versions of Excel. CosmoCalc only propagates the analytical uncertainty of the measured TCN concentrations. No uncertainty is assigned to the production rate scaling factors, radioactive halflives or other potentially illconstrained quantities. On the banana plots, the user is offered the choice between error bars or ellipses with the latter being the default. Banana plots are graphs of the type N_{1}/N_{2} vs. N_{2} which are always associated with some degree of “spurious correlation” (Chayes, 1949). This causes the error ellipses to be rotated according to the following correlation coefficient:
If N_{1} stands for ^{26}Al or ^{21}Ne and N_{2} for ^{10}Be, then N_{1} and N_{2} are the measured concentrations of these respective nuclides while σ_{N1} and σ_{N2} are the corresponding measurement uncertainties. 5 Age/erosion rate calculationsEquation 4 has three unknowns: t (exposure age), ϵ (erosion rate) and τ (burial age). If only one nuclide was measured, we must assume values for two of these quantities in order to solve for the third. If two nuclides were analysed (of which at least one is radioactive), only one assumption is needed. CosmoCalc is capable of both approaches. In this section, we will first discuss how to solve for ϵ (assuming infinite exposure age and zero burial) and t (assuming zero erosion and burial) using a single nuclide (Section 5.1). Then, numerical methods will be presented to simultaneously solve for t and ϵ (assuming zero burial), t and τ (assuming zero erosion) or ϵ and τ (assuming infinite exposure age), using two nuclides (Section 5.2). Note that in the case of two nuclides (^{26}Al or ^{21}Ne combined with ^{10}Be), the assumption of zero burial can be verified on the banana plot.
5.1 Calculations using a single nuclideCosmoCalc requires three pieces of information to calculate an exposure age or erosion rate: the TCN concentration (corrected for topography), its analytical uncertainty and a composite correction factor for production rate scaling with latitude/elevation and shielding (Equation 6). We somehow need to incorporate this scaling factor into the ingrowth equation (Equation 4). This poses a problem because the scaling factor is a single number whereas Equation 4 explicitly makes the distinction between neutrons, slow and fast muons. Granger and Smith (2000) avoid this problem by separately scaling the different production mechanisms:
Instead of one scaling factor, Equation 8 has four, one for neutrons (S_{0}), two for slow muons (S_{1} and S_{2}) and one for fast muons (S_{3}). Granger et al. (2001) separately calculate each of these four scaling factors. Thus, the original method of Granger and Smith (2000) is incompatible with the common practice of lumping all production mechanisms into a single latitude/elevation scaling factor (Section 2). To ensure optimal flexibility and userfriendliness, CosmoCalc uses a slightly different approach. S_{0},...,S_{3} are calculated from the composite correction factor S, by approximating the total scaling by a single attenuation factor caused by a virtual layer of matter of thickness x (in g/cm^{2}):
so that
with F_{i} and Λ_{i} as in Equation 4 and S as defined in Equation 6. CosmoCalc
solves Equation 10 iteratively using Newton’s method. As said before, some assumptions are needed to solve Equation 8. An exposure age (t) can be calculated under the assumption of zero erosion and burial (ϵ = 0, τ = 0). For a radionuclide with decay constant λ, this yields:
whereas for stable nuclides (^{3}He and ^{21}Ne):
Alternatively, the erosion rate (ϵ) can be calculated under the assumption of steady state and zero burial (t = ∞, τ=0):
CosmoCalc solves this equation iteratively using Newton’s method. Statistical uncertainties are estimated by standard error propagation:
These error estimates do not include any uncertainties in production rates and scaling factors, which are difficult to quantify, but can be evaluated by using a range of input parameters.
5.2 Calculations with two nuclidesEquation 8 has three unknowns (t, ϵ and τ). If two nuclides have been measured (with concentrations N_{1} and N_{2}, say), only one value must be assumed in order to solve for the remaining two. By assuming zero erosion (ϵ = 0), CosmoCalc simultaneously calculates the exposure age and burial age (Section 5.2.1); by assuming steadystate erosion (t = ∞), the erosion rate and burial age are calculated; and by assuming zeroburial (τ = 0), the erosion rate and exposure age can be computed (Section 5.2.2).
5.2.1 Burial datingIf a rock surface gets buried by sediments or covered by ice, it is shielded from cosmic rays and the concentration of cosmogenic radionuclides decays with time. Such samples plot outside the steadystate erosion island of the banana plot, in the socalled field of “complex exposure history”, a feature which is considered undesirable by most studies. Other studies, however, intentionally target complex exposure histories, using radionuclides to date preexposure and burial (e.g., Bierman et al., 1999; Fabel et al., 2002; Partridge et al., 2003). CosmoCalc calculates burial ages, either by assuming negligible erosion or steady state erosion (ϵ = 0 or t = ∞, respectively). It does not handle postdepositional nuclide production. Burial  Exposure dating
The easiest case of twonuclide dating is that of simultaneously calculating exposure age (t) and burial age (τ) with one radionuclide and one stable nuclide. Because the stable nuclide is not affected by burial, it can be used to calculate the preexposure age, using Equation 12. This age can then be used to calculate the burial age:
In the case of two radionuclides, CosmoCalc finds t by iteratively solving the following equation using Newton’s method:
With λ_{1} > λ_{2}. The solution is then plugged into Equation 17, using nuclide 1. Burial  Erosion dating
These equations are easy to solve since the variables τ and ϵ are separated. If nuclide 1 has the shortest halflife (largest decay constant λ), the burial age is written as a function of the erosion rate ϵ:
The erosion rate is given implicitly by:
CosmoCalc solves Equation 21 for ϵ using Newton’s method and then plugs this value into Equation 20. If λ_{2} < λ_{1}, then N_{2} is used instead of N_{1} in Equation 20.
5.2.2 Ageerosion rate calculationsAssuming zero burial (τ = 0) yields the following system of equations f_{1} and f_{2}:
It is impossible to solve these equations for exposure age (t) and erosion rate (ϵ) separately. Instead, CosmoCalc implements the twodimensional version of the NewtonRaphson algorithm:
With J(ϵ,t) = the Jacobian matrix, which is also used for error propagation.
5.2.3 Error propagationError propagation is less straightforward in the twodimensional case than in the single nuclide case (Section 5.1). The bijection from (N_{1},N_{2})space to (ϵ,t)space is not orthogonal, particularly in the case of ageerosion dating (Section 5.2.2). For this reason, it is only possible to analytically compute upper bounds for σ(ϵ), σ(τ) and σ(t):
with x and y placeholders for ϵ and t or τ, respectively, and ∥⋅∥ the absolute
values of the matrix . In the case of ageerosion dating, the confidence intervals for
t and ϵ are very wide, often too wide to be useful. Therefore, it can be more
productive to solve each quantity separately instead of simultaneously. Thus, using
equations 11 and 13, it is possible to estimate minimum exposure ages and
maximum erosion rates (e.g., Nishiizumi et al., 1991). However, for burial
dating there is no choice and we must simultaneously solve for τ and ϵ or t.
In addition to Newton’s method, CosmoCalc offers a second way of solving equations 19 and 22 by means of Monte Carlo simulation, implementing the MetropolisHastings algorithm (Metropolis et al., 1953; Tarantola, 2004). The MetropolisHastings algorithm is a socalled Bayesian MCMC (Markov Chain Monte Carlo) method. It not only finds the best solution to the system of nonlinear equations, but actually explores the entire solution space. If the “Metropolis” option of the AgeErosion function is selected, CosmoCalc generates 1000 “acceptable” solutions to Equation 19 or 22, where “acceptable” is defined by the bivariate normal likelihood of the forwardmodeled TCN concentrations (Figure 3). The last 900 of these solutions are then ranked according to their likelihood. For a 95% confidence interval, those solutions with the lowest 5% likelihoods of the 900 results are discarded, leaving 855 values for ϵ, t or τ. The minimum and maximum values of these 855 numbers are the lower and upper bounds, respectively, of the simultaneous 95% confidence intervals. In contrast with the symmetric confidence intervals given by equation 24, the MCMC confidence limits are always greater than or equal to zero. However, as said before, the 95% confidence intervals can be very wide especially in the case of ageerosion dating.
5.2.4 An a posteriori modification of the banana plotsSection 4 discussed the construction of ^{26}Al^{10}Be and ^{21}Ne^{10}Be banana plots. To plot samples from different field locations (with different latitude, elevation and shielding conditions) together on the same banana plot, it is necessary to scale the TCN concentrations to SLHL. In other words, each TCN concentration must be divided by an appropriate scaling factor, the socalled “effective scaling factor” S_{e} (Equation 5):
With N_{meas} the measured TCN concentration, and N_{SLHL} the equivalent TCN concentration which would be measured had the sample been collected from SLHL. In the case of zero erosion, N_{SLHL} = N_{meas} / (S_{p} × S_{t} × S_{s} × S_{c}). This is no longer true when ϵ > 0, because the relative contributions of neutron spallation, slow and fast muons change below the surface. This is the “fractionation” effect that was discussed in Section 4 and quantified by Equation 10. For example, consider the case of two high latitude, high elevation samples, one with negligible erosion (ϵ=0) and one with nonzero erosion (ϵ>0). Because the relative importance of neutron spallation increases with decreasing erosion rate, and neutrons are more important (relative to muons) at higher elevations, S_{e} will be greater for the zeroerosion than for the nonzero erosion case. CosmoCalc first solves Equation 19 or 22 for ϵ and τ or t, whichever fits the measured TCN concentrations best. Plugging these solutions into Equation 4 yields the equivalent TCN concentration at SLHL. S_{e} is then given by Equation 25.
6 ConvertersSection 2 discussed four different models to scale TCN production rates from SLHL
to any other location on Earth. All these models have in common that they require
two columns of data in CosmoCalc: “latitude” and “elevation”. They differ in
how they quantify these two pieces of information. The scaling factors of
Lal (1991) are the only ones that use the actual geographical latitude (in
degrees) and elevation (in meters). Stone (2000) also uses the geographical
latitude for estimating the latitude effect, but uses atmospheric pressure (in
mbar) for modeling the elevation effect. Dunai (2000) uses the geomagnetic
inclination (in degrees) instead of latitude, and atmospheric depth (in g/cm^{2})
instead of elevation. Finally, Desilets et al. (2003, 2006) use cutoff rigidity (in
GV) for the latitude effect and atmospheric depth for the elevation effect.
All these different measures of “latitude” and “elevation” are related to each other and can be converted into each other. To facilitate the comparison of the different methods and, for example, reinterpret published literature data, CosmoCalc provides a series of easytouse conversion tools.
6.1 Converting different measures of “elevation”To convert elevation (z, in m) to atmospheric pressure (p, in mbar) (Iribane and Godson, 1992):
With p_{0} the pressure at sea level, β_{0} the adiabatic lapse rate, T_{0} the temperature at sea level, g_{0} the gravitational constant and R_{d} the universal gas constant. In the standard atmospheric model, β_{0} = 6.5 K/km, g_{0} = 9.80665 m/s^{2}, p_{0} = 1013.25 mbar and T_{0} = 288.15 K. However, these values are not valid for Antarctica, where p_{0} ≈ 989.1 mbar and T_{0} ≈ 250 K. The modified Equation 26 can be rewritten as (Stone, 2000):
Atmospheric pressure is converted to atmospheric depth (g/cm^{2}) by:
The reverse conversions are trivial inversions of these equations.
6.2 Converting different measures of “latitude”Converting latitude (L, in degrees) to geomagnetic inclination (I, in degrees) and back:
Converting latitude to geomagnetic cutoff rigidity (in GV) for a geomagnetic field strength M, compared to the 1945 reference value (M_{0} = 8.085×10^{2} A m^{2}):
The default value for M/M_{0} = 1, but can be changed by clicking the “Option” button of the CosmoCalc conversion form. e_{1},...,e_{6} and f_{1},...,f_{6} are defined in Table 8 of Desilets and Zreda (2003). The reverse operation of Equation 30 does not have an analytical solution and is solved iteratively with Newton’s method.
7 Customizing CosmoCalcThe interface of CosmoCalc is very simple because default values are set for most of the parameters that occur in the various equations discussed in this paper (Table 1). This greatly reduces the chance that novice users make mistakes when reducing their TCN data. For more advanced users, the program allows nearly all the parameters to be changed.
7.1 Specifying the production rate calibration sitesAs mentioned in Section 2, it is very important to use a consistent set of scaling factors for the unknown sample and the production rate calibration site. Failing to do so can cause significant systematic errors. To avoid this, CosmoCalc defines the SLHL production rates implicitly, by specifying a set of calibration sites and their measured TCN concentrations. Using the “Calibration sites” form of the “Settings” menu, these concentrations are scaled to SLHL and an average production rate is calculated using one of the five scaling models of Section 2. CosmoCalc comes with a default set of published production rate calibrations, some of which (^{10}Be and ^{26}Al) were borrowed from Balco and Stone (2007; http://hess.ess.washington.edu/math). The published data come from a variety of latitudes and elevations, yielding a presumably reliable estimate of the globally averaged production rates. This, however, is not always the best approach. For example, if a TCN study is carried out in the vicinity of one particular calibration site, then it makes more sense to use only this site to estimate the local production rate. Therefore, CosmoCalc offers the user the flexibility to delete or add calibration sites at will.
7.2 Changing the relative contributions of different production pathwaysBeing based on the equation of Granger and Smith (2000), the TCN production
equation consists of four exponentials: one for neutrons, two for slow neutrons, and
one for fast neutrons (see Equation 8 and Section 5). These exponentials are governed
by two sets of parameters: the fractions F_{i} and the attenuation lengths Λ_{i} (for
i=0,...,3). Default values for Λ_{i}, F_{i}(^{10}Be) and F_{i}(^{26}Al) were taken from Granger and
Smith (2000) (Table 1), but these values can be changed in the “Settings” form.
In addition to Equation 8, several alternative, but similar looking TCN ingrowth equations exist in the literature. Schaller et al. (2001, 2002) use not four but eight exponentials (two for neutrons, and three for each slow and fast muons), whereas others use three (one for each production mechanism) (e.g., Braucher et al., 2003; Miller et al., 2006) or only one exponential (neglecting muon production). CosmoCalc provides a separate set of default parameters for each of these alternatives. For example, the ingrowth equation of Schaller et al. (2002) was recast in the parameterization of Granger and Smith (2000) by a least squares fit of a virtual depth profile (Figure 4).
8 How to report data reduced with CosmoCalcCosmoCalc was designed to be a userfriendly program for both novice and advanced
users of TCN methods. Nearly every function comes with a set of options which allow
the user to change the default values of various parameters (Section 7). Although
CosmoCalc’s flexibility should be considered a positive feature, a danger exists that
his may lead to confusion. Therefore, it is important that the datareduction
method is well documented when reporting results. Here is an example:
“Ages were calculated with CosmoCalc 1.0 (Vermeesch, 2007), using Dunai
(2000) scaling factors and default values for all parameters with the following
exceptions: ρ = 3.5 g/cm^{3} and λ(^{10}Be) = 5.17×10^{7}yr^{1}.” Clearly indicating the version number is important because CosmoCalc
will be updated in the future to keep track of new developments in TCN
geochronology. Older versions will always be available from the CosmoCalc website
(http://cosmocalc.googlepages.com). Users with useful suggestions for improvements
should feel free contact the author by emailing to cosmocalc@gmail.com.
Acknowledgments. The author is financially supported by a Marie Curie Fellowship of the European Union (CRONUSEU network, RTN project reference 511927). Many thanks to Dr. Naki Akçar for testing betaversions of CosmoCalc and making useful suggestions which simplified the user interface and made the program more userfriendly, and to two anonymous reviewers for constructive comments. References
