Abstract

A 250-km-long broadband and long-period (0.01–20,000 s periods) magnetotelluric (MT) line along latitude 40½°N spans from the California-Nevada border at longitude 120°W across the buried extension of the northern Sierra Nevada, the southern Cascades arc, and the subducting Gorda plate to longitude 123°W. The resulting resistivity cross section reveals conductors at the locations of dewatering of the subducting metasomatized crust at 100 km depth, partial melting of the mantle wedge at 40–60 km depth, and melting in the crust shallower than 40 km. The conductor at 100 km is too conductive (8 S/m) to result solely from magma; very conductive fluids must be present either separately or incorporated into hydrous melts. A melt fraction of 7% is estimated for the mantle wedge (average conductivity is 0.06 S/m). Sites of Pleistocene and Holocene volcanism at Lassen Peak, Hat Creek, and Poison Lake are closely associated with crustal conductors inferred to be shallow magma. The MT study also reveals presumably basaltic magma at 40–50 km depth in the upper mantle where the eastern end of the profile extends into the Basin and Range province.

INTRODUCTION

The southern Sierra Nevada is the site of recent lithospheric removal, with crustal rocks directly underlain by mantle asthenosphere (Wernicke et al., 1996; Ducea and Saleeby, 1996). Seismic refraction data (Ruppert et al., 1998), seismic receiver functions (Jones and Phinney, 1998), and magnetotelluric (MT) measurements (Park et al., 1996) showed the absence of mantle lithosphere at depths as shallow as 35 km and partial melting in the uppermost asthenosphere beneath the Southern Sierra Continental Dynamics (SSCD) line (Fig. 1). The presence of garnet pyroxenite in xenoliths in 11 Ma basaltic flows led Ducea and Saleeby (1998) to conclude that a mantle lithosphere at least 100 km thick was present in the middle Miocene. The disappearance of garnet pyroxenite and an increase in temperatures estimated from thermobarometry of xenoliths in 3.5 Ma flows was cited by these authors as evidence that this lithosphere had been replaced by asthenosphere sometime between 11 Ma and 3.5 Ma. LePourhiet et al. (2006) suggested that the lithosphere was weakened by heating associated with Basin and Range extension to the east and had delaminated. This lithosphere flowed west and formed the high wave-speed Isabella anomaly (Jones et al., 1994). However, it is unclear whether thin modern lithosphere is restricted to the southern end or is a common feature of the entire range.

Recent work farther north with seismic receiver functions (Frassetto et al., 2011) and MT data (Ostos and Park, 2012) shows a similar picture for the lithosphere beneath Yosemite National Park (YNP, Fig. 1) in the central Sierra Nevada. A thick crustal root is absent beneath the range, with Mohorovičić discontinuity (Moho) at depths of 35–40 km beneath the highest elevations and deepening to 50–55 km beneath the western foothills (Frassetto et al., 2011). Ostos and Park (2012) concluded from low resistivities (10–30 ohm-m) below Moho that buoyant asthenosphere extends west of longitude 120°W (Fig. 1).

The Lassen MT line (LASS) is the northernmost line in the study of the Sierra Nevada lithosphere (Fig. 1). Unlike the SSCD and YNP profiles, the LASS profile lies in an area of great geologic complexity. Three tectonic plates (North America, Pacific, and Gorda) meet at the Mendocino triple junction ∼150 km west of this profile (Fig. 2). Active subduction of the Gorda plate results in magmatic activity comprising the southernmost Cascade arc, and our profile lies north of the southern edge of this plate (Verdonck and Zandt, 1994). The LASS profile passes just north of the southernmost active volcano (Lassen Peak) in the Cascade volcanic arc and crosses Pliocene and Pleistocene volcanic fields to the east (Fig. 3). Walker Lane (Fig. 1), a dextral shear zone that accommodates ∼20% of the Pacific–North American plate motion, may extend northward into the eastern portion of the LASS profile, but its associated faults apparently terminate southeast of the profile (Bennett et al., 2003; Wesnousky, 2005). However, the northern California shear zone (Wesnousky, 2005) east of Lassen Peak may represent a zone of transtensional shear. The profile crosses a region where the pre-Cenozoic structures of the northern Sierra Nevada appear to be offset westward in the Klamath Mountains to the north (Fig. 1). Dickinson (2008) proposes that this offset is the result of northwestward displacement of the Klamath Mountains along a forearc normal fault system at the northern end of the Great Valley. Finally, the profile also crosses the northern end of the Redding seismic anomaly (Benz et al., 1992), a high wave-speed body that is oriented vertically beneath the northern end of the Great Valley (Fig. 1). We thus expect the resistivity image along this profile to be more complicated than the southern profiles because of the interplay between subduction, volcanism, and possible lithospheric delamination.

MT METHOD AND DATA

The Lassen MT profile spans 250 km from the southern end of the Klamath Mountains to the Nevada border, passing 20 km north of Lassen Peak (Fig. 2). The profile consists of 14 stations with collocated short-period (0.001 s–300 s) MT-24 and long-period (30 s–10,000 s) Narod Intelligent Magnetotelluric Systems (NIMS) instruments from the Electromagnetic Studies of Continents (EMSOC) instrument pool. The MT method consists of measuring horizontal electric fields and both horizontal and vertical magnetic fields that result from natural oscillations of the Earth’s primary magnetic field (Vozoff, 1991). Because these passive measurements can contain both signal and noise, data are recorded at two sites simultaneously. We then use cross-correlation to compare signals between the sites in order to separate coherent signal from incoherent noise. The results of this analysis are frequency-dependent transfer functions between the horizontal electric and magnetic fields and between the vertical and horizontal magnetic fields. Larsen et al.’s (1996) robust processing technique was used for this analysis.

The transfer function between electric and magnetic fields is a complex 2 × 2 tensor (Z) called the impedance tensor (Vozoff, 1991). This tensor’s two principal axes are usually orthogonal and have both magnitude and phase. The magnitude is scaled by frequency and other physical constants to an apparent resistivity that would be the true resistivity if the Earth were homogeneous. The apparent resistivity and phase are then plotted versus period (the inverse of frequency) at a single site to show how the transfer function varies with depth and/or distance from station (Fig. 4). Sites with little noise have small error bars and smooth curves (site 208), while noisier sites show more scattered results and larger errors (site 213, Fig. 4). Because of the nature of electromagnetic induction, fields with longer periods penetrate deeper and/or farther away. A multi-period plot of transfer functions (called a sounding curve) thus shows the variation of resistivity with depth and/or distance.

Data were recorded with two different sets of instruments sensitive to different frequency bands, and then merged into a single set of responses. Because we used the same set of electric field dipoles for both the MT24 and the NIMS, we had no MT static shifts between the two sets of responses. Examples of merged station responses are shown in Figure 4.

The two principal axes decouple completely into two orthogonal modes if the fields are measured over one-dimensional (1-D) or two-dimensional (2-D) structures (Vozoff, 1991). Over 2-D structures, these modes may be found by rotating the tensor coordinates to a direction either perpendicular or parallel to structural strike. The transverse magnetic (TM) mode links the magnetic field parallel to strike (Hx) with the electric field perpendicular to strike (Ey) and the vertical electric field (Ez). The TM principal axis is perpendicular to strike. The transverse electric (TE) mode couples Ex, Hy, and Hz, and has a principal axis that is parallel to strike. We rotated our data to principal directions that maximized these two modes and minimized coupling between them. Principal directions between stations may still differ by 90° because the rotation could make Zxy ( = Ex/Hy) or Zyx( = Ey/Hx) larger. Principal axes at periods longer than 10 s ranged between N45W and due north (Fig. 2). We did not see the gradual shift in orientation from NNW to north observed by Park et al. (1996) and Ostos and Park (2012) from west to east. Instead, site groups usually had one site aligned due north and the others averaging N20W (Fig. 2).

The impedance tensor relates horizontal electric fields to horizontal magnetic fields, but part of the Earth’s response includes a vertical magnetic field. The transfer function relating vertical to horizontal magnetic fields is called the tipper (Vozoff, 1991), and it is also complex and frequency dependent. We present the tipper transfer function as a set of induction arrows for each site (Fig. 2). Because an induction arrow of unit length means a vertical magnetic field as large as the horizontal ones, arrows usually have lengths of less than 1. In a layered Earth, there are no induction arrows (all fields are horizontal), so lengths greater than 0 imply multidimensional structure. If all of the induction arrows are parallel or antiparallel, then the structure is ∼2-D and oriented orthogonal to the directions of the induction arrows. Our induction arrows are computed so that the real components point away from good conductors (Wiese, 1962). Real induction arrows point perpendicular to the average principal axis orientation and away from the electrically conductive Pacific Ocean in the west (Fig. 2). They gradually rotate southeast at the eastern end of our profile, indicating a conductor north of the profile. Most of the arrows are aligned approximately parallel to the profile and all of them are small, however (Fig. 2), so the fields behave as if the structure were 2-D.

Induction arrows are affected by shallower structure at shorter periods. Our data consisted of measurements made with two different MT systems, and only the NIMS system recorded the vertical magnetic fields needed to compute induction arrows. Thus, we only have induction arrows for periods longer than 20 s. We show induction arrows for 10,000 s (Fig. 2) because they are representative of periods longer than 300 s. Induction arrows for periods between 30 s and 300 s show more directional variability but still tend to be smaller than 0.1. Presumably, this is because of the more variable volcanic crustal structure.

Regional structure in California south of the LASS profile is ∼2-D with the Coast Ranges, Great Valley, and Sierra Nevada roughly parallel to one another through most of the state (Fig. 1). Most of the young faults east of site 209 also trend parallel to the regional structure (Fig. 3). While this regional structure changes north of the profile, directions of principal axes and induction arrows remain approximately parallel and perpendicular, respectively, to the regional geological strike of N25W found to the south (Fig. 2). This uniformity of geoelectrical strike is consistent with Dickinson’s (2008) inference that the Klamath Mountains are offset down to the northwest from the pre-Cenozoic Sierran arc structure. The pre-Cenozoic rocks may continue northward under the volcanic rocks from their exposures in the northernmost Sierra Nevada, or the fault structures (Fig. 3) could be dominating the MT response. The induction arrows aligned southeast are small, indicating that 3-D effects are minimal. For these reasons, 2-D modeling and inversion are appropriate for this study. While a thorough analysis of our data would theoretically require 3-D inversion, Wannamaker (1999) has shown that field behavior from 3-D bodies is often sufficiently similar to that from 2-D structures with the same cross section. With this restriction to 2-D modeling, all responses need to be rotated to a common strike for 2-D modeling. We chose N20W as an average between tensor rotations and directions of real induction arrows. The modeling also requires that we identify the TE and TM modes, so the component parallel to strike (Zxy) is the TE mode, and the component perpendicular to strike (Zyx) is the TM mode.

Presentations of sounding curves for each site along the profile are cumbersome, so data are plotted instead in pseudosections (Vozoff, 1991). One component (apparent resistivity or phase for one of the two modes) is chosen from each site and a section of that component is made by plotting those values directly beneath the site versus period. Values of one component for all periods and all sites are then contoured to make a section of horizontal position versus period. These are called pseudosections because the vertical axis is period and not depth (Fig. 5).

Pseudosections of TE and TM mode apparent resistivities (Fig. 5) show a shallow, moderately resistive (∼300 ohm-m) surface layer presumably related to Cenozoic volcanic flows (Figs. 1 and 3). East of site 209, apparent resistivities decrease to minima at ∼1 s periods and then rise to maxima between 10 s and 100 s before decreasing again at the longest periods. The TE mode apparent resistivity section has a prominent low centered beneath site 206 at longer periods, perhaps related to electrically conductive magma, fluids, or to the 3-D effect of truncation of a conductor along strike. West of site 209, regions of higher apparent resistivity are seen between sites 210 and 213 (Fig. 5). Unlike the other pseudosections across the Sierra Nevada (Park et al., 1996; Ostos and Park, 2012), we see no prominent high apparent resistivity zone where the profile crosses the northward projection of the Sierra Nevada batholith. The lower apparent resistivities seen at these stations are comparable values in the pre-batholithic metamorphic rocks west of Yosemite, however (Ostos and Park, 2012).

Qualitative examination of pseudosections can reveal structure, but numerical modeling and inversion are needed in order to establish depths and locations of anomalous bodies. This is because period, not depth, is plotted in a pseudosection. The depth of penetration in a layered structure is approximately proportional to √[(apparent resistivity) × (period)], so two sites at the same period may be sensing very different depths. At a period of 1 s, electromagnetic waves at site 213 may be affected by structure five times as deep as waves at site 212 (Fig. 5). Structure off to the side of a site can also affect its response. For example, Mackie et al. (1996) showed that the electrically conductive Pacific Ocean affected fields at long periods in Death Valley over 200 km from the ocean.

MT MODELING

We used 2-D inversion (Rodi and Mackie, 2001) to transform the MT responses into a cross section of electrical resistivity. This method assumes that apparent resistivities and phases from the TE and TM modes are available, as well as the tipper for the vertical magnetic field. However, numerous 3-D modeling studies (e.g., Wannamaker, 1999; Park, 1985) have shown that the TE mode apparent resistivity is particularly affected by truncation of elongated bodies along strike. While we began our inversion with the same weights for all MT components, we downweighted the TE mode apparent resistivity in order to achieve good fits to the rest of the data and for the inversion to converge. Rodi and Mackie’s 2-D inversion (2001) weights data with the inverses of the data errors but also requires minimum error floors. Error floors for the TE mode phase and the TM mode phase and apparent resistivity were set to 5%, based on the scatter between values at adjacent frequency points. Error floors for the tipper were set to 0.03, or ∼10% of the observed responses. Note that the inversion uses the larger of the error floor and the data errors. Given that the data had very small errors (Fig. 4), the error floors effectively determined the weights in the inversion. Downweighting the TE mode apparent resistivity was achieved by setting its error floor to 10,000%. This high error floor meant that the inversion did not attempt to fit the TE mode apparent resistivity.

Although the 2-D inversion can be started with an initial model that is a homogeneous half-space, we have found that beginning with a layered model that contains the essential features of the sounding curves results in faster convergence. An initial layered model has an additional advantage that any lateral structure in the resulting model has been placed by the inversion and is thus needed to fit the data. We found a layered model by graphically averaging the determinant of the impedances at all sites and then performing a 1-D layered inversion on the composite sounding curve. Methods for optimizing this composite curve and the resulting 1-D layered model were less important for the 2-D inversion because we were only picking a starting model. We used the layered model only for structure deeper than 20 km. The starting model had 1000 ohm-m from the topographic surface to 20 km depth below sea level, 300 ohm-m from 20 km to 58 km, 100 ohm-m from 58 km to 121 km, 30 ohm-m from 121 km to 200 km, and 10 ohm-m below 200 km. The full model extended 1000 km to the west of site 214 and included Pacific Ocean bathymetry. The model consisted of 99 columns and 53 rows, with 16 of the rows dedicated to topography. Widths of adjacent columns did not vary by more than a factor of 1.5, as required by the inversion.

Rodi and Mackie’s 2-D inversion used Tikhonov regularization to minimize the model roughness represented by the Laplacian of the model and the data misfit represented by the l2 norm (Aster et al., 2005). This meant that the inversion was trying to fit the data with a model that had the least variation between adjacent blocks (i.e., a smooth model). Because the layers introduced by the starting model added model roughness at the interfaces, the model also tried to smooth these abrupt changes in resistivities. Layers remained if they were required by the data, however. Inclusion of the model roughness is theoretically designed to ensure that only those features required to fit the data are included in the model (e.g., Aster et al., 2005). The importance of model smoothness relative to data misfit in the inversion is controlled by a weighting factor, τ, applied to the model roughness. A large value of τ results in a smooth model with larger data misfits, while a small value produces a rough model with smaller misfits. The optimal value of τ is determined by running test inversions with different values of τ and seeking a value that yields the smallest combination of model roughness and data misfit. A value of 0.25 was optimal for our model and data, so all models shown and discussed in the paper were generated with this same weighting factor.

A smooth model that fit the data was found after 199 iterations (Fig. 6). A root mean square (RMS) error of 1.82 was achieved, meaning that the data were fit to within 1.82 times the data error. The resulting model has no artifacts from the starting layered model, but instead shows large lateral contrasts across the entire section (Fig. 6). Pseudosections of observed mode data and synthetic responses from the model are relatively similar, except for the TE mode apparent resistivity (Fig. 5). The misfit for the TE mode is expected because these data were not fit by the inversion. This model has a complex crust that comprises areas with resistivities over 1000 ohm-m interspersed between areas with resistivities less than 3 ohm-m. As will be discussed later, many of the low resistivities lie beneath areas of Quaternary volcanism. Unlike profiles to the south (Park et al., 1996; Ostos and Park, 2012), there is no large, contiguous resistor corresponding to a body of granitic rock. No large crustal conductor lies beneath Lassen Peak (LP, Fig. 6), but there is a mantle conductor with resistivity less than 1 ohm-m centered at a depth of 150 km. Another deep conductive body located beneath site 201 lies at a similar depth (Fig. 6). This conductor may extend eastward, but the structure beyond 40 km east of site 201 was constrained (dotted region, Fig. 6) because preliminary inversions placed large, very conductive structures in an area with no data control.

Comparison of this model to seismic images is initially perplexing. While there is good agreement between Moho estimated from receiver functions (Frassetto et al., 2011) and the base of the short wavelength variations in resistivity at depths less than 40 km (Fig. 6), the P–wave-speed tomography (Reeg, 2008) does not match either the Moho or features in the resistivity section. A high wave-speed body dipping eastward from beneath the San Joaquin valley (Fig. 7) appears electrically discontinuous with a resistive top and bottom and a conductive middle (Fig. 6). Initial sensitivity testing (Park and Ostos, 2008) showed that the upper resistor and middle conductor were required. Additionally, the eastern conductor and a mantle resistor separating the two conductors were required to fit the data (Fig. 6).

A possible explanation for this mismatch is that electrical resistivity is a property strongly influenced by small amounts (<1%) of fluids, melt, and metallic minerals and only minimally controlled by bulk mineralogy (e.g., Jones, 1992). Seismic wave speed is determined primarily by bulk composition and less so by small amounts of fluids and melt. The same small amount of melt or fluids changes the resistivity of a rock much more than it changes the wave speed. The mismatch between the resistivity image (Fig. 6) and seismic tomography may simply indicate that the methods are sensitive to different aspects of the physical and chemical state of the rocks.

MT inversion is an ill-posed mathematical problem (e.g., Rodi and Mackie, 2001). The model grid has over 5400 unknown block resistivities and fewer than 4000 observations of apparent resistivity, phase, and tipper, so models derived from inversion will be nonunique. Additionally, electrical geophysical methods suffer from another source of nonuniqueness due to the principle of equivalence. This principle states that these methods are sensitive only to the conductance (inverse of resistance) of a target body (Telford et al., 1990), and conductance is a product of conductivity and geometry. It can be shown that the sensitivity at a site to a resistivity block in the subsurface depends on the amount of electrical current from a fictitious dipole placed at that site (Madden and Mackie, 1989; McGillivray and Oldenburg, 1990). That current is, in turn, dependent on the product of the block’s conductivity (inverse of resistivity) and its cross-sectional area (McGillivray and Oldenburg, 1990). Because models are nonunique, the areal conductance (cross-sectional area * conductivity = cross-sectional area/resistivity) is a more robust parameter in the inversion than an individual block’s resistivity. The form of this sensitivity means that the MT measurements will be most sensitive to conductors and least sensitive to resistive bodies. The common geophysical approach to this problem of nonuniqueness is to constrain models with other geophysical observations. For example, interpretations of gravitational data often use seismic data to constrain interfaces.

We constrained the MT model with results from the seismic tomography (Fig. 7). Because the MT method is most sensitive to conductive bodies, we constrained the resistivities within the outline of the high wave-speed body. The high wave-speed body is a combination of the Redding anomaly at depths less than 100 km and the subducting slab at greater depths (Reeg, 2008; Benz et al., 1992). Because the Redding anomaly represents downwelling lithosphere similar to the Isabella anomaly in the southern Sierra Nevada (Reeg, 2008), we assumed that it should be resistive because the Isabella high wave-speed body was resistive (Park, 2004). The subducting slab should also be resistive because it is cooler than the surrounding mantle, so we made the entire high wave-speed body resistive. In effect, we used the seismic data to constrain a part of the model to which MT data have little sensitivity. We then ran a second inversion with the resistivity of the high wave-speed body fixed at 1100 ohm-m and found a new model (Fig. 8) after an additional 75 iterations. This new model has all of the conductive features of the unconstrained model (Fig. 6), but their locations have shifted. There is still a mantle conductor beneath Lassen Peak, but it is shifted slightly eastward (body B, Fig. 8). Other major features of the model include a conductive (3 ohm-m) mantle west of site 214 (body C), and eastern conductors both shallow (10–20 ohm-m, body D) and deep (10–20 ohm-m, body E) beneath sites 201–204 (Fig. 8). The overall RMS error for the new model was also 1.82. We do not show the new fits to the observed data because they are equivalent to those in Figure 5. Additional inversions starting with the model in Figure 8 allowed greater weights (error floors of 20%–100%) for the TE mode apparent resistivity but resulted in higher errors (RMS = 2.87) and degraded fits to the rest of the data. We thus rejected these models.

We favor the model in Figure 8 because it provides a better match to the seismic wave-speed model (Fig. 7). In addition to the constrained slab, the low wave-speed region enclosing bodies G and H is generally more conductive in Figure 8. A second, low wave-speed region beneath sites 201–204 is equally conductive in both models (D, Fig. 8). Instead of conductor B falling in the middle of the high wave-speed body, it now is centered at a depth of 100 km at the upper edge of the slab (Fig. 8).

SENSITIVITY TO MODEL FEATURES

The two models (Figs. 6 and 8) contain a number of conductive bodies in the crust and upper mantle, but we have no bounds on how conductive these must be in order to fit the data. While a simple way to find bounds would be to change the resistivity of a feature and then see how it changes the model response, MT responses to features are coupled, and a perturbation in one could lead to changes in another. To determine bounds, we instead fix the resistivity of a feature and then run the inversion again. We reject that resistivity value if no alternate model is found that fits the data equally well. We repeat this process with multiple values of resistivity for each target conductor. In this way, we can establish the existence of and range of resistivities for features in the crust and mantle.

ROBUSTNESS OF CRUSTAL STRUCTURE

Two models (Figs. 6 and 8) with different mantle structures fit the data equally well, but both have very similar features above 50 km. The unconstrained model has a conductive channel in the uppermost mantle (F, Fig. 6); this channel is shifted upward to the lowermost crust in the constrained model (F enclosed by dashed line at top of Gorda plate, Fig. 8). Both models have resistivities less than 30 ohm-m in the crust beneath Lassen volcanic center, Hat Creek, and Poison Lake (G, Figs. 6 and 8). Both have a resistive (>10,000 ohm-m) body beneath sites 206 and 207 (H, Figs. 6 and 8), and a pattern of three conductors between sites 201 and 206. The crust east of site 201 is resistive in both models.

Despite these similarities, it is still possible that some of these crustal conductors are not required to fit the data because the inversion has many degrees of freedom. As an additional test of the crustal structure, we fixed the crustal conductor beneath site 205 (Fig. 8) at 100 ohm-m and ran the inversion. The final model had misfits at both sites 205 and 206 in the TE phase at periods of 10–100 s and TM apparent resistivity at periods shorter than 10 s. RMS errors for site 205 were much higher (1.84 versus 1.64) and slightly higher for site 206 (1.57 versus 1.52). Significantly, the alternate model still had a conductor of similar size at the same depth but adjacent to the west side of the fixed region.

We compare the areal conductance of these two alternative conductors because it is a measure of how much they influence the responses at sites. The areal conductance can be computed for a complicated body made up of multiple blocks by summing up the contributions from individual blocks (Σ [ΔxΔz/ρ]). Additionally, the areal conductance is not strongly affected by arbitrary boundaries for the body because resistive blocks contribute little or nothing to the sum. The areal conductance of the crustal conductor in the preferred model (Fig. 8) is 46 S/m·km2, while that in the alternate model is 41 S/m·km2. We conclude from the larger errors and the fact that the inversion still required comparable areal conductance between sites 205 and 206 that the conductor beneath site 205 is required by the data. We infer from this test that the other crustal conductors are also needed. Based on the similarities of the crustal models in Figures 6 and 8 and on this test, we conclude that crustal features are robust and are minimally influenced by changes in the mantle structure.

SENSITIVITY TO MANTLE STRUCTURE

We next examine sensitivity in order to provide bounds on the resistivities of features in the mantle (Fig. 8). While the MT method is sensitive primarily to conductive features, we still wanted to set bounds on the resistivity of the high wave-speed body (body A, Fig. 8). By running multiple inversions with different fixed resistivities for the high wave-speed region, we determined that any resistivity greater than 300 ohm-m provides an equally good fit to the data. This value is comparable to that found for the Isabella high wave-speed anomaly in the southern Sierra Nevada (Park, 2004).

Conductor B has a significant effect on MT fields at the stations because it is so conductive. Providing an average conductivity is complicated because the cumulative effect of all the conductive blocks governs field behavior. However, the areal conductance also provides a way to estimate the average resistivity of a region composed of multiple blocks. Conductor B has an areal conductance of 8142 S/m·km2, with most of that (8020) from the blocks at a depth range of 94 to 104 km between sites 206 and 207 (Fig. 8). Dividing the areal conductance by the area of the blocks in this depth range yields an average resistivity of only 0.05 ohm-m. This value is extremely low, so our sensitivity testing determined by how much the areal conductance could be reduced. Reducing the areal conductance to 1456 S/m·km2 (average resistivity of 0.25 ohm-m) produced a model with conductors added west of site 214 at 100 km depth and a larger RMS error of 1.83. Misfits in the TE phase at sites 205, 207, 208, and 209 for periods longer than 1000 s were larger than the data uncertainties. We thus rejected the possibility that the areal conductance of body B could be as low as 1456 S/m·km2. We then reduced the areal conductance to 3267 S/m·km2 (average resistivity of 0.12 ohm-m) and obtained a model that fits as well as that in Figure 8. Thus, we conclude that conductor B can have an average resistivity no larger than 0.12 ohm-m.

While the preferred models in both Figures 6 and 8 placed conductors (resistivity <2 ohm-m) in the mantle west of site 214 (conductor C), we found that a model with this area replaced with 100 ohm-m fit the data with an overall RMS error of 1.83. No other changes were made to the model, so we conclude that our data have little sensitivity to mantle structure beneath the Pacific Ocean.

Both models place a moderately conductive (average resistivity of 16 ohm-m) body (conductor D) just below Moho in the eastern portion of the profile (Fig. 8). Conductor D has an areal conductance of 76 S/m·km2. Establishing bounds on its average resistivity was more complicated than the other tests because we found both upper and lower limits. We first reduced the conductance to 23 S/m·km2 (average resistivity of 50 ohm-m), and the resulting model had a larger RMS error of 1.83, as well as misfits that exceeded data errors in the TM mode at site 201 for periods between 3 s and 300 s. Decreasing the areal conductance further to 12 S/m·km2 (average resistivity of 100 ohm-m) led the inversion to replace it with a slightly deeper conductor. Conversely, increasing the areal conductance to 145 S/m·km2 (average resistivity of 8 ohm-m) again resulted in an RMS error of 1.83 and misfits at site 201 in excess of data errors. We conclude from this test that an average resistivity of 8–50 ohm-m is preferred for conductor D. Finally, conductor E is a small, deep conductor with an areal conductance of only 39 S/m·km2 (Fig. 8). Replacing conductor E with the background resistivity of 100 ohm-m resulted in an RMS error identical to the preferred model (Fig. 8), so we conclude that this conductor is not necessary to fit the model.

DISCUSSION

Joint interpretation of seismic and MT data led to a much better subsurface model than analysis of either data set alone. With the simple constraint that the high wave-speed body is also resistive, we have a model (Fig. 8) that makes more sense than did our initial result (Fig. 6) because the conductor B is now at the edge of the subducting slab rather than in its interior. The original goal of our project was to determine whether thin lithosphere produced by delamination was present this far north in the Sierra Nevada, but our results in Figure 8 are equivocal. Park (2004) and Ostos and Park (2012) both interpreted large contrasts between resistive crust and shallow conductive mantle as evidence of delamination. Our model shows shallow conductors just below Moho, but the only mantle feature similar to those in the previous studies is the broad conductor at 200 km extending from the east beneath the profile (Fig. 8). We conclude that our MT data do not show clear evidence of delamination at this latitude.

Our electrical data are dominated instead by magma generated by subduction. The model shows the magmatic system for Lassen volcanic center (LP, Fig. 8) and surroundings with greater detail than is seen in the tomographic section (Fig. 7). Tomographic images from P–wave-speed data show a low-velocity zone (LVZ) that is as much as 7% slow (Reeg, 2008; Benz et al., 1992). Benz et al. (1992) interpret this LVZ as a crustal magma source for Lassen volcanic center that extends into the uppermost mantle; the Reeg study cannot discern whether the LVZ has a crustal component because of the wide station spacing used in Reeg’s study (Jones, 2012, personal commun.). The next step is to compare the geophysical interpretation (Fig. 8) to geologic processes so the model can help constrain those processes, and vice versa.

The LASS profile crosses a broad swath of Quaternary volcanic rocks spanning from site 206 to site 210 (Fig. 3). This swath includes Lassen volcanic center, Hat Creek, and Poison Lake—all locations of Quaternary and/or Holocene volcanism (Feeley et al., 2008; Turrin et al., 2007; Muffler et al., 2011). Our electrical image shows that the LVZ encompasses three smaller crustal conductors (a, b, and c in Fig. 9), a crustal resistor (H), and a more diffuse upper mantle conductor beneath this broad swath. Two additional crustal conductors (d and e in Fig. 9) are not associated with a crustal LVZ but fall within an eastern region of Quaternary eruptions (Grose et al., 1992) and lie above a mantle LVZ. Resistor H exceeds 30,000 ohm-m and is similar to values seen in Sierran granitic rocks in Yosemite (Ostos and Park, 2012). We speculate that this resistor may be a small Mesozoic intrusive body beneath the surficial volcanic rocks.

The Quaternary volcanic rocks are quite varied, with calc-alkaline basalts and andesites from arc magmatism dominating but small amounts of low-potassium tholeiitic (LKOT) basalt scattered throughout (Clynne and Muffler, 2010). The LKOT basalts appear to be related to the Basin and Range province to the east. Lassen volcanic center is a long-lived site erupting compositions from basaltic andesite to rhyolite over the past 825 Ka (Clynne and Muffler, 2010). Based on petrology and geochemistry, rising mantle basalts heat the lower crust, melting and assimilating that crust to form andesites (Clynne and Muffler, 2010; Feeley et al., 2008). Additionally, minor amounts of felsic ash and rhyolite found at Lassen volcanic center (LP, Fig. 9) are inferred to result from melting of the lower crust with little to no mixing (Borg and Clynne, 1998). The resistivity image for the LVZ beneath Lassen volcanic center shows discrete conductive foci in the crust (a and b Fig. 9) as well as a diffuse upper mantle conductor spanning the region between conductors B and G (Fig. 8). We speculate that the conductive foci are melt reservoirs at depths of 15–30 km where mantle basalts reside, heating and melting the lower crust.

Multiple calc-alkaline basaltic eruptions between 100 Ka and 110 Ka at Poison Lake (PL, Fig. 9) in the Caribou volcanic field represent some of the most mafic arc magmas in the Lassen region (Clynne and Muffler, 2010). Muffler et al. (2011) argue that these basalts are not related to a common parent through crystal fractionation but instead represent different pulses of magma from the mantle. Exhibiting little or no crustal contamination, these basalts were stored only briefly in the lower crust before eruption (Muffler et al., 2011). Crustal conductors beneath Lassen volcanic center and Poison Lake (a and b in Fig. 9, respectively) must be discrete bodies in which magmatic processes differ. In the former (a in Fig. 9), melting of the lower crust is much more extensive than in the latter (b in Fig 9). A possible explanation for this difference is that the Caribou volcanic field may be a developing volcanic center with insufficient basaltic magma input yet to trigger crustal melting and/or assimilation (Clynne and Muffler, 2010).

Tholeiitic (LKOT) basalt forms the 24 Ka eruptions at Hat Creek (HC, Fig. 9; Turrin et al., 2007). These basalts also show little evidence of crustal melting and have reached the surface along normal faults (Clynne and Muffler, 2010). They are inferred to be related to Basin and Range extension, rather than arc magmatism (Clynne and Muffler, 2010). Basalts similar to these from Medicine Lake, also related to Basin and Range extension, appear to be generated at depths of 30–45 km (Bartels et al., 1991) and most likely below Moho at ∼40 km in our section (Fig. 8). We note that the LVZ lies beneath three active loci of volcanic activity (Lassen volcanic center, Poison Lake, and Hat Creek), but we see only two crustal conductors (a and b in Fig. 9) beneath these loci. A possible explanation is that the crustal conductors are only associated with the volcanic centers because the residence time of the basaltic melt for Hat Creek is too short.

The broader conductive region just below Moho is the proximal source of the basaltic magma rising into the crust (“basalt magma,” Fig. 9). These magmas have been generated from melting at the top of the Gorda plate and/or from melting throughout the lower part of the mantle wedge (Grove et al., 2006; Till et al., 2012). A deeper region of very high conductivity at 100 km depth and just above the seismically defined slab underlies all of the melt zones but cannot be due to melt alone. The average resistivity of conductor B (Fig. 8) is 0.12 ohm-m. Basaltic magma has a resistivity of ∼0.5 ohm-m at 1200 °C and ∼0.2 ohm-m at 1400 °C (Tyburczy and Waff, 1983). Even a region with 100% melt would be more resistive than observed, so other components must be making it so conductive. We propose that this region must have a substantial fluid content because saline fluids can be much more conductive (Reynard et al., 2011).

Conductor B at a depth of 100 km (Fig. 8) is so conductive that we must assess whether it is even plausible that fluids are responsible. While good conductors have been seen beneath volcanic arcs in other studies (e.g., Soyer and Unsworth, 2006; Brasse and Eydam, 2008; Jodicke et al., 2006), none have had as high an average conductivity of ∼8 S/m (from the inverse of average resistivity of 0.12 ohm-m for conductor B). Laboratory measurements and model calculations of the conductivity of dry rocks of the subducting plate and mantle wedge show values of less than 10–4 S/m (Reynard et al., 2011; Zhu et al., 1999), so the matrix does not contribute to this high conductivity. A fluid must be capable of forming an interconnected phase and must be conductive enough to explain the high average conductivity. Mibe et al. (1998) showed that water was capable of wetting forsterite at pressures from 3 to 5 GPa in the mantle. Our depth of 100 km for conductor B places it within this range of pressures, so it is possible for aqueous fluids to form an interconnected phase at this depth and channel electrical current.

The effect of a very conductive fluid phase (σfluid) on the average conductivity (σave) can be predicted from Hashin and Shtrikman’s (1962) upper bound given by: 
graphic
where ϕ is the fraction of fluid in the rock. For 1% fluid, an average conductivity of 8 S/m requires a fluid conductivity of 1200 S/m. If the fluid fraction is increased to 10%, then a fluid conductivity of 116 S/m is required. These values seem unrealistically high, but Reynard et al. (2011) have fluid conductivities as high as 100 S/m in their models of slab melts. Gaillard (2004) reported values near 100 S/m for hydrous rhyolitic melts and attributed the high conductivity to enhanced mobility of sodium due to the presence of water in the melt. We suggest that several percent of a very saline fluid or hydrous melt could explain conductor B.

A qualitative explanation for the distribution of conductors beneath the Lassen region involves dehydration of the subducting slab releasing fluids into the overlying mantle wedge (Grove et al., 2002; Grove et al., 2006; Till et al., 2012). The depth range of 94 to 104 km for conductor B falls at the lower end of the range for metamorphic dehydration of chlorite (Till et al., 2012). The first basaltic melt produced at the lower edge of the mantle wedge is hydrous (up to 28% water, Grove et al., 2006) but is only 2.5%–3% by volume. As this melt rises, it encounters progressively warmer mantle wedge material. Fluids released by the rising melt lower the melting temperature of the material, triggering additional melting of the wedge. Melting occurs under vapor undersaturated conditions (Till et al., 2012), so there is never any aqueous fluid, and all of the enhanced electrical conductivity is due to the basaltic melt. At depths of 50–70 km, melt fractions can exceed 10% in the hottest region of the mantle wedge (Grove et al., 2006). At shallower depths in the mantle wedge, temperatures decrease and no further melting occurs. Underplating of the basaltic magma occurs (“basaltic magma,” Fig. 9), with some of this magma erupting or rising to melt and/or assimilate the lower crust (crustal conductors, a and b in Fig. 9).

The average resistivity of the mantle wedge above conductor B and below the Moho is ∼10 ohm-m (Fig. 9). Temperatures of the basalt melts are estimated to be 1225–1275 °C (Anderson et al., 1982; Bartels et al., 1991), so a melt conductivity of ∼2 S/m is probable (Tyburczy and Waff, 1983). Using this melt conductivity, a melt fraction of 7% in the mantle wedge above conductor B is predicted from Equation (1). This value of melt fraction is within the range predicted by Grove et al. (2006) for the mantle wedge. Prediction of melt fractions in the crust is not attempted here because the rock, fluid, and magma conductivities are much more variable.

Our interpretation so far has focused on the Lassen region, but the resistivity section also contains crustal and mantle conductors beneath and east of Susanville (S, Fig. 9). Conductor D (Fig. 8) has an average resistivity of 16 ohm-m, and a melt fraction of 4% is predicted for this region. It is unlikely that aqueous fluids contribute much to conductor D because fluids would not form an interconnected phase at depths of 40–60 km (Mibe et al., 1998). Conductor D is likely the source of melt for the Pleistocene basalts (Grose et al., 1992) in the vicinity of sites 201–204 (c, d, and e in Fig. 9). These melts are most likely associated with Basin and Range extension (Guffanti et al., 1990) but could result from decompression melting of rising asthenosphere driven by counterflow in the mantle wedge (Grove et al., 2002). In every unconstrained inversion, conductive material was added at depths below 150 km and east of site 201. If the Basin and Range structure was constrained, then conductive material was added below and just east of site 201 (Fig. 6). A model with a conductive Basin and Range asthenosphere (Fig. 9) provided a better fit to the vertical magnetic field transfer function at site 201 than did the one in Figure 6. We conclude from these tests that an electrically conductive asthenosphere east of the profile is preferred by the MT data. Conductor D thus appears to be a western edge of the generally conductive mantle beneath the Basin and Range, so we prefer the model of Basin and Range extension instead of the counterflow model.

All models we ran placed conductors west of site 214 in the mantle at depths of 100 km. However, sensitivity tests showed that these conductors were not required by the models. Seismic refraction-reflection studies (Beaudoin et al., 1996) place the top of the Gorda plate at 30 km depth at a position 80 km west of site 214 and dipping 11° east (Fig. 8). The top of the Gorda plate thus lies at the base of the crustal resistor west of site 214; this crustal resistor corresponds to arc-related rocks of the Klamath terrane (Beaudoin et al., 1996). When we started a model with conductive sediments inserted at the top of the Gorda plate and the deeper mantle structure constrained to 100 ohm-m, the inversion removed the conductors. We therefore conclude that conductive rocks (sediments or metasomatized crust) are not required west of the profile and that such conductors are probably artifacts of the inversion.

CONCLUSIONS

A 14-site magnetotelluric profile has imaged the magmatic system associated with the subduction of the Gorda plate beneath the northern Sierra Nevada. This system is complex, with fluids from dehydration metamorphic reactions fluxing the mantle wedge and generating basaltic magma from the wedge. An average conductivity of 8 S/m for a narrow region just above the subducting plate at 100 km requires several volume percent of very conductive (>100 S/m) fluids. Melt fractions of less than 10% are needed to explain the conductivity of the mantle wedge beneath Lassen volcanic center. Unlike previous studies of subduction zones in the Cascades (Soyer and Unsworth, 2006) and the Andes (Brasse and Eydam, 2008) that did not find conductive anomalies beneath the volcanic arcs, our study shows that the conductors are directly associated with the foci of active volcanism. Within the active volcanic arc, we suggest that there is an arc-related source of basaltic magmatism. A second source of basalts is at the eastern end of the profile and likely results from the encroaching Basin and Range asthenosphere.

This project was supported by grant EAR0607349 from the National Science Foundation’s Continental Dynamics Program. Any opinions, findings, and conclusions or recommendations are those of the author and do not necessarily reflect the views of the National Science Foundation. Magnetotelluric instruments from the EMSOC instrument pool were used for this study. We thank the Bureau of Land Management, Lassen National Forest, W.M. Beatty and Associates, and numerous private landowners for permission to install sites. Thanks to Paula Chojnacki, Bryan Smith, and Adrian Sears for helping under the difficult field conditions during the 2008 California forest fire season. Special thanks to Craig Jones and Patrick Muffler for thorough reviews of earlier versions of this manuscript.