Multi-scale, open-system magmatic and sub-solidus processes contribute to the chemical and isotopic characteristics of the Jurassic Guadalupe Igneous Complex, Sierra Nevada, California, USA

Disentangling the processes that create the chemical and isotopic characteristics of plutonic complexes (e.g., fractional crystallization, magma mixing, assimilation, sub-solidus alteration) and investigating the extent that these processes take place at the emplacement level is important for understanding continental crust formation and the volcanic-plutonic connection (e.g., Ducea et al., 2015; Bachmann and Huber, 2016). Sustained elevated temperatures supplied by basaltic, mantle-derived melts to the lower crust support the formation of lower-crustal MASH zones (i.e., zones of intensive melting, assimilation, storage, and homogenization: Hildreth and Moorbath, 1988; or deep crustal hot zones: Annen et al., 2006) as primary magma storage and differentiation sites in arcs. However, recent geophysical, petrological, and numerical modeling studies indicate that magma storage and differentiation can occur throughout the arc crust (e.g., Cashman et al., 2017; Delph et al., 2021; Barnes et al., 2021; Hildreth, 2021). Considerable debate still exists as to the extent to which differentiation processes may operate in the shallow crust (e.g., Bachmann and Huber, 2016). Cold ambient temperatures and the rarity of mafic cumulates in the shallow crust of many arcs argue against extensive differentiation in the upper crust (e.g., Keller et al., 2015). Still, large-volume eruptions equilibrated at upper-crustal pressures attest to the presence of large upper-crustal magma reservoirs (e.g., Bachmann and Huber, 2016).

Solidified shallow-crustal plutons provide opportunities for studying the extent of differentiation and magma–host rock interaction in the upper crust. These plutons reflect the time-integrated record of multiple magmatic and sub-solidus events that generate compositionally and isotopically diverse intrusive complexes. Plutonic bodies may be constructed incrementally with multiple recharging events and experience remelting (“defrosting”) of solidified units, magma mixing, fractional crystallization, and assimilation or recycling at multiple scales during their active magmatic lifespan as well as sub-solidus alteration after solidification (e.g., Harper et al., 2004; Bindeman et al., 2008; Coint et al., 2013; Barth et al., 2016; Barnes and Werts, 2022). Geochemical tools such as radiogenic and stable isotope analysis, mineral and whole-rock chemistry, and numerical modeling are used to unravel these processes (e.g., Wotzlaw et al., 2012; Fiedrich et al., 2017; Barnes et al., 2021; Payacán et al., 2023). Each tool tracks different processes at different scales; for example, zircon geochemistry records the melt composition at zircon saturation and crystallization, which can occur early or late with respect to the crystallization of other mineral phases either in the source region of a magma or at the emplacement level (e.g., Broderick et al., 2015; Samperton et al., 2015). Whole-rock geochemistry, on the other hand, integrates across a wide range of mineral compositions and thus may or may not reflect liquid compositions (Barnes et al., 2019; Cornet et al., 2022). A combination of various geochemical tools, recording the evolution of an intrusive complex at different time and length scales, allows disentanglement of the contributions of multiple processes to the observed chemical and isotopic diversity.

In this study, we combine new whole-rock and mineral isotope compositions and major and trace element mineral geochemistry to investigate the magmatic and sub-solidus processes affecting the Late Jurassic Guadalupe Igneous Complex (GIC) in the Sierra Nevada, California, USA. The GIC is a tilted, shallow-crustal, and compositionally bimodal (gabbro and granite) pluton emplaced in marine metasedimentary host rocks of similar age (Best, 1963; Haeussler and Paterson, 1993; Putirka et al., 2014). Previous work using high-precision U-Pb zircon isotope dilution thermal ionization mass spectrometry (ID-TIMS) geochronology has demonstrated that mafic and felsic magmas are coeval in the GIC (Ratschbacher et al., 2018) and syn-magmatic deformation occurred in the Bear Mountains fault zone (BMFZ) at the exposed base of the pluton (Paterson et al., 1987; Saleeby et al., 1989; Vernon et al., 1989). We aim to understand how much of the chemical and isotopic characteristics in the GIC can be attributed to (1) processes in the deeper, unexposed parts of the magmatic system; (2) in situ differentiation at the emplacement level; or (3) sub-solidus processes affecting the pluton after solidification. By combining geochemical tools capable of tracing magmatic processes on various time and length scales, we show that the chemical and isotopic characteristics of the GIC require initial magmatic differentiation below the emplacement level, during ascent, and at the emplacement level.

The GIC, located in the western foothills of the Sierra Nevada, California (Fig. 1), is one of several Jurassic arc plutons emplaced into fault-bounded, Triassic–Jurassic accreted terranes of metavolcanic and metasedimentary rocks, which collectively form the Western Metamorphic Belt (e.g., Clark, 1964; Kistler and Peterman, 1973; Tobisch et al., 1989; Snow and Scherer, 2006). The Foothills terrane, in which the GIC is located, forms the westernmost terrane of the Western Metamorphic Belt. It is bound by the Melones fault zone in the east and overlain by the Great Valley Sequence in the west (Fig. 1). The GIC intrudes metavolcanic and metasedimentary formations of oceanic affinity as well as ophiolite slivers, which experienced prehnite-pumpellyite-facies regional metamorphism (e.g., Clark, 1964; Tobisch et al., 1989; Paterson et al., 1991; Saleeby et al., 1992; Snow and Scherer, 2006). At the top of the GIC, fine-grained rhyolites are exposed (unnamed rhyolite unit in Mariposa Formation; Schweickert, 2015, unit Jmv), which yield U-Pb zircon ID-TIMS ages ~3–4 m.y. older than the GIC (Ratschbacher et al., 2018; orange unit in Figs. 1 and 2) and thus are not part of the GIC magmatic system. The GIC intrudes the Late Jurassic marine Mariposa Formation (called the Salt Spring Slate west of the BMFZ) consisting dominantly of slate and greywacke (e.g., Clark, 1964; Paterson et al., 1987; Ernst et al., 2009). The Mariposa Formation overlies the Jurassic–Triassic volcano-sedimentary Peñon Blanco Formation, which is exposed to the north of the GIC (e.g., Bogen, 1984), and the Late Jurassic Gopher Ridge island-arc volcanics exposed to the west (Clark, 1964; Tobisch et al., 1989; Schweickert, 2015). The Jurassic arc during the time of formation of the GIC is interpreted as a fringing oceanic arc built on slightly older dominantly marine sedimentary and volcanic units and local ophiolites (Schweickert and Cowan, 1975; Attia et al., 2021). The exposed base of the GIC intersects with the BMFZ, a zone of Middle Jurassic to Early Cretaceous multi-stage ductile deformation and metamorphism (Fig. 1; Paterson et al., 1987; Saleeby et al., 1989; Miller and Paterson, 1991). Metamorphic grades as high as upper amphibolite facies in the BMFZ led to partial melting (Vernon et al., 1989; Paterson et al., 1991; Lafrance and Vernon, 1993). BMFZ activity was synchronous with magmatism in the GIC and the contemporaneous Hornitos pluton to the north of the GIC (Fig. 1), indicated by outcrop-scale magmatic to solid-state transitions, sub-magmatic to high-temperature microstructures (Vernon et al., 1989; Tobisch et al., 1991), and a U-Pb TIMS multi-grain zircon age (151 ± 1 Ma; Fig. 2; Saleeby et al., 1989) of a migmatite in the BMFZ (Fig. 3A).

Best (1963), Best and Mercy (1967), and Putirka et al. (2014) showed that the GIC is compositionally bimodal and stratified, composed of fine- to medium-grained and partly layered gabbros (49–55 wt% SiO₂; Putirka et al., 2014) at the exposed base of the complex (Fig. 2). The gabbro unit grades into a compositionally heterogenous meladiorite unit (55–75 wt% SiO₂) and fine-grained granites to granophyres (67–78 wt% SiO₂; Putirka et al., 2014) at the exposed top of the complex (Fig. 2). A mingling zone separates the dominantly mafic lower part (gabbro and meladiorite units) from the upper felsic part (granites and granophyres; Fig. 2), in which mafic and felsic magmas show limited mixing (e.g., Putirka et al., 2014; “agmatite” zone of Best, 1963). Structural and paleomagnetic studies indicate an ~28° SW-side-up tilting of the GIC, exposing an ~7-km-thick plutonic section from ~10 km intrusion depth at the base to near-surface levels (Paterson et al., 1991; Haeussler and Paterson, 1993). Best and Mercy (1967) and Putirka et al. (2014) proposed an in situ differentiation model for the formation of the GIC. Using whole-rock compositional data, Putirka et al. (2014) presented a petrogenetic model involving crystal-liquid separation and mixing of differentiates. Repeated intrusion of hydrous basaltic parental magma formed the lower part of the complex, which differentiated along multiple fractionation paths to generate intermediate compositions via continuous crystal-liquid fractionation. Intermediate melts further fractionated by a discontinuous process to form mafic residues and highly evolved liquids. Local examples of such processes are felsic segregations, abundant in the upper gabbro and meladiorite units. These centimeter- to meter-sized, high-silica (68–75 wt% SiO₂; Putirka et al., 2014; Fig. 3) pods and segregations were interpreted as locally collected interstitial melt pooled from the surrounding gabbro mush (Putirka et al., 2014; Ratschbacher et al., 2018; see field description in section 4.1). Felsic differentiates were added to a growing felsic cap and mixed with underlying intermediate magmas to generate the observed compositional range in the granite and granophyre units. In situ differentiation halted when the mafic input ended.

Saleeby et al. (1989) obtained the first (bulk) zircon U-Pb ages of the GIC using the TIMS method (Fig. 1), showing that the mafic and felsic rocks have overlapping ages at 150–151 Ma. A 153 ± 2 Ma sensitive high-resolution ion microprobe (SHRIMP) zircon U-Pb age by Ernst et al. (2009) supported this result. Ratschbacher et al. (2018) presented the first high-precision chemical abrasion ID-TIMS zircon U-Pb ages, which constrained the crystallization ages of the mafic and felsic rocks between 149 and 150 Ma (Fig. 2). Calculated maximum zircon crystallization temperatures start at ~870 °C for both the mafic and felsic samples and last to temperatures of ~750 °C for the mafic and ~660 °C for the felsic rocks (see Ratschbacher et al., 2018, for details on calculations). The overall construction duration of the GIC was estimated at 295 ± 110 k.y. Two-dimensional finite-difference thermal modeling indicates that melt was present for ~200 k.y. after emplacement, resulting in total melt-present time scales of ~550 k.y. (Ratschbacher et al., 2018). Ratschbacher et al. (2018) further showed that zircons in one sample from a felsic segregation have the same ages, within uncertainty, as the zircons in the surrounding gabbro.

Saleeby et al. (1989) published seven ⁴⁰Ar/³⁹Ar biotite and hornblende cooling ages from plutons in the vicinity of the GIC, pluton aureoles, and host rocks within or near the BMFZ (Fig. 1). Ages from plutonic rocks yield 141.2 ± 0.9 Ma (biotite) and 146.6 ± 0.3 Ma (hornblende; both Santa Cruz Mountain pluton) and 140 ± 0.3 Ma (hornblende, Hornitos pluton; Fig. 1), and they are interpreted to represent cooling of these intrusive suites. Biotite ⁴⁰Ar/³⁹Ar ages of deformed host rocks were interpreted as metamorphic cooling ages and range from 144.2 ± 0.9 to 137 ± 0.5 Ma within the BMFZ (Fig. 1), suggesting that ductile deformation and metamorphism persisted until at least 137 Ma (Saleeby et al., 1989). One 114.1 ± 0.1 Ma age from a host rock outside the BMFZ and within the aureole of the adjacent Cretaceous Fine Gold Intrusive Suite (biotite; Bass Lake Tonalite; Fig. 1) reflects cooling of this intrusive suite. One multi-grain zircon U-Pb TIMS age of 123 ± 2 Ma from the undeformed White Rock pluton provides a minimum deformation age for the BMFZ (Fig. 1).

Detailed descriptions of analytical procedures for whole-rock Sr-Nd isotope analyses, whole-rock oxygen and major element geochemistry, and zircon geochemistry are in Text S1 in the Supplemental Material.¹ Brief summaries of these analytical methods are given below. Figures 1 and 2 show the locations of analyzed samples (white symbols are GIC samples; black, green, and brown samples are host rock samples).

Twenty-two (22) whole-rock samples were analyzed for Sr and Nd isotopes using a VG Sector multi-collector instrument at the University of Arizona (Tucson, Arizona, USA) following the analytical protocol detailed in Otamendi et al. (2009). The estimated standard errors for the samples analyzed are ⁸⁷Sr/⁸⁶Sr = 0.000006–0.000019 and ¹⁴³Nd/¹⁴⁴Nd = 0.000005–0.000011. Table S1 (see ^{footnote 1}) lists the data, and Text S1 describes the analytical protocol.

Twenty-nine (29) samples from the GIC were analyzed for whole-rock δ¹⁸O values: 18 from the lower gabbro unit (Shultz Mountain area; “sample transect” in Fig. 2), five from the upper gabbro unit, one from a felsic segregation at the contact of the upper gabbro and meladiorite unit, three from the meladiorite unit, and two from the granophyre unit (Fig. 2). Additionally, we analyzed one slate from the Mariposa Formation north of the GIC and three strongly deformed gneissic Mariposa Formation samples from the exposed base of the GIC within the BMFZ (Figs. 1 and 2). Whole-rock oxygen isotope analyses were carried out at Indiana University (Bloomington, Indiana, USA) using Finnigan Delta-E and Finnigan MAT 252 mass spectrometers. Results for whole-rock oxygen isotope analyses are reported in the delta (δ) notation relative to Vienna standard mean ocean water (VSMOW; Baertschi, 1976). Analytical precision is better than 0.06‰ (2σ). Table S2 (see ^{footnote 1}) lists the δ¹⁸O values in parts per mil VSMOW. Whole-rock major and selected trace elements were measured for the same samples by solution inductively coupled plasma source emission spectroscopy using a Jarrall Ash Atom Comp II at Indiana University and by instrumental neutron activation analysis (INAA) at the Phoenix Memorial Laboratory at the University of Michigan (Ann Arbor, Michigan, USA) following the procedures described in Webster and Wintsch (1987). Table S2 lists the data.

Major and minor element abundances of pyroxene, amphibole, feldspar, and titanite were analyzed using an electron microprobe at the University of California, Los Angeles (UCLA) and the U.S. Geological Survey (USGS) Menlo Park (California). Analysis conditions at UCLA and USGS have been documented in Ratschbacher et al. (2018). Table S3 (see ^{footnote 1}) lists the data.

Zircon trace-element compositions were analyzed with a laser ablation–inductively coupled plasma mass spectrometer (LA-ICP-MS) at the University of California, Santa Barbara (UCSB), following the protocol described in Kylander-Clark et al. (2013). Figure S1 (see ^{footnote 1}) shows the location of each laser spot on cathodoluminescence (CL) images of polished zircon grains. Additionally, surfaces of unpolished zircon grains from the lower and upper gabbro units have been analyzed. Text S1 describes analytical protocols and standards, and Table S4 lists results including 2σ uncertainties.

A total of 28 in situ zircon oxygen isotope analyses were performed on 23 individual grains from five samples (Fig. 2). Oxygen isotope analyses in zircon were performed at UCLA using a CAMECA IMS 1270 in multi-collection mode with a cesium cation primary beam with a spot size of 20–25 μm and a beam current of ~15 nA, following the protocol of Trail et al. (2007). Zircons were coated with gold, and raw ¹⁸O/¹⁶O values were corrected using zircon standard AS3 (δ¹⁸O = 5.34‰; Trail et al., 2007), which yielded an in-run reproducibility of 0.48‰ at the 2σ level. Values of δ¹⁸O are reported in the delta (δ) notation relative to Vienna standard mean ocean water (VSMOW; Baertschi, 1976) and are listed in Table S5 (see ^{footnote 1}).

Early-crystallizing mineral–whole-rock equilibrium tests and mineral-melt chemometric calculations were applied to mineral analyses to test whether whole-rock compositions represent frozen melts or whether they represent cumulates that have lost melt (e.g., Barnes et al., 2016; Barnes and Werts, 2022; Cornet et al., 2022). The Fe-Mg exchange coefficient [^Fe/MgK_D = (Fe/Mg)_mineral/(Fe/Mg)_rock] between clinopyroxene and host rocks in equilibrium is 0.28 ± 0.08, and ^Fe/MgK_D = 0.29 ± 0.06 between orthopyroxene and host rocks (Putirka, 2008). The anorthite-albite (An-Ab) equilibrium test of Putirka (2008) was similarly applied to plagioclase (pl) analyses (; mineral and liquids components, X, are calculated as cation fractions). An ^An-AbK_D value of 0.1 ± 0.05 indicates equilibrium conditions at temperatures <1050 °C, and an ^An-AbK_D value of 0.27 ± 0.11 is suitable for temperatures >1050 °C.

Chemometric approaches leverage experimental data sets to quantify the dependence of a mineral composition on the temperature (± pressure) and composition of the melt it crystallized from (e.g., Ridolfi and Renzulli, 2012; Putirka, 2016; Zhang et al., 2017; Scruggs and Putirka, 2018). Higgins et al. (2022) used machine learning techniques to calibrate these relationships for arc magmas in the Lesser Antilles and predict the pressure, temperature, and anhydrous melt composition (SiO₂, Al₂O₃, CaO, Na₂O, K₂O, MgO, FeO, TiO₂) for a given clinopyroxene or amphibole analysis. Standard error estimate values for the amphibole and clinopyroxene equilibrium melt models range from 0.2 to 4.7 wt% (element and mineral dependent; see Higgins et al., 2022). The equations 2 and 3 of Scruggs and Putirka (2018) were applied to plagioclase analyses to calculate equilibrium melt SiO₂ and CaO contents. SiO₂ and CaO are estimated to ±3.9 wt% and ±1.2 wt%, respectively (Scruggs and Putirka, 2018).

We use multiple mineral thermometers to constrain crystallization temperatures of mafic and felsic rocks in the GIC. Application of mineral thermometers may rely on (1) knowledge of the melt composition in equilibrium with the mineral composition (e.g., clinopyroxene-liquid thermometer [Putirka et al., 2003]; plagioclase-liquid thermometer [Putirka, 2005]). A melt composition could be either a whole-rock, glass, melt inclusion, or a calculated equilibrium melt composition; (2) pressure (e.g., Putirka, 2016, equation 6); or (3) mineral compositions alone (e.g., Putirka, 2016, equation 5; Higgins et al., 2022). Where thermometers required a melt input, the whole-rock composition was used only if it passed the mineral-whole rock equilibrium tests described in the previous section. Given the low number of analyses that fulfill this requirement (see section 4), the calculated equilibrium melt composition was used if available (e.g., clinopyroxene and amphibole; Higgins et al., 2022). The following thermometers were applied to the data set:

• Clinopyroxene (cpx): Higgins et al. (2022) (cpx-only); Putirka et al. (1996, 2003) (cpx-liquid; whole-rock composition was used as a liquid composition if in equilibrium with cpx; valid for only one gabbro from the upper gabbro unit);

• Orthopyroxene (opx): Kohler and Brey (1990) (opx-only; assuming a pressure of 2 kbar);

• Plagioclase (pl): Putirka (2005, equation 2; pl-liquid; whole-rock composition was used as a liquid composition if in equilibrium with pl; water contents of 3.5 wt% in the lower and upper gabbro and 4 wt% in the mingling zone and granophyre units); rhyolite-MELTS (Gualda et al., 2012) modelling was used to determine the onset of plagioclase crystallization in each unit; the lowest-SiO₂ equilibrium melt calculated from pyroxene and amphibole compositions (Higgins et al., 2022) was used as the starting melt composition;

• Amphibole (amp): Putirka (2016, equation 5; amp-only; pressure-independent); Higgins et al. (2022; amp-only);

• Titanite (ttn): Hayden et al. (2008, equation 7; ttn-only; activities of TiO₂ and SiO₂ of 1);

• Zircon (zrn): Ferry and Watson (2007, equation 15; zrn-only; TiO₂ activity of 0.5 (presence of ilmenite), and SiO₂ activity of 1 for quartz-bearing and 0.8 for quartz-free samples; zrn crystallization temperatures were also estimated by calculating zircon saturation distributions; see Ratschbacher et al., 2018, for details).

Varying mineral proportions define modal layering in the gabbro, commonly by the presence or absence of amphibole and/or varying grain size. Compositional layering is more pronounced in the lower gabbro than the upper gabbro, but monomineralic rocks rarely occur. Subtle layering between amphibole- and pyroxene-dominated gabbro occurs in the upper gabbro (Fig. 3B). Igneous layering dips to the east, and dip magnitude ranges from sub-vertical dips to 40° (Best, 1963) toward the top of the complex. Numerous centimeter- to meter-sized felsic segregations occur in the upper gabbro, exhibiting elongated tear shapes with high aspect ratios (usually >5:1; Figs. 3C–3E). On average, their long axes point to the top of the complex, commonly at a high angle to the sub-vertical layering in the gabbros. The contact between the felsic segregation and the surrounding gabbro is gradational in small segregations (Fig. 3C) but sharp in large ones and characterized by vesiculation (Figs. 3D and 3E). The tips of segregations are round and/or irregular (Fig. 3C and 3E). The segregations have a larger grain size (as large as several centimeters) than the surrounding gabbros and abundant hydrous mineral phases (amphibole, biotite). The meladiorite unit is compositionally heterogenous (gabbro to granite), but limited exposure does not allow the estimation of proportions of mafic and felsic domains. Like in the upper gabbro unit, felsic segregations occur, but they tend to be larger in size (as large as multiple tens of meters). The mingling zone separating dominantly mafic and felsic units shows evidence for mutual intrusions of mafic and felsic magmas. Field observations indicate emplacement of sub-parallel gabbroic dikes into granitic hosts (Putirka et al., 2014, their figure 8A) and net veining of felsic magma breaking up mafic blocks (Fig. 3F). This indicates that both felsic and mafic magmas were present during the formation of the mingling zone. Although exposure is limited, a gradual disappearance of mafic dikes and enclaves marks the transition into the felsic units (granite and granophyre) at the top of the GIC. The felsic units were pervasively altered in parts as evidenced by red staining on surfaces and along fractures (Fig. 3G).

Petrographic descriptions are based on thin section observations of samples collected in this study and observations from Best (1963), Best and Mercy (1967), Putirka et al. (2014), and Canchola (2016). Petrographic observations from each unit are summarized in Table 1 and described in detail in Text S1 (see ^{footnote 1}). Representative thin section images are shown in Figure 4.

Six samples were collected for textural investigations and geochemistry of zircon (Fig. 2): (1) amphibole (meta)hornblendite (sample BRGIC13) from the lower gabbro unit within the BMFZ; (2) amphibole-pyroxene gabbro (BRGIC10) from the upper gabbro unit, which was collected ~1 m away from (3) felsic segregation sample 10D, collected from the upper gabbro unit; (4) pyroxene gabbro (BRGIC3) from the meladiorite unit; (5) granite (BRGIC2) from the granite unit; and (6) a rhyolite (BRGIC5; no zircon oxygen isotope data). These samples are described in detail in Ratschbacher et al. (2018) and have been dated with zircon U-Pb CA-ID-TIMS geochronology. Zircon textures in the six investigated samples are simple, lacking obvious inherited cores and showing either no zoning or one zoning type (Fig. S1 [see ^{footnote 1}]). The cathodoluminescence (CL) images show concentric, sector, and patchy zonation, defined as follows: (1) Concentric zoning exhibits repeated, oscillatory and parallel zones of various thicknesses; (2) sector zoning shows bright and dark zones separated by straight boundaries; and (3) patchy zoning has irregularly shaped and chaotic bright and dark domains, which do not show a repeated pattern.

Mg# [molar Mg/(Mg + Fe_total)] of clinopyroxene ranges from 0.54 to 0.80 and is higher in amphibole gabbros (0.73–0.80; n = 3) and lower in pyroxene-gabbros (0.54–0.65; n = 2; Fig. S2 [see ^{footnote 1}]). Compositional zonation is absent in clinopyroxenes. Calculated Fe-Mg exchange coefficients (K_D) vary between 0.34 and 0.75, above the equilibrium value of 0.28 ± 0.08 (Putirka, 2008; Fig. S2). Equilibrium melts calculated following Higgins et al. (2022) yielded melt compositions of 57–68 wt% SiO₂ for amphibole and pyroxene gabbros from the gabbro and meladiorite units (Table 2).

Mg# ranges from 0.4 to 0.66 in orthopyroxene. Amphibole gabbros have higher Mg# (0.61–0.66; n = 2) than pyroxene gabbros (0.41–0.50; n = 2; Fig. S2). Calculated Fe-Mg exchange coefficients (K_D) are higher than equilibrium K_D values, ranging between 0.56 and 1.22 (Fig. S2). Several transects through individual grains show no zonation in orthopyroxene.

Primary magmatic green amphiboles are magnesio- and ferri-hornblendes. Secondary brown, magmatic amphiboles replacing pyroxenes are dominantly tschermakites. Actinolite occurs in samples that show widespread alteration and are interpreted as deuteric (e.g., Werts et al., 2020). Amphibole Mg# ranges from 0.22 to 0.78 (Fig. S2). Amphiboles from the lower gabbro sample deformed in the BMFZ (sample BRGIC13) yield lower Mg# (0.58–0.66) compared to upper gabbro samples (0.61–0.78). Amphiboles from an upper gabbro adjacent to a felsic segregation have Mg# between that of the felsic segregation and upper gabbro amphiboles, corresponding to the transition from green primary magmatic amphiboles (segregation) to brown secondary magmatic amphiboles (gabbro) observed in thin section. Amphibole grains from the felsic hosts in the mingling zone have lower Mg# (0.22–0.63) compared to those of mafic enclaves (0.51–0.69; Fig. S2). Most amphiboles from the granite and granophyre samples are actinolitic. Secondary magmatic amphiboles, replacing clinopyroxenes, from the lower and upper gabbro units yielded equilibrium melt SiO₂ contents between 62 and 77 wt%, using the Higgins et al. (2022) chemometer (Table 2). Felsic segregations containing primary magmatic amphiboles yielded melt compositions between 70 and 78 wt% SiO₂. The meladiorite unit yielded amphibole melt compositions between 73 and 78 wt% SiO₂. Upper GIC units, including the mingling zone, granite, and granophyre, produced calculated amphibole melt compositions with SiO₂ contents between 72 and 78 wt% (Table 2).

Plagioclase is predominantly found in the gabbro, felsic segregations, and mingling zone. Alkali feldspar characterizes the upper felsic units and felsic segregations. Anorthite contents progressively decrease up-section through the GIC, from An₇₃ (anorthite end member component in plagioclase) in the upper gabbro unit, to An₆₃ in the felsic segregations within the upper gabbro unit, to An₃₂ in the meladiorite unit (Fig. S2). The mingling-zone plagioclase compositions are distinct, with the mafic enclaves containing cores and rims with An₇₂–An₁₆ and the felsic host containing An₃₆–An₁₂ (Fig. S2). Individual grains span a wide range in anorthite content, indicating compositional zoning within all units. Although rims generally have lower anorthite contents than cores, no simple core-rim zoning pattern occurs within most grains. All plagioclase data from the meladiorite and mafic mingling zone plot well above the lower-temperature An-Ab equilibrium boundaries, at K_D >0.3. Some felsic segregation plagioclases and felsic mingling zone plagioclases plot within the equilibrium zone. Equilibrium melt calculations (Scruggs and Putirka, 2018, equations 2 and 3) estimate SiO₂ melt compositions between 58 and >80 wt% (Table 2) and CaO contents <6 wt%.

Crystallization temperatures were calculated for pyroxene, amphibole, plagioclase, and titanite and compared to zircon crystallization temperatures reported in Ratschbacher et al. (2018). Results are shown in Figure 5 and reported in Table 2.

Clinopyroxene crystallization temperatures agree across methods in the upper gabbro unit, clustering at 1000 °C. The meladiorite unit shows clinopyroxene crystallization temperatures clustering at 960 °C. The Ca-in-orthopyroxene thermometer yielded crystallization temperatures >900 °C in both the upper gabbro and meladiorite units, corresponding to petrographic observations that orthopyroxene is an early crystallization phase.

The amphibole-only thermometer of Putirka (2016) yielded crystallization temperatures ranging from 950 to 750 °C in lower mafic units and 780–620 °C temperatures in intermediate and felsic upper units. These temperatures show general agreement with estimates using the machine learning calibration method of Higgins et al. (2022) for lower mafic and intermediate GIC units, but the Higgins et al. (2022) thermometer predicts significantly higher temperatures (>750 °C) than Putirka (2016) for the upper GIC units.

Calculations using Putirka (2005, equation 2) produced plagioclase crystallization temperatures that decrease upwards in the GIC but at unusually high temperatures relative to other minerals in the assemblage (1100–1200 °C). These estimates were compared to rhyolite-MELTS models. In the upper gabbro, models produced An_73–67 plagioclase at 1020–990 °C, which corresponds to the highest anorthite content found in the upper gabbro samples, and the model recreated the approximate (high-temperature, near liquidus) crystallization sequence observed in thin sections. Modeled plagioclase crystallization temperatures in the felsic segregations are slightly lower at 940–910 °C with good agreement between clinopyroxene- and amphibole-liquid models. The meladiorite results are comparable to those from the upper gabbro, with plagioclase crystallization temperatures modeled at 950–990 °C, producing An_68–66 plagioclase. This is more calcic than found in the meladiorite samples, but thin section observations suggest plagioclase cores have experienced alteration. The mingling zone, granite, and granophyre units produced similar plagioclase crystallization temperatures in rhyolite-MELTS, between 870 and 820 °C. The models in this case produced quartz at high temperatures and could not reproduce high-temperature phases. This may be due to the reliance on amphibole-equilibrium melts (no clinopyroxene found) and the poor reproducibility of amphibole in MELTS models. Anorthite contents predicted by the model also showed some mismatch compared to natural samples. For example, the mingling zone model produced An₅₀ plagioclase, but the most calcic plagioclase measured is An₃₆. The model could not reproduce the calcic plagioclase grains (An₇₁) in the mafic mingling zone end member, and we suspect this is due to inheritance. Rhyolite-MELTS models were also used to estimate liquidus temperatures for each unit, showing a general decrease in temperature between the lower and the upper units.

Titanite crystallization temperatures yielded temperatures from 780 °C to the solidus in the upper gabbro, mingling zone, and granite units.

The Ti-in-zircon thermometer of Ferry and Watson (2007) was used to calculate zircon crystallization temperatures in addition to the calculation of zircon saturation distributions (see Ratschbacher et al., 2018, for details). Temperatures cluster across an interquartile range of 30–70 °C but span a much larger range from 890 to 650 °C. Figure 5 illustrates the relative shift in crystallization temperature between the major and accessory minerals studied. Initially, zircon crystallized at lower temperatures (i.e., later) relative to early plagioclase, pyroxene, and amphibole (where available) within the gabbro units and the felsic segregations. In contrast, zircon is an early phase relative to amphibole and plagioclase in the meladiorite and granite units.

Whole-rock initial ε_{Nd(149 Ma)} (ε_Ndi; where the subscript _{(149 Ma)} indicates age-corrected values) and initial ⁸⁷Sr/⁸⁶Sr _{(149 Ma)} (Sr_i; where the subscript _{(149 Ma)} indicates age-corrected values) values of 13 samples from the GIC vary between +7.02 and +2.13 and between 0.7026 and 0.7069, respectively (Fig. 6). There is no systematic variation in ε_Ndi and Sri_i values from the base to the top of the GIC. The upper gabbro and felsic segregation samples (BRGIC10 and 10D), collected in close proximity, differ in ε_Ndi values by 4.90 and in Sri_i values by 0.001 (“Adjacent samples” in Fig. 6). All samples except two samples from the granophyre unit show values of Sri_i <0.706, consistent with the GIC’s location west of the ⁸⁷Sr/⁸⁶Sr 0.706 line (Kistler and Peterman, 1973; Kistler, 1990). Three samples from the granophyre unit show a trend toward increasing Sri_i with slightly decreasing ε_Ndi values similar to trends observed in mid-ocean-ridge basalt (MORB) affected by seawater alteration (Hart et al., 1974; Fig. 6). One 152 Ma rhyolite sample collected directly east of the granophyre unit (Fig. 2) has a ε_{Nd(152 Ma)} value of –11.33 and a Sri_{(152 Ma)} value of 0.7078 (Fig. 6). Five samples collected from several units of the Foothills terrane and Calaveras Complex analyzed in this study (Fig. 1) yielded ε_{Nd(152 Ma)} values from +9.63 to +3.55 and Sri_{(152 Ma)} values from 0.7029 to 0.7047 (Fig. 6). For comparison, Figure 6 shows data from the Jurassic Hornitos pluton (this study) and published data from Cretaceous plutonic rocks (pink outline), Cretaceous mafic ring dike complexes (black outline), two samples from the Calaveras Formation (brown squares), and one sample from a volcanic rock west of the Melones fault (light green square; see figure captions for references).

Whole-rock δ¹⁸O [δ¹⁸O_(whole-rock)] values in the GIC vary from +5.47‰ to +8.57‰ (Figs. 7 and 8) and show an increase from the lower gabbro to the meladiorite unit (Fig. 8). Granophyre δ¹⁸O_(WR) values depart from the generally upward increase in δ¹⁸O with increasing stratigraphic position, having values that are lower (+5.74‰ and +7.01‰) than those of the units immediately below (Fig. 8).

A detailed transect of 21 samples, collected perpendicular to compositional layering in the lower to upper gabbros (Shultz Mountain area; Fig. 7A; “sample transect” in Fig. 2), shows a variation in δ¹⁸O_(WR) values from +5.47‰ to +7.24‰ (Fig. 7B). Although scatter is observed, the lowest δ¹⁸O_(WR) values are from the base of the transect and increase by ~2‰ up stratigraphy (Fig. 7B). This correlates with a decrease in whole-rock Mg# values from ~75 to 65 in the lower gabbros and an abrupt shift to lower values in the upper gabbros (Fig. 7C).

Mariposa Formation host rock δ¹⁸O_(WR) values at the base of the GIC in the BMFZ are low (+4.31‰ to +5.95‰) compared to the host rock sample collected ~1 km north of the GIC (+12.12‰; Figs. 2 and 8). The δ¹⁸O_(WR) values of this Mariposa Formation sample is similar to δ¹⁸O_(WR) averages of ~13‰ in the Mariposa Formation from farther north (Boehlke and Kistler, 1986).

Zircons from the lower and upper gabbro units show overlapping oxygen isotope compositions (+3.8‰–+4.6‰), below the value for mantle-derived zircons (+5.3‰ ± +0.3‰; Valley et al., 1998; Fig. 9). Zircons from the granite and meladiorite samples show values above the mantle value, between 5.6‰ and 6.6‰, while values from the felsic segregation zircons (+4.8‰–+5.4‰) mostly overlap with the mantle values (Fig. 9). Zircon δ¹⁸O trends upward with stratigraphic position: At the base of the GIC (gabbros), the zircons show sub-mantle isotope compositions, whereas at the top (granite), the zircons record elevated oxygen isotope compositions. This upward increase in zircon δ¹⁸O values mirrors the upward increase in whole-rock δ¹⁸O values (Fig. 8). Zircon grains show an intra-sample range in δ¹⁸O of as much as 0.5‰ (Fig. 9).

A comparison to zircon δ¹⁸O values from Late Jurassic gabbros of the Eastern Sierra Nevada (Gevedon, 2013; Gevedon and Clemens-Knott, 2013; green arrow in Fig. 8) show similarly low zircon δ¹⁸O values to GIC gabbros. A compilation of δ¹⁸O_(zircon) values of Cretaceous plutons in the Western Sierra Nevada (Lackey et al., 2008; light pink arrow in Fig. 8) reveals that the GIC gabbros have lower zircon δ¹⁸O values. In addition, GIC δ¹⁸O_(zircon) values are lower than typical mantle values obtained in Cretaceous gabbros in the Western Sierra Nevada (Fu et al., 2008; Lackey et al., 2008, 2012).

Zircon trace element variations occur at the grain, hand sample, and pluton scales. At the grain scale, zircon rims from the lower gabbro and meladiorite units have higher Yb/Dy but lower Th/U ratios and show higher Hf concentrations compared to grain interiors (Figs. 10A–10B and 10G–10H). Rim compositions of unpolished zircons from the lower and upper gabbro correspond to rims on polished grains (Fig. 10). Zircons from the gabbro of the meladiorite unit show a positive correlation between Hf and Yb/Dy and a negative correlation between Hf and Th/U and between Hf and Ti (Figs. 10G–10I). Zircons from the upper and lower gabbro show a weak positive correlation between Hf and Yb/Dy (Figs. 10A and 10D). Zircons from the felsic segregation (Figs. 10J–10L) do not show systematic core-rim compositional relationships, while zircons for the granite sample show slightly higher Yb/Dy in the rims compared to cores (Fig. 10M). The intra-grain compositional range observed in the mafic unit and felsic segregation samples is larger than in the zircons from the granite. Figure S3 (see ^{footnote 1}) shows the difference in the minimum and maximum values of Hf concentrations and Yb/Dy within single zircons analyzed from each sample. These variations in intra-grain composition do not correlate with a zircon CL-specific zonation type.

Although compositional overlap exists at the hand-sample scale, the zircons from the upper gabbro unit can be distinguished from the zircons in the lower gabbro unit by slightly higher Ti contents and Lu/Hf and lower Hf concentrations and Yb/Dy (Figs. 11A–11C). Zircons from the felsic segregation sample have higher Hf concentrations and Yb/Dy and lower Ti contents than the zircons of the adjacent upper gabbro sample (Figs. 11A–11C). Zircons from the gabbro of the meladiorite unit are dissimilar to the other gabbro zircons in having generally higher Yb/Dy and lower Th/U (Figs. 10 and 11). Zircons from the granite sample have the lowest Eu/Eu* ([Eu_N/(Sm_N × Gd_N)^0.5]; N indicates normalization to chondritic values; Fig. 11D); their Hf content is similar to that of the rhyolite zircons (Fig. 11). The latter have higher Eu/Eu* (Fig. 11D), but the Ti contents and element ratios overlap with those of the granite unit zircons (Figs. 11A–11C). The zircons of the felsic segregation sample have distinctively higher Yb/Dy than the majority of zircons from the other samples (Fig. 11C). Figure 11 also shows that the mafic zircons and the zircons from the felsic segregation sample span a wider compositional range than the zircons from the granite and rhyolite samples.

On the pluton scale, the Eu anomalies (Eu/Eu*) are increasingly more pronounced from the lower gabbro to the granite unit (Fig. 11D); this corresponds with decreasing ages and correlates with the whole-rock composition trend: More mafic and older samples have higher Eu/Eu* compared to younger and more felsic samples.

To better understand the observed geochemical range in zircon trace elements, we conducted coupled major and trace element modeling integrated with accessory phase saturation calculations using Perple_X (Connolly 2005, 2009) and StatGeochem.jl (Keller, 2021). This allows estimation of the composition of zircon expected to crystallize from the interstitial melts resulting from equilibrium bulk crystallization of each whole-rock composition. As demonstrated in previous sections, the majority of early-crystallizing mineral phases are not in equilibrium with their whole-rock composition, thus whole-rock compositions cannot be equated to melt compositions. However, by using respective whole-rock compositions, we can predict what the zircon compositions would be given closed-system equilibrium crystallization and compare that with observed compositions. We first conducted equilibrium batch crystallization calculations from liquidus to solidus using the thermodynamic modeling program Perple_X (Connolly 2005, 2009) and the following melt and solid solution models: melt (Holland et al., 2018); clinoamphibole and clinopyroxene (Green et al., 2016); olivine, orthopyroxene, spinel, and garnet (Jennings and Holland, 2015); and mica (White et al., 2014), with the internally consistent end-member thermodynamic database of Holland and Powell (2011) at a pressure of 3 kbar and 4 wt% H₂O. Using the resulting modal abundances of each solid phase and melt, we then conducted trace element partitioning calculations using the composition-dependent partition coefficient data set of Keller (2021), compiled via the Geochemical Earth Reference Model (GERM) Partition Coefficients Database (https://kdd.earthref.org/KdD/). With these melt trace element compositions and temperatures, we then applied the zircon saturation model of Boehnke et al. (2013), the apatite saturation model of Harrison and Watson (1984), the sphene model of Ayers et al. (2018), and the monazite saturation model of Montel (1993). Because the saturation of several of these accessory phases depends on trace elements that also partition into other accessory phases, bulk solid-melt partition coefficients and melt trace element compositions were updated iteratively after adding each accessory phase. Trace element concentrations in model zircon were calculated using the final melt trace element compositions and temperatures and the temperature-dependent partition coefficients of Claiborne et al. (2018). Finally, Ti in model zircon was calculated using the Ti-in-zircon relation of Ferry and Watson (2007) with melt Ti activities checked against the rutile saturation model of Hayden and Watson (2007). Results are shown as modeled paths in Figure 10.

In this study, we present a combination of (1) textural observations, (2) mineral-whole rock equilibrium calculations, (3) mineral chemometric calculations, (4) zircon trace element geochemistry and modeling, and (5) pluton-scale Sr-Nd, δ¹⁸O_(WR), and δ¹⁸O_(zircon) isotope data that together indicate that open-system magmatic and sub-solidus processes operated at all scales to create the observed chemical and isotopic characteristics of the GIC. In the following discussion, we provide evidence for this interpretation by first discussing the isotopic and chemical signatures that were acquired below the emplacement level and then the signatures acquired at the emplacement level from the outcrop scale to the pluton scale. Finally, we integrate the outcomes of this study with published geochronology to present a new conceptual model of the evolution of the GIC.

Mineral chemometric calculations of early-crystallizing phases (e.g., pyroxene) indicate that melts arriving in the GIC were at least andesitic in composition (≥58 wt% SiO₂; Table 2) and did not represent direct mantle-derived primary magmas. An additional differentiation step would have been required below the emplacement level to generate intermediate melts of ≥58 wt% SiO₂.

GIC Sr-Nd isotope compositions are regionally consistent with their spatial position in the western foothills of the Sierra Nevada arc. Plutons and volcanic rocks in this region, west of the Foothills suture and the ⁸⁷Sr/⁸⁶Sr_i 0.706 line, display a juvenile, largely mantle-derived isotopic signature, modified by intrusion into a pre-existing arc basement of oceanic and distal continental affinity sediments (Kistler and Peterman, 1973) that include hydrothermally altered materials (Lackey et al., 2008). A comparison to other mafic complexes in the foothills, such as the Cretaceous mafic ring dike complexes (Clemens-Knott, 1992) further south, shows that the GIC Sr-Nd isotope compositions are similar (Fig. 6). However, the GIC shows slightly lower δ¹⁸O_(WR) values at low-SiO₂ (gabbro) compositions compared to the Cretaceous mafic ring dike complexes and lower δ¹⁸O_(WR) values compared to other Cretaceous plutonic rocks (δ¹⁸O_(WR) panel in Fig. 8). Additionally, GIC gabbro δ¹⁸O_(zircon) values are lower than typical mantle values and δ¹⁸O_(zircon) values obtained from Cretaceous gabbros in the Western Sierra Nevada (Fu et al., 2008; Lackey et al., 2008, 2012; δ¹⁸O_(zircon) panel in Fig. 8). The low δ¹⁸O values observed in whole rocks and zircons from GIC gabbros can be explained by either (1) an enriched source with unusually low δ¹⁸O values present in the Jurassic but not in the Cretaceous Sierra Nevada arc, and/or (2) assimilation of high-temperature, hydrothermally altered material (either bulk assimilation or assimilation of partial melts of altered material). Late Jurassic to Early Cretaceous gabbros from the Eastern Sierra Nevada show similarly low ¹⁸O_(zircon) values (Gevedon, 2013; Gevedon and Clemens-Knott, 2013; Fig. 8). Thus, mafic magmas in the Late Jurassic arc might have re-worked earlier Jurassic marine sediments that had been hydrothermally altered along arc-wide shear zones, such as the BMFZ. Studies specifically targeting Late Jurassic mafic rocks are required to further test the possible presence of a low-δ¹⁸O source and the extent of alteration in the Late Jurassic arc. Assimilation of pre-existing hydrothermally alerted crust is commonly invoked to explain low-δ¹⁸O magmas (e.g., Bindeman et al., 2008; Troch et al., 2020). Bulk assimilation of host rock units is not observed at the emplacement level (e.g., there is an absence of stoped blocks and host rock xenoliths), thus if present, such assimilation must have occurred below the emplacement level and/or during ascent. Contamination of rocks from the Calaveras Complex, east of the Melones fault zone, is excluded due to their distinctively more radiogenic isotope composition (i.e., Calaveras Formation; Fig. 6) and their fault-separated location in another terrane (Fig. 1). The Mariposa Formation wall rocks have a similar Sr_i isotopic composition (0.7047 and 0.7046, respectively; Lackey et al., 2012) to the GIC magmatic Sr_i values (0.7026–0.7046) but higher δ¹⁸O_(WR) values (Fig. 8). However, at the base of the GIC in the BMFZ, Mariposa Formation δ¹⁸O_(WR) values are as low as +4.3‰ (Fig. 8), indicating the influence of high-temperature exchange with meteoric water and/or seawater. Given the proximity of the BMFZ to the gabbro unit (Fig. 1) and the presence of migmatites in the BMFZ at the base of the GIC (Fig. 3A), bulk assimilation or assimilation of partial melts of hydrothermally altered Mariposa Formation might account for the low δ¹⁸O_(WR) and δ¹⁸O_(zircon) of GIC gabbros.

Equilibrium partitioning of trace elements between zircon and melt is a function of temperature and melt composition (e.g., Claiborne et al., 2018). Because the start of zircon crystallization in mafic and felsic samples occurred at similar temperatures, ~880 and ~865 °C, respectively (Fig. 5; Ratschbacher et al., 2018), we interpret the differences in trace element composition of zircon from the mafic and felsic samples as a function of changing melt composition and not temperature. Zircon Hf concentration can be used as an indicator of magmatic differentiation and melt evolution (Hoskin and Schaltegger, 2003; Claiborne et al., 2006). Zircons from the lower gabbro and meladiorite units show a wide range in Hf concentrations that are generally greater than those of the granite and rhyolite samples (Fig. 11). This is counterintuitive, given that zircon from felsic rocks should record more evolved and thus higher Hf concentrations (e.g., Hoskin and Schaltegger, 2003). This can be reconciled if zircon in the mafic samples crystallized later relative to the major phases (e.g., pyroxene, feldspar, amphibole; Fig. 5). None of these major phases incorporate significant amounts of Hf, thus the high Hf concentrations in the mafic samples can be explained by late-stage zircon saturation and growth in interstitial melts that were highly enriched in Hf. Ratschbacher et al. (2018) determined that >80% bulk crystallization occurred prior to zircon saturation in the mafic samples (as opposed to ~10% bulk crystallization in felsic samples). Thus, at 80% crystallinity, zircon in the mafic samples likely crystallized from residual, interstitial melt pools recording the latest stages of crystallization.

Zircon from mafic samples that crystallized from late-stage melts should display a small compositional range, but the opposite is observed (Figs. 10 and 11): For example, zircons from the granite sample show a much narrower compositional range in Hf concentrations and Eu/Eu* than the zircons from mafic samples (e.g., Fig. 11D). Perple_X modeling in combination with zircon trace element partitioning calculations shows that closed-system crystallization in mafic samples, taking the whole-rock composition as a parental liquid approximation (although this is not a good approximation—see section 5.2.3. below), produces a much narrower range in zircon trace element compositions than is observed in the data (Fig. 10). Consequently, we infer that the observed large range in zircon trace element compositions is not the result of a single differentiation path but rather the consequence of zircon crystallization from diverse interstitial melts, collected from an area within the gabbro crystal mush larger than the felsic segregation–upper gabbro sample distance (~1 m). This implies that the gabbro crystal mush was constructed from multiple magma batches with slightly different melt compositions and degrees of differentiation and that the interstitial melt was probably able to percolate through the mush, seen at the mesoscale through felsic segregations.

Previous work using whole-rock chemistry (Putirka et al., 2014) and high-precision U-Pb zircon geochronology (Ratschbacher et al., 2018) indicates that the felsic segregations in the upper gabbro unit formed through fractional crystallization processes from the surrounding gabbro crystal mush. However, the felsic segregations and surrounding upper gabbro have distinct whole-rock Sr-Nd isotopic compositions and δ¹⁸O_(zircon) values (Figs. 6 and 9), precluding closed-system fractionation of a single melt from the gabbro crystal mush. Models to explain the whole-rock Sr_i and ε_Ndi composition of the gabbro and adjacent felsic segregation include: (1) fractionation of an isotopically hybrid melt from an isotopically heterogeneous crystal mush (un-mixing trend; Beard, 2008); (2) selective accumulation of mineral phases that record different magmatic histories and isotopic ratios (e.g., Tepley and Davidson, 2003, their scenario 2), including isotopic zonation within phases (e.g., Tepley and Davidson, 2003, their scenario 3); (3) sub-solidus metasomatism or re-equilibration of isotopic composition in either the felsic segregation or the adjacent gabbros; or (4) localized and variable assimilation of crust. The samples of the upper gabbro and felsic segregation units are not located near any external contacts (Fig. 2), and field and petrographic observations indicating limited to absent sub-solidus alteration (Figs. 3 and 4) suggest that metasomatism (hypothesis 3 above) and crustal assimilation (hypothesis 4) are not important in driving the isotopic differences. Additionally, neither the gabbro nor the felsic segregation samples plot within the Sr_i seawater-alteration trend (Fig. 6; Hart et al., 1974). All these observations suggest that the observed isotopic differences record magmatic histories.

The higher Yb/Dy and Hf concentrations and lower Ti contents and Eu/Eu* (Fig. 11) in the felsic segregation zircons support the interpretation that the felsic segregation melt was more differentiated than the melt from which the upper gabbro zircons crystallized. Despite both rock types being sampled in close vicinity, the isotopic discrepancy supports the conclusion that the felsic segregation melt was not simply the interstitial melt from the immediately surrounding gabbro. We thus infer that hypotheses 1 and 2 above more likely explain the difference in isotopic composition between the gabbro and the felsic segregation samples: Pre-existing differences in the crystals and melts entering the mush system led to whole-rock isotopic heterogeneity, which was further amplified during late-stage melt extraction from the mush to form felsic segregations. To discriminate between hypotheses 1 and 2, in situ mineral isotopic data should be collected from the major phases. This would also be effective at evaluating the extent (if any) of hydrothermal alteration on Sr isotopes (Tepley and Davidson, 2003).

Several physical mechanisms have been proposed to drive the separation of crystals and melt in a mush at high crystallinity, with compaction involving crystal repacking and/or crystal deformation considered the most efficient (Bachmann and Bergantz, 2006; Holness et al., 2017; Bachmann and Huber, 2019). Elongated plagioclase crystals in the upper gabbro, which define a planar foliation, show deformation twinning and undulose extinction (Figs. 4A–4C), indicating crystal-plastic deformation of a crystal framework during igneous compaction rather than gravitational settling in a melt-rich environment (e.g., Holness, 2018). Field observations also indicate the importance of a fluid phase in the formation of the felsic segregations, which may aid in crystal-melt segregation processes (Sisson and Bacon, 1999): (1) the presence of vesicular textures around many segregations (Figs. 3D and 3E), and (2) large grain sizes in the segregations compared to the gabbros, indicating that fluids likely aided crystal growth. We thus propose that interstitial melts in the gabbros were able to percolate through the crystal mush at high crystallinity and collect to form felsic segregations by a combination of crystal-mush compaction processes and gas-driven filter pressing (e.g., Sisson and Bacon, 1999). The tear shapes of the segregations with rounded to irregular tips indicate that ductile rather than brittle fracturing operated in the crystal mush, which is consistent with the observed visco-plastic behavior of a crystal mush upon intrusion of melt (Bergantz et al., 2015). The structures described here resemble melt segregation pools from granitic plutons described by Weinberg (2006).

The gabbro unit is compositionally and texturally heterogenous, as evident by different modal abundances of olivine, pyroxene, amphibole, and fine-grained and coarse-grained varieties of these minerals (Best, 1963; Putirka et al., 2014; Table 1). For example, field observations in the upper gabbro unit show centimeter- to meter-scale layering of amphibole- and pyroxene-rich layers (Fig. 3B). The lack of a single, well-defined differentiation trend in the whole-rock compositions from the gabbro and meladiorite unit samples (TiO₂, Al₂O₃, and P₂O₅ versus SiO₂; Putirka et al., 2014, their figure 7) indicates either (1) different degrees of crystal accumulation or melt loss; (2) emplacement of multiple compositionally distinct magma batches; or (3) a combination of both. The previous discussion presented evidence for both processes. Mineral chemistry further supports crystal accumulation processes: Most pyroxene and plagioclase compositions from the gabbro and meladiorite samples are not in equilibrium with their respective host rock compositions (Fig. S2 [see ^{footnote 1}]), indicating that whole-rock compositions are too mafic to have been liquid compositions (e.g., Cornet et al., 2022). Mineral chemometers of Higgins et al. (2022) indicate that pyroxene equilibrium melts are andesitic to dacitic (57–68 wt% SiO₂; Table 2) in the lower GIC units, supported by plagioclase chemometers of Scruggs and Putirka (2018) (Table 2). Variation in whole-rock isotope chemistry, on the other hand, provides evidence for the emplacement of individual compositionally distinct magma batches. We thus argue that both mineral accumulation and the emplacement of compositional distinct magma batches occurred. Batch-specific whole-rock isotope signatures could be preserved if most mineral phases crystallized prior to mineral accumulation and melt loss (i.e., adcumulate formation). A rigid mineral framework at these relative high crystallinities would then prevent large-scale mixing between incoming magma batches, but interstitial liquid would still be capable of percolating through a crystal mush to collect in segregations.

The value of δ¹⁸O_(melt) can increase during closed-system fractionation due to crystallization of low-δ¹⁸O minerals (e.g., olivine, pyroxene, Fe-Ti oxides; e.g., Bindeman et al., 2004). The extent of this increase depends on the exact differentiation pathway (e.g., crystallization assemblage, temperature, melt SiO₂) and is ≤1‰ during fractionation from basalt to rhyolite (Bindeman et al., 2004; Bucholz et al., 2017). Zircon δ¹⁸O values would vary only by tenths of a per mil due to closed-system fractionation (Bucholz et al., 2017). The variations in δ¹⁸O_(WR) and δ¹⁸O_(zircon) in all GIC samples (Fig. 8) are larger than can be explained by closed-system fractionation, thus open-system behavior would have been required to explain the observed isotopic compositions.

The low δ¹⁸O_(WR) and δ¹⁸O_(zircon) values in the gabbro units contrast with higher values observed in the meladiorite and granite units (Fig. 8). Low δ¹⁸O_(WR) values in the GIC are similar to δ¹⁸O_(WR) values of deformed Mariposa Formation rocks at the base of the GIC (BMFZ; Fig. 8). Textural observations and the presence of melt leucosomes (Fig. 3A) indicate that the BMFZ was active prior, during, and after the emplacement of the GIC under upper amphibolite–facies conditions (Paterson et al., 1987; Vernon et al., 1989). The BMFZ likely provided a pathway for low-δ¹⁸O seawater and/or meteoric waters (δ¹⁸O_seawater: 0‰; Valley et al., 2005) to reach deeper parts of the crust during the marine Jurassic arc conditions, or seawater could also have been liberated from marine sediments of the Mariposa Formation during deformation (Paterson and Vernon, 2001). High deformation temperatures (>500 °C; Vernon et al., 1989) in the BMFZ would have allowed for water-rock oxygen isotope exchange to produce low-δ¹⁸O rocks in the vicinity of the lowermost units of the GIC. Tobisch et al. (1991) documented low δ¹⁸O_(WR) values (~+5‰) in mylonitic Santa Cruz Mountain tonalites (Figs. 1 and 8), interpreted as a result of fluid-enhanced deformation. In addition, a combination of contemporaneous partial melting of previously high-temperature, hydrothermally altered rocks of the Mariposa Formation in the BMFZ and GIC magmatism, as documented by field evidence (Fig. 3A) and zircon U-Pb ages (Fig. 1; Saleeby et al., 1989), would have produced high-silica, low-δ¹⁸O melts (Montel and Vielzeuf, 1997) that could have mixed with the gabbro crystal mush at the base of the GIC. To test such a scenario, we conducted assimilation–fractional crystallization (AFC)–mixing calculations using Sr_i and δ¹⁸O isotopes (Fig. S4; DePaolo, 1981). The composition of a GIC primary magma is unknown but is taken here to be similar to the “western primary magma” for the Sierra Nevada arc from Lackey et al. (2012). We assumed a Mariposa Formation sample as the assimilant with Sr_i of 0.7046 and 184 ppm Sr (Lackey et al., 2012) and δ¹⁸O_(WR) of +4.3‰ (lowest value of Mariposa Formation). AFC-mixing calculations show that the low gabbro δ¹⁸O_(WR) values can be generated using a ratio of assimilation to fractionation of 0.05–0.4. If the gabbro mush assimilated partial melts from the BMFZ, the amount was likely minor, such that the δ¹⁸O_(WR) values in the gabbros did not develop a low-δ¹⁸O magmatic signature (e.g., Troch et al., 2020). However, if the partial melts entered the gabbro mush (potentially through percolation processes) after most Fe-Mg phases had crystallized but before zircon saturation, the δ¹⁸O_(zircon) values should record the δ¹⁸O composition of a contaminated melt. The gabbro unit zircon rims show lower δ¹⁸O_(zircon) values than their cores (Fig. 9); this might reflect increasing assimilation of low-δ¹⁸O melts during zircon crystallization.

Contrary to the pluton-wide trend of increasing δ¹⁸O_(WR) from the exposed base to the top, the uppermost granophyre unit shows lower δ¹⁸O_(WR) values than the meladiorite unit below (Fig. 8). The Sr_i values for the granophyres are also distinct from those of the other GIC units: They show a trend toward higher Sr_i at similar ε_Ndi (Fig. 6). The relatively low δ¹⁸O_(WR) and higher Sr_i values of the granophyres are best explained by high-temperature isotopic exchange with seawater (δ¹⁸O: 0‰; ⁸⁷Sr/⁸⁶Sr: 0.7091; Hart et al., 1974; Valley et al., 2005). The widespread granophyric textures (Fig. 4O) attest to large degrees of undercooling at the top of the GIC, i.e., at shallow emplacement levels. The GIC intruded into slightly older marine Mariposa Formation greywacke and slates (Schweickert and Cowan, 1975). Exchange with seawater would have increased the Sr-isotope ratio of altered rocks but have a limited impact on ε_Nd (Hart et al., 1974). To acquire the observed low δ¹⁸O_(WR) values in the granophyres, alteration must have occurred during elevated temperatures, when isotope mass fractionation is small and altered rocks may approach the low δ¹⁸O values of seawater (Valley et al., 2005; Troch et al., 2020). Hydrothermal alteration likely took place shortly after solidification, when ambient temperatures were still high (~450 °C) and a hydrothermal system developed at the top of the GIC. This interpretation is supported by the pervasive alteration of the granophyres observed in outcrop (Fig. 3G), thin sections (e.g., sericitization in feldspars, chloritization of mafic minerals), and whole-rock compositions that indicates that some granite and granophyre samples experienced albitization (high Na index; Putirka et al., 2014). Assimilation of Mariposa Formation rocks as an explanation for the high Sr_i and low δ¹⁸O_(WR) values in the granophyres is unlikely because the typical Mariposa Formation rocks in the region have high δ¹⁸O_(WR) (as high as 13‰; Fig. 8; Lackey et al., 2012) and relative low Sr_i (~0.7046; Lackey et al., 2012). Equally, assimilation of rhyolites at the top of the GIC is also unlikely, given that they have a lower radiogenic Nd isotope composition than the GIC units (Fig. 6).

In the following, we summarize the main observations and interpretations resulting from this study and place them into a temporal context by combining them with geochronology results from Ratschbacher et al. (2018) and Saleeby et al. (1989) to arrive at a new conceptual model for the GIC:

• The shallow-marine Mariposa Formation was deposited on top of the Peñon Blanco and Gopher Ridge volcano-clastic formations. Silicic volcanism was active ~2–3 m.y. prior to the onset of GIC construction, as evident by zircon U-Pb crystallization ages of a rhyolitic unit interbedded with the Mariposa Formation (Ratschbacher et al., 2018; Attia et al., 2021; Fig. 12A).

• Construction of the GIC started with the emplacement of multiple isotopically and chemically distinct mafic to intermediate magma batches at ca. 148.5 Ma, originating from a deeper-seated magma reservoir undergoing differentiation, which formed the lower section of the pluton (Fig. 12B). Distinct compositions are manifested in different mineralogy (e.g., hornblende poor versus hornblende rich) and slightly different δ¹⁸O_(WR) and Sr-Nd isotope values. The low δ¹⁸O_(WR) and δ¹⁸O_(zircon) in GIC gabbros were likely produced through contact with the BMFZ at the base of the GIC, which modified the oxygen isotope composition of the melt by bulk assimilation and/or assimilation of partial melts derived from high-temperature, hydrothermally altered Mariposa Formation rocks (Fig. 12B).

• Construction continued for ~300 k.y. (Ratschbacher et al., 2018) with the emplacement of felsic magmas in the upper part of the GIC and coeval intrusions of mafic and felsic magmas in the mingling zone (Fig. 12C). The range of isotopic values preclude a simple closed-system, in situ differentiation origin for the entire pluton, but open-system differentiation in a crystal mush operated on a smaller scale (e.g., crystal-melt separation). Thermal modeling (Ratschbacher et al., 2018) indicates that melt was present in the central part of the GIC for an extended amount of time (maximum ~500 k.y. including construction and cooling time; pink ellipse in Fig. 12C) from 1000–900 °C liquidus temperatures to ~750–650 °C solidus temperatures (Fig. 5). This facilitated crystal-melt separation and melt migration of chemically and isotopically diverse interstitial melt into felsic segregations and pods (upper gabbro and meladiorite units). The presence of cumulates is also supported by textural observations and pyroxene, amphibole, and plagioclase mineral chemistry.

• Shortly after solidification, sub-solidus, high-temperature hydrothermal alteration (~450 °C) involving seawater took place at the roof of the complex, affecting δ¹⁸O_(WR) and Sr_i values of granite and granophyric rocks (Fig. 12D). At the base of the GIC, high temperatures persisted for 4–5 m.y., as indicated by re-crystallization textures in the lower gabbros (Fig. 4F) and ca. 144 Ma ⁴⁰Ar/³⁹Ar biotite cooling ages from the BMFZ adjacent to the GIC (closure temperatures 300–400 °C). Deformation continued along the BMFZ until ca. 137 Ma.

Integrating zircon U-Pb ages and thermal modeling with ⁴⁰Ar/³⁹Ar cooling ages (Saleeby et al., 1989) and mineral thermometry indicates that (1) the host rocks to the GIC were deposited at least ~2–3 m.y. before magmatism started; (2) construction of the GIC lasted for ~300 k.y., followed by ~200 k.y. of cooling to temperatures of ~750–650 °C; and (3) after solidification, a hydrothermal system developed at the top of the GIC likely caused rapid cooling, whereas the base of the GIC (BMFZ) remained at elevated temperatures (>300–400 °C) for ~4–5 m.y. despite the shallow emplacement level.

This study highlights the benefit of combining different geochemical tools to probe processes at various time (e.g., zircon versus rock-forming mineral chemistry) and length (e.g., detailed transects versus pluton-wide trends) scales to decipher the complex, multi-scale, and open-system magmatic and sub-solidus processes that ultimately result in the integrated chemical and isotopic characteristics of plutonic bodies. In particular, whole-rock data can be misleading if processes, such as crystal accumulation and open-system behavior, operate to modify compositions, which can only be traced by mineral chemistry (e.g., zircon oxygen isotopes and trace elements). However, using mineral chemistry without a context about host magma composition can be misleading because zircon would have crystallized relatively late in mafic samples compared to felsic samples and thus records different stages of melt evolution. This has implications for detrital zircon studies inferring processes such as crustal thickness variations from zircon geochemistry alone (e.g., Tang et al., 2021).

The authors thank Andrew Kylander-Clark at UCSB for help obtaining zircon trace element data, Leslie Hayden at the USGS for help with electron microprobe analyses, and Elizabeth Bell at UCLA for help with secondary ion mass spectrometry (SIMS) analyses. Snir Attia provided a basemap for Figure 1. We also acknowledge the feedback and comments provided by two reviewers and Guest Associate Editor Jonathan Miller that helped to improve the manuscript.