Late Jurassic paleogeography of the U.S. Cordillera from detrital zircon age and hafnium analysis of the Galice Formation, Klamath Mountains, Oregon and California, USA Open Access

The Late Jurassic was a time of tectonic transition in the North American Cordilleran orogen, as the North American plate shifted northward and increased in velocity (e.g., Seton et al., 2012). However, agreement remains elusive regarding the Late Jurassic tectonic evolution of the Oregon–California, USA, segment of the continental margin during the Nevadan orogeny. Here, we use the term Nevadan orogeny to refer to the relatively brief (ca. 157–145 Ma) deformational event in the Klamath Mountains province and Western Sierra Nevada metamorphic province (e.g., Harper and Wright, 1984; Schweickert et al., 1984; Harper et al., 1994; Schweickert, 2015), rather than the longer period of deformation (>20 m.y.) hypothesized by others (e.g., Tobisch et al., 1987, 1989; Saleeby et al., 1989; Hacker and Ernst, 1993).

Tectonic models for the Nevadan orogeny range from double-sided subduction of an intervening oceanic plate that resulted in collision of an east-facing island arc with the west-facing continental margin arc (e.g., Schweickert and Cowan, 1975; Ingersoll and Schweickert, 1986; Ingersoll, 2008; Schweickert, 2015) to transtensional forearc and intra-arc extension, followed by contraction along the continental margin above east-directed subduction due to changes in relative plate-motion vectors (e.g., Davis et al., 1978; Saleeby, 1981; Burchfiel et al., 1992; Harper et al., 1994) and possible Middle Jurassic accretion and northward migration of the Insular superterrane (Saleeby and Busby-Spera, 1992; Saleeby and Dunne, 2015). Furthermore, a global-scale tectonic model based on geophysical interpretations postulates west-directed subduction of the North American plate beneath a stationary east-facing island archipelago, which caused North America to migrate westward until Late Jurassic–Cretaceous collision with the archipelago (Sigloch and Mihalynuk, 2013, 2017, 2020; Clennett et al., 2020).

The Upper Jurassic Galice Formation in the Western Klamath Mountains province was deposited prior to and during the Nevadan orogeny (Garcia, 1979, 1982; Wyld and Wright, 1988; Harper et al., 1994). These strata comprise the youngest rocks deformed during Late Jurassic Nevadan deformation, so details of their depositional age and sediment provenance can inform proposed models of the Late Jurassic tectonic evolution of this region. Here, we present provenance interpretations to identify possible source-to-sink transport pathways based on new U-Pb detrital zircon ages (n = 2792) and ε_Hf values (n = 603) from the Galice Formation. We use these results to evaluate tectonic models of the Nevadan orogeny for development of a paleogeographic reconstruction of the U.S. segment of the Cordillera during Late Jurassic time. Combined with new ε_Hf data from eight Late Jurassic plutons of the Klamath Mountains, detrital zircon laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) U-Pb age and ε_Hf data from 11 metasandstone samples collected along 200 km of strike length permit detailed characterization of the depositional age and provenance of the upper turbiditic strata within the Galice Formation.

The Klamath Mountain province is a system of fault-bounded, imbricated lithotectonic units of broadly oceanic affinity that regionally dip eastward (Fig. 1). The lithotectonic units generally decrease in age to the west and structurally downward, with ages ranging from Neoproterozoic remnants in the Eastern Klamath Mountains (Wallin et al., 1995; Mankinen et al., 2002; Lindsley-Griffin et al., 2003, 2006; Grove et al., 2008) to Upper Jurassic in the Western Klamath Mountains (Diller, 1903; Irwin, 1960, 1994). Although recent interpretation of mantle tomography data suggests that amalgamation of the Klamath terranes may have occurred significantly west of the continental margin prior to latest Jurassic to Cretaceous accretion (e.g., Sigloch and Mihalynuk, 2013, 2017; Clennett et al., 2020), whole-rock Sr and Nd and U-Pb detrital zircon results from the Western Klamath terrane demonstrate connection to the North American continent during Jurassic time (e.g., Frost et al., 2006; LaMaskin et al., 2022).

Irwin (1960) divided the Klamath Mountains province into four lithotectonic units (his “belts”), or tectonostratigraphic terranes in modern terminology. From east to west, the composite terranes are: Eastern Klamath, Central Metamorphic, Western Paleozoic and Triassic, and Western Jurassic (Fig. 1). Below, we summarize the major lithologies, depositional setting, and published geochronology of the Klamath Mountains province.

The Eastern Klamath composite terrane consists of the Trinity, Yreka, and Redding subterranes. The Trinity subterrane consists of Neoproterozoic to lower Paleozoic ultramafic rocks intruded by mainly Silurian–Devonian gabbroic plutons (e.g., Quick, 1981; Lindsley-Griffin, 1991; Lindsley-Griffin et al., 2008; Wallin et al., 1995; Wallin and Metcalf, 1998). The Yreka subterrane structurally overlies the Northeastern Trinity terrane and consists of a sequence of tectonically juxtaposed and disrupted units that range in age from Neoproterozoic to Devonian (Grove et al., 2008, and references therein). In contrast, the Redding subterrane is a broadly homoclinal-dipping sequence of volcanic and sedimentary strata that overlie the Eastern Trinity subterrane and range in age from Devonian to Jurassic (e.g., Watkins, 1985, 1993; Renne and Scott, 1988; Noble and Renne, 1990; Wallin and Metcalf, 1998). Detrital zircon U-Pb, biogeographic data, and paleomagnetic data suggest a Baltic rather than Laurentian origin for Neoproterozoic through Paleozoic rocks of the Eastern Klamath terranes (Wright and Wyld, 2006; Grove et al., 2008).

The Central Metamorphic terrane is in fault contact beneath the Eastern Klamath terrane along the Trinity fault and consists of the structurally lower Salmon Formation and structurally higher Grouse Ridge Formation (Davis, 1968; Barrow and Metcalf, 2006). The Salmon Formation consists of amphibolitic metabasite, and the Grouse Ridge Formation is primarily metasedimentary, comprising calc-silicates and metacarbonates, with some metabasite (Holdaway, 1965; Peacock and Norris, 1989). Lanphere et al. (1968) reported Devonian Rb/Sr whole-rock ages for the Grouse Ridge Formation and a wide range of K-Ar ages for the Salmon and Grouse Ridge formations. However, Barrow and Metcalf (2006) determined Early Permian ⁴⁰Ar/³⁹Ar cooling ages (274 ± 2 Ma) on hornblende from the Salmon Formation.

The Western Paleozoic and Triassic (now determined as Triassic–Jurassic; see below) composite terrane occupies the central part of the Klamath Mountains province (Fig. 1) and structurally underlies the Central Metamorphic terrane along the Siskiyou fault. It consists of a stack of east-dipping, fault-bounded, Triassic and Jurassic terranes. Structurally downward (westward), these are the Stuart Fork Formation (= Fort Jones terrane; Snoke and Barnes, 2006), North Fork, Eastern Hayfork, Western Hayfork, and Rattlesnake Creek terranes. The Stuart Fork Formation is characterized by blueschist-facies metasedimentary and metabasic rocks with local eclogite blocks (Hotz et al., 1977; Goodge, 1990). K-Ar dating of white mica yielded Middle Triassic ages (Hotz et al., 1977). The North Fork terrane is an assemblage of metabasites and overlying clastic metasedimentary rocks that typically display greenschist-facies assemblages (Ando et al., 1983; Ernst, 1990, 1998; Scherer and Ernst, 2008). U-Pb (zircon), ⁴⁰Ar/³⁹Ar, and fossil assemblages indicate depositional ages from Permian to Middle Jurassic (Ando et al., 1983; Irwin and Blome, 2004; Hacker et al., 1993). Detrital zircon data indicate deposition of metasandstone units in Early to Middle Jurassic time, depending on whether or not two anomalously young Middle Jurassic grains are included in maximum depositional age considerations (Scherer and Ernst, 2008).

The Eastern and Western Hayfork terranes were identified by Wright (1982) as two distinct parts of Irwin’s (1972) Hayfork terrane. The Eastern Hayfork terrane is a variably chert–argillite- and feldspathic wacke-matrix mélange and broken formation (Irwin, 1972; Wright, 1982; Ernst et al., 2017) with blocks of chert, limestone, metasandstone, metabasite, metagabbro, and metaserpentinite. Detrital zircons from mélange matrix display sensitive high-resolution ion microprobe–reverse geometry (SHRIMP-RG) U-Pb ages ranging from late Archean to Late Triassic (Ernst et al., 2017; Barnes et al., 2021). These younger ages are consistent with fossil ages of mainly Permian to Late Triassic (Irwin, 1972; Irwin and Galanis, 1976; Stevens et al., 1987). In contrast, detrital zircon assemblages from metasandstone blocks yield Proterozoic and Archean ages (Scherer et al., 2010).

The Western Hayfork terrane consists of arc-related metasandstone and volcaniclastic rocks with intercalated siliceous argillite (Wright, 1982; Wright and Fahan, 1988; Donato et al., 1996; Barnes and Barnes, 2020). It rests with faulted depositional contact on the underlying Rattlesnake Creek terrane (Wright and Fahan, 1988; Donato et al., 1996). ⁴⁰Ar/³⁹Ar and K/Ar dates from igneous hornblende of the Western Hayfork terrane indicate deposition from 177 Ma to 168 Ma (Wright and Fahan, 1988; Hacker et al., 1995; Donato et al., 1996); thus, this terrane is best interpreted as a sedimentary apron adjacent to the Middle Jurassic arc (Barnes and Barnes, 2020).

The Rattlesnake Creek terrane forms the depositional basement to the Western Hayfork terrane and consists of block-on-block ophiolitic mélange (Irwin, 1972; Donato, 1987, 1989; Wright and Wyld, 1994) that is locally overlain by layered volcanic, volcaniclastic, and crystal-lithic arenite strata that Wright and Wyld (1994) termed the “cover sequence” (also see Gray, 1986). Scant fossil evidence from limestone and radiolarian chert blocks indicates that deposition of mélange blocks was as young as Late Triassic to Early or possibly Middle Jurassic (Silberling and Irwin, 1962; Irwin and Galanis, 1976; Irwin et al., 1982, 1983, 1985; Irwin and Blome, 2004). Wright and Wyld (1994) reported Late Triassic to Early Jurassic ages (ca. 207–193 Ma) of plutons that crosscut the cover sequence. However, LaMaskin et al. (2022) reported detrital zircon LA-ICP-MS U-Pb ages as young as Middle to earliest Late Jurassic (ca. 170–161 Ma) from cover-sequence samples.

The Western Jurassic terrane was subdivided by Blake (1984) into subterranes; however, Harper (2006) pointed out that the terrane consists of two tectonic elements: the Rogue-Chetco arc complex and the Josephine ophiolite. The Josephine ophiolite is largely interpreted to have resulted from suprasubduction-zone rifting from 164 Ma to 162 Ma (Harper, 1984; Harper et al., 1994), whereas the Rogue-Chetco arc formed outboard of the ophiolite. The Josephine ophiolite and volcanic strata of the Rogue Formation are overlain by hemipelagic and flysch deposits of the Galice Formation. Additional details about development of the Western Jurassic terrane are presented in the following sections.

Jurassic tectonism in the Klamath Mountains province involved late Early to early Middle Jurassic contractional deformation (Siskiyou orogeny; Coleman et al., 1988; Sullivan, 2009); Middle Jurassic oblique, margin-parallel rifting; and Late Jurassic shortening and marginal basin collapse (Nevadan orogeny; after Coleman et al., 1988; Harper et al., 1994). The Siskiyou orogeny juxtaposed the Stuart Fork Formation, North Fork terrane, Eastern Hayfork terrane, and composite Rattlesnake Creek and Western Hayfork terrane assemblage along a series of east-dipping thrust faults (Fig. 1; e.g., Coleman et al., 1988; Wright and Fahan, 1988; Barnes et al., 2006; Barnes and Barnes, 2020). Timing of this contractional event is constrained to ca. 170 Ma by the age of the Ironside Mountain pluton (multi-grain thermal ionization mass spectrometry [TIMS] U-Pb age of zircon; Wright and Fahan, 1988), which intrudes the Wilson Point thrust fault that juxtaposes the Rattlesnake Creek/Western Hayfork assemblage beneath the Eastern Hayfork terrane.

Following Siskiyou deformation, intra-arc rifting ruptured the Rattlesnake Creek terrane and resulted in the formation of the Josephine ophiolite 164–162 Ma (Harper et al., 1994; Wright and Wyld, 1986). Evidence for rifting consists of rift-edge assemblages on both the eastern (Snoke, 1977; Saleeby et al., 1982) and western (Yule et al., 2006) sides of the Josephine ophiolite basin (Fig. 1). Reconstructions of extension directions and ophiolite pseudo-stratigraphy indicate north–south extension (in present coordinates) along a series of well-developed, en echelon spreading ridges linked by long transform faults akin to the Andaman Sea of the Northeastern Indian Ocean (Harper, 1982; Harper et al., 1986; Yoshinobu and Harper, 2004). At the same time, broadly calc-alkaline arc magmatism was active to the east and south of the Josephine basin (168–156 Ma Wooley Creek suites; Wright and Fahan, 1988; Allen and Barnes, 2006) and more tholeiitic magmatism to the northwest (166–157 Ma Rogue-Chetco arc; Yule, 1996; this study).

Opening of the Josephine rift basin provided accommodation space for deposition of hemipelagic sediments and turbidites of the Galice Formation and volcanogenic rocks of the Rogue Formation beginning at ca. 162 Ma. However, by 157–155 Ma (Harper et al., 1994; Dailey and Barnes, 2020), under-thrusting of the Josephine-Galice basin marked initiation of the Nevadan orogeny (Blackwelder, 1914; see discussions in Hacker et al., 1995; Snoke and Barnes, 2006; Dickinson, 2008). Contractional deformation is interpreted to have continued through ca. 150 Ma (Saleeby and Harper, 1993; Harper et al., 1994), with juxtaposition of the Rogue-Chetco arc (possibly including the distal interfingering Galice Formation; e.g., Yule, 1996) and its basement beneath the Josephine ophiolite + Galice assemblage, and of the Josephine ophiolite + Galice assemblage beneath the Rattlesnake Creek and overlying terranes along the Orleans thrust (Hershey, 1906, 1911; Harper et al., 1994).

Voluminous Middle Jurassic plutonism began in the Southern Klamath Mountains province with emplacement of the Ironside Mountain batholith from ca. 170 Ma to 168 Ma (Fig. 1; Wright and Fahan, 1988). This activity occurred immediately after regional thrusting that placed the Eastern Hayfork terrane and overlying units above the Western Hayfork/Rattlesnake Creek terranes along the east-dipping Wilson Point thrust (Barnes et al., 2006). The two largest plutons of the Ironside Mountain batholith (Ironside and Wildwood; Fig. 1) consist primarily of quartz monzodiorite to quartz monzonite characterized by mafic assemblages of biotite + 2- or 3-pyroxene or biotite–hornblende–pyroxene assemblages. These plutons display a K₂O enrichment trend unlike most younger Mesozoic plutons in the Klamath Mountains province (Lanphere et al., 1968; Charlton, 1979; Barnes et al., 2006), and the most evolved rocks display the highest Zr contents (>200 ppm; Barnes et al., 2006; Angulo, 2022). Initial ε_Nd and ⁸⁷Sr/⁸⁶Sr values are approximately uniform at 5.2 and 0.7037, respectively (Barnes et al., 2006). Two small satellite plutons vary from olivine pyroxenite to quartz diorite and are zircon-poor.

Starting at ca. 168 Ma and ending at ca. 156 Ma, arc plutonism shifted northward from the Ironside Mountain batholith to the Central and Northern Klamath Mountains province to form the Wooley Creek suite (Fig. 1; Allen and Barnes, 2006; Coint et al., 2013; Barnes et al., 2021; this study). Wooley Creek suite plutons are broadly calc-alkalic and display a wide range of compositions, from olivine pyroxenite to two-mica granite. Quartz diorite, tonalite, and granodiorite are the most common rock types. These plutons display variable Sr, Nd, oxygen, and Hf isotopic values, including the highest δ¹⁸O and initial ⁸⁷Sr/⁸⁶Sr, and the lowest ε_Nd of any Klamath Mountain plutons, and in the Wooley Creek–Slinkard system, a range of ε_Hf (zircon) values from 13.3 to −3.1 (Barnes et al., 1990, 2021; Allen and Barnes, 2006, and references therein).

West of the Josephine ophiolite, Middle to Late Jurassic arc magmatism resulted in the deposition of volcanic rocks and volcanogenic sediments of the Rogue Formation (Garcia, 1979, 1982) and coeval gabbroic through tonalitic plutons of the Chetco batholith (aka Chetco complex, Illinois River gabbro; Jorgenson, 1970; Yule, 1996; McLachlin, 2011). Multi-grain TIMS U-Pb (zircon) ages of Chetco plutonic rocks range from 166.6 ± 2.2 Ma to 157.4 ± 1.3 Ma (Yule, 1996; Yule et al., 2006). Plutons of the Chetco batholith are mainly mafic to intermediate in composition (olivine gabbro to hornblende diorite) and display arc tholeiite geochemical affinities with zirconium contents <100 ppm (Yule, 1996; McLachlin, 2011; Weiss, 2014). The batholith also includes smaller-volume tonalitic bodies.

Plutonic activity resumed at ca. 150 Ma, with emplacement of the Western Klamath plutonic suite (Allen and Barnes, 2006). Western Klamath plutonism spans 150–143 Ma and is marked by plutons that intrude the Galice Formation (e.g., Summit Valley, Pony Peak, Bear Peak, and the early stage of Grants Pass; Fig. 1; Saleeby and Harper, 1993; Harper, 2006). The Summit Valley and Pony Peak plutons intrude the Orleans thrust (Harper, 2006, and references therein). Western Klamath suite plutons range from olivine pyroxenite to biotite–hornblende tonalite, with one example of a potassic (monzodioritic) unit in the Bear Mountain complex (Snoke et al., 1981). Although this suite of plutons encompasses the same SiO₂ range as the Wooley Creek suite, it displays lower initial ⁸⁷Sr/⁸⁶Sr and higher ε_Nd than the Wooley Creek suite (Allen and Barnes, 2006).

The Galice basin is interpreted to have formed as a result of suprasubduction zone extension that produced the 164–162 Ma Josephine ophiolite and was coeval with the Rogue-Chetco arc (Harper, 1984; Harper and Wright, 1984; Harper et al., 1994; MacDonald et al., 2006). The Galice Formation was deposited on and interbedded with Rogue volcanic rocks in the northern Klamath Mountains, and the Josephine ophiolite farther south (Figs. 1 and 2). Although primitive boninitic lava compositions and the spatial distribution of the ophiolite across and trenchward of the contemporaneous Wooley Creek magmatic suite is consistent with generation in a forearc environment (Harper, 2003), a number of authors suggest a backarc or intra-arc basin setting (e.g., Snoke, 1977; Saleeby et al., 1982; Harper et al., 1994; Yule et al., 2006).

The Galice Formation includes a lower hemipelagic sequence, a middle transition zone, and an upper turbidite sequence (Fig. 2; Pinto-Auso and Harper, 1985; Harper, 1994; MacDonald et al., 2006). The hemipelagic sequence is ~50 m thick, and at its base it is intercalated with Josephine ophiolitic lavas (Harper, 1994; MacDonald et al., 2006). The hemipelagic sequence is mainly green to black slaty argillite and lesser radiolarian chert, with components of volcanic detritus (e.g., bipyramidal quartz; Pinto-Auso and Harper, 1985; MacDonald et al., 2006; Pessagno, 2006). Rare thin, graded sandstone beds within the hemipelagic sequence are compositionally similar to the upper turbidite sandstone, and one shows evidence of scour into underlying chert (MacDonald et al., 2006). Radiolarian biostratigraphy from the hemipelagic sequence indicates late Callovian to middle Oxfordian depositional age (Pessagno et al., 1993; Pessagno, 2006), but the 162 Ma age of the underlying Josephine ophiolite provides a maximum depositional age of the hemipelagic sequence, resulting in an age of ca. 162–154.8 Ma (time scale of Gradstein et al., 2020). Correlation of radiolarian tuffs within the hemipelagic sequence with 157 ± 2 Ma volcanic rocks of the Rogue Formation (Saleeby, 1984) suggests that the uppermost hemipelagic sequence is missing, which may represent a disconformity (Pessagno and Blome, 1990), or disrupted bedding due to normal faulting or submarine landsliding (Harper, 2006; MacDonald et al., 2006).

Where present, the hemipelagic sequence is overlain by a transition zone consisting of ~55 m of radiolarian-bearing argillite and minor sandstone (Fig. 2; Harper, 1994). MacDonald et al. (2006) interpreted a gradational contact between the hemipelagic sequence and the transition zone that is characterized by an up-section increase in sandstone beds and a decrease in the proportion of radiolarians in argillites. The transition zone is interpreted to have formed between 157 Ma and 153 Ma, based on age constraints from the hemipelagic sequence and reported detrital zircon ages of ca. 153 Ma from near the base of the upper turbidite sequence (Miller et al., 2003; MacDonald et al., 2006).

The transition zone is overlain by a thick sequence of turbidites with feldspathic to lithic wacke metasandstone and scant conglomerate (Fig. 2; Snoke et al., 1977; Harper, 1984; Wyld, 1985; MacDonald et al., 2006). In its northern outcrop area, the Galice Formation lacks the hemipelagic sequence and transition zone, and consists of turbidites that interfinger with volcanic members of the Rogue Formation (Fig. 1; MacDonald et al., 2006). Slaty units in the turbidites are typically thin, but reach at least 20 m thick locally (Frost et al., 2006) and contain trace fossils Chondrites, Cosmophorae, and Spirophycus that indicate abyssal to bathyal depths (MacDonald et al., 2006). Sandstone units include mud rip-up clasts, partial Bouma sequences, and scoured bases. Extensive deformation of the turbidite sequence precludes accurate determination of its thickness, but the unit was likely several kilometers thick, given its significant outcrop area (MacDonald et al., 2006).

The turbidite sequence contains the bivalve Buchia concentrica (Sowerby), which indicates late Oxfordian to middle Kimmeridgian depositional age (Imlay, 1952; Imlay et al., 1959). Two sills intruding the turbidite sequence northeast of Gasquet were dated by ⁴⁰Ar/³⁹Ar on hornblende, yielding ages of 150.5 ± 1.4 Ma and 146.2 ± 1.0 Ma (1-sigma uncertainty; Harper et al., 1994). The oldest pluton that intrudes the roof thrust (Orleans thrust) in outcrop is the 150 Ma Summit Valley pluton (multi-grain TIMS age; Harper et al., 1994). However, gabbro-diorite of the Bear Mountain complex (150–148 Ma; Snoke et al., 1981; Chamberlain et al., 2006) intrudes a klippe of the Rattlesnake Creek terrane (part of the eastern rift margin of the Josephine ophiolite basin; Snoke et al., 1977, 1981) that overlies the Galice Formation (Fig. 1). A 40 mgal positive gravity anomaly centered on the Bear Mountain complex (Roberts et al., 1981; Jachens et al., 1986) is best explained by the Bear Mountain complex cutting the roof thrust. Thus, according to biostratigraphy and known crosscutting relationships, the age of Galice turbidite deposition is constrained to between ca. 157 Ma and 150–148 Ma.

Deposition of the Galice Formation was originally interpreted to precede contractional deformation of the Nevadan orogeny (Diller, 1903). However, Harper et al. (1994) inferred synorogenic deposition, because both roof and basal thrusts were active from ca. 155 Ma to 150 Ma, synchronous with Kimmeridgian (ca. 154.8–149.2 Ma; Gradstein et al., 2020) deposition of Galice greywacke. Recent dating of migmatites formed above the roof thrust indicates that deformation began by 157 Ma (Dailey and Barnes, 2020), which suggests that Nevadan tectonism involved eastward thrusting of low-grade metasedimentary rocks beneath the Rattlesnake Creek terrane by 157 Ma, even as Galice deposition continued. Harper et al. (1994) suggested that deposition of Galice turbidites beginning at ca. 157 Ma represented a rapid increase in terrigenous sedimentation that indicated uplift in the sediment source area, linking turbidite deposition with Nevadan mountain building. Hacker et al. (1995) suggested that the Galice flysch formed in a “trench-like” setting related to the under-thrusting of the Josephine ophiolite to the west, with sediment derived from Middle Jurassic arc rocks and older units in the advancing Orleans thrust sheet to the east.

Nevadan orogenesis resulted in the development of slaty cleavage, folding, flattening of sand and pebble grains, and stretching lineations within the Galice Formation, with strain increasing southward (Harper, 1984, 2006; Cashman, 1988; Jones, 1988). Metamorphic grade likewise increases from prehnite-pumpellyite in the north to lower greenschist facies to the south (Harper, 1984; Harper et al., 1988). Paleomagnetic data suggest clockwise rotation of the Klamath Mountains by as much as 100° (Bogen, 1986; Harper and Park, 1986), which may have occurred as a single block (Mankinen and Irwin, 1982; Bogen, 1986) or through oroclinal bending (Renne and Scott, 1988; Saleeby and Harper, 1993). MacDonald et al. (2006) assumed clockwise rotation of 65°, which restores the trend of the Western Klamath terrane to parallel the structural trend of the Sierra Nevada foothills, which have not been rotated (Bogen et al., 1985; Frei, 1986). Paleoflow was variably west-directed in the northern Galice outcrop area and north-directed in the Smith River drainage (modern coordinates; Harper, 1984; Park-Jones, 1988; MacDonald et al., 2006), which restore to south-directed and west-directed, respectively, when corrected for 65° of rotation.

Galice turbidite detrital components are mostly siliceous argillite, chert, plagioclase, and volcanic rock fragments, and both monocrystalline and polycrystalline quartz (MacDonald et al., 2006). MacDonald et al. (2006) noted that detrital modes reported from the Galice Formation are semiquantitative, given matrix content upward of 20%, and often 30%, which suggests abundant pseudomatrix resulting from the alteration of original unstable lithic grains. Nevertheless, Galice samples plot largely within arc provenance fields on QFL and Q_mFL_t ternary plots (Harper, 1980; Norman, 1984; Wyld, 1985; MacDonald et al., 2006). Heavy mineral assemblages include zircon, tourmaline, apatite, biotite, muscovite, and Cr-spinel, and less commonly garnet, epidote, staurolite, and glaucophane. Volcanic-rich basal turbidites also include clinopyroxene and hornblende (Snoke, 1972, 1977; Harper, 1980, 1984; Wyld, 1985; MacDonald et al., 2006). Nd and Sr isotopic data from Galice argillites and metagreywackes indicate a significant continental isotopic signature, which suggests derivation from both arc rocks and previously accreted Eastern Klamath terranes and/or the North American continent (Frost et al., 2006). Detrital zircon age spectra published to date include one sample of primarily Mesozoic ages with age modes of 153 Ma and 227 Ma (Miller et al., 2003) and four samples (n = 341) with significant age modes at 158–157 Ma as well as abundant (15%–55%) Precambrian zircon with characteristic North American age modes (LaMaskin et al., 2022).

These compositional and age data suggest that the Galice Formation was sourced by both an active Mesozoic arc and continentally derived sediment that was perhaps recycled through previously accreted Klamath terranes and/or sources farther east on the continent. MacDonald et al. (2006) speculated that ca. 230–225 Ma detrital zircon in the Galice Formation was derived from the Pine Nut terrane and basal Luning assemblage of Western Nevada, which may have been east of the Klamath Mountains during Late Jurassic time (Wyld and Wright, 2001). LaMaskin et al. (2022) concluded that Galice sediment sources included previously accreted terranes of the Klamath Mountains and Sierra Nevada, recycled transcontinental sand, and primary and/or recycled sources in the Southwestern United States. Our robust detrital zircon age and ε_Hf data set, systematically sampled along 200 km of strike-length in the Klamath Mountains, combined with our new ε_Hf data from Klamath Mountains province Jurassic plutons and an extensive compilation of published detrital and igneous zircon age and ε_Hf data, permit more detailed assessment of the depositional age and provenance of Galice strata, which in turn improve evaluation of tectonic models and Late Jurassic paleogeography.

We collected 11 metasandstone samples from the siliciclastic turbidite section of the Galice Formation along ~200 km of strike length within the Klamath Mountains of California and Oregon (Fig. 1 and Table 1). Pervasive deformation and structural truncation precluded accurate determination of the stratigraphic level within the sections sampled. We collected samples with fine- to medium-sand–sized grains and visible quartz. Thin sections of each sample were manufactured at Texas Tech University and stained for potassium feldspar.

Samples of the Chetco complex (= Illinois River gabbro) and the marginal facies of the Grants Pass pluton were collected by J.D. Yule (Yule, 1996). Samples of the Grayback, Thompson Ridge, and Ashland plutons were collected by C.G. Barnes and R.F. Gribble (Gribble et al., 1990; Barnes et al., 1995). Samples were collected from homogeneous in situ bedrock outcrops. See Text S1 in the Supplemental Material¹ for sample locations.

We followed standard mineral separation methods to isolate zircon grains (e.g., DeGraaff-Surpless et al., 2002; Fedo et al., 2003), with slope Frantz current settings at 0.5 A, 0.8 A, 1.0 A, and 1.2 A to capture slightly magnetic, metamict zircon or zircon with inclusions in the final splits. The final zircon splits were sent to the Arizona LaserChron Center for isotopic analysis, where they were mounted in epoxy, polished to expose crystal interiors, and documented with backscattered electron images to guide spot selection in order to avoid cracks and inclusions.

Laser ablation–multicollector–inductively coupled plasma-mass spectrometry (LA-MC-ICPMS) was used to analyze ~300 randomly selected zircon grains per sample, or all zircon grains if the sample contained fewer than 300 (e.g., 21-GF-20 and 21-GF-27). We used a 30-μm-diameter spot and followed the analytical methods described in Pullen et al. (2018). Isotopic data were reduced using AgeCalcML version 1.42 (Sundell et al., 2021). Age data were filtered to remove from further consideration grains >500 Ma that were more than 20% discordant or 5% reverse discordant and grains <500 Ma with 2-sigma error ellipses that plotted off concordia and/or exceeded 5% analytical uncertainty at 1 sigma. Grains younger than 1000 Ma are reported as ²³⁸U-²⁰⁶Pb ages, and grains older than 1000 Ma are reported as ²⁰⁷Pb-²⁰⁶Pb ages, following Gehrels et al. (2006; Table S1). U-Pb age spectra are shown in histograms and kernel density estimates (KDE), which show the continuous relative probability of ages (plots generated using IsoplotR; Vermeesch, 2018).

In addition to visual comparison of KDE plots, we compare samples using nonmetric multi-dimensional scaling (MDS), a technique based on quantified pairwise comparisons of zircon ages in different samples (Vermeesch, 2013; Wissink et al., 2018); we use the Kolmogorov-Smirnov K-S statistic as the measure of dissimilarity in MDS comparisons. MDS describes the variation between sample-age frequency distributions by transforming pairwise similarities among age distributions into 2-D or 3-D Cartesian coordinates, such that the distance between samples represents the degree of dissimilarity (Vermeesch, 2013; Saylor et al., 2018; Wissink et al., 2018). Thus, the closer two samples are in MDS space, the more similar their age distributions are. Nonmetric MDS does not quantify the dissimilarities between samples, but instead approximates the relative ranks of the dissimilarities; the goodness-of-fit of the MDS solution is evaluated by the minimum stress parameter, with values of <0.1 considered fair, < 0.05 good, and <0.025 excellent (Vermeesch, 2013). MDS plots and stress values were generated using DZmds (Saylor et al., 2018).

Complex post-depositional deformation of Galice strata precludes determining the stratigraphic level of our samples to estimate their true depositional age (TDA). However, because the Galice Formation formed within an active convergent margin system, the calculated maximum depositional age (MDA) determined from the ages of the youngest zircon grains within each sample is likely to be close to the TDA because the time from crystallization to deposition may be negligible (Cawood et al., 2012; Coutts et al., 2019). Accurate determination of MDA has become a controversial topic as the use of detrital zircon MDAs has increased, with various methods for calculating MDA proposed (e.g., Dickinson and Gehrels, 2009; Spencer et al., 2016; Coutts et al., 2019; Herriott et al., 2019; Sharman and Malkowski, 2020; Vermeesch, 2021). We report the results of several of the common methods used to calculate the MDA: youngest single grain (YSG), weighted mean of youngest cluster of at least three grains that overlap at 1 sigma (Youngest Grain Cluster [YGC] 1σ) or 2 sigma (YGC 2σ; Dickinson and Gehrels, 2009), the youngest statistical population with a mean square weighted deviation close to 1.00 (YSP; Coutts et al., 2019), and the maximum likelihood age (MLA; Vermeesch, 2021; Table 1; Fig. 3). Following Coutts et al. (2019), the YGC 1σ and YGC 2σ methods limit the sub-sample used to calculate the MDA by the age of the youngest upper limit of uncertainty, and means are weighted by date uncertainty. All MDA ages are reported with 2σ uncertainties, and MDAs calculated from pooled analyses include both analytical (internal) and systematic (external) uncertainties; YSG ages include only analytical uncertainty.

Some of our samples have negatively skewed age populations in the youngest age peaks that could indicate Pb-loss (e.g., Spencer et al., 2016), and we cannot rule out Pb-loss in the youngest single grains. We therefore removed from MDA calculations any very young outliers that plot away from the youngest cluster of grains to avoid biasing MDA estimators (0–4 grains per sample; Fig. 3). All of our samples have abundant near-depositional–age grains, which increases the risk that MDAs determined only from grains within the young tail of an age distribution (e.g., YGC 1σ and YGC 2σ) will be younger than TDA because of the increased effect of statistical outliers in these calculations (Herriott et al., 2019). Therefore, we consider YSP values to be the best determinants of MDA in this study, although we acknowledge that the YSP values may be older than the TDA values if young grains did not experience Pb-loss. Following Herriott et al. (2019), we used the YSP to represent coeval zircon crystallization, without ascribing that crystallization to a singular geologic event (e.g., volcanic eruption).

Zircon from the Chetco complex and Wooley Creek suite plutons was originally dated by isotope dilution–thermal ionization mass spectrometry (ID-TIMS) on multi-grain aliquots of individual size and magnetic fractions at the California Institute of Technology (Yule, 1996). Zircon and titanite from the Bear Mountain complex were dated using ID-TIMS methods by Chamberlain et al. (2006). Reported ages were determined by a range of methods, including single-crystal and multi-grain analyses, some of which used titanite data. A summary of reported ages is presented in Text S1, and the reader is referred to Chamberlain et al. (2006) for further details. Splits of these same zircon separates were mounted in epoxy, imaged by cathodoluminescence at Texas Tech University, and analyzed for U-Pb by LA-MC-ICP-MS at the University of California, Santa Barbara. Details of the analytical methods may be found in Barnes et al. (2021). Sample descriptions, reported U-Pb ages from Chamberlain et al. (2006) and Yule (1996), and results of LA-MC-ICP-MS dating (Tera-Wasserburg concordia diagrams and weighted mean plots) are given in Text S1.

Lu-Hf isotopic data were measured from 603 selected Triassic and Jurassic detrital zircon grains from all 11 metasandstone samples. We selected Mesozoic zircon grains that were large enough to accommodate a 40-μm spot in addition to the spot used for the U-Pb analysis. We followed the analytical methods described in Gehrels and Pecha (2014), and used HfcalcML to complete data reduction. We converted Lu-Hf data to epsilon (ε) units, and removed from further consideration any grains with 1-sigma uncertainties of >2 epsilon units. For Triassic and Jurassic zircon, ε_Hf values >9 are considered juvenile or highly radiogenic, between 3 and 9 are moderately radiogenic, and <3 are considered to be evolved (Bahlburg et al., 2011; Gehrels and Pecha, 2014). We used Hf Plotter (Sundell et al., 2019) to produce contour plots for the new ε_Hf data from this study as well as compilations of published data, and DZstats2D (Sundell and Saylor, 2021) to generate MDS comparisons of combined U-Pb age and ε_Hf zircon data.

Lu-Hf isotopes on plutonic zircons were measured by LA-MC-ICP-MS at the University of California, Santa Barbara (see Barnes et al., 2021). Unlike samples from the Wooley Creek batholith (Barnes et al., 2021), six of the nine samples analyzed for this project yielded ε_Hf values higher than those of modern mid-oceanic-ridge basalts. These samples were reanalyzed using different standards (Table S2). This means that for the three samples that were not reanalyzed, age and Hf isotope spots are coincident and individual spot ages were used to calculate ε_Hf, whereas for the remaining six samples, age and Hf isotope spots are not coincident. In the latter case, the weighted mean average of grains whose ²³⁸U/²⁰⁶Pb–²⁰⁷Pb/²⁰⁶Pb ages are less than 10% discordant were used to calculate ε_Hf.

Metasandstone samples are medium- to fine-grained, and most display phyllitic to schistose textures (Fig. 4). Metamorphic grade varies from subgreenschist to greenschist, as indicated by the presence of fine-grained biotite and white mica. The least-deformed samples display weak foliation and consist of variable proportions of quartz, plagioclase, and volcanic clasts ± fossil fragments (e.g., 21-GF-20) ± rip-up clasts of very fine-grained sediment. With increasing deformation and metamorphism, as determined from the degree of fabric development and mineral assemblages, the proportions of quartz and plagioclase vary slightly, but the abundance of “matrix” material increases. This “matrix” material consists of fine-grained white mica + quartz ± biotite ± unidentified phases. The increasing abundance of “matrix” with deformation suggests that much of the “matrix” is extensively metamorphosed and recrystallized lithic fragments, i.e., pseudomatrix (Dickinson, 1970), rather than a primary mud matrix. Because more than 20% of each sample is pseudomatrix, samples likely have more abundant quartz and fewer lithic grains than their original detrital composition, rendering quantitative point counts inaccurate. Accessory minerals include pyrite, hematite (oxidized pyrite?), tourmaline, scant apatite, and zircon. None of the samples contains alkali feldspar.

The KDE plots for all samples are shown stacked by MDA (YSP method), with oldest on the bottom (Fig. 5); MDAs range from 160.1 ± 1.6 Ma to 150.8 ± 1.9 Ma. Because the MDA is likely close to the true depositional age and the Galice Formation is too complexly deformed to determine the relative location of the samples in the stratigraphy, we use MDA as a proxy for stratigraphic level. All KDE plots show a significant Middle to Late Jurassic age mode, with younging of the Late Jurassic age mode corresponding to younger MDAs, which is consistent with deposition within an active convergent margin (Fig. 5; Cawood et al., 2012). Most samples have few to no Triassic and Early Jurassic grains, with the exception of two samples (21-GF-25 and 21-GF-27), which both include a second Mesozoic age mode at ca. 199 Ma, and a range of grains between 250 Ma and 230 Ma in age. Sample GASD-18-09 also contains Early Jurassic grains, but these form part of a broader distribution of Jurassic grains in the sample rather than a distinct age mode. Most samples include abundant Paleozoic and Precambrian grains, although samples display different age modes and different relative proportions of similar age modes. The percentage of grains that are Mesozoic in age varies significantly between samples, from 25% to 90% (Fig. 5). Overall, Galice Formation detrital zircon age spectra display considerable sample-to-sample variability. MDS comparison demonstrates minimal clustering of samples and no systematic correlation between location along strike and dissimilarity (Fig. 6).

With the exception of sample AP-10, data reported here represent new LA-MC-ICP-MS analyses of zircon grains originally dated by Yule (1996) and Chamberlain et al. (2006). Redating of these samples was conducted as part of our analyses of Hf isotopes. Analytical data are included in Supplemental Data S2. All results reported below are ²⁰⁶Pb/²³⁸U ages with 2σ uncertainties. Comparisons of the original results with our new LA-MC-ICP-MS ages are provided in Text S1.

Four samples from the Chetco complex were dated. Sample IR-3 is a hornblende diorite. Fifteen concordant analyses yielded a weighted mean age of 158.2 ± 1.3 Ma (MSWD = 2.0). Sample IR-4 is an agmatitic diorite from which eight concordant grains gave an age of 163.2 ± 2.0 Ma (MSWD = 1.19). The remaining zircons from this sample displayed discordance on either side of concordia curve. Two tonalite bodies from the Chetco complex were dated, IR-1 and DY-PP-64. Sample IR-1 displayed considerable reverse discordance, but five concordant grains yielded an age of 159.3 ± 2.4 Ma (MSWD = 1.04). Fourteen concordant zircons from DY-PP-64 yielded an age of 158.4 ± 1.2 Ma (MSWD = 1.02).

Three samples from the Wooley Creek suite (one sample each from the Grayback, Thompson Ridge, and Ashland plutons) and one sample from the Grants Pass pluton were analyzed. Thirteen grains from sample GBP-1688, a hornblende diorite of the main stage of the Grayback pluton (Barnes et al., 1995), yielded an age of 161.9 ± 1.8 Ma (MSWD = 1.6). Seventeen concordant U-Pb dates of zircons from Thompson Ridge pluton sample TR-4, a biotite–hornblende–pyroxene gabbro, gave a weighted mean age of 159.6 ± 1.6 Ma (MSWD = 1.9). Six concordant zircons from biotite–hornblende quartz monzodiorite sample AP-79 from the Ashland pluton (Gribble et al., 1990) yielded a weighted mean average of 156.1 ± 1.8 Ma (MSWD = 0. 81). Sample AP10 is a foliated hornblende diorite collected near the intrusive contact between the Grants Pass pluton and the andalusite slate of the Galice Formation. Twelve concordant zircons from this sample yielded a weighted mean average age of 148.1 ± 1.2 Ma (MSWD = 1.15).

Zircon from three samples from the Bear Mountain complex (Snoke et al., 1981; Barnes et al., 2006; Chamberlain et al., 2006) was analyzed. Ten concordant grains from hornblende–pyroxene gabbro sample DP4-8 gave a weighted mean average age of 146.3 ± 1.3 (MSWD = 0.78). Five concordant grains from biotite–pyroxene monzodiorite sample DP1-13 gave a weighted mean age of 151.6 ± 2.3 Ma (MSWD = 0.013). Ten concordant grains from biotite–hornblende tonalite sample DP2-20, which is intrusive into the Blue Ridge pluton (part of the Bear Mountain complex), yielded a weighted mean age of 147.8 ± 0.9 Ma (MSWD = 0.78).

Lu-Hf isotopic analyses reveal a wide range of ε_Hf values for 603 Mesozoic detrital zircon grains in the Galice Formation, with most grains falling between +14 and −15 in all samples (Fig. 7A). The majority of the grains analyzed in all samples are Middle to Late Jurassic (<174.7 Ma; n = 416; 69%), of which 36% are juvenile (ε_Hf > 9), 36% are moderately radiogenic (9 > ε_Hf > 3), and 28% are evolved (ε_Hf < 3). The remaining 31% (n = 187) of grains analyzed are Triassic and Early Jurassic, of which 26% are juvenile, 33% are moderately radiogenic, and 42% are evolved.

The ε_Hf values of zircon from most Middle–Late Jurassic plutons in the Klamath province vary significantly (Fig. 7B). For example, ε_Hf values in Chetco complex sample 91DY-PP-64 range from +16.6 to +8.6 (Fig. 7B). Nevertheless, averaged ε_Hf values display distinct differences on the basis of age and location in the province. The four samples from the Chetco complex have average ε_Hf ranging from +17.4 to +13.0. In contrast, three samples from the Wooley Creek suite (Grayback, Thomson Ridge, and Ashland plutons) have average ε_Hf ranging from +9.9 to +8.4. These values are consistent with data from the Wooley Creek batholith and Slinkard pluton, which range from +10.8 to +2.0 in main-stage mafic–intermediate rocks and are as low as −3.1 in granite (Barnes et al., 2021). Gabbro and monzodiorite from the Bear Mountain complex yield average ε_Hf values of +11.2 and +11.5, respectively, and the tonalite intrusive into the Blue Ridge pluton gives average ε_Hf of +10.5. The ε_Hf value of the early-stage sample of the Grants Pass pluton averages +10.9.

Because we compare calculated geological ages (MDAs) with the chronostratigraphic time scale to place our results in a larger tectonic context, the selection of time scale is critical. The Geologic Time Scale 2020 (Gradstein et al., 2020) includes revision of the Middle and Late Jurassic stage boundaries by Hesselbo et al. (2020), and these updates were incorporated in the International Chronostratigraphic Chart v2022/10 (Cohen et al., 2013) at https://stratigraphy.org. Bounding ages for the Oxfordian, Kimmeridgian, and Tithonian stages are 2.0–2.9 m.y. younger than on previous time scales, and they are now calibrated to magnetostratigraphy (Hesselbo et al., 2020). The Geological Society of America Geologic Time Scale v. 6.0 (Walker and Geissman, 2022) does not yet reflect these updates, so we follow the Gradstein et al. (2020) time scale herein.

Our detrital zircon MDA data suggest that turbidite deposition began as early as ca. 160.1 Ma during Oxfordian time and continued through Kimmeridgian time, which is consistent with Oxfordian to Kimmeridgian (Gradstein et al., 2020) biostratigraphic age constraints on Galice deposition based on the pelecypod Buchia concentrica (Imlay et al., 1959), and nearly synchronous with the 150.5 ± 1.4 Ma sill in the hemipelagic sequence (Harper et al., 1994). In addition, LaMaskin et al. (2022) report detrital zircon age data from four Galice Formation turbidite samples. YSP MDAs calculated from their data following the same MDA methods used here result in MDAs of 159.5–154.5 Ma, which is consistent with our results. These depositional ages document a 9 m.y. duration of turbidite deposition in the Galice basin and suggest either more rapid deposition of the ca. 162–154.8 Ma underlying hemipelagic and transition zone strata or diachronous onset of turbidite deposition across the Galice basin. Furthermore, age constraints on the formation of the Josephine ophiolite from 164 Ma to 162 Ma, followed by deposition of the overlying hemipelagic and transition zone from ca. 162 Ma to 155(?) Ma, and the beginning of turbidite deposition as early as ca. 160 Ma, suggest that the possible disconformity at the top of the hemipelagic sequence (see above; Pessagno and Blome, 1990) may be a local feature, and/or related to normal faulting or submarine sliding (e.g., Harper, 2006) rather than a true depositional hiatus. Our results show that deposition of the upper turbidite sequence likely was established by 160 Ma, during Oxfordian time (ca. 161.5–154.8 Ma; Gradstein et al., 2020), if not earlier.

Our results are consistent with interpretations of the Galice Formation as a synorogenic deposit (Harper et al., 1994; Hacker et al., 1995), with deposition continuing even as older strata were underthrust eastward beneath the Rattlesnake Creek terrane during Nevadan orogenesis. Thus, the youngest strata of the Galice Formation were deposited after the initiation of thrusting along the Orleans thrust and nearly synchronous with crosscutting of that roof thrust by the 150 Ma Summit Valley pluton (Harper, 2006) and the 150–148 Ma Bear Mountain complex (Snoke et al., 1981; Chamberlain et al., 2006).

Our Galice Formation detrital zircon age results display considerable inter-sample variability and no discernable trends along almost 200 km of strike length (Fig. 6). In addition, post-depositional deformation makes determining stratigraphic level and accurately correlating among sample locations a challenge. To avoid over-interpreting variability between individual samples that may result from natural hydrodynamic sorting processes rather than distinct provenance (Ibañez-Mejia et al., 2018), and lacking means to accurately correlate among sample localities, we herein compile our samples into one KDE (Fig. 8). The compiled Galice age spectrum includes 50% Mesozoic grains, 10% Paleozoic grains, and 40% Precambrian grains (Fig. 8). The Mesozoic age spectrum includes a large Middle–Late Jurassic mode (ca. 156 Ma), and smaller Early Jurassic (ca. 199 Ma) and Triassic (ca. 240 Ma) modes (Fig. 8). The polymodal Paleozoic age spectrum spans the entire era, with age modes at ca. 263 Ma, 315 Ma, and 430–400 Ma (Fig. 8). The polymodal Precambrian age spectrum includes modes at ca. 607 Ma, 1200–980 Ma, 1435 Ma, 1850–1610 Ma, and 2680 Ma (Fig. 8). ε_Hf data from 603 Mesozoic detrital zircon grains reveal a wide range of ε_Hf values, from +16.9 to −26.2 (Fig. 7A). These data suggest mixed provenance in an active Middle to Late Jurassic magmatic arc with abundant input from Precambrian detrital zircon presumably recycled from older sources.

To document the provenance of the Galice Formation, we compare our Galice zircon age and ε_Hf data with compilations of previously published and new U-Pb age and ε_Hf data for rocks older than 148 Ma from potential sources in the U.S. Cordillera. Following LaMaskin et al. (2022), our compilations present igneous spot data wherever possible (i.e., for all data generated by SIMS) rather than interpreted intrusive ages or ε_Hf values, to better match our detrital zircon age distributions with similar distributions from plutonic sources. We acknowledge that our data compilations are limited by sampling density, sampling bias, and the availability of published data, and may require revision as more data become available. However, these compilations represent known ages and ε_Hf ranges within each of these source regions, and are therefore useful for first-order provenance comparisons with our compiled data from the Galice Formation.

Plausible potential sediment sources for the Galice Formation are shown in Figure 1A and include, from north to south: (1) arc rocks of the Northwest Cascades System and the Methow Terrane in Washington State (Sauer et al., 2017); (2) arc rocks and sedimentary cover strata of the Blue Mountains province in Oregon (Schwartz et al., 2011a, 2011b, 2014; LaMaskin et al., 2011; Anderson, 2013; Johnson et al., 2015; Kurz et al., 2017 [whole-rock ε_Nd data converted to zircon ε_Hf values using the terrestrial array equation of Vervoort et al., 2011]); (3) arc rocks and older terranes within the Klamath Mountains province (Wallin and Metcalf, 1998; Gehrels and Miller, 2000; Wallin et al., 2000; Allen and Barnes, 2006; Barnes et al., 2006; Chamberlain et al., 2006; Johnson and Barnes, 2006; Yule et al., 2006; Scherer and Ernst, 2008; Scherer et al., 2010; Ernst et al., 2017); (4) the Sierra Nevada arc and pre-batholithic framework (Harding et al., 2000; Spurlin et al., 2000; Barth et al., 2011; Cecil et al., 2012, 2018; Attia et al., 2018, 2020); (5) Sierra Nevada retroarc plutons, Nevada–Utah backarc basin, and passive margin strata (Darby et al., 2000; Manuszak et al., 2000; Riley et al., 2000; Workman, 2012; Gehrels and Pecha, 2014; Colby, 2017; Holm-Denoma et al., 2017; Chapman et al., 2018); and (6) Jurassic Eolianites of Nevada and Utah, USA (Dickinson and Gehrels, 2008, 2009). In addition, arc rocks and older strata of the Insular superterrane (White et al., 2016, and references therein; Alberts et al., 2021, and references therein) represent a potential source region if the Insular superterrane accreted to North America during Middle Jurassic time at the approximate latitude of the Klamath Mountains (e.g., Saleeby and Dunne, 2015) rather than remaining far offshore until Cretaceous accretion (e.g., Balgord et al., 2021; Tikoff et al., 2023).

We first compare pre-Mesozoic age spectra (grains >250 Ma) of the Galice Formation with compilations of potential source region data (Fig. 9). A visual comparison of KDE plots suggests that the Galice Formation is very similar to the Sierra Nevada pre-batholithic framework, Blue Mountains Overlap, and Jurassic Eolianites (Fig. 9A). The older Klamath terranes that are most proximal to the Galice Formation do not provide as good a match for Galice age spectra because they include a large late Paleozoic mode that is not well-represented in the Galice Formation and lack the 1200–980 Ma and 1435 Ma modes that characterize the Galice samples (Fig. 9A).

MDS analysis confirms the visual assessment of KDEs (Fig. 9B). The Blue Mountains Overlap and Sierra Nevada pre-batholithic framework plot closest to the Galice Formation in MDS space, and the Jurassic Eolianites plot very close to the Blue Mountains Overlap (Fig. 9B). The MDS plot shows that the older Klamath terranes plot farther from the Galice Formation than either the Blue Mountains Overlap or Sierra Nevada pre-batholithic framework sources, but also that the Galice Formation is the closest neighbor of the older Klamath terranes (Fig. 9B). The Nevada–Utah backarc basin and passive margin region is most closely related to the Sierra Nevada pre-batholithic framework, but it lacks the prominent Paleozoic age modes that characterize the Galice Formation and therefore plots farther away. The Insular superterrane plots farthest from the Galice Formation, with the older Klamath terranes identified as its nearest neighbor (Fig. 9B).

These provenance comparisons indicate that distal sources were important contributors of pre-Mesozoic zircon to the Galice basin, in addition to proximal sources in the older Klamath terranes. Thus, the pre-Mesozoic age spectrum compiled from the Galice Formation suggests mixed provenance in continental sources.

We combine zircon U-Pb age and ε_Hf data of Mesozoic grains to further assess potential source regions, because the Mesozoic arc extended the length of the Cordillera and therefore age alone is nondiagnostic. We present compiled zircon age and ε_Hf data for Mesozoic igneous rocks in the potential source regions, including our new Hf data from Jurassic plutons in the Klamath Mountains (Fig. 10A), and use MDS to compare these compilations with data from the Galice Formation (Fig. 10B). ε_Hf data compilations are shown as shaded fields that represent 95% of the compiled data, and we consider all of the Galice ε_Hf data in aggregate (Fig. 10A). We note that the density of published ε_Hf data varies significantly by location, with the Sierra Nevada better characterized than other possible source regions; therefore, as more data become available, the relationship between these potential source regions and the Galice Formation may need to be reconsidered.

Arc sources north of the Galice Formation in the Northwest Cascades System, Methow Terrane, and the Blue Mountains, as well as the proximal Klamath arc, lack evolved zircon and therefore cannot fully account for the Galice Formation ε_Hf signature. However, the juvenile Middle–Late Jurassic zircon in the Klamath and Blue Mountains arcs matches the cluster of juvenile zircon grains in the Galice Formation (Fig. 7), and these proximal arcs likely represent a component of the sediment sources of the Galice basin. These three potential sources plot as a cluster in MDS space, with the Blue Mountains and Klamath arcs plotting nearly on top of one another, which suggests identical arc ε_Hf signatures for Middle–Late Jurassic plutons in these two regions (Figs. 10B–10E).

The Sierra Nevada arc can account for some of the Middle–Late Jurassic juvenile grains, although not the highest ε_Hf values, as well as many of the evolved grains of Triassic, Early Jurassic, and Middle–Late Jurassic age (Fig. 10A). The Sierra Nevada retroarc plutons lack abundant juvenile grains but can partially account for the most evolved Middle–Late Jurassic grains in the Galice Formation. However, unlike the Sierra Nevada arc and retroarc sources, the Insular superterrane source region can account for much of the ε_Hf data from the Galice Formation and plots as the nearest neighbor of the Galice Formation in MDS space (Figs. 10B–10E). As noted above, the Insular superterrane is not a good match for pre-Mesozoic detrital zircon ages in the Galice Formation, but we acknowledge the possibility that the Insular superterrane was part of the mixed provenance of the Galice Formation. In that case, the Paleozoic grains that characterize the pre-Mesozoic age spectrum of the Insular superterrane were diluted in relative abundance, as other sources contributed more abundant Proterozoic detrital zircon to the Galice basin. In any case, a few of the Galice Mesozoic zircon grains are not accounted for by the combination of the Klamath arc, Sierran arc, Sierran retroarc, and Insular superterrane, which suggests that the data representing each source is incomplete and/or our potential sources are not comprehensive.

Taken together, pre-Mesozoic zircon age and Mesozoic zircon ε_Hf comparisons suggest likely sources of pre-Mesozoic zircon in the Sierra Nevada pre-batholithic framework, Blue Mountains Overlap, older Klamath terranes, and possibly Jurassic Eolianites, and sources of Mesozoic zircon in the Blue Mountains, Klamath, and Sierran arcs, and perhaps retroarc plutons. However, in the Blue Mountains province, LaMaskin et al. (2015) used detrital zircon age data to revise the Middle Jurassic depositional age of the fluvial to marine Coon Hollow Formation to Late Jurassic, recognizing that Middle Jurassic corals were transported into the younger Coon Hollow Formation during deposition that occurred from 160 Ma to 150 Ma. Thus, Middle Jurassic strata of the Blue Mountains Overlap may have been submerged during Late Jurassic time and were therefore unlikely to be a significant source of sediment of the Galice Formation. Although the Insular superterrane is a poor match for the pre-Mesozoic Galice age spectra because it lacks abundant Proterozoic grains (Fig. 9), the Insular superterrane arc is a close match for the arc-derived grains of the Galice Formation (Fig. 10) and may have been an additional source of sediment.

Our provenance interpretations for the Galice Formation have implications for Late Jurassic tectonics in the U.S. segment of the North American Cordillera. We reconstructed Late Jurassic paleogeography by restoring Cenozoic extension in the Basin and Range (following Wyld et al., 2006); Cretaceous dextral displacement within the Sierra Nevada arc, Western Nevada, and Western Idaho (e.g., Wyld and Wright, 2001); clockwise rotation of the Blue Mountains province (Wilson and Cox, 1980; Housen and Dorsey, 2005); and clockwise rotation (Bogen, 1986; Harper and Park, 1986) and post-Jurassic westward motion of the Klamath Mountains (Wyld and Wright, 2001; Ernst, 2013). Although our data do not require it, we show the Insular superterrane located offshore of North America and north of the Klamath Mountains province in Late Jurassic time, following the model of Saleeby and Dunne (2015). East of the Cordilleran arc, late-phase deformation in the Eastern Luning-Fencemaker thrust belt during Middle to Late Jurassic time led to uplift and erosion of Triassic backarc basin strata (Oldow, 1983; Wyld, 2002; LaMaskin et al., 2011). The resulting Late Jurassic paleogeography (Fig. 11) is similar to the restorations of Saleeby and Dunne (2015), Yonkee and Weil (2015), and Balgord et al. (2021).

Changing plate motions during Late Jurassic time led to increased contractional deformation in the region as well as probable sinistral transpression, as North American plate motion shifted from west to northwest and plate velocity increased (e.g., Seton et al., 2012; Saleeby and Dunne, 2015). This contractional deformation marked the beginning of Nevadan orogenesis in the Klamath province. Closure of the Galice basin began by ca. 155 Ma (Harper and Wright, 1984; Harper et al., 1994), and perhaps as early as 157 Ma (Dailey and Barnes, 2020), as rocks of the Galice/Josephine assemblage and perhaps the outer Condrey Mountain Schist (Saleeby and Harper, 1993) were thrust beneath inboard terranes along the Orleans thrust, even while younger Galice strata continued to be deposited (Harper et al., 1994; this study). To the west–northwest of the Galice basin, arc activity continued in the Rogue-Chetco complex until 157 Ma, and then these rocks and their basement assemblage were thrust beneath the Galice/Josephine assemblage along the Madstone Cabin thrust (Harper et al., 1996). In addition, Saleeby and Dunne (2015) suggested that sinistral motion associated with increased northwestward motion of the North American plate during this time brought the Insular superterrane to a position north of the Klamath Mountains province and west of the Blue Mountains province and Northwest Cascades system (Fig. 11).

The abundance of pre-Mesozoic zircon derived from sources in the Sierran arc pre-batholithic framework and Sierran retroarc regions suggests that sediment transport systems crossed the Cordilleran magmatic arc throughout Late Jurassic time. Thus, during Oxfordian–Kimmeridgian time, fluvial systems with headwaters well into the Sierran retroarc (e.g., the Luning-Fencemaker thrust belt) may have traversed the Klamath-Sierran arc, delivering sediment into the shallow marine system, with turbidity flows funneled down submarine canyons into the bathyal to abyssal Galice basin (Fig. 11). Although active arc volcanoes may have been emergent (Young, 1978; Garcia, 1979, 1982), the older Klamath terranes east of the Galice basin and the Blue Mountains province may have been at least partially submerged, perhaps maintaining some of the interconnected marine system across the Klamath-Sierran arc that characterized the Middle Jurassic (Attia et al., 2021), with deposition of the Coon Hollow Formation overlapping some older Blue Mountains province rocks during Late Jurassic time (LaMaskin et al., 2015). Mesozoic arc magmatism within the proximal Klamath terranes, Sierra Nevada arc, and possibly the Insular superterrane may have contributed Triassic, Early, and Middle–Late Jurassic zircon with wide-ranging ε_Hf values. Uplift of the Klamath Mountains province terranes and the embedded Wooley Creek suite arc system east of the Galice basin during Nevadan orogenesis, as well as arc activity in the Rogue-Chetco complex, likely contributed to Late Jurassic arc-derived zircon with abundant juvenile ε_Hf values in the Galice basin.

Our tectonic model for Late Jurassic time postulates fluvial systems with headwaters east of the Klamath-Sierran arc that transported sediment westward into the forearc region during Oxfordian–Kimmeridgian time (Fig. 11). Our model predicts that Oxfordian–Kimmeridgian strata west of the Sierran arc share provenance with the Galice Formation, while forearc strata north of the Klamath Mountains, in the Blue Mountains and Northwest Cascades System and Methow Terrane, have distinct provenance.

The Mariposa Formation crops out in the Western Sierra Nevada foothills region (Fig. 11) and has been correlated with the Galice Formation by numerous authors (e.g., Diller, 1907; Taliaferro, 1942; Imlay, 1952; Irwin, 2003; Snow and Scherer, 2006; Ernst et al., 2008; Ernst, 2013). Like the Galice Formation, the Mariposa Formation has been recrystallized to subgreenschist grade (Snow and Ernst, 2008), but protolith lithologies and sedimentary structures remain recognizable. Bogen (1984) documented sole markings, mudstone rip-up clasts, conglomerate lenses, and partial to nearly complete Bouma sequences in sandstone beds with consistent southeasterly flow directions (modern coordinates) in the southern outcrop area. Snow and Ernst (2008) reported SHRIMP-RG detrital zircon ages from five samples, which all contain Early to Late Jurassic zircon and abundant pre-Mesozoic zircon (Fig. 12A). The youngest 10 detrital zircon grains from a composite of all five samples yield a weighted mean age of 152.1 ± 1.1 Ma, and individual reported MDAs range from 155 Ma to 151 Ma (Snow and Ernst, 2008), which suggests Kimmeridgian deposition (time scale of Gradstein et al., 2020).

The Great Valley forearc basin developed in latest Jurassic time between the Franciscan subduction zone to the west and the Sierra Nevada magmatic arc to the east (e.g., Dickinson, 1995). Recent detrital zircon results from the basal Great Valley Group confirm Late Jurassic deposition for limited regions, with much of the mapped Upper Jurassic strata deposited in Early Cretaceous time (Surpless et al., 2006; Orme and Surpless, 2019), and suggest that transtension within the forearc region resulted in isolated, fault-bounded depocenters during latest Jurassic time. Accordingly, we limit our detrital zircon age compilation to only samples of the Great Valley Group that have Kimmeridgian MDAs (Fig. 12A; time scale of Gradstein et al., 2020; Surpless et al., 2006; Orme and Surpless, 2019).

In the Blue Mountains Province, LaMaskin et al. (2015) presented detrital zircon age data from the fluvial to marine Coon Hollow Formation. Based on similarities in depositional age, lithology, and limited provenance data, LaMaskin et al. (2015) suggested that the Coon Hollow Formation and Galice Formation were part of the same belt of suprasubduction-zone extensional back-arc basins during Late Jurassic time.

Farther north, Late Jurassic sedimentary units include the Lower Newby Group (Twisp Formation) of the Southern Methow subterrane in North-Central Washington (Sauer et al., 2017) and sedimentary units within the fault-bounded nappes of the Northwest Cascades System. These include the Fidalgo complex, Lummi Formation, Constitution Formation, and Easton Metamorphic suite, which each have nearly unimodal, Mesozoic detrital zircon age spectra (Brown and Gehrels, 2007), and the Yellow Aster complex, with Paleozoic and Precambrian detrital zircon age spectra (Brown and Gehrels, 2007). Here, we compare the pre-Mesozoic detrital zircon age spectra of the Galice Formation with ages compiled from Late Jurassic units of the Northwest Cascades System (Yellow Aster complex; Brown and Gehrels, 2007) and the Methow Terrane (Twisp Formation; Sauer et al., 2017).

Visual comparison of KDE plots of pre-Mesozoic detrital zircon suggests that the Galice and Mariposa formations share very similar provenance, as does the basal Great Valley Group (Fig. 12A). In contrast, the samples compiled from the Coon Hollow Formation, Northwest Cascades System, and Methow Terrane look very different from any other samples (Fig. 12A). In MDS space, the Galice and Mariposa formations plot nearly on top of one another, with the basal Great Valley Group in close proximity (Fig. 12B). Samples compiled from the Northwest Cascades System, Methow Terrane, and the Coon Hollow Formation plot much farther away from any other samples, which reflects their distinctive detrital zircon age distributions (Fig. 12B). Taken together, detrital zircon age comparisons from coeval strata are consistent with our proposed tectonic model for the Late Jurassic margin.

One aspect of long-standing debate about the tectonic development of the Mesozoic convergent margin of North America focuses on whether much of Western North America is the product of west-dipping, intraoceanic subduction during Jurassic–Cretaceous time (Johnston, 2001, 2008; Hildebrand, 2009, 2013, 2015; Sigloch and Mihalynuk, 2013, 2017; Spencer et al., 2019; Clennett et al., 2020), or east-dipping subduction of oceanic plates beneath Western North America (e.g., Monger and Price, 1979; Burchfiel et al., 1992; Dickinson, 2008; Ingersoll, 2008; Saleeby and Dunne, 2015; Boschman et al., 2018; Pavlis et al., 2019). For example, Sigloch and Mihalynuk (2013, 2017) and Clennett et al. (2020) interpret seismically imaged, linear features that now reside in the lower mantle beneath Eastern North America and the Atlantic Ocean as near-vertical slab walls associated with subduction at stationary, intraoceanic trenches that were located far west of Jurassic North America. In this interpretation, the Western Klamath terranes formed above an east-facing subduction zone and are considered exotic to North America and were later accreted to North America as the westward-migrating continent overrode the oceanic trench and collided with the archipelago during the opening of the Atlantic Ocean (e.g., Sigloch and Mihalynuk, 2017). In contrast, endemic models postulate that long-standing east-dipping subduction beneath the continental margin produced the outboard arcs and associated terranes, and therefore these endemic terranes should be clearly linked to the continent (e.g., Gray, 1986; LaMaskin et al., 2022). Based on their interpretations of Galice Formation and Rattlesnake Creek Terrane detrital zircon provenance, LaMaskin et al. (2022) argued an endemic origin for these outboard terranes. Our provenance results expand this initial data set and confirm that Galice Formation provenance was linked to the North American continent throughout deposition. Furthermore, we document the development of a composite marginal basin system that included the Galice and Mariposa formations and the basal Great Valley Group and place this basin system within a Late Jurassic tectonic reconstruction.

Deposition of the Galice Formation occurred during the Late Jurassic transition from transtensional to contractional (transpressional?) tectonics in the U.S. Cordillera (Harper et al., 1986, 1994; Saleeby, 1992; Saleeby and Dunne, 2015). Our results document abundant detrital zircon derived from sources within the retroarc region, in addition to magmatic arc sources. Furthermore, our data demonstrate that pre-Mesozoic detrital zircon in the Galice Formation apparently was not derived from older Klamath terranes or the Blue Mountains province (Fig. 9). These results suggest that the Klamath Mountains province, Blue Mountains, and the Western Sierra Nevada had muted topography and/or were covered by a Mesozoic volcanic carapace that was relatively zircon-poor. In the Blue Mountains, LaMaskin et al. (2015) document deposition of the transgressive, fluvial to deep-marine Coon Hollow Formation during a period of trench retreat at ca. 160–150 Ma, which suggests that the Blue Mountains Province was the locus of deposition during this time. In their stratigraphic overview of the Mesozoic Sierra Nevada, Attia et al. (2021, and references therein) documented deep-marine sedimentation in the Southern Sierra Nevada and Northern and Central Western metamorphic belt through ca. 157 Ma and 152 Ma, respectively, and shallow-marine deposition in the Eastern Sierra Nevada through ca. 148 Ma. Furthermore, Attia et al. (2021) inferred the presence of an integrated marine depositional system across the entire Sierra Nevada by Middle Jurassic time, which suggests that the arc only later became emergent through differential uplift related to Late Jurassic contraction.

We suggest that fluvial systems with headwaters in the retroarc region traversed the Klamath-Sierran arc, transporting detritus from both the retroarc region and the pre-batholithic framework of the Eastern and Central Sierra Nevada into deep-marine systems in the Galice basin. With changing plate kinematics during Late Jurassic time, west-directed transport systems continued to traverse the rising arc, even as the emergent arc provided abundant sediment to the forearc.

The Middle Jurassic Siskiyou orogeny in the Klamath Mountains province involved amalgamation of terranes that constitute the Western Paleozoic and Triassic belt, and accretion of these terranes to inboard components of the Klamath province (Coleman et al., 1988). A possible explanation for the cause of Siskiyou contractional deformation is Middle Jurassic collision of the Insular superterrane, as proposed by Saleeby and Busby-Spera (1992) and further developed by Saleeby and Dunne (2015). Following Saleeby and Dunne (2015), our preferred tectonic model places the Insular superterrane north of the Galice Formation in Late Jurassic time, with Galice sediment derived in part from Insular superterrane sources. We note that our provenance data do not require Insular superterrane sources, and that others have proposed that the Insular superterrane was offshore at this time (e.g., Balgord et al., 2021; Tikoff et al., 2023). However, the timing of Siskiyou orogenesis agrees with the Middle Jurassic timing of collision proposed by Saleeby and Dunne (2015), who attributed the subsequent transition from contraction to extension as resulting from a late Middle Jurassic shift of the Euler pole as the Insular superterrane migrated north.

Our Galice Formation provenance results indicate that deposition of Galice Formation turbidites began during early Oxfordian time and continued well into Kimmeridgian time (ICS time scale) even as older deposits were underthrust along east-dipping faults. Furthermore, our results suggest that contraction associated with Nevadan orogenesis increased local topography, resulting in detrital zircon derived from proximal sources within the Klamath arc and possibly the nearby Insular arc, even as the basin continued to receive detritus from sources in the backarc region. Although our data do not require it, our results are consistent with Middle Jurassic accretion of the Insular superterrane and a hypothesized latest Jurassic position of the Insular superterrane just north of the Klamath Mountains province.

Our detrital zircon MDA calculations document Galice deposition as late as 150.8 Ma, well into Kimmeridgian time. We describe the provenance of Galice strata that is linked to arc and retroarc sources as the continental margin shifted from transtensional to contractional. Finally, we infer sediment transport pathways that traversed the arc from the retroarc region to deliver sediment to the forearc, even as the arc gained topographic relief during contraction associated with Nevadan orogenesis. Galice Formation provenance matches well with coeval strata in the basal Great Valley Group and the Mariposa Formation in the Western Sierran foothills, which is consistent with the development of integrated westward drainage systems into the forearc region of the Oregon–California segment of the North American continental margin during latest Jurassic time.

Funding for this project was provided by National Science Foundation Division of Earth Sciences (NSF EAR) grant no. 2052111 (principal investigator Kathleen Surpless) and NSF EAR grant no. 2052255 (principal investigator Aaron Yoshinobu and co-principal investigator Calvin Barnes), and by the University of Arizona LaserChron Center with support from NSF EAR grant no. 2050246. We thank Shane Dailey for collecting the initial two samples in this study, Glenna and Dave Atwood for their hospitality at Somes Bar, and Katie Gates and Melanie Barnes for thoughtful discussions and assistance in the field. We thank the LaserChron Center staff for support with analyses of detrital zircon and data processing and the isotope lab at the University of California, Santa Barbara, for assistance with the Hf data from plutons of the Klamath Mountains province. We are grateful to Doug Yule and Kevin Chamberlain for providing zircon separates, Joshua Schwartz for providing original Hf spot data, and Richard Silver for lab preparation and support at Trinity University.