3. MAIN AGE/PRODUCTION RATE CALCULATIONS

3.A. Governing equations

The general equation relating exposure age to nuclide concentration is: 

We measure the nuclide concentration N. We then solve either for the exposure age t (if we are calculating an exposure age) or the production rate P (if we are estimating the production rate at a site where we already know t). 

The spallogenic production rate at the site is: 

Except that in LSDn case, solves for a correction factor

3.B. Measurement standardization

Measurements are standardized in various ways. Explain. 

3.C. Atmosphere model

Uses Nat Lifton's atmosphere model derived from the ERA40 reanalysis. It is not clear whether this is better for production rate estimation than the NCEP reanalysis used in 2008, but from the atmospheric-science perspective it appears that ERA40 is supposed to be generally better overall. See discussion here:

https://cosmognosis.wordpress.com/2015/10/16/elevationatmospheric-pressure-models/

Stone/Radok Antarctic atmosphere model also available. This is the same as 2008 code. 

3.D. Other correction factors

Thickness. Same as 2008 code.

Shielding. As in 2008 code, accepts a precomputed shielding factor. 

3.E. Production by muons

3.E.1. For exposure-age and production-rate calculations. 

Uses a highly simplified method for surface production by muons that is described in a not-yet-published paper that is posted here:

https://cosmognosis.wordpress.com/2016/11/04/its-all-about-mu/

Download the paper and refer to Figures 15 and 16 and the accompanying text. 

Update: As of Feb 2017, this paper is now published in Quaternary Geochronology. 

Parameters for 10Be, 26Al, and 14C are based on fitting to the Beacon Heights core as described in that paper, as described in the text associated with those figures. 

Parameters for 21Ne cross-sections are based on cross-section estimates by Fernandez-Mosquera et al. This could maybe be improved, but for exposure-dating it doesn't really matter. 

3He. Production by muons assumed zero. Unclear whether this is OK or not. 

3.E.1. For erosion rate calculations. 

Nothing here yet.

3.F. Spallogenic production

Overview:  There are three scaling methods: St, Lm, LSDn. 

3.F.1. Depth dependence for spallogenic production

Single exponential. For exposure-dating, constant e-folding length. Show calculations indicating that latitudinal variability in the e-folding length does not cause errors above reasonable threshold. 

3.F.2. Non-time-dependent scaling for spallogenic production

"St." From Stone et al., 2000, which is based on Lal (1991). The actual scaling factor computation is the same as in the 2008 code. Of course, the reference production rate associated with this scaling method for a given calibration data set will be different for v2 and v3 because of the different muon calculation method and the different atmosphere model. 

3.F.3. Time-dependent scaling for spallogenic production

Lm, which is similar to 2008, and LSDn. 

These are fast implementations that use a look-up table of precalculated scaling factors for the forward integrations. 

3.F.3.a. Magnetic field reconstruction

The LSDn and Lm scaling methods described below require a paleomagnetic field reconstruction. The magnetic field reconstruction is based on Lifton (2016) and is intended to be the simplest possible field reconstruction that is not known to be wrong or oversimplified. This paper shows that:

1. A Holocene and late-glacial field reconstruction ("SHA.DIF.14K") by Pavon-Carrasco et al. (2014), which is simpler than an alternative reconstruction ("CALS3k.4/CALS10k.1p") by Korte and Constable (various references), performs better at matching production rate calibration data, and 

2. Geocentric (not axial) dipole (GD) approximations to the full field reconstructions do not degrade agreement with calibration data. 

Thus, the simplest Holocene field reconstruction that does not appear to be either incorrect or oversimplified is the GD approximation to the SHA.DIF.14K field reconstruction. The magnetic field reconstruction used here, therefore, consists of:

1. GD approximation to the DGRF field reconstructions between 1950-2010 

2. GD approximation to the SHA.DIF.14k field reconstruction period between 1950 and 14,000 years before 1950.  

3. GAD (Geocentric axial dipole) with field strength from the GLOPIS-75/PADM2M ("GLOPAD") record of Lifton et al. (2016) and Lifton et al. (2014) between 14,000 years before 1950 and 2 million years ago. Lifton (2014) (I think) discusses the various advantages and disadvantages of this part of the field reconstruction. 

4. A GAD with the mean field strength of the GLOPAD record for any time prior to 2 Ma. 

Thus, for any time in the past the position of the sample in the magnetic field is represented by a geomagnetic latitude referred to the pole position of the reconstructed GD at that time. For any time before 14 ka, of course, this is the same as the geographic latitude of the sample. The geomagnetic cutoff rigidity needed for input to the production rate scaling method is then computed from the geomagnetic latitude and the dipole field strenght using Equation 2 of Lifton et al. (2014). 

Time steps for forward integration defined in field reconstructions

3.F.3.b. LSDn scaling

Fast-LSDn uses precalculated LUT/interpolation with linearization in S. 

Explanation of correction factor vs. reference production rate, again

3.F.3.c. Lm scaling

Reparameterized in Rc using Lifton equation. Otherwise the same as 2008 code. 

Uses LUT-interpolation for forward integration, which is the same as for LSDn except that there is no solar parameter, so can look up scaling factor directly. 

3.F.3.d. Ag scaling. 

Not functional yet. Some confusion as to how exactly cutoff rigidity is defined for this implementation in relation to the Lifton/Sato code. 

3.F.3.e. Forward integration and age estimate.

Time steps

Get vectors of Rc and S

Compute corresponding SF using LUT and linearization in S (for LSDn) 

Apply P0 and other adjustments. 

Piecewise-constant integration formula

Interpolation to derive age (when calculating exposure age)

Precision test for interpolation scheme WRT complete code

Speed gain results

3.G. Uncertainty estimates.

Linearized, approximated using simple age equation and an effective production rate, internal and external uncertainties. Basically the same as the 2008 code. 

Similar linearization for production rate uncertainties. However, those are not used in most cases. 

3.H.  Production rate calibration. 

Calibration data are sourced from the ICE-D calibration database. Not dynamically (yet). But that is coming. 

10Be and 26Al in quartz and 3He in px/ol are complete and stable at this point and can be used for production. They are calibrated from the CRONUS-Earth primary calibration data sets:

Be-10

Al-26

He-3

There are some differences in averaging/data fitting between here and the Borchers paper that have a substantive effect. Explain why this matters. 

Show calibration results. 

Discuss uncertainty estimates. Bootstrapped from scatter with respect to larger calibration data sets. 

Show this in figures. 

Compare to Borchers MSWD scheme. 

14C, 21Ne, and 3He in quartz are developmental and not yet stable. 

14C is calibrated from some measurements of the CRONUS-A sample (saturated) by Brent Goehring in the now-defunct Tulane lab. This needs work. It is also not integrated with the ICE-D database. At present, I recommend supplying your own calibration for 14C calculations. 

Update: 21Ne in quartz is now calibrated using Cassie Fenton's 'SPICE' data. There is a blog entry about this somewhere. This should be stable now. 

3He in quartz is based on Vermeesch artificial-target data, separately calibrated for each scaling scheme. However, this is not completely self-consistently implemented yet; needs work. Fortunately as a practical matter the precision of this calibration is pretty much irrelevant.