Deconvolution of fission–track length distributions and its application to dating and separating pre– and post–depositional components

To enable the separation between pre– and post–depositional components of the length distribution of (partially annealed) horizontal confined fission tracks it is corrected by deconvolution. Probabilistic least–squares inversion corrects natural track length histograms for observational biases considering the variance of data, modelization, and prior information. 10 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, 3D–histograms are introduced as an alternative to histograms of c–axis projected track lengths. Thermal history modeling of samples is not necessary for the calculation of track age distributions of corrected tracks. 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 15 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 share a common thermal history. Therefore, 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, titanite, and zircon create tracks in the crystal lattice. Track density is reduced as a function of 20 temperature and time as tracks anneal and shorten due to atom diffusivity (Li et al., 2011; Fleisher et al., 1975; Afra et al., 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 25 derived from either the length distribution of tracks within the crystal or from the length distribution of projected tracks 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 model (Keil et al., 1987): τ(λ0) = 1 ε ∫ n(λ)dλ λ0 1 , (1) 30

temperature and time as tracks anneal and shorten due to atom diffusivity (Li et al., 2011;Fleisher et al., 1975;Afra et al., 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 25 derived from either the length distribution of tracks within the crystal or from the length distribution of projected tracks 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 model (Keil et al., 1987): where 0 is the present track normalized length, is the fission-track production rate, ( ) is the measured number of fission tracks of length produced per volume. Equation (1) is understood as follows: The annealing properties of the mineral affect the final apparent age as well as the expected age of the oldest randomly oriented unetched track. However, the expected age of this track can be determined by counting the number of all tracks in a volume and divide by the track generation rate. In the simplified case of no spread of track lengths, the expected age of any given track is calculated by counting the number of 35 shorter tracks plus one. This age is determined without the use of an annealing model. The unetched tracks are not routinely measured. Etched confined horizontal tracks are measured instead. This manuscript explains how they can be used together with the surface track density to separate pre-and post-depositional components.
The temperature history is derived from the distribution of projected tracks together with an annealing model. Keil et al. (1987) 40 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 is recommended instead (Laslett et al., 1984), that is "tracks identified by the constancy of focus over their entire length and strong reflection 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 45 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 of projected tracks. Donelick (1988) introduced a fission-track inversion procedure based on principles like that of Jensen et al. (1992) and the present paper. 50 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 of a track being exposed to etching through fractures and tracks cutting the etched surface. We assume it to be proportional to 55 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 D. 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 a) fluid filling track tips, b) poorly chosen track tips, c) analyst inexperience. Variation of the number of tracks in a given length interval is dependent 60 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 corrected track ages can be derived separately from the thermal history as shown by Keil et al. (1987). The columns of the 65 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. We are not assigning an age to the bins of observed fission-track length histograms. Ages are assigned to the bins of the deconvolved (corrected) track length histograms which can then be used to age date thermal events. The procedure of deconvolution described in Jensen et al. (1992) is for measurements where the caxis angle is not measured. When they are available it is the practice to project the tracks on the crystal c-axis following the 70 procedure described by Donelick et al. (1999). As an alternative, the inversion method presented here uses 3D-histograms.
Age-length relations for corrected tracks are given explicitly in contrast to the indirectly embedded relations in the computer programs by Green et al. (1989), Lutz and Omar (1991), Ketcham (2005), and Gallagher (2012). Our age calculations are performed without sample temperature history modeling. The inversion procedure presented here is based on least-squares probabilistic inverse theory (Tarantola, 2005). The method is computationally fast. 75

Summary of age and temperature calculation
Before entering the mathematics of age calculation, a recipe for age dating of an uplift event is given including temperature calculation. First, an idealized model is considered with tracks generated continuously (not spontaneously), no length spreading during annealing, no etching, and no errors of measurements: 1. Construct a track length density histogram (numbers/volume) of the randomly oriented unetched tracks. 80 2. Convert the track length histogram into a histogram of time intervals by dividing the number of tracks in each column by the rate of track generation. Each column of this histogram represents the time it takes to generate the tracks belonging to the column.
3. Cumulate the time intervals from the youngest (longest tracks) to the oldest (shortest tracks) to obtain a histogram of track ages. 85 4. If the uplift is sufficiently pronounced a statistically significant break is visible in the age histogram. The event is hereby age-dated.
5. The age of the event can be used to extract the pre-event track histogram from the left part (shortest tracks) of the cumulated ages and the post-event histogram from the right part (longest tracks). To this point, there has been no need for a track annealing model. 90 6. The post-event temperature history can now be calculated separately which is useful if the pre-event tracks are of mixed origins. The temperature is calculated for each column of the time interval histogram, as explained in Jensen et al. (1992).
In practice, however, confined horizontal etched track length histograms are used together with the surface track density, for 95 age dating and calculation of the history of temperature. Tracks are produced spontaneously, their lengths are spread during annealing, observations are biased, and measurements are uncertain. The recipe then sounds: 1. Construct a track length histogram of the observed confined horizontal etched tracks.
2. Remove the spreading of track lengths by deconvolution (deblurring) of the observed histogram.
3. Use the deconvolved histogram for age dating and calculation of the history of temperature following the recipe from 100 point 2 for the idealized tracks.
The recipes are based on the development by Bertagnolli (1983), Keil et al. (1987), and Jensen et al. (1992). See also Donelick (1988) for a similar development. The procedure for calculating the temperature history based on the deconvolved track length histogram was presented by Jensen et al. (1992). At that time, deconvolution was performed by trial and error. Mathematically 105 simulated annealing was later applied instead (Jensen and Hansen, 2018). Age dating and identification of inherited tracks using Tarantola inversion are presented in the next sections. The calculations are direct with no use of Monte Carlo simulations.
Readers unfamiliar with the technique of deconvolution may benefit from reading Appendix A before entering the Tarantola inversion technique.

Correcting for biases of track length histograms by deconvolution
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 ∆ : (2) 115 is the number of time intervals. The exponential decay of U-238 nuclei is ignored, to begin with. A track is generated for 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 120 gradually annealed from the initial length by temperature and time leading to shorter tracks with the mean track length (Fig.   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 125 number of fissions in the time interval that is 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 135 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 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 , .., : 140 where is a proportionality constant. Equation (7) relates the partly annealed randomly oriented unetched tracks to the observed horizontal tracks. The 'th elements of is The number of observable tracks generated in the time interval ∆ is 145 Horizontal confined tracks of large area (the product of etched length and width) are likely to be etched. This means that the proportionality constants are as follows: where is the track length of the partially annealed track. The relation between the randomly oriented tracks to the number 150 of observable tracks is then Equation (8) and Eq. (10) lead to: Summation on both sides of the equation leads to 155 which means that It is expected that the shape and mean track length of the histogram of observable tracks with origin in the time interval 160 ∆ can be reproduced in a laboratory annealing experiment: where is a proportionality constant. is a track length distribution histogram derived from annealing in the laboratory, or an interpolation thereof: and have identical number of tracks: because is normalized for all : The natural observed histogram is the sum over all the histograms , each of which linked to time interval : 170 and with Eq. (15) Therefore, a natural fission-track histogram can be approximated by a weighted sum of interpolated laboratory track length histograms ( Fig. 2 and Fig. B1 in Appendix B). is a filter transforming the weights into the natural histogram : or 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. Predictions by this equation are assumed to be error-free. Instead, the Bayesian framework (Bardsley, 2018) is used to include the probability density for 1) the observations , 2) prior information for the solution, and 3) the solution . Probability density distributions are assumed to be Gaussian. The solution is then the posterior minimizing twice the misfit function (Tarantola, 2005): 190 where the notation || || = −1 is used.
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. 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 200 therefore used to describe variance-covariance. The diagonal elements of are the variances of the observed number of tracks for each length interval.: where = / is the probability of the observed data . is the number of measured natural tracks. The off-diagonal elements are the covariances 205 The modelization matrix is not an error-free physical model but based on uncertain measurements of tracks generated in the laboratory. This is accounted for by replacing with + where is the modelization covariance matrix (Tarantola, 2005). The elements of are calculated as a diagonal matrix: where is the identity matrix, =̃, and = ( ) being the variances of the smoothed random Gaussian realizations of the filter matrix assuming the standard deviations to be √ . Once is calculated it is reused with other data .

215
The variance-covariance matrix for the prior is calculated similarly to : for the diagonal where = ℎ / is the probability of the priors ℎ . is the number of prior tracks being equal to the number of measured tracks. The posteriors are assumed to be a Gaussian distribution. This allows for negative posterior values which are unphysical. To suppress this, we choose positive priors with restricted variance. The priors are distributed equally, ℎ = 220 / for all , in case there is no other prior information available. The weights are experimentally determined so that they are as large as possible while negative posteriors are avoided. Large weights decrease the influence of the priors and vice versa.
An alternative way of avoiding negative posteriors is given in the discussion section. The off-diagonal elements are the covariances Stability of the inversion of the matrixes , , and ( −1 + −1 ) in Eq. (24) is tested by the condition numbers for the matrixes (Calvetti, 2007). They prove to be stable in the examples given in this paper. The subject of stability is further dealt with in the discussion section.
The variance-covariance of the posterior is calculated following Tarantola (2005): 230 An approximation to the observed data is calculated forwardly using the posterior: ̃=̃ .
The deconvolved histogram ̃ is essential for the calculation of track age as discussed below.

Equations for the age distribution of corrected tracks 235
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 C): 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 240 to have the same reduced surface track density as laboratory annealed tracks with the same reduced mean track length. Equation (32) 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. (32) and Eq. (33) leads to the surface track density assigned to the time 245 interval ∆ : 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 250 track length . The natural surface track density is composed of contributions from all the time intervals Inserting the right side of Eq. (34) instead of in Eq. (36) leads to Remembering that ℎ = , Eq. (17), and that the result of the inversion ̃ approximates ℎ we get, together with Eq. (14), 255 Inserting Eq. (14) for ∆ into Eq. (34) leads to

265
The formation time (age) for the oldest corrected track in each column of the histogram ̃ is the cumulation of the time intervals ∆ corresponding to the column and younger columns.
That is 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. (43).
where is the total decay constant. The corresponding surface track density is 275 and using Eq. (40) Equation (43) and Eq. (46) are valid for tracks selected within a given angle from the horizontal. Equations for the case when 280 tracks are selected following the requirement that both ends are in focus at the same time are given in Appendix D.
Inversion following Eq. (24) 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 E.

Variance of ages 285
The oldest age of corrected tracks in column given by Eq. (43) and repeated here: and approximately The variance of is derived by the linear approximation 290 where 0 is the vector with the elements , , , , ∑ ( ) . is the Jacobian calculated numerically, and = − 0 .
The variance-covariance matrix of is then where is the variance-covariance matrix of . The variances of the oldest ages of the corrected tracks are the diagonal elements of .

Testing the deconvolution method 300
Deconvolution is a mathematical tool capable of reducing noise and correct for biases in time series such as seismic signals.
We, therefore, expect that deconvolution can reduce the length spread of confined horizontal fission tracks. This is shown in the following example. Assume that a sample starts generating tracks at a temperature of 63 °C, 150 Ma back in time. At 50 Ma it is suddenly uplifted to a temperature of 20 ºC until the present. The corresponding track length histogram is calculated forwardly resulting in the histogram seen in Fig. 3a. The track length annealing model by Stephenson et al. (2006)

7 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. 4 by a simplified model for track annealing and observation. Grains from two source areas are considered. There are more tracks in the post-depositional columns of Grain 2 than in Grain 1 because it is assumed that the 330 uranium content of Grain 2 is twice that of Grain 1. The two post-depositional columns of Grains 1 and 2 represent the same two-time intervals. For both Grains 1 and 2, there is one column where pre-and post-depositional tracks are mixed. This column and the pre-depositional columns of both grain types do not necessarily represent the same time intervals because their different thermal histories create individual track length distributions (histograms). The pre-depositional tracks of Grain 1 compared with those of Grain 2 have experienced different thermal histories and therefore different degree of annealing. Post-335 depositional tracks are ordered with decreasing length as a function of time in contrast to the pre-depositional tracks which may be disordered. Fig. 4C sums Grains 1 and 2 histograms. It is observed that the post-depositional time intervals for the bulk histogram are identical to those for the two separate grain types. The age equation, Equation (44), is valid for postdepositional corrected tracks but not for mixed or pre-depositional tracks. The post-depositional part of the deconvolved histogram ̃ are the columns having track ages less than the deposition age 345 where is the column number. Together with Eq. (44), 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 355 modeling procedure described by Jensen et al. (1992). Alternatively, the post-depositional de-convolved histogram can be convolved by forward calculation by Eq. (31) and used by random or guided-random search algorithms (Lutz and Omar, 1991;Willett, 1997;Ketcham, 2005;Gallagher, 1995;2012).

Identification of post-sedimentary fission-tracks is exemplified by applying apatites from the sample GGU103113 (Jameson 360
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. 5a and 5b.  (Hansen, 1996).

370
This histogram is the basis for direct calculation of the temperature history (Jensen et al., 1992); (e) The convolved post-sedimentary histogram, resembling a measured histogram, can be used by Monte Carlo type simulation models (Ketcham, 2005;Gallagher, 1995).
The parameter in the age equation, Eq. (44) is adjusted to 0.752 (with b = 0.784) to obtain simulated apparent age equal to the apparent age reported by Hansen et al. (2001). The apparent age is calculated using a track length distribution of an 375 unannealed track length distribution in Eq. (44). The post-depositional columns are identified as being the rightmost columns with accumulated ages of corrected tracks less than the deposition age, Fig. 5c and Fig. 5d. The deconvolved histogram Fig.   5d can be used to calculate the post-sedimentary thermal history. The histogram in Fig. 5e is the convolution of the histogram in Fig. 5d.

Discussion 380
The equations for accumulated fission-tracks derived here are based on the practice of selecting horizontal tracks for length 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 D). An alternative set of equations are given in Appendix D for this situation. 385 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 390 both cases.
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 because of inversion due to unrealistic data. An improved prior model based on other information is then necessary. 395 The problem can also be managed by smoothing input data and/or covariance matrices. Another approach to be examined is the Gaussian approximation to the Poisson distribution (Calvetti and Somersalo, 2007).
The stability of the inversions in Eq. (24) depends on the condition numbers for , , and ( −1 + −1 ) which should be less than 10 5 … 10 6 . Application of the Octave function cond (Eaton et. al, 2020) for the 2D and 3D examples in this paper 400 shows numbers less than 10 5 … 10 6 indicating stability. Refined binning of the histograms may lead to large and possible unstable inversions. In those cases and in Eq. (23)-(24) and ( −1 + −1 ) in Eq. (24) should be replaced by the Cholesky factorized matrixes, see e.g., Cheney and Kincaid (2012). It is then required that the matrixes are symmetric and positive definite. It is seen from Eq. (25)-(29) that the matrixes and are symmetric. The matrix ( −1 + −1 ) in Eq.
(24) is symmetric as well. Application of the Octave function chol (Eaton et. al, 2020) shows that the matrixes are positive 405 definite. The three matrixes can therefore be Cholesky factorized.
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 during the accumulation along the projection path. Instead, one can use 3D-track length histograms as described in Appendix E. 3D-histograms do not involve projections and therefore retain both 410 length and angle information. The inversion procedure presented here is valid for both 2D-and 3D-histograms (Appendix E).

Conclusions
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. A new procedure for deriving ages of corrected fission-tracks as a function of track length is suggested. The natural fission-track length histogram is considered linearly composed of histograms of induced tracks 415 partially annealed in the laboratory. Inversion is performed by a probabilistic inverse theory. The resulting ages of corrected tracks are given as centers of posterior with variance. The equations are valid for old-type measurements where the track angle to the c-axis is not measured as well as for recent data including measurements of the angle to the c-axis. Data with both track lengths and angles are organized in 3D-histograms presented as images. Two types of track length measurements are discussed.
If tracks with both ends in focus are selected there is a tendency to off-select longer tracks. If all tracks within a given angle 420 are selected this bias disappears. The equations presented here are prepared for both cases. The calculation of the ages of corrected tracks does not require the calculation of the sample thermal history. Instead, they are the basis for temperature history calculation.

Appendix A 425
After the formation of the unetched randomly oriented tracks, the spreading of lengths increases during the process of annealing leading to a considerable mixing of lengths of the observed track length histogram. There is not a unique relationship between time and the length of individual tracks. However, the relationship is not completely blurred. There is a tendency that old tracks appear in the short track length part of the histogram and young tracks in the opposite part. To some extent, deconvolution can reduce the spread. Deconvolution is used extensively in seismic processing to improve the signal-to-noise ratio and remove 430 biases. This requires the character of the noise to be known. The "noise" in connection with fission tracks is the spread observed in laboratory annealing experiments. This is used to reduce the spread of the observed track length histograms. That is, the observed histogram can be deblurred by deconvolution. A simplified demonstration of deconvolution of a bimodal track length histogram is shown in Fig. A1. The present paper applies a more advanced least-squares deconvolution technique.

named filters. Histogram A is histogram B multiplied by the number of tracks in histogram C. Histogram D is histogram E multiplied by the number of tracks in histogram F. Histogram G is histogram A plus histogram D. Histogram G is then the result of convolution of histogram I with the filters B and E. The synthetic histogram G is compared with the observed
histogram H however, the left part is dissimilar. Therefore, the suggested deconvolved histogram I is not successful.

445
A successful deconvolution is shown in Fig. A2. The spread of tracks in the bimodal histogram H is assembled in two columns corresponding to two thermal time periods. Note that two time periods can be separated despite the track length overlap of the histograms A and D. The deconvolved track length histogram I is arranged according to time. Each column can be converted to an equivalent time interval. The age of the oldest track of each column of deconvolved histograms is calculated by adding its associated time interval and younger time intervals. In this process there has been no use of a track annealing model. 450

Appendix B
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 track lengths are normalized by division with the number of tracks of each histogram. They represent the probability of an 460 induced track appearing in a certain length interval. Interpolated filters are obtained using a mesh of 0.2 in both directions (bins and mean track lengths), Fig. B1.

Appendix C 470
Plots of mean track length versus surface track density for tracks generated by neutron radiation and annealed in the laboratory are given by Green et al. (1986). A logarithmic expression approximates data, Fig. C1: where 0 ⁄ is the surface track density of partially annealed tracks relatively that of unannealed tracks.
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 490 microscope receives light. The criteria by Gleadow et al (2019) imply that a short track, 6 µm long, is accepted for counting if the angle to the horizontal is less than 10°. A long track, 18 µm, must 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 495 selection criteria (Gleadow et al., 2019) the set of age equations is then: The age of the oldest track is obtained summing all the time intervals: (D4) 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 E 505
The inversion principle presented in the main text is for 2D-track length histograms ignoring track angles to the c-axis.
However, the inversion principle is also applicable to data that includes the angles. The measured data vector in Eq. (24) 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. E1 with a synthetic 3D-histogram as data obtained from two annealing experiments. 510 Figure E1: 3D-histogram of track lengths versus angles derived by combining two annealing experiments (Barbarand et al., 2003). Grayscale is the number of tracks in each bin.

515
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 Barbarand et al. (2003). Each column of is a filter. Two of them are shown in Fig. E2.