Influence of pre-existing structure on pluton emplacement and geomorphology: The Merrimac plutons, northern Sierra Nevada, California, USA

The relationship between magmatism and tectonics is key to understanding processes related to crustal growth. Plutons record major inputs of new material into the crust, and on a batholith scale, pre-existing structure appears to have influenced the composition and geometry of Mesozoic arcs. For example, compositional, isotopic, and geophysical data in the Peninsular Ranges and Sierra Nevada of California (USA) mark suture zones (Fig. 1) within the batholiths that indicate emplacement across accreted oceanic arc and transitional continental crust (Silver and Chappell, 1988; Langenheim et al., 2014; Kistler, 1990; Oliver, 1977). Here we examine the influence of pre-existing structure at a more local scale on selected plutons in the northern Sierra Nevada (Fig. 2) not only on the composition and geometry of the plutons, but also its resulting influence on the geomorphic expression of these plutons.

Although much study has focused on the magmatic evolution of the central and southern Sierra Nevada, fewer studies have concentrated on the plutons of the northern Sierra Nevada (i.e., north of latitude 39°N). In the northern Sierra Nevada, several Late Jurassic to Early Cretaceous, syn- to post-tectonic plutons intrude fault-bounded belts of late Paleozoic to early Mesozoic rocks with oceanic affinities in the Sierra Nevada metamorphic belt. We chose one set of these plutons, the Merrimac plutons (Fig. 3), for in-depth study because their foliation and compositional patterns had been previously studied by one of us (Guglielmo). The Merrimac pluton was originally described as a single, compositionally zoned pluton (Hietanen, 1951, 1973, 1976, 1977) because limited outcrops curtailed recognition of internal contacts. Subsequently, Guglielmo (1993a) used foliation and compositional data to interpret two nested zoned plutons, with a southern, mostly granodioritic to tonalitic pluton and a northern, mostly tonalitic pluton (Fig. 3E). However, the relationship of these plutons to each other remained hampered by poor exposure (Fig. 2H) and inconclusive unpublished U-Pb dates for zircons of 142 ± 3 Ma reported in Guglielmo (1993a). This paper updates and refines the work of Guglielmo (1993a) with the addition of analysis of gravity and magnetic data and new ²⁰⁶Pb/²³⁸U dates for delimiting the crystallization age of zircons in order to refine the number, shape, and contact relations between the Merrimac plutons and estimate the depth extent of the plutons.

In contrast to the central and southern Sierras, plutons in the north are dispersed and separated by wide tracts of structurally disrupted metamorphic country rock. The wide variety of rock types leads to large variations in density and magnetic properties that are amenable to geophysical study by potential-field methods to investigate the relationship of plutons to regional structural trends. This study integrates geophysical, geochronologic, and geomorphic data to examine the role of pre-existing structure in the geometry and composition of the Merrimac and adjacent plutons. We present these results within a regional context of the northwestern Sierra Nevada, supplemented by dates on adjacent plutons.

The Sierra Nevada batholith is a north-northwest–trending plutonic complex that extends for >500 km east of the Great Valley (Fig. 1). In contrast to the dominance of plutonic rocks in the central and southern parts of the batholith, plutons of the northern Sierra Nevada (i.e., north of latitude 39°N) are widely dispersed, exposing range-parallel belts of deformed and metamorphosed country rock of both continental and oceanic affinity that compose the metamorphic belts of the Sierra Nevada foothills (Saleeby et al., 1989a; Cecil et al., 2012; Rousseau, 2016). At the northwestern tip of the greater Sierra Nevada batholith, the Merrimac plutons are part of a suite of Late Jurassic and Early Cretaceous plutons that intruded the Western and Central metamorphic belts of the Sierra Nevada foothills (Figs. 1 and 2). South of these plutons, the metamorphic belts and structural features within the belts trend north-south; structural trends bend sharply toward the west in the vicinity of the plutons (Hietanen, 1976; Fig. 2A). These plutons were emplaced during an apparent intrusive slowdown between the Middle Jurassic and Late Cretaceous magmatic flare-ups in the Sierra Nevada (Paterson and Ducea, 2015). Aluminum-in-hornblende geobarometry indicates pressures of 2–3 kbar for the emplacement of the Early Cretaceous Merrimac plutons, 3 kbar for the roughly coeval Grizzly pluton and Bald Rock plutons, 3 kbar for the Late Jurassic Yuba Rivers pluton, and 4 kbar for the Late Jurassic Cascade pluton (McLeod, 1990). These pressures are consistent with that predicted from the presence of staurolite with andalusite in the outer aureole of the Merrimac pluton, which indicates metamorphism at upper temperatures of 500–600 °C and pressures of ∼4 kbar (Hietanen, 1973).

The Merrimac plutons intruded the Central belt, a structurally complex, disrupted mixture of ophiolitic, volcanic, and sedimentary rocks divided into the Fiddle Creek and Slate Creek Complexes (Day and Bickford, 2004; Fig. 2A). The Big Bend fault zone (BBF in Fig. 2) separates these rocks from the Western belt, a relatively coherent assemblage of Middle and Late Jurassic volcanic and volcanoclastic rocks (the Smartville Complex) overlying and intruding the Fiddle Creek Complex and the Slate Creek Complex to the south of the area shown in Figure 2. Day and Bickford (2004) inferred that the Big Bend fault zone has generated only minor displacements since ca. 160 Ma because U-Pb geochronology of cross-cutting dikes to the south of the study area implies that the Western and Central belts have been part of the same terrane since the Middle Jurassic. The Bald Rock pluton intrudes across the Big Bend fault zone, thus precluding displacement since the time of intrusion. Note that Hietanen (1973) mapped the Big Bend fault zone as curving around the Bald Rock pluton but cut by a tongue of Bald Rock trondhjemite that extends east into the Cascade pluton to the east (Fig. 3A). The Big Bend fault zone dips steeply to the northeast and is several kilometers wide where it is best exposed and mapped along its western extent (Dilek, 1989) near where it projects beneath younger cover. Dilek (1989) suggested that this fault accommodated sinistral displacement between similar rocks in the Klamath Mountains and northern Sierra Nevada. Others, such as Compton (1955), suggested that the significant change in strike in both the fault and older rocks was caused by intrusion of plutons, specifically the Merrimac, Bald Rock, and other Early Cretaceous plutons.

Early studies of the Merrimac and Bald Rock plutons inferred a combination of emplacement processes. Both Compton (1955) and Hietanen (1951) argued for forceful emplacement of magma based on doming and buckling of previously folded country rock around the pluton bodies. They also suggested that stoping and assimilation contributed to pluton size based on observations of brecciation and migmatites along the pluton contacts with the country rock. Compton (1955) suggested that the compositional zoning of the Bald Rock pluton from trondhjemite in the center to tonalite along the pluton margin was due to contamination. Larsen and Poldervaart (1961) argued for little wall-rock contamination based on the crystal habits of zircon microphenocrysts. Subsequently Sawka (1985) used geochemical and initial strontium isotopic data to rule out wall-rock contamination as a source of compositional zoning, instead arguing for sidewall crystallization and subsequent magma mixing. He noted that variations in initial strontium ratios could reflect multiple stages of intrusion rather than wall-rock contamination.

In a regional geomorphic sense, the Merrimac plutons reside within the relatively stable Cenozoic Sierra Nevada block. This block is characterized by a gentle westward slope and a steep eastern escarpment marked by normal and dextral faults. The Merrimac region slopes westward at a gradient of ∼38 m/km over a distance of ∼50 km (Wakabayashi, 2013). The timing of uplift and incision in the Sierra Nevada block is controversial. According to Wakabayashi (2013), uplift, tilting, and stream incision began ca. 3 Ma, superposed on a relic pre-Eocene surface; minimal incision occurred from Eocene to Miocene time based on relief between the elevation of basement topographic highs and the elevation of the base of nearby Cenozoic strata. Others argue that the Sierra Nevada block has been a high-standing block since the Eocene or even the Cretaceous, based on paleoaltimetry, sedimentology, and geomorphologic analysis (e.g., Cassel et al., 2012; Gabet, 2014). As Gabet (2014) discussed, the ability to use of marker horizons, such as volcanic rocks deposited upon older bedrock, to estimate incision rates strongly depends on the ability to accurately identify the original river elevation. Additionally, Gabet (2014) argued that present gradients of Miocene volcanic rocks in the area are within the range expected for primary depositional slopes, consistent with little post-depositional westward tilting. Hence, it is possible that that deep canyons of the Feather River thought to have been recently incised were instead formed as early as the Eocene–Oligocene.

This study integrates geophysical, geochronologic, and geomorphic data to examine the role of pre-existing structure in the geometry and composition of the Merrimac plutons. The geophysical data, in particular the aeromagnetic and physical property data, were used to refine the geologic map in Guglielmo (1993a). We used gravity data (Langenheim et al., 2019) to produce isostatic gravity maps (Figs. 2B, 4A), which reflect density variations in the shallow to middle crust (Simpson et al., 1986). Accuracy is estimated to be 0.5 mGal or less.

Aeromagnetic data reflect crustal magnetization variations. The aeromagnetic map in Figure 2C is based on several surveys (Langenheim, 2015; U.S. Geological Survey, 1984, 2002) that were placed in a common datum and merged. Most of the area, including over the Merrimac plutons (Fig. 5A), was flown at a nominal height of 305 m (1000 ft) above ground along flight lines spaced 0.8 km (0.5 mi) apart. Data over the extreme southeastern corner is of considerably lower resolution because of flight height or spacing. Accuracy of the modern surveys is 1 nT or better. In limited areas, some artifacts are present in the data; for example, northeast-trending anomalies northeast of the Cohasset Ridge fault are likely flight-line artifacts because the anomalies are parallel to the flight-line direction, extend across multiple rock types, do not follow topographic grain, and have a width commensurate with the flight line spacing.

Density and magnetic properties of rock samples are crucial parameters for analyzing and modeling gravity and magnetic data. To augment density measurements of plutons by Hietanen (1973, 1976) and Larsen and Poldervaart (1961), we measured density and magnetic susceptibility of rock slabs collected for compositional analysis by Guglielmo (1993a) and of additional hand samples of not only the Merrimac plutons (Figs. 4C, 5C; Langenheim et al., 2019) but also of the Bald Rock pluton and country rock, summarized in Table 1. Densities were determined using a precision Sartorius electronic balance. Magnetic susceptibilities were measured with a Geophysica KT-5 susceptibility meter and are reported to a precision of 0.01 × 10⁻³ to 0.1 × 10⁻³ SI units.

To help delineate structural trends and gradients expressed in the gravity and magnetic fields, a computer algorithm was used to locate the maximum horizontal gradient (Cordell and Grauch, 1985; Blakely and Simpson, 1986). Contacts between various units of contrasting density and magnetization were mapped with the help of horizontal gradients in the gravity and magnetic fields. We calculated magnetization boundaries on the residual magnetic potential field data. The residual magnetic potential field data were obtained by first calculating a residual magnetic field by subtracting a regional field from the observed magnetic field. The regional field was calculated by upward continuation of the magnetic field by 200 m, producing a magnetic field as if it were collected 200 m higher and thus minimizing the effect of magnetic sources in the upper 1–2 km of the crust. This residual magnetic field (Fig. 2G) was then converted into the magnetic potential (also known as pseudogravity) field, which centers the anomaly over the magnetic source and produces the equivalent “gravity” field as if all magnetic material were replaced by proportionately dense material (Baranov, 1957). The horizontal gradient of the residual magnetic potential field was then calculated. Gradient maxima occur approximately over steeply dipping contacts that separate rocks of contrasting densities or magnetizations. For moderate to steep dips (45° to vertical), the horizontal displacement of a gradient maximum from the top edge of an offset horizontal layer is always less than or equal to the depth to the top of the source (Grauch and Cordell, 1987). We also produced a residual gravity field by subtracting a regional field that was upward continued 1 km to highlight short-wavelength, near-surface density features that would be more directly comparable to mapped contacts at the ground surface (Fig. 2F). We calculated density boundaries from the isostatic gravity anomaly field rather than the residual gravity anomaly field because of the nonuniform distribution of gravity measurements.

To determine the depth extents and geometry of the Merrimac plutons, we used the Geosoft GM-SYS program (https://www.geosoft.com) to model gravity and magnetic data (Fig. 6). We chose profile A-A′ (Figs. 3–5) to maximize nearby gravity measurements and to be collocated with the cross section in Guglielmo (1993a, his figure 7b). Because gravity and magnetic modeling is non-unique, we used geology and physical property measurements to provide an initial estimate of model parameters (depth, shape, and density of suspected sources). The amplitude of the anomaly is not the only attribute to match; matching its gradients provides important information on source depth and shape.

We collected samples for U-Pb dating from the two main Merrimac plutons and adjacent plutons (Cascade, Granite Basin, Bald Rock, and Bucks Lake; Fig. 2A), where the only previously published ages were obtained by thermal-ionization mass spectrometry of multi-crystal fractions. These older data are associated with somewhat large uncertainties and yield a wide range of crystallization ages for the same pluton. Single zircon crystals were separated from rock samples at the U.S. Geological Survey (USGS) in Menlo Park (California) by standard methods of grinding, heavy liquids, and magnetic separation. These zircon grains, including fragments of the Temora-2 zircon standard with a ²⁰⁶Pb/²³⁸U date of 418.4 Ma (Mattinson, 2010), were mounted in epoxy, sectioned to reveal crystal interiors, and polished. The grains were then imaged at Stanford University (Stanford, California) with a scanning electron microscope operated in secondary electron and cathodoluminescence mode in order to identify compositional zoning, cores and overgrowths, and fractures.

The zircons were analyzed using the sensitive high-resolution ion microprobe–reverse geometry (SHRIMP-RG) managed jointly by the USGS and Stanford University. Data were gathered for U-Pb geochronology and trace element concentrations (Table S1¹) with a data acquisition protocol and analytical conditions as described in Matthews et al. (2015). Raw data were reduced using the SQUID2 software version 2.51 of Ludwig (2009). Trace element concentrations were calculated from secondary ion yields normalized to ⁹⁰Zr₂¹⁶O⁺ calibrated to the yields from concurrent analyses of the MAD-559 (Coble et al., 2018) zircon concentration standard. Titanium-in-zircon temperatures were derived from the formulation of Ferry and Watson (2007) assuming a SiO₂ activity of 1.0 and TiO₂ activity of 0.7, values that are appropriate for quartz-saturated and titanite-saturated silicic magmas (Claiborne et al., 2006; Hayden et al., 2007), as is the case for the studied plutons. Reported ²⁰⁶Pb/²³⁸U dates and their uncertainties were calculated with Isoplot (Ludwig, 2003) version 3.76 using a Pb/U versus UO/U calibration (Williams, 1998) to the Temora 2 age standard and correction for common Pb using age-appropriate values from the Stacey and Kramers (1975) crustal evolution model. Of the 269 analyses (Table S1 [^{footnote 1}]), 235 were from the rims of individual grains; the remainder were paired core-rim analyses from 17 grains.

We characterized the geomorphology of the Merrimac plutons and the surrounding country rock as an indicator for relative erosion in response to late Neogene (ca. 5 Ma) accelerated uplift (Wakabayashi and Sawyer, 2000, 2001) using the best available digital elevation model (DEM) at a cell size of 10 m from The National Map (U.S. Geological Survey, 2019). We generated a range of topographic indices to explore differences in landscape and channel morphology within and between different geologic units. First we used ESRI ArcGIS (version 10.6.1) with the TauDEM toolbox (http://hydrology.usu.edu/taudem/taudem5 index.html) to derive a hydrologically conditioned (pit removed) DEM in order to calculate a continuous stream network and associated derivatives such as drainage density (D_d), Strahler stream order, and stream power. The drainage delineation approach of TauDEM uses the steepest downslope flow direction and a Peucker Douglas stream definition, as described in Tarboton (1997) and Peucker and Douglas (1975). Drainage density D_d is defined as the ratio between the total stream lengths (l) and drainage area (A), where higher values of D_d generally correlate with areas of more evolved erosion or less competent rock. Stream order classification was adopted from Strahler (1952). To examine local topographic characteristics at geologic unit boundaries, we extracted longitudinal profiles of the primary throughgoing drainages that traverse both the northern or eastern and southern Merrimac plutons to identify possible knickpoints, as expressed by changes in channel slope arising from variable rock resistance (or uplift).

To represent the long-term, spatially integrated influence of mass wasting and channel incision on the landscape morphology, ESRI ArcGIS focal mean statistics and Whitebox Geospatial Analysis Tools (https://www.uoguelph.ca/∼hydrogeo/Whitebox/) (Lindsay, 2016) were used to generate estimates of local relief (100 m circular search radii), ruggedness, and multiscale maximum elevation deviation from mean elevation of not only the Merrimac plutons themselves, but also the broader region as a means to evaluate influence of variable rock resistance and valley incision on landscape form. Ruggedness, a measure of local topographic relief, was calculated following the procedure of Riley et al. (1999) such that the root-mean-square deviation for each grid cell in a DEM denotes the residuals (i.e., elevation differences) between a grid cell and its eight neighbors. Given that elevation residual attributes, such as ruggedness, are inherently scale dependent because they are defined in the context of a local neighborhood search of grid cells, we also calculated the maximum deviation from mean elevation (DEV_max) across a range of spatial scales in accordance with Lindsay et al. (2015). The DEV_max integral image-based approach generates a unitless measure of topographic position that is scaled by the local ruggedness. DEV_max images were calculated using a minimum search window radius of 30 m and a maximum radius of 10,230 m with a step between searches of 100 m. The resulting raster images were combined into a single color-composite image providing a synthesized representation of multiscale relief properties with output values ranging from −3.0 to 3.0 where negative values represent topographic valleys. The geomorphic indices are intended to represent spatially varying functions reflecting the interplay between the underlying bedrock properties, spatial distributions of erosive transport processes, and crustal deformation.

Gravity data reveal a prominent gravity low centered on the Bald Rock pluton, with low values extending northward across a septa of country rock and the southern part of the Merrimac plutons and northeastward across country rock and the northern part of the Cascade pluton, coincident with a tongue of Bald Rock trondhjemite (Figs. 2A, 2B, and 2F). Gravity highs coincide with Western and Central belt rocks, most notably to the south and east of the Bald Rock gravity low. The highest gravity values, surprisingly, are not located over exposed basement rock, but over Quaternary deposits southwest of the Chico monocline; this high is part of a much larger inferred basement feature that extends south-southeastward for ∼250 km concealed beneath the Sacramento Valley (thin blue lines in Fig. 1; Cady, 1975) and connotes a very large body of very dense rock beneath the Cenozoic and relatively thin Cretaceous strata. Although the largest and most prominent gravity low is associated with the Bald Rock pluton, all of the other early Cretaceous plutons, excluding the Bucks Lake pluton, coincide with gravity lows relative to values measured on adjacent country rock; these lows are accentuated in the residual gravity field (Fig. 2F).

Aeromagnetic anomalies generally correspond with mafic and ultramafic rock types, most notably the Feather River ophiolite in the northeastern corner of Figure 2. Magnetic highs associated with serpentinite, metagabbro, and peridotite partially encircle many of the Early Cretaceous plutons, particularly the Merrimac, Concow, Bald Rock, and Granite Basin plutons, and the eastern side of the Late Jurassic Cascade pluton. Although variations in basement magnetizations produce most of the anomalies in Figure 2, some of the short-wavelength anomalies arise from the overlying Cenozoic volcanic rocks. The sinuous magnetic lows in the western part of the map area are caused by Miocene Lovejoy Basalt, exposed at Table Mountain (TM in Fig. 2G) but mostly concealed beneath the valley floor. The Lovejoy Basalt flows are characterized by a strong remanent magnetization of negative polarity (Coe et al., 2005).

Both gravity and aeromagnetic anomalies reflect the regional structural grain exposed along the foothills (Fig. 2), from north-south orientations in the south, curving to east-west in the central map area, and curving to northwest-southeast in the north. The Big Bend fault west of the Bald Rock and Merrimac plutons is well expressed in the regional potential-field data (Fig. 2). Southwest of the Concow pluton, the gravity and magnetic highs bound by the Big Bend fault on their southern margin appear westward to merge with, or be truncated by, a feature marked by fairly continuous, north-northwest–striking magnetic highs and gravity lows, a view enhanced in the residual magnetic and gravity anomalies in Figures 2D and 2G. These north-northwest–striking anomalies appear to coincide with outcrops of serpentinized ultramafic rocks that are sporadically exposed protruding from beneath Cenozoic volcanic cover (Fig. 2G). These anomalies project toward the Cohasset Ridge and Beaver Creek faults (red lines, Fig. 2A), which are Quaternary faults with vertical slip rates of 0.02 and 0.05 mm/yr, respectively, based on offset of a ca. 1 Ma basalt flow (Wakabayashi and Sawyer, 2000); note that the existence of strike-slip displacement is unknown, but possible. The general persistence of these north-northwest–striking geophysical anomalies does not support direct projection of the Big Bend fault across the Great Valley, as suggested by Dilek (1989). If the Big Bend fault did continue westward, its westward trace would be offset some distance to the north by normal offset on the younger cross-cutting faults; however, no such feature is obvious in the geophysical data within the area of Figure 2.

The Merrimac plutons are bisected by gravity and aeromagnetic gradients that support the presence of at least two plutons, confirming and refining the original analysis by Guglielmo (1993a). A significant gravity gradient trends N60°W, with higher values to the northeast (Fig. 4). This gravity gradient coincides with a magnetic gradient, also with higher values to the northeast (Fig. 5). The potential-field anomalies indicate that the rocks northeast of the gradients are denser and more magnetic relative to rocks southwest of the gradients, which is borne out by density and magnetic susceptibility measurements (Figs. 4 and 5) and is consistent with composition data of Guglielmo (1993a) which show that the northern pluton has somewhat lower potassium feldspar (Fig. 3D) and higher mafic contents. The gravity and magnetic gradients are more linear than the originally inferred contact between the two plutons (Figs. 3A, 3E) and indicate that the diorite and quartz diorite mapped in the eastern part of the plutons originally included as part of the southern pluton in Guglielmo (1993a) are more closely related to the northern pluton.

The new geophysical and physical property data, along with foliation and compositional data, identify additional contacts within the Merrimac plutons. Both lithology and magnetic susceptibility exhibit considerable variability at outcrop scale, but gravity and magnetic fields are less sensitive at that scale. Integration of these two different scales of mapping with foliations defined by alignment of biotite and hornblende crystals and interpreted to be magmatic (Guglielmo, 1993a) reveals a third pluton on the east (Fig. 3A, 3B). A northeast-trending magnetic gradient in the northern pluton appears to separate a more mafic eastern pluton from the rest of the northern pluton as suggested by the composition data. This eastern pluton includes a body of diorite (Fig. 3C), the most mafic rock type in the Merrimac plutons, which is parallel to and elongated along the southern margin of the eastern pluton. The contact between the eastern and northern plutons is refined by foliation measurements, where those in the eastern pluton are nearly orthogonal to those in the adjacent northern pluton (Fig. 3B). Where the foliation patterns are more ambiguous, gradients in magnetic anomalies help guide the location of the northwestern contact of the eastern pluton (Fig. 5). The magnetic gradients that mark the boundary between the northern and eastern plutons bracket a similarly oriented drainage (Fig. 5E). The drainage may reflect a zone of less-resistant rock along the pluton contact. Another difference between the northern and eastern plutons is that foliation in the eastern pluton is weaker and includes compositional cumulate textures. The southern margin of the eastern pluton is defined by a steep magnetic gradient that is aligned with the magnetic gradient that separates the southern and northern plutons.

The location of the maximum horizontal gravity and magnetic gradients 1–2 km northeast of the inferred pluton contacts between the grouped northern and eastern plutons and the southern pluton suggests that this boundary dips steeply to the northeast; such an orientation is borne out by the joint gravity and magnetic model (Fig. 6). The N60°W strike and steep northeast dip of the contact between the southern and grouped northern and eastern plutons is parallel to the dominant metamorphic fabric in the country rock as well as to the trend of the Big Bend fault 7–8 km to the southwest. The observed anomalies can be fit assuming average density and magnetic susceptibility values of the northern and southern plutons and that they extend to depths of ∼5 km. In contrast, lower gravity values over the Bald Rock pluton indicate the pluton extends to depths of >12 km. Note that if the true pluton density is greater, the density contrast is less, and thus its depth extent is greater.

New ²⁰⁶Pb/²³⁸U zircon dates indicate that the northern and southern Merrimac plutons were assembled at different times and refine zircon crystallization ages of nearby plutons (Table 2; Table S1 [^{footnote 1}]). In some of these samples, a minority of zircons yield dates that are younger than the dominant population, likely reflecting sampling of crystal domains that have experienced Pb loss (e.g., Watts et al., 2016). Accordingly, these outliers are excluded from the apparent populations that are taken to reflect the last interval of crystallization and in turn the age of pluton emplacement. In addition, several zircons yield Mesozoic dates that are “old” outliers relative to the dominant population and are interpreted to represent antecrysts (e.g., Samperton et al., 2015).

The southern Merrimac pluton yields a weighted mean crystallization age of 133.9 ± 1.3 Ma, whereas the northern Merrimac pluton is statistically older, with a mean crystallization age of 140.2 ± 1.6 Ma (Fig. 7A). This pair of mean crystallization ages differs from the 142 ± 3 Ma U-Pb dates for both plutons cited by Guglielmo (1993a) based on thermal ionization mass spectrometry (TIMS) analyses of multi-crystal fractions of zircon, as well as a SHRIMP-derived crystallization age of 142 Ma reported in Irwin and Wooden (2001; no uncertainty given) for a sample of the northern pluton near the contact with the eastern pluton. However, recalculation of the original U-Pb multi-crystal data reveals discordance for some of the sample aliquots, likely reflecting Pb loss, with the oldest concordant fractions yielding ²⁰⁶Pb/²³⁸U dates of 140.7 ± 0.7 Ma and 134.1 ± 0.6 Ma for the northern and southern plutons, respectively (Table S2; Fig. S1 [^{footnote 1}]). Hence, the ²⁰⁶Pb/²³⁸U dates and interpretation of Pb loss from this study are supported by the previous analyses of bulk zircons. The recalculated TIMS date for the southern Merrimac pluton is nearly identical to one of our Bald Rock pluton U-Pb zircon dates (133.7 ± 2.2 Ma; Fig. 7B) and a recently published U-Pb date of 133.4 ± 3.3 Ma for zircons from the Bald Rock pluton (Rousseau, 2016). These dates, combined with the large Bald Rock gravity low that extends across the southern Merrimac pluton, suggest that the southern Merrimac pluton and the Bald Rock pluton are closely related in terms of intrusion age and geophysical signature. However, our crystallization age of a sample close to the eastern margin of the Bald Rock pluton (Fig. 2A) is considerably older (142.0 ± 2.0 Ma; Fig. 7B), suggesting that the Bald Rock pluton may have grown over a period of several million years or, more likely, reflects multiple intrusions like the Merrimac suite.

New ²⁰⁶Pb/²³⁸U dates for zircons from other plutons (Table 2) indicate that the Bucks Lake pluton, both the rim and center portions, intruded at 139.0 ± 2.0 Ma, i.e., effectively during the same time span as the northern Merrimac pluton. This new Bucks Lake date is indistinguishable from a previously published date of 140 Ma based on TIMS ²⁰⁶Pb/²³⁸U geochronology of multi-zircon-crystal fractions (Saleeby et al., 1989a; date S140 in Fig. 2A), as well as reported dates of 138 Ma and 141 Ma generated via the same technique for the Grizzly and Concow plutons, respectively (date S140, D<141, and I141 in Fig. 2A; Saleeby et al., 1989a; Irwin and Wooden, 2001). New crystallization ages for zircons from the center and rim portions of the Granite Basin pluton (Fig. 7C) are virtually indistinguishable from each other (141.6 ± 1.4 Ma versus 141.8 ± 1.0 Ma, respectively). All of these are within the uncertainty of the older northern Merrimac and Bald Rock pluton dates, indicating a period of increased intrusive activity. However, geophysical differences between the plutons suggest that they are not all distributed parts of a single intrusion.

Our crystallization ages suggest that the Cascade pluton intruded at ca. 155 Ma (Fig. 7C) and is somewhat older than the nearby plutons. This apparent age agrees with the mean U-Pb date for zircons determined by Rousseau (2016; 151.5 ± 5.9 Ma; mean-square weighted deviation [MSWD] = 11.4). An ⁴⁰Ar/³⁹Ar date of 149.8 ± 0.2 Ma for Cascade pluton hornblende (Fagan et al., 2001), which would reflect pluton cooling below ∼500 °C (Harrison, 1982), is consistent with an intrusion age of ca. 155 Ma. It is noteworthy that Day and Bickford (2004) derived a concordia intercept ²⁰⁶Pb/²³⁸U date of 143 ± 2 Ma for Cascade pluton zircons, which they inferred to reflect the presence of a younger domain within the pluton boundaries based on the discordance of their date with the aforementioned ⁴⁰Ar/³⁹Ar cooling age. One fraction of their bulk zircon analyses yielded somewhat elevated ²⁰⁷Pb/²⁰⁶Pb, which they interpreted to reflect a component of Paleoproterozoic zircon in their multi-crystal separates. Rousseau (2016) noted two distinct Cascade zircon morphologies: (1) older prismatic zircons with oscillatory zoning and truncated overgrowths, and (2) younger zircons with inherited cores. This observation, along with the considerably high MSWD (i.e., scatter of dates), led Rousseau to model his zircon age data as three populations of crystallization ages at 165.7 ± 3.9 Ma, 145.0 ± 4.7 Ma, and 125.1 ± 8.7 Ma. We interpret these complex age distributions for the Cascade pluton to result from both Pb loss and recycling of antecrysts. Despite the complex age distributions from this pluton, the interpreted age of ca. 155 Ma is supported by the ca. 150 Ma ⁴⁰Ar/³⁹Ar cooling age reported by Fagan et al. (2001) and is similar to that obtained from the nearby Hartman Bar and Lumpkin plutons (155.8 ± 3.0 Ma and 150.8 ± 4.7 Ma, respectively; Rousseau, 2016).

Trace element concentrations in the dated zircons (Fig. 8) provide insight into petrologic differences within and between the plutons, which in turn can be related to magmatic evolution. Titanium concentrations in zircon are positively correlated with crystallization temperature (Watson and Harrison, 2005) such that cooling-induced fractionation of granitic magmas is typically associated with decreasing Ti concentrations in zircon (Barth and Wooden, 2010). Uranium and Hf typically increase in zircon with increasing degree of fractionation in granitic magmas (e.g., Barth and Wooden, 2010; Claiborne et al., 2006), whereas divalent Eu is preferentially partitioned into plagioclase and orthoclase, leaving melts with lowered Eu/Eu* (europium anomaly, where Eu* is the expected Eu value for a smoothed rare earth element pattern normalized relative to CI chondrite) that is imparted to the rare earth element (REE) distribution of crystallizing zircons (e.g., Hoskin and Schaltegger, 2003). The positive correlation between europium anomaly (Eu/Eu*) and Hf as well as U, and the negative correlation between Eu/Eu* and Ti, in zircons from the northern and southern Merrimac granitoids are opposite to the trends expected if the zircons grew from silicic melts reflecting fractional crystallization (e.g., Claiborne et al., 2006), but could reflect parental melts related by magma mixing of high Hf-Ti and low Hf-Ti end members.

Trace elements in zircons from the Granite Basin pluton roughly follow the trend expected for crystallization from a series of melts related by fractional crystallization, apparently without distinction between the crystals from the pluton center or rim. The negative correlation between Hf and Nd/Yb, Ti, and Eu/Eu* for the Granite Basin zircons (Fig. 8) suggests crystallization from melts related by cooling and fractionation of silicic magma (e.g., Samperton et al., 2015). Apparent temperatures of crystallization are mostly in the ∼750–600 °C interval, indicating the zircons formed at the near-solidus conditions of this granitoid magma (cf. Boettcher and Wyllie, 1968; Johannes, 1984).

Zircons from the rim and center of Bucks Lake pluton are nearly identical in some of their trace element concentrations (Fig. 8). However, zircons from the pluton center have higher Hf concentrations and lower Eu/Eu* than those from the pluton rim, suggesting crystallization from somewhat more evolved melts (e.g., Claiborne et al. 2006) than those responsible for zircons from the pluton rim, although the concentration of U is generally constant and the same for both center and rim. The range of apparent crystallization temperatures (∼750–800 °C) recorded by the Bucks Lake zircons as indicated by the Ti-in-zircon geothermometer is somewhat higher than for zircons from the other plutons, which is consistent with the relatively unevolved quartz diorite to pyroxene diorite (Saleeby et al., 1989a) composition of the Bucks Lake pluton rim and center, respectively. The broad range of Nd/Yb in Bucks Lake zircons (Fig. 8) may indicate that light REE–rich minerals, such as allanite or monazite, were not saturated or fractionated during the crystallization interval recorded by the zircons (cf. Samperton et al., 2015). In contrast, lower-temperature zircons from the more-evolved Granite Basin pluton have Nd/Yb ratios that form coherent trends when coupled to an index of differentiation such as Hf concentration (Fig. 8).

Zircons from the rim and interior portions of the Cascade pluton have similar trace element compositions but record generally different Th/U ratios. Temperature-dependent changes in Th/U partitioning can be ruled out as a cause for this difference because the zircons from both portions of the pluton yield indistinguishable crystallization temperatures using the Ti-in-zircon geothermometer (Fig. 8). Differences in Th/U could reflect melt composition and extent of melt evolution, with a decrease in the ratio typically correlating with an increase in degree of melt fractionation (e.g., Samperton et al., 2015) and/or crystallization of U- or Th-rich accessory minerals. Hafnium concentrations, Eu/Eu*, and U versus Eu/Eu*, which also mirror degree of melt evolution in granitoids (Claiborne et al., 2006; Barth and Wooden, 2010), are similar between the Cascade zircons from the pluton rim and interior, suggesting that the magmas responsible for the rim and interior of the Cascade pluton evolved similarly but differed in their crystallization of accessory minerals. Zircons from the ca. 142 Ma eastern margin of the Bald Rock pluton are similar to those from the Cascade pluton in terms of Hf, Eu/Eu*, Th/U, and apparent crystallization temperatures of ∼600–700 °C (Fig. 8). However, Bald Rock zircons from the ca. 134 Ma pluton interior yield a larger range of apparent crystallization temperatures (∼600–850 °C) and U concentrations. Zircons from the Bald Rock eastern margin have relatively constant low U concentrations and lack a negatively correlated U versus Eu/Eu* that is trend indicative of fractionation, despite their variability in terms of Eu/Eu* and Nd/Y and in turn fractionation of feldspars and accessory minerals.

Color shaded-relief topography of the Merrimac pluton area, derived from the 10 m DEM, reveals the large drainage area of the North Fork and Middle Fork Feather River, but also numerous, isolated smaller watersheds (Fig. 9A). Topographic relief calculated within a 100-m radius search window (hereafter “100 m relief”), ruggedness, and the DEV_max images highlight deeply incised, steep-walled valleys of the Feather River surrounding lower-relief upland drainages within the boundaries of the Merrimac plutons (Figs. 9B, 9C, 9D). Forks of the Feather River appear to spatially bracket the Merrimac plutons. The 100 m relief, ruggedness, and DEV_max indices all similarly depict Jurassic and Cretaceous plutons that are generally characterized by substantially lower relief relative to adjacent metamorphic country rock.

Surface elevations of the northern and eastern Merrimac plutons are on average 400 m higher than those of the southern Merrimac pluton (1132 ± 185 and 762 ± 111 m, respectively; Fig. 10A). Because the northern and eastern plutons reside higher within the drainage system with a lower stream order (Fig. 10B), the change in mean elevation may simply be due to proximity to the drainage divide, not to differing erodibility or structural offset. However, the abrupt change in slope at the contact between the northern and southern plutons lies at the base of generally steep south-facing topography within the northern pluton and runs parallel to and essentially coincident with the pluton contact (Figs. 9, 11).

To explore whether the slope difference between the northern and southern Merrimac plutons results from compositional influence on erodibility (e.g., Eggleton, 2017) or from structure, we examined differences in drainage density, D_d. Drainage density can be influenced by a range of climate, vegetation, time, basin morphometry, sediment transport process, and underlying bedrock resistance or infiltration capacity influencing erosion rates (Oguchi, 1997; Sangireddy et al., 2016). It is slightly higher for the area encompassed by the southern pluton (2.59 km⁻¹) versus that of the northern and eastern plutons (2.39 km⁻¹). Assuming all factors but sediment transport processes and rock permeability and resistance to erosion are generally constant across the landscape of the adjacent plutons, the slightly lower value of the northern and eastern plutons is consistent with either (1) the headwater areas with smaller drainage basins being more dominated by threshold landsliding rather than advective fluvial processes (Oguchi, 1997) or (2) more permeable or erosionally resistant rocks. Alternatively, the higher D_d values of the southern pluton are consistent with highest rates of erosion in steep threshold channels (Howard, 1997). Similarly, Hurst et al. (2013) concluded that shorter hillslope lengths, inversely related to D_d, correlate with higher erosion rates along the Middle Fork Feather River. The southern pluton has a much higher stream order (Fig. 10B), a proxy for stream discharge, despite roughly similar values of stream power (Equation 1; area multiplied by slope squared; Fig. 10C), considered to represent the potential for erosion of a river into rock. Although the potential for erosion via a higher stream order in the southern pluton is higher, the ruggedness is lower and the 100 m landscape relief is roughly similar (Figs. 9B, 9C), also perhaps indicating that the southern pluton rocks are more easily eroded and less able to support steep channel walls needed for deeply incised channels than the northern and eastern plutons.

If main-stem and tributary channels are in steady-state conditions of balanced uplift and erosion, all channels should be co-linear in the context of a chi plot depicting χ and channel bed elevation z. To identify potential transient signals of uplift and erosion, we compare chi values (Equation 2) of the channel network within and across the Merrimac pluton boundaries (Fig. 12). Similar to the single profile depicted in Figure 11, the distribution of upstream distance as a function of bed elevation for all channels in the watershed reveals pronounced inflections (Fig. 12A). Viewing this channel network as a chi plot (Fig. 12B) highlights these inflections that deviate from the linear expectation at steady-state conditions where uplift (U) equals erosion (K). The z-intercept of the line is the elevation at x_b (the downstream terminus of the profile), and the dimensionless slope is (where A_o is the reference drainage area of 250 m², and m and n are positive exponents as described for Equation 1). Because (chi steepness index) provides a means of constraining changes proportional to U/K, it is evident that changes to U and/or K correlate with pluton boundaries (Fig. 12B). Segments of the transformed profiles are systematically steeper within the lower-elevation exposures of the southern pluton and within metamorphic rocks downstream (south) of the southern pluton contact (350-m elevation). This steepness indicates either (1) an increase in uplift rate, (2) a headward-migrating wave of erosion, or (3) a reduction in erodibility in the metamorphic rocks and lower-elevation rock of the southern pluton. In contrast, rock of the southern pluton at higher elevations (∼550–700 m) correlates with very shallow inclinations of the chi steepness index, consistent with either lower uplift rates or more readily erodible material. At higher channel bed elevations occupied by the northern and eastern plutons (∼700–1200 m elevation), the steeper- segment may correlate with less-erodible granitic rock, higher uplift rates, or erosion more dominated by mass-wasting transport processes near the drainage divide. We assume that discharge and sediment supply are secondary influences at this watershed scale. Furthermore, knickpoints, located at sharp convexities in channel profiles identified by TopoToolbox, tend to cluster with inflections of chi steepness index (Figs. 9E, 9F, 12A, 12B). Although transient erosional features, such as knickpoints, originating from a common source should plot at the same value of χ in all affected channels, Figure 12 reveals that they are concentrated within the northern pluton and within tributaries, not the main-stem channel. Although the watershed that encompasses the Merrimac plutons is small (Figs. 9 and 12), the uplift field may not be spatially uniform. Over the 15 km length of the basin, the influence of regional westward tilting may be present (Fig. 11), but the wavelength of the chi-space channel deviations does not match what would be expected from a fixed regional tilt. Rather, the accelerated erosion downstream, adjacent to the Feather River, has migrated, or is migrating, into dissected uplands, resulting in broad concave-up channel profiles downstream of knickpoints. The largest knickpoints (dz >80 m) (Fig. 9F) are located near the lower elevation contact between the southern pluton and metamorphic rocks upstream of the basin outlet and at elevations above the contact between the two plutons, consistent with a wave of erosion emanating from the southern pluton but slowing or stalling within potentially more resistant rock of the northern pluton.

On a regional scale, large fluvial canyons dominate the ruggedness and DEV_max images (Fig. 13). Using cosmogenic ¹⁰Be radionuclide concentrations in fluvial sediment, Riebe et al. (2000) documented higher erosion rates for catchments traversed by canyons, such as the Bald Rock and Grizzly plutons. Hurst et al. (2012) arrived at similar denudation rates of generally hundreds of millimeters per thousand years for samples from the northeastern part of the Bald Rock pluton. Away from these canyons, the plutons tend to be of low relief and lower ruggedness and characterized by an order of magnitude lower denudation rates (Riebe et al., 2000; Hurst et al., 2012). In particular, the Merrimac, Concow, and Granite Basin plutons, as well as parts of the larger Grizzly and Bald Rock plutons, are ringed by higher ruggedness and DEV_max values. This may reflect differing erodibility of the plutons relative to that of the generally fine-grained contact aureole and/or metamorphic country rock.

Compositional, geophysical, geochronologic, and geomorphic analyses support the presence of at least two Merrimac plutons. Despite the absence of exposed contacts, a third eastern pluton is consistent with magnetic and foliation patterns, the presence of a topographic break (drainage) that could reflect less-resistant rocks between the northern and eastern plutons, and higher rock densities along the southern margin of the eastern pluton. Although the eastern pluton has not been directly dated, the 142 Ma age of Irwin and Wooden (2001; Fig. 3A) for a sample of the northern pluton near the contact and its similar geophysical expression suggest that the eastern pluton may be the same age as or slighter older than the northern pluton.

Density measurements within the eastern pluton (Fig. 4C) show a gradient from high densities to lower densities from southwest to northeast, consistent with compositional data that indicate diorite to the southwest and tonalite to the northeast. These density and compositional patterns, along with weak foliation and cumulate textures (Fig. S2 [^{footnote 1}]), can be interpreted as a differentiation sequence from diorite at the bottom (or southern) part of the eastern pluton to tonalite at the top (northern part of eastern pluton), or perhaps as emplacement of compositionally diverse magmas adjacent to each other. Such an orientation of compositions suggests a regional top orientation to the northeast, as has been documented by Guglielmo (1993b) based on compositional asymmetry and a slight east to northeast plunge of the poles of magmatic foliation planes for the other two plutons.

Geophysical data define a distinct northeast-dipping internal contact between the more felsic southern pluton and the more mafic northern and eastern plutons; this contact is aligned with regional structural trends, especially a fault mapped to the southeast that is nearly a continuation of the contact (Fig. 3). The contact between the northern and southern plutons is somewhat more linear than might be expected from the interpolated magmatic foliation patterns (Fig. 3B); its orientation parallel to the Big Bend fault suggests that structure within the country rock influenced pluton emplacement. Guglielmo (1993a) posited that the curved contact based on his interpolated foliation pattern would suggest the northern pluton intruded the southern pluton, and thus would be younger. However, in light of the geochronologic data and refinement of contacts between the most mafic eastern pluton and slightly less mafic northern pluton, we posit instead that the southern pluton partially molded itself to the already crystallized contacts of the northern and eastern plutons. The absence of solid-state deformation within the plutons and along the contact between the three main plutons (Guglielmo, 1993a, 1993b), combined with wall rocks wrapping around the plutons (Hietanen, 1976) and brecciation and stoping along discordant contacts between the plutons and the wall rocks, still suggests forceful rather than passive emplacement. However, regional structure likely provided the pathway for that emplacement, particularly where dikes extend ∼5 km into the country rock along and parallel to the westward projection of the contact between the northern and southern plutons (Hietanen, 1951, p. 593), which in turn is parallel to the structures in the country rock.

Interestingly, the positive correlation between Eu/Eu* and Hf or U for zircons from the northern and southern Merrimac plutons (Fig. 8) suggests the zircons grew from different granitic melts from a compositional spectrum related by mixing. Given the apparent several-million-year difference in emplacement age between the northern and southern Merrimac plutons and that the crystallization history of individual plutonic intrusions is typically several hundreds of thousands of years (e.g., Matzel et al., 2006; Schoene et al., 2012; Ratschbacher et al., 2018), this spectrum of melts parental to the zircons likely reflects the mixing within a deeper-crustal source zone rather than shallow-level mixing between the magmas responsible for the plutons (e.g., Gray et al., 2008). Granitoids composing the Sierra Nevada batholith are considered to be generated in the lower crust (Ratajeski et al., 2005; Sisson et al., 2005), reflecting mixtures between mantle-derived and crustal components (e.g., Gray et al., 2008; Lackey et al., 2008). Correlations between trace element ratios, e.g., Th/U versus Nd/Yb, point to crystallization of accessory minerals such as allanite or titanite that have high Th/U and light-REE/heavy-REE ratios, likely within a granitic crystal mush (e.g., Deering et al., 2016; Reid and Vazquez, 2017). Based on the Ti-in-zircon geothermometer (Ferry and Watson, 2007), Ti concentrations suggest crystallization temperatures of ∼750–650 °C, with slightly lower Ti and generally slightly higher Hf concentrations in zircons from the southern pluton (Fig. 8). Thus, our interpretation is that the granitoid magmas that produced the two main plutons reflect mixing at lower crustal levels, but final emplacement of these magmas was at depths of 5–7 km (McLeod, 1990) and was influenced by pre-existing crustal structure forming a boundary between the two sets of plutons. A biotite K-Ar date (Saucedo, 1992, recalculated from Evernden and Kistler, 1970) indicates that the northern pluton cooled to temperatures of ∼300 °C by 129 Ma, which is reasonable given the size and emplacement depth of the Merrimac plutons and cooling rates for plutons elsewhere in the Sierra Nevada (Saleeby et al., 1989b).

Although pre-existing structure or map-scale lithologic differences may have localized the site of intrusion and the position of the contact between the main southern and northern Merrimac plutons, the nature of such controlling structure is cryptic. The magnetic anomalies aligned along the contact extend away from the plutons in either direction, indicative perhaps that the modest compositional difference between the plutons may be influenced by the country rock, even though the country rock surrounding both is quite similar (Fig. 3A). Hietanen (1976) mapped a fault between the Merrimac and Cascade plutons (Fig. 3) that is along the southeastward projection of the Merrimac plutons contact and appears to separate higher magnetic values to the northeast from lower values to the southwest (Fig. 5). This magnetic gradient is well expressed, especially in the residual magnetic field, but is subtly shown in the residual gravity field (Fig. 2). Within the Cascade pluton, the expression of this structure weakens considerably but could curve southward within the pluton or extend farther east. This structure likely does not extend beyond the eastern margin of the pluton, where prominent magnetic anomalies are sourced by ultramafic rocks at the base of the Slate Creek terrane and bounded by a thrust. To the west-northwest of the Merrimac plutons, projection of the contact between the southern and northern plutons coincides with the northeastern margin of a belt of magnetic highs, merging with the Camel Peak fault which curves toward the Quaternary Cohasset Ridge and Beaver Creek faults.

Wakabayashi and Sawyer (2000) showed a north-striking fault (the Cascade fault) within the Cascade pluton with an estimated slip rate of 0.007 mm/yr and no appreciable offset of the pluton margin (very little long-term offset), as well as one that lies along the eastern margin of the pluton that may be a strand of the Grass Valley fault zone with a slip rate of 0.05 mm/yr. It is possible that the Mesozoic structure marked by the Merrimac plutons boundary is serving to connect the north-striking late Cenozoic faults in and near the Cascade pluton to the Quaternary Cohasset Ridge and Beaver Creek faults in the northwestern part of the map area. The timing of this proposed fault connection is not known, given the temporal gap in the rock record between the Cretaceous and Miocene and slow rates of faulting during the Neogene.

The geophysical boundaries also coincide with topographic features. A south-facing abrupt slope break separates the Merrimac plutons, with the northern and eastern plutons being of higher elevation and steeper slope than the southern pluton. The east-southeast trend of the slope break within the pluton and along the contact continues beyond the pluton margins in the topography, subtly expressed in the ruggedness and DEV_max grids, as short, linear valley-and-ridge pairs that align with faults (Fig. 13). Trends of these segments curve to a more northerly orientation toward the Cohasset Ridge and Beaver Creek faults. To the east-southeast, similar-trending segments extend midway across the Cascade pluton before north-northeast–trending ridges and valleys dominate the terrain, as expressed by the DEV_max grid, to the east of the north-south–trending Cascade fault. The ridges are capped by Lovejoy Basalt, which is characterized by smooth topography, such as exemplified at Table Mountain. We interpret the general colocation of geophysical and topographic features that link the Quaternary faults to the northwest with late Cenozoic faults in the southeast as reflecting differences in erodibility along the pre-existing structure. This pre-existing structure likely influenced the emplacement of the Merrimac plutons, leading to a linear contact between the northern and southern plutons. The compositional differences between the plutons also affected landscape erosion as evidenced by knickpoints associated with headward erosion of longitudinal channel profiles clustering north of the pluton boundary. We speculate that crustal stresses associated with regional late Cenozoic tilting may have possibly reactivated this contact and/or fault, generating differential vertical uplift along this zone. Additionally, regional tilting is consistent with the generation of a wave of valley incision resulting in channel profile inflections and large-magnitude knickpoints north of the contact.

The geomorphic expression of deep valleys of the Feather River incising into a west-tilted relict upland surface is consistent with a wave of propagating incision initiated by accelerated uplift ca. 3 Ma (Wakabayashi, 2013). Of note is the bracketing of the deeply incised Middle Fork and North Fork Feather River around the boundaries of the Merrimac plutons. Gabet (2014) recognized this drainage network pattern, with high-relief, steep reaches through resistant metamorphic rock down-valley of relatively low-relief, higher-elevation granitic reaches, arguing that it could also be consistent with multiple generations of incision evolved over this lithological template. Because this drainage system clearly post-dates pluton emplacement, we speculate that the course of the canyons around these plutons (highly visible in the DEV_max of Fig. 13) may result from lithologic influence on the development of the drainage system; we suggest that country rock located outside of contact metamorphism associated with pluton emplacement is more easily eroded than plutonic rock and its metamorphic halo. However, it must be noted that the Middle Fork and North Fork Feather River extend across the Bald Rock and Grizzly plutons, respectively, although these forks may be along internal pluton contacts, as suggested by variations in the residual gravity and magnetic fields (Figs. 2F, 2G). The southeastern lobe of the Grizzly pluton and southeastern lobe of the Bald Rock pluton are denser than the rest of the respective plutons, as highlighted by residual gravity variations (Fig. 2F), suggesting that these areas are underlain by separate, more mafic plutons, with the contact between plutons being more erodible than the surrounding, more felsic plutonic rock. The presence of longitudinal-profile knickpoints identified in the chi analysis, however, indicates that the landscape has not yet fully equilibrated to the accelerated uplift starting in the late Neogene (Wakabayashi, 2013), with areas up-valley from the highest knickpoints likely decoupled from erosion in the deeply incised Feather River valleys. For context, cosmogenic radionuclide studies constrain basin-averaged erosion rates of ∼20 mm/k.y. for the low-relief uplands to ∼250 mm/k.y. in the incised canyons (Riebe et al. 2000; Hurst et al., 2012).

In contrast to the Merrimac plutons, the Bucks Lake pluton is a single-age composite pluton with a pyroxene diorite core surrounded by quartz diorite in a clear concentric fashion (Hietanen, 1973). This relationship is supported by the dense and magnetic core imaged by the geophysical data (Fig. 2), although at face value the zircon trace element chemistry data indicate the pluton rim is less evolved than its center (Fig. 8, Hf and Eu/Eu* values). We note that the shape of the pluton and the internal compositional changes does not appear to be aligned with the regional northwest-trending structural grain. The Bucks Lake pluton intrudes Calaveras Complex accretionary rocks; perhaps this sequence is too finely interlayered, homogeneous, or lacking in significant throughgoing faults to influence pluton emplacement, unlike the Merrimac plutons.

An interesting outcome of the geochronologic data is that the Granite Basin pluton (mapped area of 18 km²) may have been constructed within a relatively short period of time, perhaps within 200,000 yr, whereas the larger Cascade and Bucks Lake plutons (exposed areas of 130 and 124 km², respectively) intruded during a timespan of ∼1 m.y. The Bald Rock pluton, with an area of 190 km² and whose depth extent is at least twice that of the Merrimac plutons, may have been constructed over the span of 8 m.y. or consist of multiple plutons intruded over a long time period. Taken at face value, these data support the hypothesis that smaller plutons are constructed during short time frames (de Saint Blanquat et al., 2011). More precise dating of these plutons using isotope-dilution TIMS would resolve timing of pluton construction and its relationship to pluton size.

Integration of geologic, geophysical, geochronologic, and geomorphic analyses confirms and refines division of the Early Cretaceous Merrimac pluton (originally mapped a single large pluton) into at least two main plutons. The northern pluton can be separated into a third pluton along its eastern margin on the basis of foliation, rock property, and geophysical data. This third eastern pluton is of unknown age, but likely older.

Trace element compositions of zircons suggest that the northern and southern Merrimac plutons are the product of mixing of two source magmas rather than fractionation. Because the plutons were intruded several million years apart, this mixing likely occurred at depth in the lower crust. The eastern pluton shows the textural and compositional signs of differentiation after emplacement. Emplacement of the plutons in the shallow to middle crust appears to have been influenced by earlier Mesozoic structure, given the alignment of the contact between the southern pluton and the northern and eastern plutons with Jurassic regional structural grain. Projection of regional structure from geophysical anomalies allows linking of Quaternary faults to the northwest and southeast through the contact between the northern and southern Merrimac plutons. Furthermore, the contact between the two main plutons lies at the base of a nearly 400 m difference in elevation, as supported by the identification of large knickpoints in longitudinal channel profiles. A chi analysis reveals that the southern pluton has a generally lower steepness index than the northern pluton, indicating that the northern pluton may be more resistant to erosion with large-magnitude, longitudinal channel profile knickpoints just upstream of the contact between the plutons. Although much of this difference in elevation can be attributed to differences in erosional susceptibility caused by lithologic variations, we speculate that there may also be a Quaternary component of the elevation difference.

We thank Russell Graymer (USGS) for a very thorough review of an earlier version of the manuscript and acknowledge the National Geologic Mapping Program for financial support. John Wakabayashi and an anonymous reviewer helped sharpen the presentation. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

ERRATUM: Influence of pre-existing structure on pluton emplacement and geomorphology: The Merrimac plutons, northern Sierra Nevada, California, USA

V.E. Langenheim, J.A. Vazquez, K.M. Schmidt, G. Guglielmo, Jr., and D.S. Sweetkind

ORIGINAL ARTICLE: 2021, v. 17, no. 2, p. 455–478, https://doi.org/10.1130/GES02281.1. First published 19 January 2021.

ERRATUM PUBLICATION: 2021, v. 17, no. 2, p. 478. First published 12 March 2021.

When the pre-issue publication version of this paper was first published online (19 January 2021), Figure 3 appeared incorrectly. The foliation layer did not display when the corresponding radio button, “B) Foliations,” was selected. The figure and radio button have now been corrected, and the correct Figure 3 now appears in the final published version of the paper to which this erratum is appended.