ACE : Age Calculation Engine

The recent ³⁶Cl calibration dataset of Licciardi et al 2008 is included with ACE, but this dataset is designed to calibrate Ca production rates, so that if the standard ACE multiple pathway calibration is invoked strange things happen:

As the figure above shows, with the dataset calibrated Ca spallation rates are high (100 atoms / g/ yr), K spallation rates are negative (-1345 atoms / g / yr) and low-energy neutron capture rates are about half their typical values (343 atoms / g /yr).

This happens because the Licciardi calibration dataset was never intended to calibration all three production pathways, just Ca spallation. Therefore in order to use this calibration dataset for the purpose with which it is intended, we need to modify ACE to calibrate only for Ca spallation production rates. A patch is supplied below which provides a new 36Cl workflow, but this post outlines ACE was modified to allow this so that similar changes can be made by users.

Modifying ACE to constrain calibrations

First we convert the ³⁶Cl/³⁷Cl ratios reported by Licciardi et al 2008 to ³⁶Cl/Cl ratios to be input into ACE as clRatio. The procedure is documented here.  Next we create a new workflow called ’36Cl Constrained Calibration’. To do this, in the ACE/repo/workflow directory copy the original workflow (’36Cl Time Stepping Calibration.txt’) to a new file called ’36Cl Licciardi Calibration.txt’. Then open this file and change the last line from:

OutputFor36ClAges => Output36ClCalibration

to

OutputFor36ClAges => OutputConstrainedCalibration

We are telling ACE that instead of the standard calibration procedure, the workflow ’36Cl Licciardi Calibration’  has a new routine (called OutputConstrainedCalibration.py) to calibrate production rates. Now we need to create this new routine. In the folder ACE/src/ACE/Components copy the file ‘Output36ClCalibration.py’ to ‘OutputConstrainedCalibration.py’. As always ACE will find this module from the header information so don’t forget to rename the header instances on lines 43 and 45:

class OutputConstrainedCalibration(Component):
    def __init__(self, collections, workflow):
        super(OutputConstrainedCalibration, self).__init__(collections, workflow)

In the original code the calibration is performed using least squares regression (by  QR decomposition to minimise round-off errors). In particular ACE calibrates observed vs modelled inventories, with the modelled inventory sets up as a coefficient of HLSL production rates.  In this way, the calibration coefficient regressed for is the HLSL production rate. In matrix form, if A is the production rate coefficients, x is the production rates, the modelled inventory is Ax. We know the observed inventory of the calibration samples, and call this B. In the multiple linear regression:

Ax = B

To solve for x is the production rates, we decompose A into matrices Q and R such that A = Q^TR. As a result:

x = R^-1QB

The python code for this procedure is given below:

numData = i # Number of samples in database
numParams = 3 # Number of production mechanisms to constrain
(Q,R) = Sci.linalg.qr(LHS_weighted)
R = R[0:numParams,:];
Q = Q[:,0:numParams];
matmult = dot(transpose(Q),RHS_weighted)
rates = Sci.linalg.solve(R,matmult)

Where LHS_weighted is the weighted coefficients of production rates (A in the above sense) and RHS_weighted is the weighted, observed inventory of the calibration samples minus any contributions from muons (B in the above sense).

Our problem at hand is that the method used in ACe by default simultaneously calibrates production rates for Ca and K spallation together with low-energy neutron capture, whereas the Licciardi et al dataset is designed to contain only information regarding Ca spallation. To correct this we need to change x to be only one production pathway (Ca spallation), and assume some production rates for the other pathways. We  adopt the Licciardi et al values from Table 2, that K spallation is 161 ± 9 atoms / g / yr and low energy neutrons is 626 ± 40 atoms / g / yr. Now the code has to be updated to reflect these changes (fully modified code supplied below, the below is for reference):

psi_k_0  = 161. / 1.54E22
Pf_0     = 626.
psi_k_0_one_sigma = 9 / 1.54E22
Pf_0_one_sigma = 40

As the atomic weight of K is 39.10, there are 1.54E22 atoms per gram of K .

The next step is to modify the matrix of production rate coefficients (A in the above sense) to only represent coefficients for the spallation of Ca:

LHS_weighted[i][0] = Inv_err[i] * LHS[i][0];

As a result we are calibrating Ca spallation rates only against that part of the inventory which was generated by Ca spallation. Therefore in the code we also need to subtract the ³⁶Cl production due to K spallation and low-energy neutron capture from the observed inventory (B in the above sense):

#Remove K and N contributions from Total Inventory
RHS[i] = RHS[i] - LHS[i][1] * self.KConstrained *1E20 / 1.54E22 - LHS[i][2] * self.NConstrained

In addition, the Licciardi production rates contain components for Fe and Ti spallation.  ACE does not consider these, so they must be introduced.  As for K spallation and low-energy neutron capture, they must be removed from the inventory component which is calibrated against (B in the above sense).  As ACE does not track these inventory contributions at each timestep, they must be changed in AddConstrainedCalibrationCoefficients.py, the module where contributions are added:

# Remove inventory contribution from Fe and Ti spallation here
self.FeSpallRate = 1.9 #Spallation rate from Stone (2005) Table 2 doi:10.1016/j.epsl.2007.11.036
self.TiSpallRate = 6 #Spallation rate from Masarik (2002) Table 2 doi:10.1016/j.epsl.2007.11.036
s["RHS_1"] -= mul1 * self.shielding_factor * self.Q_s * self.FeCalc * self.FeSpallRate / norm_factor
s["RHS_1"] -= mul1 * self.shielding_factor * self.Q_s * self.TiCalc * self.TiSpallRate / norm_factor

Results

Now we are ready to calibrate for Ca spallation using this new workflow. Save the changes introduced above, and rerun the calibration using the new code:

When the K spallation and low-energy neutron capture rates are constrained in ACE, the Licciardi et al dataset estimates Ca spallation rates of 68 a / g / yr.  This is about 20% higher than the esimate of Licciardi et al (57 a / g / yr).

There are several reasons why the ACE calculated value for Ca spallation rates could be different to the Licciardi et al calculation:

• The conversion of ³⁶Cl/³⁷Cl ratios reported by Licciardi et al 2008 to ³⁶Cl/Cl could be erroneous
• Scalings could be different
• The calibration technique could be different

Different Calibration Techniques

Using the code above, we can force ACE to use a Ca spallation rate of 57 a / g / yr.  The results are shown below:

With a reduced chi-square of 7.41, higher than for a Ca spallation rate of 68 a / g / yr (4.67 shown above).  From the inventory plots, it would appear that the higher spallation rate generates inventories closer to observed for large inventory, older samples (ie closer to y = x on the calibration curves).  This can be verified by using both experiments to estimate ages of the samples in the calibration database:

By clicking on these samples in the sample browser and selecting Copy in the ACE menu toolbar, they can be pasted into Excel or OpenOffice for further analysis:

sample	experiment	computed	experiment	computed	independent	independent
id	type	age	type	age	age	age uncertainty
IC02-16 4	Constrained	3960	Forced	4500	4040	250
IC02-17 4	Constrained	4270	Forced	4850	4040	250
IC02-19 4	Constrained	3660	Forced	4170	4040	250
IC02-20 4	Constrained	3540	Forced	4030	4040	250
IC03-17 4	Constrained	4500	Forced	5120	4040	250
LEIT-1 4	Constrained	4940	Forced	5620	5210	110
LEIT-2A 4	Constrained	4570	Forced	5200	5210	110
LEIT-2B 4	Constrained	5110	Forced	5820	5210	110
LEIT-3 4	Constrained	5240	Forced	5960	5210	110
LEIT-4 4	Constrained	4560	Forced	5190	5210	110
LEIT-5 4	Constrained	4680	Forced	5330	5210	110
BUR-1 4	Constrained	8190	Forced	8410	8060	120
BUR-4 4	Constrained	6990	Forced	9730	8060	120
BUR-2 4	Constrained	7450	Forced	7920	8060	120
BUR-3 4	Constrained	8600	Forced	7540	8060	120
BUR-6 4	Constrained	8620	Forced	9750	8060	120
BUR-5 4	Constrained	6670	Forced	9300	8060	120
IC02-11 4	Constrained	10050	Forced	11650	10330	80
IC02-10 4	Constrained	10100	Forced	11590	10330	80
IC03-9 4	Constrained	11650	Forced	15800	10330	80
IC02-1A 4	Constrained	13720	Forced	14200	10330	80
IC02-1B 4	Constrained	12330	Forced	9950	10330	80
IC02-7 4	Constrained	8720	Forced	13410	10330	80
	Weighted RMS error	152.60		421.435869

By weighted RMS, the constrained (68 a / g / yr) experiment appears to have a smaller error than the forced (57 a / g / yr) experiment.  This would suggest that the calibration technique is not the source of difference between ACE and Licciardi et al calibrated production rates.

Different Scalings

The Licciardi et al calibration uses Lal and Stone scalings, which are not optimal for ³⁶Cl as they do not specify scaling for low-energy neutrons, just spallation neutrons and muons (in the case of Stone). We can check for similarities in scaling numbers between ACE and Licciardi et al 2008.  The below compares ACE output with Table 4 of Licciardi et al 2008 for the Lambahraun group of samples:

So the scalings appear to be very similar, for example for sample IC02-16 Licciardi et al report spallation scalings of 1.586, where ACE reports 1.59.  For muons Licciardi et al report a scaling of 1.249, for ACE 1.24. It is difficult to conceive how differences in these numbers could generate such a large difference in calibrated Ca spallation production rates. However scalings for low-energy neutrons could be different, as they are not formally stated with the Lal or Stone scaling.  However, as contributions of low-energy neutron capture are small (Figure 3 of Licciardi et al) this is unlikely.

The calibration dataset is provided with two different atmospheric specifications, the standard atmosphere:

• default sea-level pressure = 1013.2 mb,
• default sea-level temperature = 15 °C,
• default lapse rate = 6.5 °C / km.

And also a low-pressure atmospheric specification, with default sea-level pressure = 1005.8 mb to account for the Icelandic low pressure system.  Both of these datasets are available in ACE format here. The following calibration plot shows the Ca spallation rate for this low-pressure dataset together with Stone scaling, which accounts for pressure instead of elevation:

With a lower atmospheric pressure and Stone scaling the calibrated Ca spallation rate is 63 a / g / yr. This corresponds closely to Licciardi et al’s results, as they also found a similar reduction in Ca spallation for lower atmospheric pressure.

The question remains: Why does ACE and Licciardi et al generate different Ca spallation rates when using the same dataset and scaling?

Constrained Calibration Patch

This entry was posted on Friday, January 9th, 2009 at 5:27 pm and is filed under Software News. You can follow any responses to this entry through the RSS 2.0 feed. Both comments and pings are currently closed.