Age distribution of horizontal fission tracks

Equations for the distribution of age versus length of partially annealed horizontal fission tracks in apatite is presented. Probabilistic least–squares inversion corrects natural track length histograms for observational biases considering the variance of data, modelization, and prior information. The corrected histogram is validated by its variance–covariance matrix. It is considered that horizontal track data can be with or without measurements of angles to the c–axis. In the last case, 10 3D–histograms are introduced as an alternative to histograms of c–axis projected track lengths. Thermal history modeling of samples is not necessary for track age distribution calculation. In an example the age equations are applied to apatites with pre–depositional (inherited) tracks, to extract the post–depositional track length histogram. Fission tracks generated before deposition in detrital apatite crystals are mixed with post–depositional tracks. This complicates the calculation of the post– sedimentary thermal history as the grains have experienced different thermal histories until deposition. Thereafter the grains 15 share a common thermal history. The extracted post–depositional histogram without inherited tracks may be used for thermal history calculation. 1 Overview of the formation of fission tracks Fission of U–238 in apatite, sphene, and zircon create tracks in the crystal lattice. Track density is reduced as a function of temperature and time as tracks anneal and shorten due to atom diffusivity (Li et al., 2011; Fleisher et al., 1975; Afra et al., 2

3D-histograms are introduced as an alternative to histograms of c-axis projected track lengths. Thermal history modeling of samples is not necessary for track age distribution calculation. In an example the age equations are applied to apatites with pre-depositional (inherited) tracks, to extract the post-depositional track length histogram. Fission tracks generated before deposition in detrital apatite crystals are mixed with post-depositional tracks. This complicates the calculation of the postsedimentary thermal history as the grains have experienced different thermal histories until deposition. Thereafter the grains 15 share a common thermal history. The extracted post-depositional histogram without inherited tracks may be used for thermal history calculation.

Overview of the formation of fission tracks
Fission of U-238 in apatite, sphene, and zircon create tracks in the crystal lattice. Track density is reduced as a function of temperature and time as tracks anneal and shorten due to atom diffusivity (Li et al., 2011;Fleisher et al., 1975;Afra et al., 20 2011). The ongoing track generation through time and simultaneous annealing is used to derive the thermal history of the sample. Bertagnolli et al. (1983) derived the differential equation describing the track length distribution within a mineral due to annealing through time. The surface density of tracks is also described. Forward calculation examples are given by Bertagnolli et al. (1983). On this basis, Keil et al. (1987) developed an inversion procedure, where the temperature history is derived from either the length distribution of tracks within the crystal or from the length distribution of projected tracks 25 intersecting the surface. The calculation procedure is direct and without the use of a Monte Carlo type optimization. The time of track generation can be derived from the observed track length distribution independent of any annealing law (Keil et al., 1987): with an annealing model. Keil et al. (1987) showed, by a synthetic example for projected tracks, that the thermal history can be calculated numerically stepwise backward in time starting with the present temperature. The large uncertainties of projection of tracks mean that the approach by Keil et al. (1987) is dubious. The selection of horizontal, confined tracks are recommended instead (Laslett et al., 1984), that is "tracks identified by the constancy of focus over their entire length and strong reflection 35 in incident illumination" (Gleadow et al., 1986b;Gleadow et al., 2019). The model by Keil et al. (1987) does not include blurring of track length histograms caused by the initial distribution of fission fragment energy (Jungerman and Wright, 1949), annealing and etching anisotropy, mineral composition, the uncertainty of measurement (Ketcham, 2003), and track selection biases Jensen et al. (1992). Proportionality of the number of tracks in a track length histogram column and time assumed by Keil et al. (1987) is disturbed by the blurring. Jensen et al. (1992;1993) extended Eq. (1) to confined horizontal tracks instead 40 of projected tracks. Belton and Raab (2010) also used a backward cumulative method to estimate the track ages.
Three major biases appear to be important when deriving the equation for track age: 1) Surface track density bias reflecting the likelihood of a track to be exposed to etching on the surface. We use an exponential approximation to surface track density versus mean track length for induced tracks annealed in the laboratory (Green, 1988). 2) Track length bias due to the likelihood 45 of a track being exposed to etching through fractures and tracks cutting the etched surface. We assume it to be proportional to the track length. 3) Selection bias due to the likelihood of an etched track being accepted as horizontal. We assume that tracks within a given angle from the horizontal are counted. The alternative that all tracks in focus are accepted is discussed in Appendix C. Besides biases, there is a range of effects contributing to large observational deviations of track lengths as has been documented in interlaboratory comparisons (Ketcham et al., 2009). The effects are caused by personal etching and 50 observational practice, and sampling statistics. Variation of the number of tracks in a given length interval is dependent on the spontaneous character of fission. The disordering of the unique relationship between track length and time caused by various biases and variances is essentially restored by deconvolution of natural track length histograms (Jensen et al., 1992).
Deconvolution is performed by mathematically simulated annealing (Kirkpatrick et al., 1983) and with the use of filters based on annealing of induced tracks in the laboratory. Jensen et al. (1992;1993) and Jensen and Hansen (2018) used that the track 55 ages can be derived separately from the thermal history as shown by Keil et al. (1987). The columns of the deconvolved track length histogram are then converted to equivalent time intervals. The track ages are obtained by backward cumulation from long tracks toward short tracks. The deconvolution procedure described in Jensen et al. (1992) is for measurements where the c-axis angle is not measured. When they are available it is the practice to project the tracks on the crystal c-axis following the procedure described by Donelick et al. (1999). As an alternative, the inversion method presented here uses 3D-histograms. 60 Age-track length relations are given explicitly in contrast to the indirectly embedded relations in the computer programes by Green et al. (1989), Lutz and Omar (1991), Ketcham (2005), and Gallagher (2012). Our age calculations are calibrated without sample temperature history modeling. The inversion procedure presented here is based on least-squares probabilistic inverse theory (Tarantola, 2005). The method is computational fast. https://doi.org/10.5194/gchron-2021-8 Preprint. Discussion started: 2 March 2021 c Author(s) 2021. CC BY 4.0 License.

Equation for the age distribution 65
Consider a mineral as apatite and a time interval ∆ in which randomly oriented fission tracks are generated with the initial track length . The number of tracks per unit volume generated in the time interval is proportional to the fission decay frequency , the uranium U-238 concentration , and the length of the time interval ∆ : is the number of time intervals. The exponential decay of U-238 nuclei is ignored, to begin with. A track is generated for 70 each fission; therefore, is also the number of initially randomly oriented tracks per unit volume. The track length histogram of the randomly oriented tracks per unit volume generated in a single time interval ∆ is initially Vectors are bold and the transpose transforms rows into columns. is the number of track length bins. The tracks are gradually annealed from the initial length by temperature and time leading to shorter tracks with the mean track length (Fig.  75 1). The tracks generated in the time interval ∆ are after partial annealing distributed into various track length columns. At present, after partial annealing, the track length histogram is For the tracks that are not annealed below the detection limit the number of tracks per unit volume (not bold) is equal to the number of fissions in the time interval that is 80 The surface track density is initially 0 and at present , Fig. 1.  A fraction of the randomly oriented tracks cut other tracks connected to the surface providing paths for etchants making them observable in the light microscope. The near-horizontal tracks are selected for length measurements. The track length histogram of the observed near-horizontal tracks generated in the time interval ∆ is 90 is expected to be linearly dependent on the randomly orientated unexposed tracks which means that the histogram is derived by multiplying the histogram by a set of constants 1 , .., : where is a proportionality constant. Equation (7) relates the partly annealed randomly oriented unetched tracks to the 95 observed horizontal tracks. The 'th elements of is The number of observable tracks generated in the time interval ∆ is It is expected that long tracks are more likely to be etched and observed than short tracks when using the angle selection criteria (Ketcham, 2019). This means that the proportionality constants where is the track length of the partially annealed track. The relation between the randomly oriented tracks to the number of observable tracks is then 105 Equation (8) and Eq. (10) lead to: Summation on both sides of the equation leads to 110 which means that https://doi.org/10.5194/gchron-2021-8 Preprint. It is expected that the shape and mean track length of the histogram of observable tracks with origin in the time interval ∆ can be reproduced in a laboratory annealing experiment: 115 where is a proportionality constant. is a track length distribution histogram derived from annealing in the laboratory, or an interpolation thereof: and has identical mean track lengths. 120 because is normalized for all : The natural observed histogram is the sum over all the histograms , each of which linked to time interval : Therefore, a natural fission-track histogram can be approximated by a weighted sum of interpolated laboratory track length histograms ( Fig. 2   is a filter transforming the weights into the natural histogram : or https://doi.org/10.5194/gchron-2021-8 Preprint. Each column of the matrix is a filter with increasing mean track length from column 1 to . Solving Eq. (21) for (the model) is named deconvolution. The equation does not take data variance into account. Therefore, a probabilistic inverse method is used (Tarantola, 2005) assuming Gaussian distributions of variance-covariance of data, prior information, and modelization. The solution (the model) is named the posterior ̃ (Tarantola, 2005): Prior information is derived from other independent information or, if not available, simply homogeneous with the elements equal to the mean value of the number of measured tracks. is the variance-covariance matrix of the natural observed data . The measurements of a given number of track lengths are statistically considered as several possible outcomes of a trial greater than two. The multinomial statistical distribution is therefore used to describe variance-covariance.
The diagonal elements of are the variances of the observed number of tracks for each length interval.: 150 where = / is the probability of the observed data . is the number of measured natural tracks. The off-diagonal elements are the covariances

155
Modelization variance caused by the variance of , Eq. (21), is calculated by forward modelling using Eq. (21) and added to (Tarantola, 2005). This variance is calculated by random Gaussian realizations of the filters assuming the standard deviations to be √ . The variance-covariance matrix for the prior is calculated similarly to . When the influence of prior information is unwanted large variance values are given as diagonal elements of the matrix . The variance-covariance of the posterior is calculated following Tarantola (2005): 160 An approximation to the observed data is calculated forwardly using the posterior: ̃=̃ .
The deconvolved histogram ̃ is essential for track age calculation as discussed below.

165
The density of induced tracks, generated in the time interval ∆ , intersecting a polished surface plane of the mineral is related to the mean track length (Green et al., 1986). It is here approximated by a logarithmic expression (Appendix B): where 0 is the reduced surface track density, 0 is the reduced mean track length, and is a calibration constant dependent on the mineral composition. Natural tracks generated in the time interval ∆ and with a present mean track length is expected 170 to have the same reduced surface track density as laboratory annealed tracks with the same reduced mean track length. Equation (28) is therefore also valid for the natural tracks. The initial surface track density 0 generated in the time interval ∆ is expected to be proportional to the initial mean track length and time where is a calibration constant. Combining Eq. (28) and Eq. (29) leads to the surface track density assigned to the time 175 where The value of ℒ , having the unit of length, is less than the mean track length and is considered as a correction to the mean 180 track length . The natural surface track density is composed of contributions from all the time intervals Inserting the right side of Eq. (30) instead of in Eq. (32) leads to Remembering that ℎ = , Eq. (17), and that the result of the inversion ̃ is an approximation to ℎ we get together with 185 Eq.(14) Inserting Eq. (14) for ∆ into Eq. (30) leads to The ratio between and is 190 .
(37) 195 The formation time (age) for the oldest track in each column of the histogram ̃ is the cumulation of the time intervals ∆ corresponding to the column and younger columns. (38) That is 200 The unit of ̃ is number/ 3 but since it appears both in the nominator and denominator ̃ can as well be considered dimensionless. Therefore the volume need not be measured. The decreasing uranium concentration through time is considered by introducing the logarithm in Eq. (39).
where is the total decay constant. The corresponding surface track density is and using Eq. (36) 210 Eq(39) and Eq. (42) are valid for tracks selected within a given angle from the horizontal. Equations for the case when tracks are selected following the requirement that both ends are in focus at the same time are given in Appendix C.
Inversion following Eq. (23) applies both for data measured in the old-fashioned way where track angle to the c-axis is not measured as well as for new measurements which include the angle, Appendix D. 215 and approximately The variance of is derived through the variance of the three parts 1 , 2 , and 3 of the equation: The age of the oldest track in column p is then 225 The variance of each part is calculated using standard textbook statistics. The variance of 1 is The second part 2 is a cumulation of correlated histogram columns . The variance of 2 therefore includes the covariance: ( 2 ( )) = ∑ ( ) + 2 ∑ ( , The variance of 3 is where (ℒ ) = ( 0 ln ( 0 )) 2 1 4 ( ).
(53) 235 The equation for the age, Eq. (48) repeated, is written in terms of the three parts The variance of is then

Inherited tracks
The deconvolution technique can be used to extract the post-depositional part of track length histograms with inherited tracks. This is illustrated in Fig. 3   The post-depositional part of the deconvolved histogram ̃ are the columns having track ages less than the deposition age : ≤ , (56)  280 where is the column number. Together with Eq. (40), we obtain Summations are performed for all columns representing the post-sedimentary history. Together with the surface track density the post-depositional histogram may be used to calculate the post-depositional thermal history applying the backward modeling procedure described by Jensen et al. (1992). Alternatively, the post-depositional de-convolved histogram can be convolved by forward calculation by Eq. (27) and used by random or guided-random search algorithms (Lutz and Omar, 290 1991;Willett, 1997;Ketcham, 2005;Gallagher, 1995;2012).
Identification of post-sedimentary fission-tracks is exemplified by applying apatites from the sample GGU103113 (Jameson Land, East Greenland). The Middle Jurassic sandstone sample (170 Ma) has the apparent fission-track age 245 Ma (Hansen et al., 2001). The measured and the deconvolved track length distributions are shown in Fig. 4a and 4b.   (Hansen, 1996). with accumulated ages less than the deposition age, Fig. 4c and Fig. 4d. The deconvolved histogram Fig. 4d can be used to calculate the post-sedimentary thermal history. The histogram in Fig. 4e is the convolution of the histogram in Fig. 4d.

Discussion
The equations for accumulated fission-tracks derived here are based on the practice of selecting horizontal tracks for length 310 measurements. We follow the recommendation by Ketcham (2019) who selects all tracks within ±10° from the horizontal.
Alternatively, Gleadow et al. (1986b) and Galbraith (2005) select tracks observed in reflected light being within ±10° from the horizontal and with both ends in focus in transmitted light. Fewer long tracks (> 8.5 µm) are then selected relatively to shorter tracks (Appendix C). An alternative set of equations are given in Appendix C for this situation.

315
Standard deviations of the ages are calculated based on the variance of input data, prior model, and modelization. The calculated deviations are high but they agree with deviations of histograms reported in an interlaboratory comparison (Ketcham et al., 2009). It is required that the laboratory annealed tracks are measured in the same way as the natural tracks. E.g., in the Jameson Land example, the recommendations by Gleadow et al., (1986b) and Gleadow et al., (2019) are used in both cases.

320
The inversion procedure by Tarantola (2005) requires a prior model being independent of data. As the default, a prior model histogram with an equal number of tracks in the columns is chosen. In some cases, negative track length histogram columns appear as a result of inversion because due to unrealistic data. An improved prior model based on other information is then necessary. The problem can be managed by smoothing of input data and/or of covariance matrices.

325
When track angles to the c-axis is available it is common practice to calculate the histogram of the c-axis projected tracks.
However, detailed information on track density is lost along the accumulation along the projection path. Instead, one can use 3D-track length histograms as described in Appendix D. 3D-histograms do not involve projections and therefore retain both length and angle information. The inversion procedure presented here is valid for both 2D-and 3D-histograms (Appendix D). is a tendency to off-select longer tracks. If all tracks within a given angle are selected this bias disappears. The equations presented here are prepared for both cases. The calculation of track ages does not require the calculation of the sample thermal history. Instead, they are the basis for temperature history calculation. As an example, the deconvolution method presented here is used to identify the post-depositional part of a track length histogram from Jameson Land, Greenland. 340

Appendix A
The filters used for deconvolution are the columns of matrix , Eq. (22). They are based on measurements of track lengths after annealing in the laboratory (Gleadow et al., 1986a;Green et al., 1986;Barbarand et al., 2003). The histograms of these

Appendix B
Plots of mean track length versus surface track density for tracks generated by neutron radiation and annealed in the laboratory 355 are given by Green et al. (1986). A logarithmic expression approximates data, Fig. B1:

Appendix C 365
Ketcham (2019) recommends the selection of all etched tracks being within ±10° from the horizontal. It is expected that the likelihood of tracks being etched is proportional to their length. This view has been adopted in the main text of this paper.

370
The axial (vertical) resolution of a light microscope is With typical values for wavelength of = 550 nm (green light), refraction index of oil = 1.51, and numerical aperture = 1.25, the axial resolution = 0.74 µm. The numerical aperture is a measure of the maximum angle for which the microscope receives light. The criteria by Gleadow et al (2019) imply that a short track, 6 µm long, is accepted for counting if 375 the angle to the horizontal is less than 10°. A long track, 18 µm, has to be within 2.5° to the horizontal to be accepted. The tendency is therefore that long tracks are off-selected from counting. The likelihood of selection is then inverse proportional to the track length for tracks longer than 6 µm. At the same time, the likelihood of a horizontal track being etched is expected to be proportional to its length. These two biases essentially cancel each other (Jensen et al., 1992). With the focus window selection criteria (Gleadow et al., 2019) the set of age equations is then: 380 The age of the oldest track is obtained summing all the time intervals: Introducing the mean track length of the inverted histogram ̃ : The age of the oldest observable track is equivalent to the derivation by Jensen et al. (1992) when ℒ is replaced by the mean track length.

Appendix D
The inversion principle presented in the main text is for 2D-track length histograms ignoring track angles to the c-axis. 390 However, the inversion principle is also applicable to data that includes the angles. The measured data vector in Eq. (22) is then derived from a 3D-histogram (length, angle, and number). The elements of the vector is obtained by sequentially numbering each bin (Tarantola, 2005). An example is shown in Fig. D1 with a synthetic 3D-histogram as data obtained from two annealing experiments. They are the result of heating to 275ºC and 320ºC respectively (Barbarand et al., 2003). Each bin of the 3D-histogram is 5º times 0.5 µm. The matrix , Eq. (22) consists of twenty 3D-filters interpolated between the six annealing experiments by 400 Barbarand et al. (2003). Each column of is a filter. Two of them are shown in Fig. D2.

410
The standard deviations are the diagonal values of the posterior covariance matrix. The result of the inversion is evaluated by comparing forward modeling using Eq. (21) with the inverted histogram as input data (Fig. D4). Compliance is good.

Code availability
Octave codes for reproducing the results shown in Fig. A1 and D1-4 are available online on Zenodo