Exploring the law of detrital zircon: LA-ICP-MS and CA-TIMS geochronology of Jurassic forearc strata, Cook Inlet, Alaska, USA

Detrital zircon (DZ) U-Pb geochronology is a staple of modern stratigraphic research that proliferated with increasingly widespread use of laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) (e.g., Gehrels, 2014). Rapid data acquisition renders LA-ICP-MS well suited for the DZ analyses that are extensively used in provenance work and maximum depositional age (MDA) assessments (e.g., Schaltegger et al., 2015). This study examined MDAs, which are based on a logical premise that Gehrels (2014) referred to as the law of DZ: A sedimentary rock cannot be older than the youngest zircon crystal it contains (Houston and Murphy, 1965).

The validity of a DZ MDA is always complicated by uncertainty, including analytical, systematic, and geologic sources. Laboratory-reported confidence intervals, however, principally reflect analytical precision and reproducibility of standard materials, and repeat measurements do not mitigate sample-specific systematic uncertainty (Schoene, 2014). Inter-element fractionation during laser ablation requires frequent within-session analyses of reference zircon and is a significant source of systematic uncertainty in LA-ICP-MS geochronology (e.g., Schaltegger et al., 2015). Well-characterized zircon yield LA-ICP-MS dates that typically coincide with associated chemical abrasion–thermal ionization mass spectrometry (CA-TIMS) dates, but systematic offsets, likely reflecting matrix effects, are observed (Schoene, 2014). In fact, LA-ICP-MS dates of relatively young (i.e., Mesozoic–Cenozoic) zircon are prone to incorporating fractionation-associated matrix effects, imparting too-young biases of as much as ∼5% (Allen and Campbell, 2012). Mesozoic–Cenozoic strata are common in basin analysis, and MDAs that are younger than existing stratal age constraints may have considerable implications (e.g., Surpless et al., 2006).

Lead loss is largely unconstrained by single LA-ICP-MS analyses of Mesozoic–Cenozoic zircon (Spencer et al., 2016) but is often cited to account for DZ dates that are ostensibly too young. Additionally, the impact of material properties on ablation behavior—the above-noted matrix effects—and the statistical nature of overlapping dates within youthful (i.e., near stratal age) DZ populations are rarely discussed (Coutts et al., 2019). Fortunately, total uncertainty can be reduced with complementary CA-TIMS geochronology, which mitigates and assesses Pb loss for Mesozoic–Cenozoic zircon, is not subject to laser-induced matrix effects, and yields dates commonly ∼50× more precise than LA-ICP-MS. Recent DZ studies have combined LA-ICP-MS and CA-TIMS to determine MDAs (e.g., Wainman et al., 2018), but experiments that explicitly explore the law of DZ and compare dates and MDAs from these two methods are lacking.

Within this context, we conducted LA-ICP-MS and CA-TIMS geochronology of DZ in Jurassic forearc strata of southern Alaska (United States). An oceanic island arc provenance for the sampled sandstones, which have large proportions of youthful zircon, renders these strata an excellent case for evaluating best practices for establishing MDAs. Enhanced LA-ICP-MS analytical protocols are aimed to improving accuracy. The youngest DZ LA-ICP-MS dates are complemented by CA-TIMS dates from the same crystals. We compare the LA-ICP-MS and CA-TIMS dates and age constraints, provide recommendations in light of these new data, and briefly discuss their potential to refine the Jurassic geologic time scale.

Cook Inlet forearc basin (Fig. 1) hosts an ∼18-km-thick Mesozoic–Cenozoic stratigraphic record (e.g., LePain et al., 2013). The Jurassic manifestation of this long-lived basin was coupled with the Talkeetna arc, which migrated northward in the paleo–Pacific Ocean above a north-dipping subduction zone (e.g., Clift et al., 2005). The Talkeetna arc’s forearc stratigraphy comprises the Middle Jurassic Tuxedni Group and Chinitna Formation and the Upper Jurassic Naknek Formation (e.g., LePain et al., 2013; Fig. 1). Accretionary tectonics ultimately extinguished Talkeetna arc magmatism in the Late Jurassic (e.g., Clift et al., 2005), although the exact geometry and timing of terrane amalgamation and final collision with North America is debated (e.g., Stevens Goddard et al., 2018).

The Jurassic forearc stratigraphy is well exposed between Iniskin and Tuxedni Bays of the lower Cook Inlet (Detterman and Hartsock, 1966; Fig. 1). In this area, the Chinitna Formation comprises an ∼700-m-thick succession of principally shallow-marine strata of the Tonnie Siltstone and Paveloff Siltstone Members (Herriott et al., 2016; Herriott and Wartes, 2017; Fig. 1). The overlying Naknek Formation in the area is an ∼1500-m-thick interval of shallow- and deep-marine strata of the Chisik Conglomerate, lower sandstone (informal), Snug Harbor Siltstone, and Pomeroy Arkose members (Herriott et al., 2017; Fig. 1).

Detterman and Westermann (1992) summarized the largely ammonite-based Jurassic biostratigraphy of southern Alaska, reporting that the Chinitna Formation is Callovian; Chisik, lower sandstone, and Snug Harbor are Oxfordian; and Pomeroy is Kimmeridgian (Fig. 1). However, faunal assemblages noted by Detterman and Westermann (1992) also suggest that lowermost Tonnie may be latest Bathonian, and Pomeroy is potentially as old as Oxfordian. The Chisik is not fossiliferous, but stratigraphic and biostratigraphic relations indicate that this unit is probably Oxfordian, and may be associated with Callovian–Oxfordian transition climate change (Herriott et al., 2017). Recent DZ LA-ICP-MS studies in southern Alaska yielded Chinitna and Naknek constraints that are notably younger than biostratigraphic correlations suggest (Finzel and Ridgway, 2017; Reid et al., 2018; Stevens Goddard et al., 2018; Herriott et al., 2019).

We sampled sandstone beds from the base of each member in the Chinitna and Naknek Formations, bracketing ∼1400 m of stratigraphy that extends across the Middle–Late Jurassic boundary. All samples were prepared and analyzed at Boise State University’s Isotope Geology Laboratory (Boise, Idaho, USA) (see the GSA Data Repository¹).

We conducted three rounds of U-Pb geochronology: (1) LA-ICP-MS of ∼115 zircon grains per sample, with one date per grain; (2) two additional LA-ICP-MS dates per zircon for the ∼15 youngest grains per sample identified in round 1; and (3) CA-TIMS of the ∼5 youngest grains (labeled z1–z5) per sample based on sorting LA-ICP-MS multiple-analysis results by weighted mean (WM) date where n = 3 and probability of fit (PoF) is >0.05. Round 3 zircons were plucked from their mounts and broken into fragments; selected fragments were chemically abraded and analyzed by CA-TIMS (Mattinson, 2005).

We assessed round 1 youthful DZ via three approaches: (1) youngest single grain (YSG; Dickinson and Gehrels, 2009); (2) youngest mode of the kernel density estimation (KDE) (YMKDE; cf. YPP of Dickinson and Gehrels, 2009); and (3) youngest statistical population with a mean square weighted deviation of ∼1.00 (YSP; Coutts et al., 2019). Round 2 data enabled a fourth LA-ICP-MS-based determination: youngest single grain with multiple analyses (YSGMA; WM, n = 3 [rounds 1 + 2], PoF >0.05; e.g., Spencer et al., 2014). WMs of the youngest CA-TIMS dates that overlap at 2σ uncertainty and have PoF >0.05 provide round 3 maximum constraints.

Round 1 yielded nearly entirely Jurassic dates, with KDEs revealing unimodal date distributions (Fig. 2). Dates from round 1 single analyses and rounds 1 + 2 WMs for zircon crystals z1–z5 per sample are younger than their associated (i.e., same zircon) round 3 CA-TIMS dates, although ∼60% of the date pairs overlap at 2σ (Fig. 2). All CA-TIMS dates are concordant. Follow-up CA-TIMS of an archived grain fragment from each of three critically young DZ reproduced results within 2σ for each pair (Fig. 2). Average per-sample U values are low (77–112 ppm), and U versus date plots indicate no correlation or subtly increasing U toward older dates. Maximum depositional dates (MDDs) and MDAs are discussed below, adapting the date-versus-age distinction noted by Schoene (2014).

CA-TIMS dates establish the MDAs of this study. All single-grain MDDs (YSG and YSGMA) are younger than corresponding CA-TIMS MDAs and rarely overlap at 2σ (Fig. 3). YMKDE MDDs are ∼1–4.5 m.y. older than the CA-TIMS MDAs, whereas YSP MDDs commonly overlap at 2σ with the CA-TIMS MDAs (Fig. 3). Mitigating laser-induced matrix effects and evaluating the distribution of dates within youthful DZ populations are keys to determining accurate LA-ICP-MS–based MDAs.

Matrix effects–related uncertainty principally reflects varying degrees of radiation damage among unknowns and references, rendering variable ablation rates and concomitantly variable inter-element fractionation (e.g., Sliwinski et al., 2017). These effects are potentially more problematic for young zircon with low U (and Th) values that are dated relative to older references and/or references with higher U (and Th) content (Allen and Campbell, 2012). Thermal histories and other matrix factors further complicate these relations (Marillo-Sialer et al., 2016), but structural homogenization by thermal annealing reduces ablation rate variability and mitigates this systematic uncertainty (e.g., Allen and Campbell, 2012; Sliwinski et al., 2017; see also Mattinson, 2005).

We addressed uncertainty due to matrix effects by thermally annealing zircon prior to LA-ICP-MS and carefully selecting secondary references. Comparisons of our LA-ICP-MS dates with published DZ LA-ICP-MS dates from Jurassic strata of southern Alaska are hampered by stratigraphic, biostratigraphic, and analytical complexities. However, we collected samples 09BG010-14.5A (this study) and 09BG010-14.5C (Herriott et al., 2019) from the same bed, and KDEs of these LA-ICP-MS results indicate a modest systematic offset (−2.2% at mode for 09BG010-14.5C; Fig. 2) that may reflect matrix effects.

Establishing accurate MDAs also requires consideration of how densely distributed DZ dates impact interpretations. YSG and YSGMA MDDs of this study, if regarded as MDAs, would impart a too-young bias on chronostratigraphic interpretations. Two factors, which are not mutually exclusive, likely account for this bias: (1) selectively sampling the low-probability tail of a normal distribution of data resulting from random statistical fluctuations during an analytical session, and (2) Pb loss. Round 2 experiments, with multiple laser ablation spots placed on the same youngest grains, aimed to distinguish these two factors. YSGMA determinations are older than the YSG results (Fig. 3), indicating a low-probability tail bias for the YSGs. However, a consistent residual young bias—attributable to Pb loss—remains in the YSGMA MDDs relative to the CA-TIMS MDAs (Fig. 2).

Older components of closely clustered DZ dates present additional challenges. The YMKDE MDDs are older than the CA-TIMS MDAs because the full probability distribution incorporates a range of truly older dates (PoF = 0.00 for WMs of all dates per sample). The YSP approach selects the youngest subset of dates with scatter that can be explained by the uncertainties, extracting a normally distributed subsample from the youngest tail of the distribution (e.g., Fig. 2, Tonnie Siltstone Member). YSP MDDs are our preferred LA-ICP-MS constraints due to their explicit tie to the aforementioned statistical fluctuations during analysis, high n, and best overall coincidence with the CA-TIMS MDAs.

Spencer et al. (2016) and Coutts et al. (2019) also noted that the (normal) distribution of dates within youthful DZ populations can undermine the accuracy of LA-ICP-MS MDAs, with single-grain assessments prone to underestimating stratal age. Nevertheless, Spencer et al. (2016) favored single-grain, multiple-analysis MDAs (e.g., YSGMA) for in situ techniques, understandably citing the inability to determine with certainty that multi-grain detrital population samples (e.g., YSP) record truly cogenetic crystallization of zircon. Geologic context should always be considered, but the YSP MDDs of our study do not presume a narrowly defined genesis for the selected zircon; rather, we simply assert that these subsamples reflect coeval zircon crystallization as resolved by LA-ICP-MS. Our results—benchmarked by CA-TIMS—suggest that conducting statistical assessments of peaks, clusters, and/or tails of DZ LA-ICP-MS date distributions will consistently render more reliable results than potentially problematic single-grain determinations. CA-TIMS of the youngest DZ in a sample circumvents the need to favor either single- or multi-grain MDAs for LA-ICP-MS.

We recommend that DZ MDA studies of Mesozoic–Cenozoic strata include thermally annealing zircon prior to analysis and employ the YSP method, although alternative WM approaches (e.g., YC2σ of Dickinson and Gehrels, 2009) may be suitable for lower n youthful population samples. Focusing on the single youngest LA-ICP-MS date in a densely sampled DZ population is not well suited to characterizing the age of that population, although multiple LA-ICP-MS analyses on the same youthful grain(s) would improve results. CA-TIMS of the youngest grains should be conducted for critical applications, including chronostratigraphy. The most robust MDAs are derived from equivalent CA-TIMS dates from multiple grains; ideally, multiple fragments per grain would be dated to test intra-grain reproducibility and minimize geochronologic uncertainty.

Our recommendations aim to diminish or eliminate too-young biases in DZ MDAs. However, there will be cases where MDAs truly are younger than previous constraints suggest. In this study, the Tonnie MDA indicates that uppermost Bathonian(?)–Callovian strata are not older than 159.57 ± 0.11 Ma, which is ∼4 m.y. younger than even the Callovian–Oxfordian boundary (Fig. 1). Furthermore, the base of the Naknek Formation yielded a 157.33 ± 0.11 Ma MDA, coinciding with the Oxfordian–Kimmeridgian boundary (Fig. 1), yet the entire Oxfordian stratigraphy overlies the sampled stratum. These relations may reveal discrepancies with faunal correlations or time-scale calibration. Currently, the Middle–Late Jurassic time scale has few radioisotopic constraints (Gradstein et al., 2012) and would benefit from additional CA-TIMS dates. Further time-scale refinements could consider CA-TIMS MDAs for fossiliferous strata along Jurassic convergent margins.

Discussions with M. Rioux, J.V. Jones III, P.B. O’Sullivan, R.O. Lease, A.L. Willingham, D.L. LePain, and M.R. Guthrie improved this study. We thank G.E. Gehrels, C.J. Spencer, and N.M. McLean for thoughtful reviews; thanks also to science editor C. Clark for handling the manuscript. D. Pierce and K. Grosswiler prepared samples; M. Gavel digitized laser ablation spot locations for Figure DR1. The State of Alaska and U.S. Geological Survey’s National Cooperative Geologic Mapping Program (award G17AC00211) funded this work. Funding for the analytical infrastructure of the Boise State University Isotope Geology Laboratory (Boise, Idaho, USA) was provided by the National Science Foundation (grants EAR-0521221, EAR-0824974, EAR-1337887, EAR-1735889).