Statistical modelling of annealing kinetics of fission tracks in zircon; Reassessment of laboratory experiments

Chemical Geology 236 (2007) 75 – 91 www.elsevier.com/locate/chemgeo

Statistical modelling of annealing kinetics of fission tracks in zircon; Reassessment of laboratory experiments Ryuji Yamada a,⁎, Masaki Murakami b , Takahiro Tagami c a Department of Earth Science, University College London, WC1 6BT, UK Department of Earth and Planetary Science, University of Tokyo, Tokyo 113-0033, Japan Division of Earth and Planetary Sciences, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan b

c

Received 5 December 2005; received in revised form 6 September 2006; accepted 6 September 2006

Abstract We remodel the annealing kinetics of fission tracks (FTs) in zircon to fit data from 43 laboratory experiments in three earlier publications using spontaneous FT lengths in zircon from the Nisatai Dacite. Heating temperature and time are in the ranges of 350–912 °C and 10− 3–104 h, respectively. The data sets illustrate a transition in annealing process at a certain temperature zone. In order to express the annealing behavior involving the possible transitional process with improved goodness of fit, two types of models are constructed; the hybrid-linear model that consists of the fanning-linear model at high temperature and the parallel-linear model at low temperature, connected with a transitional temperature zone, and the parallel-curvilinear model that approximates smoothly the transition in the annealing process. The validity of the extrapolation of the model to geological time scale is confirmed by comparison with some deep borehole core data. The length estimates by the hybrid-linear model agree well with the borehole data, whilst the estimates by the parallel-curvilinear model show a significant disagreement with highly annealed borehole samples although this model provides good agreement with the experimental data. For the hybrid-linear model, estimates of the partial annealing zone corresponding to isothermal heating for 106 and 107 years are 281–352 °C and 262–330 °C where the top and bottom of the partial annealing zone are defined by the mean length reduction ratio of 0.4 and 0.8, respectively. © 2006 Elsevier B.V. All rights reserved. Keywords: Fission track; Annealing kinetics; NST zircon; Laboratory experiment; Thermal history analysis

1. Introduction Zircon fission-track (FT) analysis is applied in various geological settings to reveal their thermotectonic evolution combined with apatite FT analysis and/or other methods of low-temperature thermochro⁎ Corresponding author. Permanent address: National Research Institute for Earth Science and Disaster Prevention, Earthquake Research Department, Tennodai 3-1, Tsukuba 305-0006, Japan. Tel.: +81 29 863 7615; fax: +81 29 863 7610. E-mail address: [email protected] (R. Yamada). 0009-2541/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.chemgeo.2006.09.002

nometry. In order to obtain quantitative thermal history constraints from zircon FT data, a rigorous quantitative understanding of the annealing behavior of FT in zircon both in geological and laboratory conditions is required. A series of laboratory studies of thermal annealing characteristics of zircon were performed using confined spontaneous track length (Tagami et al., 1990; Hasebe et al., 1994; Yamada et al., 1995a,b; Yamada et al., 2003). Yamada et al. (1995b) first published a practical annealing model of zircon FT based on laboratory experiments using spontaneous tracks from Nisatai Dacite (known as NST; Tagami et al., 1995) zircons

76

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

annealed for ca. 4.5 min–1000 h. Galbraith and Laslett (1997) re-analysed the data of Yamada et al. (1995b) with an improved statistical modelling scheme developed by them earlier (Laslett and Galbraith, 1996). Tagami et al. (1998) subsequently revised their annealing kinetics after adding two 10 000 h annealing data. Recently, Murakami et al. (2006) performed short term (i.e. several to 100 s) annealing experiments heating NST zircon samples in a graphite furnace. Their experiments were aimed at expanding the applicable field of the zircon FT method to cover short term and high temperature geological phenomena such as analyses of thermal anomaly or activity history in and around fault zones. The time and temperature ranges estimated for the formation of pseudotachylites (ca. 5 s and 1000 °C, e.g. Otsuki et al., 2003), which are believed to result from frictional heat at fault zones, correspond to those of zircon FT shortening predicted by previous kinetic models. In this paper, we present a set of new annealing kinetic models for zircon FT using data from three previous laboratory annealing experiments (Yamada et al., 1995b; Tagami et al., 1998; Murakami et al., 2006). It is reasonable to combine their data because these experiments were all performed using the same NST zircon samples and experimental criteria. However, we need to inspect the adequacy of combining the data from different experiments. Data were carefully reexamined to minimise the uncertainty associated with different experimental conditions, such as annealing temperature monitoring, track length measurement systems, and the observers. Consequently, we used data sets with temperature and time ranges of 350–912 °C and 10− 3– 104 h, respectively. Statistical modelling was performed with the scheme developed by Laslett and Galbraith (1996). We took a two-step approach for the modelling of annealing kinetics of zircon FT. First, we searched for some fixed parameters to describe the forms of each model's geometry. Both linear and curvilinear models were investigated (Crowley et al., 1991). Second, the statistical modelling method (Laslett and Galbraith, 1996; Galbraith and Laslett, 1997) was applied in order to find a set of parameters to best represent the data. As a result, we propose two types of models that fit the data well: a hybrid-linear model that consists of a fanninglinear part and a parallel-linear part, and a more simple parallel-curvilinear model. Both proposed models are chosen to take into account the possible transition of annealing process implied by Murakami et al. (2006). The validity of extrapolation of these models to geological time scales and the usage of models are discussed.

2. Data The data employed for model analysis are listed in Table 1. They are taken from three published sets of laboratory annealing experiments and spontaneous horizontal confined track length data, using the same NST zircon sample and etching procedure but different heating times (4.5 min–1000 h, Yamada et al., 1995b; 4.5 min– 10 000 h, Tagami et al., 1998; ca. 3.6–100 s, Murakami et al., 2006). NST has a K–Ar biotite age of 21.0 ± 0.3 Ma (2σ) and a spontaneous track density (ρs) of ca. 4 · 106 cm− 2 (Tagami et al., 1995). These spontaneous tracks have experienced no annealing since their formation, confirmed by the comparison with induced track length (Yamada et al., 1995b). Spontaneous tracks are used throughout annealingmodel studies with the NST zircon although induced ones are used for most annealing studies on apatite. The differences between spontaneous and induced tracks are the freshness of the tracks and the effect of accumulated Table 1 Laboratory annealing data of spontaneous tracks in NST zircon t (h)

T (°C)

Control 0.0010 0.0010 0.0011 0.0010 0.0011 0.0029 0.0030 0.0030 0.0029 0.0028 0.0029 0.028 0.028 0.028 0.028 0.028 0.028

N

M (μm)

sd (μm)

58 11.04 0.67 599 700 751 800 912 599 649 700 750 800 858 549 649 698 750 805 858

37 38 38 38 4 31 35 30 35 29 3 38 14 25 26 6 1

10.72 9.90 8.92 8.50 5.86 10.57 10.13 9.40 8.35 7.32 5.71 10.72 9.34 8.74 7.78 5.23 4.05

0.56 0.72 1.36 1.29 0.97 0.86 1.18 1.00 1.18 1.31 0.27 0.76 0.64 0.88 1.01 1.19 –

t (h)

T (°C)

N

M (μm)

sd (μm)

1 1 1 1 1 1 1 1 11 11 11 11 11 100 100 100 100 100 1000 1000 1000 1000 1000 10 000 10 000

395 446 500 550 599 650 696 748 397 499 548 598 649 350 449 501 549 599 398 450 498 535 546 398 498

31 36 32 38 37 31 8 0 28 29 30 30 0 35 32 35 26 3 33 41 40 10 0 38 25

10.85 10.59 10.62 10.15 9.47 8.16 6.28 – 10.89 10.16 9.22 7.72 – 10.79 10.13 9.22 8.10 6.73 10.27 9.47 8.59 7.11 – 9.71 7.18

0.91 0.53 0.75 0.59 0.63 0.92 1.66 – 0.65 0.59 0.69 0.70 – 0.71 0.76 0.75 0.71 1.63 0.55 0.72 0.73 1.16 – 0.88 1.51

Measurement data used for model calculation in this study are selectively cited. Original data are from Yamada et al. (1995b), Tagami et al. (1998), and Murakami et al. (2006). t = heating time; T = heating temperature; N = number of TINTs above 60° to crystallographic c-axis; M = mean length of TINTs N 60° to the crystallographic c-axis, sd = standard deviation of M.

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

radiation damage. Freshness of the track is important for apatite annealing because tracks undergo slight annealing at ambient surface temperatures. However, there is practically no annealing in zircon at such temperatures because of the increased resistance of annealing compared to apatite; thus there is no significant difference in length between spontaneous and induced tracks (e.g., Hasebe et al., 1994). The merit in using spontaneous tracks for annealing experiments on which reference kinetic models are based is that the condition of samples is closer to that of unknown natural samples in terms of accumulated radiation damage. It controls etching behavior (Gleadow et al., 1976; Gleadow, 1981), and possibly the annealing behavior (discussed below). Because NST zircon has intermediate properties with respect to age and spontaneous track density, kinetic models in this study may stay applicable as long as the same etching conditions and measurement criteria are applied. Of the measurements of confined track length, tracks-in-track (TINTs; Lal et al., 1969) were selectively used for data processing because of the large variation in the proportion of TINCLEs in each data quoted from three articles. Murakami et al. (2006) reported that many more tracks-in-cleavage (TINCLEs; Laslett et al., 1982) were found (up to ca. 50%) for their higher temperature heating experiment than for other lower temperature experiments. This may be due to the sudden change in the volume of the zircon crystals caused by the rapid

77

temperature increase and decrease. Yamada et al. (1995a) confirmed experimentally no significant difference in mean length between TINCLEs and TINTs in annealed zircon when the proportion of TINCLEs was low (up to ca. 10%). If the proportion of TINCLEs is very high, however, the effect of possible biases could be emphasized, such as anisotropic sampling due to the preferential occurrence of cleavages parallel to the caxes in zircons, or the length bias (e.g., Laslett et al., 1982) caused by many cleavages that act as large etching hosts. In order to avoid possible effects when the TINCLEs' proportion in measurement data is highly uneven, TINCLEs were excluded from all measurements in the three articles (Yamada et al., 1995b; Tagami et al., 1998; Murakami et al., 2006) and the mean lengths were recalculated in Table 1. Fig. 1 shows a plot of all measurement data to be fitted, together with fitting curves by new models documented below. The quality of measurements is ensured by the strict application of the same etching conditions and measurement criteria as follows; zircons were etched until the width of the surface tracks perpendicular to the c-axis became 2 μm using the NaOH–KOH eutectic etchant at 248 ± 1 °C, and the horizontal confined tracks with widths of 1 ± 0.5 μm were selectively measured to minimise the overetching bias (Yamada et al., 1995a). About 50 tracks were measured for each experiment, and we applied the

Fig. 1. Laboratory annealing data (from Table 1) and curves of equal experimental duration from two investigated kinetic models in a mean length versus temperature diagram. The fitted hybrid-linear (A) and parallel-curvilinear (B) models correspond to the parameters from Table 3 (#10⁎ and #16⁎ for (A), and #18 for (B)). Dashed lines in (A) represent the fanning-linear model. The measured mean lengths are plotted against temperatures for 0.001, 0.003, 0.03, 1, 11, 100, 1000 and 10 000 h experiments. A control sample is plotted at a nominal temperature of 20 °C.

78

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

orientation criteria to distinguish the tracks at an angle of 60° or higher to the c-axis for reasons given by Yamada et al. (1995a). For the highly annealed samples, this distinction effectively prevents the erroneous estimate of the length that can be caused by the anisotropy in the number (Hasebe et al., 1994) and the length (Yamada et al., 1995a) of confined tracks, although the reduction in number of analysable confined tracks is generally small (Yamada et al., 1995a). The sample heating technique, the observer and the track length measurement system including microscopes in Murakami et al. (2006) were different from the earlier studies. For sample heating and temperature monitoring, Yamada et al. (1995b) and Tagami et al. (1998) used a muffle furnace and sheathed thermocouples whereas Murakami et al. (2006) used a graphite furnace and infrared thermometers, due to the difference in target heating time and temperature. These pairs of furnaces and temperature gauges were calibrated with an industrial standard thermocouple (Murakami et al., 2006). From these experiments, Murakami et al. (2006) found that the response time of the metal sheathed thermocouples (ϕ 3.2 mm) to the change in ambient temperatures of 500 or 600 °C was over a couple of minutes, and therefore the annealing temperature variation measured with sheathed thermocouples was less reliable for heating times of several minutes. Accordingly, we excluded the 4.5-minute experiments

of Yamada et al. (1995b) from this analysis and consequently 43 data were included. The effect of the difference in track length measuring systems was assessed by comparing the measurements of FTs with various lengths by two observers using different systems. Fig. 2(A) shows the length dependence indicated by the correspondence of measurements obtained by the newer Nikon® system used in Murakami et al. (2006; observer MM) versus those by the older Hamamatsu® system used in Yamada et al. (1995b) and Tagami et al. (1998; RY). For both observations, lengths and crystallographic orientations to the c-axes were measured on the digital images of the confined tracks, with the precisions of ±∼ 0.1 μm and ∼ 1° for track lengths and orientations, respectively (Yamada et al., 1995a; Murakami et al., 2006). Fig. 2(B) shows the orientation dependence indicated by the difference (open circles) and the ratio (solid circles) between measurements obtained by the two systems. Both indicate that measurements are not biased for track length at any orientation due to the difference in measurement systems. Close agreement is seen between ∼ 4 and 13 μm, i.e., the entire range of measured FT lengths in Fig. 2(A) and (B). The maximum absolute value, the mean and the standard deviation of the difference in measurements between two systems are 0.47, − 0.01, and 0.23 μm, respectively. Because there is no evidence of bias between the two measurements and

Fig. 2. Comparison of FT measurements in zircon with various lengths and orientation using different measurement systems. (A) Comparison between measurements obtained by the new and old systems. Measurements were performed for two sets of 14 FTs in zircons heated for 1 h at 700 °C (indicated by triangles) or 600 °C (squares), and one set of 10 FTs un-annealed (circles); 38 FTs in total. A dashed line shows a 1:1 correlation between new and old measurement systems. (B) Difference values (open symbols) and the ratios (solid symbols) between measurements by the new and the old systems. The standard deviation of the difference is shown by a gray band.

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

the discrepancies are comparable to the normal uncertainty associated with measurements by a single observer (e.g., Yamada et al., 1995a), the effect of the measurements by different observers will be ignored for the analyses hereafter. 3. Analysis 3.1. Modelling approach FT annealing has been described on an empirical basis by analogy with other diffusion processes (e.g. Green et al., 1986, 1988). On an Arrhenius diagram, the fading contours that represent different mean lengths, or mean length reductions, are plotted in 1 / T − log t space, where T is absolute temperature and t is time. For a given mean length, a contour can be described by the fundamental equation logt ¼ A þ BdFðT Þ

ð1Þ

where A is the intercept, B is the slope, and F(T) is the transform function of the temperature. Crowley et al. (1991) categorized the geometry of contours based on the function of B and F(1 / T) as follows. Parallel and fanning: The fading contours are parallel if the slope B is constant (i.e., the same for all mean lengths) whereas the contours are non-parallel or fanning if B varies with mean track length. Fanning contours are presumed to intersect at a common crossover point (log tc, 1 /Tc) where Tc is called the critical temperature. Linear and curvilinear: The fading contours are linear in shape on the Arrhenius diagram if F(T) = 1 / T, whereas the fading contours are curvilinear if F(T) = log(1 / T). Fully parameterized model forms for the combination of these geometries are as follows: Fanning-linear model: gðrÞ ¼ c0 þ c1

logt−logtc ; 1=T −1=Tc

ð2:1Þ

Parallel-linear model: gðrÞ ¼ c0 þ c1 ðlogt−Bd 1=T Þ;

ð2:2Þ

Fanning-curvilinear model: gðrÞ ¼ c0 þ c1

logt−logtc ; log1=T −log1=Tc

ð2:3Þ

79

Parallel-curvilinear model: gðrÞ ¼ c0 þ c1 ðlogt−Bd log1=T Þ;

ð2:4Þ

where g(r) is an appropriate function to make a statistical choice, and the mean length reduction ratio r = μ/μmax, where μ = mean length at the annealing condition of (t, T), and μmax = maximum mean length (normally measured in an un-annealed sample). The constant slope B in a parallel-linear model is normally described in the form of E / k with a single activation energy E and k being the Boltzmann's constant. However, most of the widely used FT annealing kinetic models prefer equation forms of the fanning-linear model (e.g. Laslett et al., 1987; Yamada et al., 1995b; Galbraith and Laslett, 1997; Tagami et al., 1998). Presumably a single activation energy cannot adequately describe the complexity in FT annealing recognized as the combined process of shrinking and segmentation along latent tracks, which has been confirmed by direct observation of annealed FTs using transmission electron microscope (Paul and Fitzgerald, 1992) and a step etching experiment using annealed FTs (Green et al., 1986; Yamada et al., 1995a). For the g(r) function of Eqs. (2.1)–(2.4), previous workers (Laslett et al., 1987; Crowley et al., 1991) have used a double Box and Cox type transformation (Box and Cox, 1964) of the form gðrÞ ¼ ððð1−rb Þ=bÞa −1Þ=a

ð3Þ

where 0 b r b 1. The Box–Cox transformation can be used for approximating the data to a normal distribution with arbitrary parameters that determine the shapes of the transformation function to allow the process capability to be easily determined. The proposed models in the early literature are described as special (c.f. Laslett et al., 1987) or intermediate cases with two fitting parameters α and β. During the fitting procedure, however, the modelling results tend to be unstable because the estimates of α and β may suffer from a local minima problem with multimodal likelihood under some conditions (Laslett and Galbraith, 1996), i.e. a subtle change in these two parameters may strongly affect the estimate of the other parameters. Laslett and Galbraith (1996) suggested fixing these parameters at integer values, because there was a small difference in the best log-likelihood between the cases of α and β estimated and fixed at appropriate values. Another good reason for fixing them is to avoid over-parameterization. Laslett and Galbraith (1996) effectively fixed α = 0 (equivalent to a log transformation) and at the same time treated the maximum mean

80

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

length (μmax) as another parameter to be estimated. Information on the experimental conditions where no confined tracks were found was also taken into account. In addition, they introduced re-parameterization at two fitting points to overcome the problem of non-convergence when estimating the crossover points, and observed a more stable and faster convergence. We took a two-step approach for the modelling of the annealing kinetics of zircon FT in this study. First, we sought combinations of integers for the α and β parameters in the g(r) transformation function to describe the forms of each model's geometry by maximizing the likelihood for each combination. Second, the statistical modelling method (Laslett and Galbraith, 1996; Galbraith and Laslett, 1997) was applied after proper integers for α and β were found. In order to find a set of parameters to maximize the likelihood functions, we performed grid search computation by means of narrowing and modifying the search area stepwise until the relevant digits of all parameters did no longer change. 3.2. Annealing-model forms 3.2.1. Linear model We first applied the fanning model expressed by Eq. (2.1) because it allows to express more general contour forms, and exceeded the parallel model in terms of degree of agreement with experimental data in most previous models. The goodness of fit was evaluated by the maximum log-likelihood LL (Laslett and Galbraith, 1996),
 nþ1 SS LL ¼ − log 2 nþ1 m X þ logð−g Vðri ÞÞ−ðn þ 1Þloglmax ; ð4:1Þ i¼0 Ni N0

where SS ¼

m X

ðgðri Þ−f ðti ; Ti ÞÞ2 ;

ð4:2Þ

i¼0 Ni N0

and f ðti ; Ti Þ ¼ c0 þ c1

logti −logtc : 1=Ti −1=Tc

ð4:3Þ

Here i= 0 corresponds to the control sample, assumed to have T0 = 293.15 K and t0 = 100 h; m is the total number of experiments excluding the control; N is the number of TINTs above 60° to crystallographic c-axis, n is the number of these for which Ni N 0, and the apostrophe in Eq. (4.1) denotes the differentiation operation. Combinations of α

and β parameters proposed previously are [α= 0, β =3] for apatite and zircon (Laslett and Galbraith, 1996; Galbraith and Laslett, 1997), and [α =0, β = 1] for zircon (Tagami et al., 1998). Laslett and Galbraith (1996) documented that setting α = 0 is reasonable to avoid the potential multimodality problems of α = −1, and the identifiability problems of α = 1. Therefore, we started the analysis with the case of [α = 0, β= 3] and changed the integer b until the maximum value of LL was obtained. Combinations of [α = 0, β = − 2] and [α = 0, β = − 1] give comparable maximum scores of LL N 63 (Table 2). However, the μmax is estimated identically as 11.04 μm when the maximum LL is obtained. This is because μmax was constrained to be longer than 11.04, the mean length of the control sample (Table 1), to avoid divergence when r ≥ 1. This boundary condition for μmax was applied throughout the modelling. In order to inspect the effect of the control sample, analysis was repeated without the control. LL (Eq. (4.1)) and SS (Eq. (4.2)) in this case need to start the summation at i= 1 instead of 0. The LL score for the combination of [α =0, β = −2] is the highest with those for the sets [α = 0, β = −1] and [α= 0, β = −3] showing comparable LL values (Table 2). 3.2.2. Curvilinear model Because curvilinear models have not been commonly used to fit annealing kinetics of FTs, we need to estimate α and β parameters before finding good integer combinations. The parallel model was investigated first because it has fewer parameters to be estimated than the fanning model. The goodness of fit was evaluated by the log-likelihood LL given by Eq. (4.1) where f ðti ; Ti Þ ¼ C0p þ C1p logti þ C2p log1=Ti

ð5Þ

This is derived from a simple transformation of Eq. (2.4) where C0p = c0, C1p = c1 and C2p = − c1 · B. A combination of [α = − 0.010, β = − 0.801; LL ∼ 68] was estimated where the control was included (Table 2). This suggests that fixing α = 0 is also acceptable for the curvilinear model and integers of − 1, 0 are strong candidates for β. Therefore, we attempted three combinations of [α = 0, β = − 1], [α = 0, β = 0] and [α = 0, β = 1]. Among these three, the set [α = 0, β = 0; LL ∼ 66.5] yielded the highest LL score but lower by ∼ 1.5 from that estimated with non-integer values for α and β. In that case, the μmax was estimated as 11.04 when the maximum LL was obtained. In order to observe the effect of the control, analysis was repeated without the control, and a combination of [α = 0.098, β = 0.191; LL ∼ 56.6] was estimated. This suggests again that fixing α = 0 is acceptable and integers of 0 and 1 are

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

81

Table 2 Results of model fitting for determining g(r) functions • Fanning-linear model μmax

c0

c1

log(tc)

1 / Tc

LL

Control included 0 −3 0 −2 0 −1 0 1 0 3

α

β

11.04 11.04 11.04 11.04 11.04

− 18.94 − 20.21 − 21.85 − 26.14 − 32.85

0.000469 0.000502 0.000546 0.000666 0.000856

− 55.07 − 64.65 − 77.12 − 111.20 − 166.00

− 0.000307 − 0.000572 − 0.000920 − 0.001890 − 0.003458

62.03 63.43 63.33 56.99 45.09

Control excluded 0 −3 0 −2 0 −1 0 1 0 3

11.06 11.12 11.23 11.62 12.68

− 11.49 − 10.75 −9.36 −7.00 −4.24

0.000272 0.000251 0.000213 0.000148 0.000073

− 32.39 − 34.07 − 34.35 − 37.56 − 41.10

0.000287 0.000238 0.000226 0.000134 0.000030

58.75 59.06 58.39 53.82 47.15

• Parallel-curvilinear model α

β

μmax

C0p

C1p

C2p

LL

Control included − 0.010 0 0 0

− 0.801 −1 0 1

11.04 11.04 11.04 11.04

− 75.55 − 72.43 − 69.81 − 67.43

0.2566 0.2457 0.2375 0.2298

10.921 10.467 10.059 9.688

67.95 65.83 66.54 58.18

Control excluded 0.098 0 0 0 0

0.191 −2 −1 0 1

11.04 11.04 11.04 11.09 11.19

− 55.64 − 82.65 − 75.65 − 65.82 − 54.79

0.1879 0.2791 0.2565 0.2236 0.1859

7.989 12.000 10.943 9.476 7.837

56.55 53.86 56.46 56.03 53.96

LL

• Fanning-curvilinear model α

β

μmax

c0

c1

log(tc)

log(1 / Tc)

Control included 0 −1 0 0 0 1

11.04 11.04 11.04

160.20 85.89 84.14

− 3.804 − 2.076 − 2.047

691.1 401.3 405.3

9.482 2.741 2.880

60.46 61.62 61.53

Control excluded 0 −1 0 0 0 1

11.04 11.06 11.22

50.49 41.97 91.78

− 1.220 − 1.031 − 2.214

207.1 192.8 522.5

− 1.914 − 2.215 5.604

52.70 53.07 53.08

candidates for β, and we chose the same three parameter combinations as with the control included. Under these conditions the set [α = 0, β = − 1; LL ∼ 56.5] yielded an LL score close to that obtained when α and β were estimated as non-integers. The μmax for [α = 0, β = − 1] was estimated as 11.04. Without a given μmax, a combination of [α = 0, β = 0] yielded an LL score lower by ∼ 0.5 than that estimated, and its μmax estimate was 11.09. Based on the above results, a combination of [α = 0, β = 0] is expected to best describe the form of the g(r) function for the parallel-curvilinear model.

Next, we analysed the fanning model by fixing α and β at the values of candidates for parallel models. The goodness of fit was evaluated by the maximum loglikelihood LL given by Eq. (4.1) where f ðti ; Ti Þ ¼ c0 þ c1

logti −logtc : log1=Ti −log1=Tc

ð6Þ

Three combinations, namely [α = 0, β = − 1], [α = 0, β = 0] and [α = 0, β = 1] were chosen in analogy to the analyses of the parallel models. Based on the LL scores,

82

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

the set [α = 0, β = 0] is slightly preferred among the three combinations with the control included. In order to observe the effect of the control, analysis was repeated without the control. The sets [α = 0, β = 0] and [α = 0, β = 1] yielded comparable scores of LL ∼ 53, and their estimates of μmax were N 11.04 (Table 2). Based on the above results, we also expect the combination of [α = 0, β = 0] to determine the form of the g(r) function for the fanning-curvilinear model. Fanning and parallel models can be compared on the basis of the position of the crossover points estimated for the fanning models. Estimates for log tc and log 1 / Tc parameters indicate that the crossover points are located very far from the distribution of measurement data, and the shapes of the contours for these good models are very close to parallel. This means that the fanning model may not improve the goodness of the fit significantly compared with the parallel model whilst a large number of parameters are to be estimated. In fact, the best LL score for the parallel model was better than the fanning model. Therefore, the parallel model is preferable and the combination [α = 0, β = 0] is adequate to determine the g(r) function for the curvilinear type transformation of heating temperature. 3.3. Fitting models In order to find suitable parameters, we followed the statistical modelling procedures of Laslett and Galbraith (1996) and Galbraith and Laslett (1997). 3.3.1. Fanning-linear model In the previous section, it was indicated that the function g(r) with the combination of [α = 0, β = −2] is suitable for the fanning-linear model. Eq. (2.1) can be rewritten as

k logt−logtc lf l ðt; T Þ ¼ lmax 1 þ exp c0 þ c1 ; 1=T −1=Tc ð7:1Þ where α = 0 and β is negative, and λ = 1 /β. If β is positive, Eq. (7.1) needs to be slightly modified as

k logt−logtc : lf l ðt; T Þ ¼ lmax 1−exp c0 þ c1 1=T −1=Tc ð7:2Þ Eq. (7.2) is the same form as Eq. (5) of Laslett and Galbraith (1996). For the set [α = 0, β = 0], another transformation is required, but we do not deal with it in this study.

The five parameters μmax, c0, c1, Tc and tc are estimated by the maximum likelihood from the annealing data. In order to obtain a stable numerical procedure, they are reparameterized at two fitting points (log ta, 1 /Ta) and (log tb, 1/Tb) as logta −logtc ; 1=Ta −1=Tc

k logta −logtc la ¼ lmax 1Fexp c0 þ c1 1=Ta −1=Tc

ð8:1Þ

logtb −logtc ; 1=Tb −1=Tc

k logtb −logtc lb ¼ lmax 1Fexp c0 þ c1 1=Tb −1=Tc

ð8:2Þ

ba ¼

bb ¼

where ba and bb are the slopes and μa and μb are the mean lengths of the contour lines on the Arrhenius plots that pass through a pair of fitting points. The sign of exponential term in the equations of μa and μb depends on the sign of β parameter in the same way as Eqs. (7.1) and (7.2). We selected fitting points such that μa ∼ 7.5 and μb ∼ 10.5 for all cases; they are located within the domain of the data, and therefore ba and bb should be well estimated. We confirmed that the choice of μa and μb does not affect the results significantly by computation tests beforehand. For a single annealing experiment at temperature T and time t that produces a sample of N tracks, the observed mean length M is assumed to be normally distributed with a mean μ and a variance given by r2 þ

1 mw ðlÞ N

ð9Þ

The parameter σ represents the “between experiment” standard deviation, which is estimated along with other parameters. The second term in Eq. (9) represents “within experiment” sampling variation and is treated as a known component of variance in the main fitting procedure. We reassessed νw(μ) by linear fitting to the scatter plot of the logarithm of standard deviation sd against log M using the present data set (Fig. 3), because this simple form would be sufficient enough to fit the present data set. The function of νw(μ) is given by mw ðlÞ ¼ expð−1:40logl þ 2:93Þ:

ð10Þ

In some experiments, annealing was so advanced that no confined tracks could be observed. This indicates that μ is likely to be small under such conditions. In order to use this information in the likelihood function, we use

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

83

where δi = 1 if tracks are observed and δi = 0 if no tracks are observed. For the control sample, assumed to have T0 = 293.15 K and t0 = 100 h, we have δ0 = 1, p(μ0) = 1 and log-likelihood 
 M0 −l0 LL0 ¼ log 1−U pffiffiffiffiffi ð14  2Þ m0 where Φ is the standard normal distribution function. The log-likelihood function LL for the control and m experiments is LL ¼

m X

LLi :

ð14  3Þ

i¼0

Fig. 3. Relation between the within-mount standard deviation and mean length of confined tracks above 60° to crystallographic c-axis. The solid circles are the points (log M, log sd), one for each annealing experiment. The dashed line is a fitted regression line to the data weighted by the number of measured TINTs. The outlying datum at the lower end of the diagram is included for the fitting.

the function p(μ) to express the probability that confined tracks are observed, given by pðlÞ ¼

e2ðl−5Þ 1−e2ðl−5Þ

ð11Þ

which equals 0.01, 0.5 and 0.99 when μ = 2.7, 5.0 and 7.3, respectively. In order to fit models to data, the log-likelihood is given in the same technical scheme documented in Laslett and Galbraith (1996) and Galbraith and Laslett (1997). For an experiment i without visible tracks, the likelihood is 1−pðli Þ

ð12Þ

where μi is the calculated mean length for experiment i, given by the equation for a relevant model (e.g., Eq. (7.1)). For an experiment i with visible tracks, the likelihood is proportional to ( ) 1 ðMi −li Þ2 pðli Þ pffiffiffiffiffiffiffiffiffi exp − : ð13Þ 2mi 2pmi Hence the log-likelihood for the data from experiment i = 1, 2,…, m is LLi ¼ ð1−d( i Þlogf1−pðli Þg þ di

1 ðMi −li Þ2 logpðli Þ− logmi − 2 2mi

) ð14  1Þ

The goodness of the fit is indicated by the Schwarz's Bayesian information criterion (BIC; Schwarz, 1978) in order to take the difference in numbers of parameters estimated into account, given by BIC ¼ −2LL þ logN d p;

ð15Þ

where N is the number of experiments analysed and p is the number of parameters estimated. The better the model, the lower the BIC. Analytical results for the fanning-linear model fitting are listed in the top section of Table 3 (#1–#4), corresponding to the cases of β = − 1, − 2, 3, and 1, respectively. The case of β = 3 was tried by Galbraith and Laslett (1997). Fig. 4 illustrates the plots of standardized residuals, which are calculated as Mi −lî pffiffiffiffi mî

ð16Þ

against annealing temperature and time. Here μˆ i is the fitted mean and ˆνi = σˆ 2 + νw(μˆi) / ni is the estimated variance for experiment i. Based on the LL obtained from analysis with all data including experiments where no tracks were found, the case with λ = 1 / 3 (β = 3; #3 in Table 3) is the best (i.e. lowest BIC) but the difference according to λ is not significant. For the best model, however, the BIC (− 30.99) was higher than that from the analysis by Tagami et al. (1998) (BIC = − 44.22, calculated for reference purpose) and the residual diagram appears structured in the residual plot in Fig. 4(A), indicating the inadequate choice of the model for the present measurement data set. We tried other integers for the β parameter and confirmed that a variation in this parameter did not improve the situation. The reason why the classical linear model could not fit well is implied by the alignment of the data points on the

84

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

Table 3 Analytical results of model fitting Linear model analysis #

Data

λ

Model N

log log 1000/ 1000/ σb ta (h) tb (h) Ta (K) Tb (K)

μmax

c0, c0p

1000c1, bp c1p

log tc

1000/ Tc

LL

BIC

− 57.41 − 54.22 − 140.95 − 94.33

−0.376 −0.290 −2.596 −1.351

26.64 26.51 26.78 25.09

− 30.71 − 30.45 − 30.99 − 27.60

[1] Analysis with all data, no-track data included 1 All − 1.000 fan 43 0.0 0.0 2 All − 0.500 fan 43 0.0 0.0 3 All 0.333 fan 43 0.0 0.0 4 All 1.000 fan 43 0.0 0.0

1.07 1.07 1.07 1.07

1.30 1.30 1.30 1.30

0.127 0.127 0.164 0.140

11.10 10.98 11.69 11.33

−15.71 −16.73 −15.31 −17.56

0.378 0.424 0.390 0.423

[2] Longer time 5 All 6 N0.002 h 7 N0.02 h 8 N0.2 h

1.07 1.06 1.16 1.20

1.30 1.30 1.42 1.45

0.164 0.167 0.047 0.000

11.92 12.00 11.90 11.93

−9.33 −6.46 −4.45 −2.13

0.233 0.157 0.112 0.047

− 95.63 −1.391 31.81 − 41.48 − 66.15 −0.611 27.69 − 34.05 − 36.02 0.097 30.56 − 40.91 − 9.35 0.836 28.93 − 39.04

[3] Lower temperature domain, no-track data excluded 9 b775 °C 0.333 fan 34 5.0 − 4.0 1.20 10 b725 °C 0.333 fan 31 5.0 − 3.0 1.20 11 b675 °C 0.333 fan 27 8.5 0.0 1.30

1.20 1.20 1.30

0.143 11.87 0.067 11.76 0.000 11.70

−7.40 −6.20 −4.98

0.183 0.156 0.126

− 71.17 −0.765 32.72 − 44.28 − 51.75 −0.307 35.48 − 50.35 − 36.47 0.082 32.81 − 45.84

[4] Higher temperature domain, no-track data excluded 12 N425 °C 0.333 fan 35 1.0 − 7.0 1.10 13 N475 °C 0.333 fan 32 1.0 − 7.0 1.10 14 N525 °C 0.333 fan 27 2.0 − 6.0 1.12 15 N425 °C 0.333 par 32 1.0 − 7.0 1.10 16 N475 °C 0.333 par 32 1.0 − 7.0 1.10 17 N525 °C 0.333 par 27 2.0 − 6.0 1.12

1.10 1.10 1.12 1.10 1.10 1.12

0.155 0.102 0.104 0.159 0.105 0.123

[5] For hybrid model, no-track data included 10⁎ b725 °C 0.333 fan 33 5.0 − 3.0 1.20 16⁎ N475 °C 0.333 par 35 1.0 − 7.0 1.10

1.20 1.10

0.063 11.63 0.146 11.63

domain, no-track data excluded 0.333 fan 40 0.0 0.0 0.333 fan 35 0.0 0.0 0.333 fan 29 3.5 3.5 0.333 fan 23 4.6 4.6

12.46 −13.90 0.349 − 177.51 −3.471 27.83 − 34.33 12.82 14.37 − 0.372 213.22 6.505 28.11 − 35.43 12.29 5.17 − 0.137 65.36 2.717 21.84 − 23.91 12.10 3.38 0.088 38 900 27.56 − 37.35 13.18 2.37 0.061 39 400 27.94 − 38.54 13.13 2.42 0.061 40 150 21.11 − 25.75 −6.95 4.35

0.176 0.113 38 760

− 53.84 −0.363 30.92 − 40.87 19.67 − 21.56

[6] Curvilinear model analysis # Data Model N

bp

LL

18

43.00

31.29 − 43.78

All

par

log c1p log 1000/ 1000/ σ μmax c0p ta (h) tb (h) Ta (K) Tb (K) 43 0.0 0.0 − 6.85 − 6.65 0.051 11.17 −63.37 0.212

BIC

Numbers in “#” column are serial numbers of analyses. In “data” column, “All” denotes analyses without dividing domain, and others the sizes of domains. The shapes of models applied for each analysis are indicated in the “model” column: “fan” = fanning model; “par” = parallel model. See text for other acronyms and symbols in the column headings.

Arrhenius diagram: e.g., the solid squares on Fig. 5 are somewhat curved at a temperature near 600 °C or a time near 1 h. This curvature may be related with a change in the annealing process at a certain time or temperature that takes place within the experimental condition (Murakami et al., 2006). 3.3.2. Hybrid model In order to verify whether the model fit improves if a change in the annealing process is taken into account, we sought the time or temperature conditions that divide the analytical domain of measurement data to find the best fanning-linear model that will be combined smoothly with a counterpart model to describe the

annealing behavior in all domain. Considering the high priority in applying models to geological time scale and temperature conditions, we divided the analytical domain by excluding data for the shorter time or higher temperature heating experiments stepwise. We attempted the domain sizes of t N 0.002 h, t N 0.02 h and t N 0.2 h for the time division, and T b 775 °C, T b 725 °C and T b 625 °C for the temperature division, in order to know whether time or temperature governs the transition of the dominance in annealing process. Fitting points are selected so that μa ∼ 7.5 and μb ∼ 10.5, and λ is fixed at 1 / 3 hereafter. The comparison target is an all domain analysis with no-track data excluded (#5 in Table 3). No-track data are excluded here, because

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

85

data from T N 725 °C experiments; 2. Use some common data by gradually excluding data from lower temperature experiments until a good model is found. We took the second option to maximize the amount of data for statistical modelling, and the two models should be smoothly combined for their application to thermal history modelling. Attempted domain sizes to facilitate

Fig. 4. Standardized residuals (calculated according to Eq. (16)) plotted against temperature (left) and time (right). (A) Fanning-linear model of Eq. (7.2). Residual patterns show a structured distribution against temperature and time. (B) Hybrid-linear model of Eq. (19) comprised of the fanning-linear (solid circles) and parallel-curvilinear models (open squares). The degree of structured distribution for mixture of two components appears smaller than in (A). (C) Parallelcurvilinear model of Eq. (24), with no structured distribution.

they only lower the LL and the significance of this varies according to the number of experiments included and excluded in the calculation at different domain sizes. This is a tentative treatment for seeking the best domain; for the final parameter assessment, it was attempted to include all no-track analysis. Results are shown in the second and third sections of Table 3. By dividing the domain by time (#5–8), no other set of conditions yielded a lower BIC than the target. In contrast, by dividing the domain by temperature (#9– 11), all three conditions provided lower BIC than the target, and T b 725 °C had the lowest BIC (#10). It is thus indicated that the goodness of the fit depends more on temperature than on time, although we were first concerned with the dependency on the heating time because we used two different heating techniques and excluded the 4.5-minute data. In search for a counterpart model to fit a higher temperature domain and make the best use of the experimental results, there are two options: 1. Only use

Fig. 5. Arrhenius plots for NST zircon showing the time–temperature conditions for the investigated annealing experiment data set and contour lines of equal length reduction for the fitted linear models (A: with parameters of #10⁎ and #16⁎ in Table 3) and curvilinear model (B: #18 in Table 3). Different degrees of length reduction of experimental data are shown with different symbols. Dotted and dashed lines in (A) express the length reduction contours of the fanning- and parallel-linear models having parameters of #10⁎ and #16⁎ in Table 3, respectively. Solid lines in (A) express the contours of the hybrid model formed by weighted sum of these two models as expressed by Eq. (19). Also shown are borehole data (listed in Table 4), that were not included for model fitting.

86

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

fanning model analysis are for temperature ranges of T N 425 °C, T N 475 °C and T N 525 °C, with no-track data excluded. Results are shown in the fourth section of Table 3 (#12–14). The analysis in the T N 475 °C domain yielded the lowest BIC (#13). Because the crossover point for this analysis is located at an extremely long time and a low temperature (#13 in Table 3), the parallel model may describe the Arrhenius contours well. Subsequently, a parallel-linear model fitting was attempted (#15–17), given by rewriting Eq. (2.2) as lpl ðt; T Þ ¼ lmax ð1−expðc0p þ c1p  ðlogt−bp d1=T ÞÞÞk ;

ð17Þ

where λ is positive. The parameters μmax, c0p, c1p, bp and σ are estimated by the maximum likelihood method from the annealing data. In order to obtain a stable numerical procedure, they are re-parameterized at two fitting points (log ta, 1 / Ta) and (log tb, 1 / Tb) as la ¼ lmax ð1−expðc0p þ c1p ðlogta −bp d1=Ta ÞÞÞk

ð18:1Þ

lb ¼ lmax ð1−expðc0p þ c1p  ðlogtb −bp d1=Tb ÞÞÞk

ð18:2Þ

where μa and μb are the mean lengths of the contour lines on the Arrhenius plots that pass through the fitting points. Fitting points are selected such that μa ∼ 7.5 and μb ∼ 10.5, and λ is fixed at 1 / 3. Other conditions and procedures are the same as those for the fanning-linear model fitting. As a result, the parallel model in the T N 475 °C domain yielded a lower BIC (#16) than the fanning-linear model. Here we build a hybrid model that consists of the fanning-linear model for the T b 725 °C domain (μ¯ f l) and the parallel-linear model T N 475 °C (μ¯ pl). The two elemental models are combined smoothly with the transitional temperature zone so that the reciprocal function and its derivation can be obtained. This feature is important in order to apply the kinetic models to thermal history modelling that requires various transformations of the functional form of the models. Because the above-described fanning model for the T b 725 °C domain can be applicable to most geological problems and has a lower BIC than all domain analyses, μ¯ fl should be treated as primary and μ¯ pl secondary. Before combining the two models, re-analysis is required with a common μmax and with no-track data included. Fitting points are selected so that μa ∼ 7.5 and μb ∼ 10.5, and λ is fixed at 1 / 3. The results of this re-

analysis are shown in the fifth section of Table 3 (#10⁎ and #16⁎). For the parallel model analysis, μmax is fixed at the estimated value for the fanning model, and thus the number of estimated parameters is four. The hybrid function is expressed by the weighted sum of two components of μ ¯ fl and μ¯ pl, given by ll ðt; T Þ ¼ wðT Þ l¯f l ðt; T Þ þ ð1−wðT ÞÞ l¯pl ðt; T Þ

ð19Þ

The weighting factor varies continuously depending on the distance in temperature from cT, given by wðT Þ ¼ 1=ð1 þ expðqð1=T −cT ÞÞÞ;

ð20Þ

where cT ¼ ð1=Ta þ 1=Tb Þ=2:

ð21Þ

The two models of μ¯ fl and μ¯ pl have an overlap around the region between Ta = 998.15 K (= 725 °C) and Tb = 748.15 K (= 475 °C). The constant q in Eq. (20) can be obtained as − 11 080, by solving wðTa Þ ¼ 1=expðRÞ;

wðTb Þ ¼ 1−1=expðRÞ:

ð22Þ

The constant R defines the smoothness of the combined contours and is set at 2. The contour lines for the hybrid model are drawn on the Arrhenius diagram in Fig. 5(A), together with the μ¯ f l (dotted line) and μ ¯ pl (dashed line) models. The degree of disagreement between μ¯ fl and μ¯ pl is small at small μ (∼ 6 μm), and it becomes larger as μ increases. This diagram indicates that the choice of the weighting factor and some arbitrary parameters to control the crossover range may not introduce significant effects on the shape of the hybrid model. The substantial modification of the hybrid model is the parallelism added to μ¯ f l at the higher temperature domain. Fig. 4(B) shows the residual plot for the analysis of #10⁎ (μ¯ f l: solid circle) and #16⁎ (μ¯ pl: open square) in Table 3 that constitute this hybrid model. Compared with the fanning-linear model obtained without dividing the analytical domain (Fig. 4(A)), the temperature dependence of a structured distribution in residuals is weakened for both the μ¯ f l and μ¯ pl components. The estimates of σ for the fanning and parallel models are ∼ 0.07 and ∼ 0.15, which are smaller than and comparable pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiwith the minimum value of mw ðlÞ= N¯ ¼ 0:67= 35 ¼ 0:14, respectively (where ¯ = an average of N = 35). μ = 11 and N 3.3.3. Parallel-curvilinear model Fitting with curvilinear models can be an alternative to improve the goodness of the fit by smoothly approximating the transition in the annealing process.

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

In the previous section, it was indicated that the function g(r) with the combination of [α = 0, β = 0] is suitable for a parallel-curvilinear model. The limit of the g(r) function at [α = 0, β = 0], is expressed as gðrÞ ¼ logð−logrÞ

ð23Þ

and Eq. (2.4) can be rewritten as lpc ðt; T Þ ¼ lmax expð−expðc0p þ c1p  ðlogt−bp log1=T ÞÞÞ

ð24Þ

The four parameters μmax, c0p, c1p and bp are estimated by maximum likelihood from the annealing data. In order to obtain a stable numerical procedure, they are reparameterized at two fitting points, namely (log ta, 1 / Ta) and (log tb, 1 / Tb) as la ¼ lmax expð−expðc0p þ c1p  logta −bp dc1p log1=Ta ÞÞ

ð25:1Þ

87

with geological constraints. The models proposed in this study are compared here; with regard to the statistical aspect, the parallel-curvilinear model fits well the laboratory annealing data and has the least structured residual plot (Fig. 4C). This model is preferable in terms of the plainness of form because it has a smaller number of parameters to be estimated and good agreement was obtained without complicated data management. With regard to the agreement with geological samples, however, the hybrid-linear model is better than the curvilinear model. Zircon FT data of some ultra-deep borehole samples with subnormal geothermal gradients are compared with the lengths predicted by different models; Vienna Basin boreholes (Tagami et al., 1996), MITI boreholes (Hasebe et al., 2003), KTB HB ultra-deep borehole (Coyle and Wagner, 1996), and Kola Peninsula Superdeep Drillhole SG-3 (Green et al., 1996). Table 4 summarizes the estimated heating time and temperature, length measurements for these core samples, and the predicted data under the presumed heating conditions of each core using the hybrid-linear model and the parallelcurvilinear model in this study, together with those using the fanning-linear models in Tagami et al. (1998) and Rahn et al. (2004). The data of four borehole samples in the

lb ¼ lmax expð−expðc0p þ c1p  logtb −bp dc1p log1=Tb ÞÞ

ð25:2Þ

where μa and μb are the mean lengths of the contour lines on the Arrhenius plots that pass through the fitting points, selected so that μa ∼ 7.5 and μb ∼ 10.5. Other conditions and procedures are the same as those for the parallel-linear model fitting. Analysis was performed with all data including experiments where no tracks were found. The result is shown in the sixth section of Table 3 (#18). The BIC (− 43.78) was lower than for any of the fanning-linear model analyses with all data, and the residual diagram does not appear to be structured (Fig. 4 (C)). The estimate of σ is ∼0.06, which is significantly ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi smaller than the minimum value of mw ðlÞ=35 ¼ 0:14. The parallel-curvilinear model fits better than the fanning-linear model from a statistical viewpoint. The contour lines for the parallel-curvilinear model are drawn into the Arrhenius diagram in Fig. 5(B). 4. Discussion 4.1. Evaluation of the present models The evaluation of the models should be discussed with respect to both statistical aspect and agreement

Table 4 Comparison of mean zircon fission-track length between borehole core samples and model predictions Borehole data

Vienna [1]

MITI [2]

Kola [3]

KTB [4]

Estimated temperature (°C) Estimated time (m.y.) Measured mean length (μm)

197 5 10.6 (0.28)

205 1 9.8 (0.2)

212 250 9.8 (0.20)

255 50 8.9 (–)

10.41 (0.14) 10.66 (0.09) 10.59 (0.10) 10.67 (0.13) 10.08 (0.36)

10.46 (0.14) 10.68 (0.08) 10.62 (0.09) 10.71 (0.13) 10.26 (0.30)

9.79 (0.24) 10.19 (0.20) 10.10 (0.21) 10.07 (0.23) 8.17 (0.91)

9.10 (0.39) 9.53 (0.37) 9.45 (0.37) 9.35 (0.39) 6.90 (1.20)

Model prediction (μm) Fanning-linear (Tc ≠ ∞), NST [5] Fanning-linear (Tc = ∞), zero damage zircons [6] Fanning-linear (Tc = ∞), NST [6] Fanning-linear (Tc ≠ ∞), NST [7] Parallel-curvilinear, NST [7]

Original references for borehole data are; [1] Tagami et al. (1996), [2] Hasebe et al. (2003), [3] Green et al. (1996), [4] Coyle and Wagner (1996); and those for models are [5] Tagami et al. (1998), [6] Rahn et al. (2004), and [7] this study. In the borehole data, standard errors of measurements are given as 2σ in parentheses. In the model prediction, statistical mean lengths are calculated for the predicted lengths of the annealed tracks produced at the end of 100 time intervals of equal duration within the estimated time of each core. Standard deviations for the predicted lengths are given as 2σ in parentheses.

88

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

upper section of Table 4 are also plotted in Fig. 5. Because the contour lines in the Arrhenius plot are drawn for isothermal heating condition in laboratory experiments, they cannot be directly compared with the geological data that consist of tracks accumulated over time–temperature histories. The lower section of Table 4 shows the statistical mean and standard deviation of the lengths of the annealed tracks produced at the end of 100 time intervals of equal duration within the estimated time of each core (an analogue of 100 track length measurements), predicted by each model. In comparing the geological constraints with the model predictions, we need to take into account the estimated uncertainties of the between experiments component σ in Eq. (9) (0.06 and 0.15 μm for μ¯ f l and μ¯ pl in Eq. (19), respectively) in addition to the reported standard errors in the measurements. The prediction by the fanning-linear model in this study, which is practically identical to that by the hybrid-linear model, shows good agreement with the Vienna, Kola and KTB core samples within the errors considered (Fig. 5). Exceptionally, the discrepancy in mean length between measurement and prediction for the MITI sample is ∼ 0.9 μm, in excess of their 2σ uncertainty levels. For the parallel-curvilinear model, the estimated mean length of the MITI core sample is concordant with the prediction. The mean length predicted for other cores is, however, systematically shorter than the measured one, and the difference between measurement and prediction increases together with the uncertainties of prediction as the mean length decreases, and exceeds the 2σ uncertainty levels for Kola and KTB cores. It follows from these comparison that the difference in predicted length over geological time scale by the hybrid-linear model is comparable to the statistical uncertainty but the parallel-curvilinear model gives a shorter length estimate for the heating conditions for borehole samples and such difference may be critical for the geological samples with highly annealed FTs. It is concluded that the hybrid-linear model provides the better prediction for solving the geological problems. The length estimates by the parallel-curvilinear model show significant disagreement with geological constraints from partially annealed borehole samples despite the good statistical performance with less structured residuals (Fig. 4C). The discrepancy observed for prolonged annealing under geological conditions is presumably because FT annealing may not be extrapolated by the curvature for heating times far from laboratory conditions. The fact that linear models are superior to curvilinear models indicates that the FT annealing in zircon is fundamentally controlled by first-

order diffusion processes although somewhat affected by additional factors, such as accumulated radiation damage, and etching gaps along the annealed tracks. An appropriate use of the hybrid-linear model is recommended as follows. For most thermal history analyses, a fanning-linear model, which comprises the lower temperature part of the hybrid model, would be suitable as long as the target geological problem does not contain phenomena shorter-lived than of the order of hours. The hybrid model would be applicable to treat a wide range of heating or cooling times when applied to shorter time and higher temperature phenomena, such as frictional heating at fault movements. 4.2. Comparison with previous models The models in this study are compared with the laboratory-based annealing models using track length data, published by Tagami et al. (1998) and Rahn et al. (2004), with respect to the estimates of mean FT lengths and zircon partial annealing zone (ZPAZ). Tagami et al. (1998) provided the fanning- and parallel-linear models calculated using the data of spontaneous tracks in the NST zircon after applying the orientation criteria; the details of the difference in the data sets between this study and Tagami et al. (1998) are documented in the “Data” chapter. Their model fitting methods were fundamentally the same as those in this study proposed by Laslett and Galbraith (1996). Rahn et al. (2004) provided the fanning- and parallel-linear models using the data of “zero damage” zircons in some literatures and the same NST zircon data as Tagami et al. (1998), without adopting the orientation criteria. They used a more simplistic approach for model fitting proposed by Laslett et al. (1987), relating to a fanning model with an infinite critical temperature of the crossover point (Tc = ∞) at which fanning contours intersect on the Arrhenius plots. Fanning models (Tc = ∞) reduce the number of parameters to estimate, and predict narrower temperature ranges for zircon partial annealing zones, when compared with the fanning models with the finite critical temperature (Tc ≠ ∞). The estimates of mean FT lengths by various models are well concordant within the uncertainties (Table 4). The estimated temperatures of ZPAZ for 106 and 107 years isothermal heating by different models are recalculated (Table 5). Here we define the top and bottom of ZPAZ by the mean length reduction ratios r = 0.4 and 0.8, respectively, as discussed by Tagami et al. (1998). Errors associated with temperatures are calculated on the same base for all cases because the equational forms and the fitting procedures vary among

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

different models in the references. They are given by the uncertainty propagated from the errors involved in the estimated parameters in individual models where the error range of 1% is equally assumed for each parameter. The estimates by the parallel-curvilinear model in this study show the largest errors, mainly due to the doublelog transformation of the r parameter. In the case of fanning-linear models, the differences in ZPAZ estimates among the listed models are not significant if the errors in the calculated temperatures are taken into account. Compared with Tagami et al. (1998), the estimates of the ZPAZ bottom (lower temperature boundary) by the present μ¯ fl model show good concordance, and those of the top (higher temperature bound) are lowered by ∼ 20 °C, where the data for μ¯ fl are almost compatible with those in other studies using NST data because many of the higher temperature and accordingly the shorter-term experiments were excluded. This narrower estimate of ZPAZ by μ ¯ fl may be caused by the exclusion of the 4.5minute data that were affected by the erroneous estimate of the heating time. The fanning model of Rahn et al. (2004) that was fitted to the NST data gives the smallest differences between the top and bottom temperatures of the ZPAZ (Δ in Table 5), reflecting the nature of the fanning models with Tc = ∞ that leads to a more narrow ZPAZ. Although the model μ¯ fl gives good concordance in ZPAZ estimates with the NST-fanning model of Rahn et al. (2004), there is a distinct difference in the estimated ZPAZ between their zero damage models and μ¯ fl; both the top and bottom temperatures of the ZPAZ boundaries are systematically lower by ∼ 20–30 °C. Such difference in the estimated

89

temperature of ZPAZ could correspond to an order of magnitude of longer heating time. There are statistically insignificant differences in ZPAZ estimates among the three parallel-linear models that were fitted to the NST data, whilst the zero damage parallel model of Rahn et al. (2004) gives somewhat higher estimates. This may be because the higher thermal stability in zero damage zircon was enhanced by the parallel-linear model fitting. 4.3. Transition of annealing behavior In order to express the annealing behavior of FTs in zircon involving a transitional process with improved goodness of fit, we constructed two types of models; the hybrid-linear model (Eq. (19)) that consists of the fanning-linear model (μ¯ f l) and the parallel-linear model (μ¯ pl) with a transitional temperature zone, and the parallel-curvilinear model (Eq. (24)) that approximates smoothly the transition in the annealing process. The improvement in the residual structure obtained from these models is illustrated in Fig. 4. The fanning-linear model (Fig. 4A) has a clear structure in the residuals, that becomes fainter with the hybrid-linear model (Fig. 4B) and is nearly absent in the parallel-curvilinear model (Fig. 4C). Throughout the assessment of the dividing conditions of the analytical domain for constructing the hybrid-linear model, it was suggested that temperature governs the transition of annealing process, and that the transition takes place in a temperature range of ca. 500–700 °C. One possible factor for the transition in annealing process is the influence of the damage as produced by

Table 5 Comparison of the zircon partial annealing zones estimated for 106 and 107 years isothermal heating Model form Fanning-linear (Tc ≠ ∞) Parallel-linear Fanning-linear (Tc = ∞) Parallel-linear Fanning-linear (Tc = ∞) Parallel-linear Fanning-linear (Tc ≠ ∞) Parallel-linear Parallel-curvilinear

Data set

NST, L60 (0.1–10,000 h) Zero damage zircons, Lall NST, Lall (0.1–10,000 h) NST, L60 (b725 °C) (N475 °C) (0.001–10,000 h)

ZPAZ (106 years)

ZPAZ (107 years)

Reference

Top

Bottom

Δ

Top

Bottom

Δ

370 (7) 346 (8) 375 (8) 401 (9) 360 (8) 346 (8) 352 (12) 352 (8) 334 (58)

281 (7) 294 (7) 311 (8) 355 (9) 292 (9) 293 (9) 281 (11) 299 (8) 247 (49)

89 52 64 47 68 53 71 53 87

348 (7) 322 (8) 352 (8) 379 (9) 337 (8) 322 (8) 330 (12) 330 (8) 303 (54)

261 (7) 274 (7) 290 (8) 335 (9) 272 (8) 273 (8) 262 (11) 280 (7) 220 (46)

87 48 62 44 65 49 68 50 83

(1) (2)

(3)

Tc: critical temperature of the crossover points of annealing contours in the fanning models. In the “Data set” column, L60 and Lall indicate that tracks at an angle of N60° to the c-axis were selectively used for the model fitting, and that tracks at all orientations were used, respectively. The top and bottom of ZPAZ are recalculated by defining the mean length reduction ratio r = 0.4 and 0.8, respectively. Δ: temperature difference between the top and bottom of ZPAZ. Temperatures are given in Celsius (°C). Errors are given in parentheses; see text for the explanation of error estimation. References for the cited models are: (1) Tagami et al. (1998); (2) Rahn et al. (2004); (3) This study.

90

R. Yamada et al. / Chemical Geology 236 (2007) 75–91

alpha-decay events (alpha-damage), proposed previously (e.g., Kasuya and Naeser, 1988). This is inferred by the correspondence in the temperature range between this transition and the annealing of the alpha-damage, indicated by the change in etching behavior in zircon (e.g., Tagami et al., 1990). Tagami and Matsuura (personal communication), however, found no significant influence of the FT ages on the annealing behavior, based on the recent experimental results of isothermal annealing for 1 h using several zircons separated from rock samples with different zircon FT ages. They measured spontaneous tracks that have experienced no annealing since their formation, confirmed by the comparison with induced track length. More experiments are required for investigating the factors that control the transition in the FT annealing process. 5. Conclusion We remodelled annealing kinetics of zircon FT statistically to fit three previously published data sets of laboratory experiments using spontaneous track lengths in the NST zircon. For linear model fitting with reciprocal conversion of heating temperature, a hybrid model is obtained that consists of a fanninglinear model in the temperature domain below 725 °C and a parallel-linear model in the temperature domain above 475 °C with a transitional temperature zone in between. The constituting fanning-linear and parallellinear models have the forms

ð

ð

lf l ðt; T Þ ¼11:63 1−exp −6:95 þ0:176

logt þ 53:84 1000=T þ 0:363

1=3

and lpl ðt; T Þ ¼ 11:63ð1−expð4:35 þ 0:113ðlogt−38 760=T ÞÞÞ1=3 ; respectively. The hybrid model is expressed by the weighted sum of these two models as ll ðt; T Þ ¼ wðT Þ l¯fl ðt; T Þ þ ð1−wðT ÞÞ l¯pl ðt; T Þ; where wðT Þ ¼ 1=ð1 þ expð−11:08ð1000=T −1:169ÞÞÞ:

For the fitting of the curvilinear model (logarithmic conversion of the heating temperature), a parallel model is obtained that has the form lpc ðt; T Þ ¼ 11:17expð−expð−63:37 þ 0:212ðlogt−43:00dlog1=T ÞÞÞ: Estimates on the ZPAZ corresponding to a step-size isothermal heating event for 106 and 107 years are 281– 352 °C and 262–330 °C for the fanning-linear model, and 247–334 °C and 220–303 °C for the parallelcurvilinear model where the top and bottom of ZPAZ are defined by the mean length reduction ratios r = 0.4 and 0.8, respectively. ZPAZ estimates for the hybrid-linear model in these time ranges are equivalent to the fanninglinear model. The fanning-linear model is found to well describe annealing under geological heating conditions with time ranges longer than several hours. This is confirmed by comparison with FT length measurements of borehole core samples. The hybrid-linear model should be used if a thermal history problem contains short-time phenomena such as frictional heating at faults. Acknowledgements Professors R.F. Galbraith and A.J. Hurford are thanked for their helpful advice and comments on modelling procedures and for improving the original manuscript. Drs. A. Carter, R. Clayton and all members of the London Thermochronology Research Group are thanked for useful discussion. Dr. N. Hasebe is thanked for the helpful advice on the interpretation of a previous work. We are grateful to Dr. Paul. F. Green, Dr. Meinert Rahn, and an anonymous reviewer for giving us critical comments on the manuscript. This study was partly supported by a Grant-in Aid (no. 12440137) as well as by a Grant-in-Aid for the 21st Century COE Program (Kyoto University, G3) from the Japanese Ministry of Education, Culture, Sports, Science and Technology. [PD] References Box, G.E.P., Cox, D.R., 1964. An analysis of transformations. J. R. Stat. Soc., Ser. B 26, 211–252. Coyle, D.A., Wagner, G.A., 1996. Fission-track dating of zircon and titanite from the 9101 m deep KTB: observed fundamentals of track stability and thermal history reconstruction. Abstract of Presentation at International Workshop on Fission Track Dating, Gent.

R. Yamada et al. / Chemical Geology 236 (2007) 75–91 Crowley, K.D., Cameron, M., Schaefer, R.L., 1991. Experimental studies of annealing of etched tracks in fluorapatite. Geochim. Cosmochim. Acta 55, 1449–1465. Galbraith, R.F., Laslett, G.M., 1997. Statistical modelling of thermal annealing of fission tracks in zircon. Chem. Geol., Isot. Geosci. Sect. 140, 123–135. Gleadow, A.J.W., 1981. Fission-track dating methods: what are the real alternatives? Nucl. Tracks 5, 3–14. Gleadow, A.J.W., Hurford, A.J., Quaife, R.D., 1976. Fission track dating of zircon: improved etching techniques. Earth Planet. Sci. Lett. 33, 273–276. Green, P.F., Duddy, I.R., Gleadow, A.J.W., Laslett, G.M., Tingate, P.R., 1986. Thermal annealing of fission tracks in apatite 1. A qualitative description. Chem. Geol., Isot. Geosci. Sect. 59, 237–253. Green, P.F., Duddy, I.R., Laslett, G.M., 1988. Can fission track annealing in apatite be described by first-order kinetics? Earth Planet. Sci. Lett. 87, 216–228. Green, P.F., Hegarty, K.A., Duddy, I.R., Foland, S.S., Gorbachev, V., 1996. Geological constraints on fission track annealing in zircon. Abstract of Presentation at International Workshop on Fission Track Dating, Gent. Hasebe, N., Tagami, T., Nishimura, S., 1994. Towards zircon fission track thermochronology: reference framework for confined track length measurements. Chem. Geol., Isot. Geosci. Sect. 112, 169–178. Hasebe, N., Mori, S., Tagami, T., Matsui, R., 2003. Geological partial annealing zone of zircon fission-track system: additional constrains from the deep drilling MITI-Nishikubiki and MITIMishima. Chem. Geol. 199, 45–52. Kasuya, M., Naeser, C.W., 1988. The effect of a-damage on fissiontrack annealing in zircon. Nucl. Tracks Radiat. Meas. 14, 477–480. Lal, D., Rajan, R.S., Tamhane, A.S., 1969. Chemical composition of nuclei of ZN22 in cosmic rays using meteoritic minerals as detectors. Nature 221, 33–37. Laslett, G.M., Galbraith, R.F., 1996. Statistical modelling of thermal annealing of fission tracks in apatite. Geochim. Cosmochim. Acta 60, 5117–5131. Laslett, G.M., Kendall, W.S., Gleadow, A.J.W., Duddy, I.R., 1982. Bias in measurement of fission-track length distributions. Nucl. Tracks 6, 79–85. Laslett, G.M., Green, P.F., Duddy, I.R., Gleadow, A.J.W., 1987. Thermal annealing of fission tracks in apatite, 2. A quantitative analysis. Chem. Geol. 65, 1–13.

91

Murakami, M., Yamada, R., Tagami, T., 2006. Short-term annealing characteristics of spontaneous fission tracks in zircon. Chem. Geol. 227, 214–222. Otsuki, K., Monzawa, N., Nagase, T., 2003. Fluidization and melting of fault gouge during seismic slip: identification on the Nojima fault zone and implications for focal earthquake mechanisms. J. Geophys. Res. 108 (B4), 2192, doi:10.1029/2001JB001711. Paul, T., Fitzgerald, P.F., 1992. Transmission electron microscope investigation of fission tracks in apatite. Am. Mineral. 77, 336–344. Rahn, M.K., Brandon, M.T., Batt, G.E., Garver, J.I., 2004. A “zero damage” model for fission track annealing in zircon. Am. Mineral. 89, 473–484. Schwarz, G., 1978. Estimating the dimension of a model. Ann. Stat. 6, 461–464. Tagami, T., Ito, H., Nishimura, S., 1990. Thermal annealing characteristics of spontaneous fission tracks in zircon. Chem. Geol., Isot. Geosci. Sect. 80, 159–169. Tagami, T., Uto, K., Matsuda, T., Hasebe, N., Matsumoto, A., 1995. K–Ar biotite and fission-track zircon ages of the Nisatai Dacite, Iwate Prefecture, Japan: a candidate for Tertiary age standard. Geochem. J. 29, 207–211. Tagami, T., Carter, A., Hurford, A.J., 1996. Natural long-term annealing of the zircon fission-track system in Vienna Basin deep borehole samples: constraints upon the partial annealing zone and closure temperature. Chem. Geol. 130, 147–157. Tagami, T., Galbraith, R.F., Yamada, R., Laslett, G.M., 1998. Revised annealing kinetics of fission tracks in zircon and geological implications. In: Van den haute, P., De Corte, F. (Eds.), Advances in Fission-Track Geochronology. Kluwer Academic Publishers, Dordrecht, The Netherlands, pp. 99–112. Yamada, R., Tagami, T., Nishimura, S., 1995a. Confined fission-track length measurement of zircon: assessment of factors affecting the paleotemperature estimate. Chem. Geol., Isot. Geosci. Sect. 119, 293–306. Yamada, R., Tagami, T., Nishimura, S., Ito, H., 1995b. Annealing kinetics of fission tracks in zircon: an experimental study. Chem. Geol. 122, 249–258. Yamada, K., Tagami, T., Shimobayashi, N., 2003. Experimental study on hydrothermal annealing of fission tracks in zircon. Chem. Geol. 201, 351–357.