We resolved the architecture of the early Proterozoic Penokean orogen suture and late middle Proterozoic (Keweenawan) Midcontinent rift system magmatic overprint in east-central Minnesota and western Wisconsin through recovery and analysis of a legacy magnetotelluric (MT) data set. We digitized printed plots of off-diagonal MT and controlled-source audio (CSA) MT responses, including error intervals, to provide 22 soundings along a profile of ∼225 km length extending from north of Lake Mille Lacs, Minnesota, southeastward to Flambeau Ridge, Wisconsin. The MT data were inverted to a smoothed electrical resistivity structure using a two-dimensional finite-element, regularized Gauss-Newton algorithm emphasizing the transverse magnetic (TM) mode data subset. Our model reveals a major electrically conductive zone dipping moderately to the southeast for >50 km in the 5–35 km depth range, which marks the probable Penokean suture in easternmost Minnesota. We interpret the conductor to reflect a package of graphitized metasediments of the former Archean continental margin and near foreland zone, underthrust as the Penokean terrane collided with the Superior Province. The large-scale conductor is now hidden beneath mainly Yavapai-aged plutonic rocks of the East-Central Minnesota batholith. Below the axis of the later Midcontinent rift (St. Croix horst, subsequently), a compact resistive body ranging from 5 to 20 km deep overlies the large conductor. We interpret this resistor to be mafic volcanic and intrusive rocks of the Midcontinent rift event, which correlate spatially with a high Bouguer gravity anomaly similarly modeled. The rift here coincides with the lower-crustal reaches of the suture, but the specific influence of the suture on rift emplacement is unclear.

The northern midwestern region of the United States, and in particular the Lake Superior region, contains unique elements in the construction and evolution of North America. Here lie remains of the Penokean compressional orogeny, the oldest Proterozoic accretionary orogen along southern Laurentia (Van Schmus et al., 2007; Schulz and Cannon, 2007; Whitmeyer and Karlstrom, 2007). It was the forerunner of Laurentia’s southward growth from the Archean Superior core via island-arc and small terrane assembly (Fig. 1). In Michigan and Wisconsin to the northeast of our study, the Superior-Penokean suture is marked by the well-known Niagara fault zone. In Minnesota, where our study lies, the suture has been associated with the Malmo structural discontinuity, but it is considered to have been extensively disrupted and obscured by post-Penokean gneiss upwelling and plutonism, for example, by the McGrath gneiss dome and the East-Central Minnesota batholith (Schulz and Cannon, 2007; Holm et al., 2007; Chandler et al., 2007, 2008; Boerboom et al., 2011). The west end of our magnetotelluric data set coverage overlies these gneisses and the Isle and Warman granitic panel of the East-Central Minnesota batholith. The gneisses appear to be remobilized from the Archean Minnesota River Valley terrane, more preserved to the southwest, which in turn is separated from the Superior granitoid Wawa terrane by the Great Lakes tectonic zone (Fig. 1; Chandler et al., 2007; Bauer et al., 2011; Southwick, 2014).

The Keweenawan Midcontinent rift system was imposed upon the terrane suture in our study area at ca. 1100 Ma. This event nearly separated the ancient continent and imparted enormous volumes of mafic magma to various depths in the crust (Hinze et al., 1997; Nicholson et al., 1997; Vervoort et al., 2007; Stein et al., 2015), engendering numerous high-grade base and precious metal deposits (Miller and Nicholson, 2013). Northeastward from our transect area in Michigan, a strong clockwise rotation of Midcontinent rift system strike regionally appears to be primary and not subsequent to emplacement, suggesting influence by preexisting large-scale structures in the evolution of the Midcontinent rift system (Hnat et al., 2006). Similarly, major Penokean structures continuing eastward into Ontario, Canada, may have influenced rift orientation there (Manson and Halls, 1997; Riller et al., 1999). Recently, the Midcontinent rift system has been interpreted as the failed arm of a triple junction as seafloor spreading between Laurentia and Amazonia was established (Stein et al., 2015, 2016). The volume of magmatism appears to require mantle hotspot input beyond passive rifting (Stein et al., 2015, 2016).

Geophysical methods can illuminate deep crustal architecture in the face of event overprinting and later sedimentary cover. The magnetotelluric (MT) method is well known for imaging subsurface electrical resistivity variations related to temperature and fluid/melt content (Wannamaker et al., 2009; Chave and Jones, 2012). However, MT also has a role in delineating fossil terrane sutures based upon low resistivity created by underthrust, metamorphosed, and remobilized graphite- or sulfide-bearing metasediments (e.g., Boerner et al., 1996; Wannamaker, 2000, 2005; Jones et al., 2005; Yang et al., 2015; Bedrosian, 2016). The Paleoproterozoic Era was prodigious in creation and preservation of organic matter and sulfide in sedimentary sections along continental margins (Des Marais et al., 1992; Des Marais, 1994; Canfield, 2005). The Penokean domain has been known since the early MT data of Dowling (1970) to possess deep crustal resistivity over broad scales that is lower than that expected for purely crystalline silicate rocks of cratonic interiors.

During the middle 1980s, 14 tensor MT soundings plus a dense profile of controlled-source audio magnetotelluric (CSAMT) data were acquired in eastern Minnesota and northwestern Wisconsin over the presumed location of the Superior-Penokean suture from its near foreland area at the northwest end, across the Midcontinent rift system, and into the internal domains of the accreted Penokean arc (Fig. 1). The MT data acquisition was motivated in part by nearly contemporaneous Consortium for Continental Reflection Profiling (COCORP) and Great Lakes International Multidisciplinary Program on Crustal Evolution (GLIMPCE) regional seismic investigations of the Midcontinent rift system (Wunderman, 1988; Hinze et al., 1992), while the CSAMT data were taken for mineral resource assessment by the State of Minnesota Department of Natural Resources (EMS Inc., 1985). These older data have limitations, as discussed herein, but their quality is sufficient to reveal first-order features of the Penokean suture zone and the Midcontinent rift system overprint, and to locally focus the larger-scale resistivity structure imaged by the EarthScope Magnetotelluric Transportable Array (MT TA; Yang et al., 2015; Bedrosian, 2016). In this study, we recovered these results and applied modern methods of inversion imaging (for analysed MT data file, see GSA Data Repository Item1).

The MT method makes use of naturally occurring, low-frequency electromagnetic (EM) wave fields as sources for imaging subsurface electrical resistivity (or its inverse, conductivity) structure (Chave and Jones, 2012). Propagation of MT wave fields in the conducting media of Earth is diffusive; for example, in a uniform half-space, penetration depths are determined by an e-folding distance (skin depth), which is dependent upon wave period (inverse of frequency) and resistivity. The observed horizontal electric (E) and magnetic (H) fields in the frequency domain are related in a tensor fashion to provide apparent resistivities (ρa) and impedance phases (φ; see  Appendix), which are primary MT response functions determined solely by Earth properties. The x and y coordinate orientations of the responses commonly have x assigned to geological strike, if strike is apparent, so that ρxy and φxy in a two-dimensional (2D) situation correspond to electric current flow parallel to strike (transverse electric or TE mode), while ρyx and φyx correspond to current across strike (transverse magnetic or TM mode) (e.g., Wannamaker et al., 2009, 2014, 2017). The MT response functions become the inputs to nonlinear inversion algorithms, a form of tomography, to derive quantitative and stable resistivity structural models (Chave and Jones, 2012).

The 14 MT soundings, for which positions are shown in Figure 1, were obtained in early 1986 by Michigan Technological University (Wunderman, 1988; see  Appendix). Sites were located according to land access, avoidance of substantial overburden, and lack of urban development. Also, the general trend toward the southeast kept the MT profile within the Penokean Pembine-Wausau arc terrane in that direction rather than entering the Marshfield terrane. The x axis for all wave periods and all sites is N060E, which is approximately parallel to overall strike of both the Penokean suture and the Midcontinent rift system in this area. The wave period range of the data set is ∼0.01 to 300 s (frequencies of 100–0.0033 Hz), which we show is sensitive to resistivity structure from around 300 m to 100 km depth. The eight CSAMT soundings, also located in Figure 1, extend over a distance of ∼17.5 km and were obtained under private contract (EMS Inc., 1985; see  Appendix). Only ρxy and ρyx values were recorded (no φ) at seven discrete frequencies over a range of 1000–10 Hz. This yields a useful depth of exploration of ∼5 km. Given the CSAMT profile orientation, the x axis for these sites differs modestly from the MT data set at N080E.

The MT and CSAMT data are assembled into pseudosections in Figure 2, where the location serves as the abscissa and the wave period serves as the ordinate for contour plots of ρa and φ. As indicated in Figure 1, the CSAMT sites lie just north of the Malmo discontinuity upon metamorphosed, early Proterozoic deposits of the Mille Lacs Group. These rocks originally were continental margin sediments and volcanics, with possible early Penokean foreland sediments (Schulz and Cannon, 2007; Boerboom et al., 2011). Two small discrete conductors are evident: the ρa minimum centered on T = 0.01 s under the third CSAMT site from the profile’s northern end, and a minimum toward T = 0.1 s at the fifth and sixth site from that end (Fig. 3). However, the sites otherwise are predominantly resistive (few thousand Ω·m), showing no further evidence of low resistivity, at least in the upper crust under this area of the Mille Lacs lithology.

Close to the southeast portion, over the Snake River area, the westernmost four MT sites show both ρxy and ρyx values plunging toward long periods with high values of φxy and φyx (Fig. 2). These lie mainly upon high-grade Archean-cored gneisses of the McGrath–Little Falls domain and post-Penokean plutons such as the Isle and Warman granitic panel of the East-Central Minnesota batholith (Chandler et al., 2008; see Fig. 1). This behavior signifies a large-scale conductor at depth in the crust, which we will argue represents remnants of the Penokean suture in our study area. We note that the westernmost MT site (N03), in particular, shows rising ρyx and low φyx values for T > 10 s (Fig. 3A, feature A) relative to sites immediately to the southeast (Fig. 3B, feature C in site S01). This suggests that the large conductor is pinching out toward the surface to the northwest and is underlain by resistive material. At site S01, a maximum in ρyx and a low in φyx develop at high frequencies (feature B in Fig. 3B), which expand in magnitude for sites southeastward, as seen in the pseudosections. This reflects an increase in depth to top of the crustal conductor in that direction.

Further southeast, below St. Croix State Park and Clam Falls area, over the core of the Midcontinent rift system, ρxy and ρyx rise to moderate values (Fig. 2). In particular, there are anomalies of low φxy and φyx in the 0.1 to 2 s period range, signifying a laterally confined resistor in the middle crust, which we will equate to a cooled rift igneous volume (Fig. 3C, feature D). Only a relatively weak high in φyx around T = 10 s persists southeastward along the profile, representing a fading amplitude for the large dipping conductor. Reaching the Rice Lake and Flambeau Ridge area, near the southeast end of the MT profile, very high ρa values appear, corresponding to the Penokean internal (arc) domain. However, somewhat elevated values of φxy and φyx, spanning the 1–100 s range, reveal that resistivity in the lowermost crust may diminish slightly to moderate values (Fig. 3D, feature E). Features such as these will be correlated further in the discussion of the inversion model.

Although the older MT/CSAMT data presented here have some limitations, their reinterpretation using modern methods produces an insightful crustal model for the survey area. The MT data lack the on-diagonal elements Zxx and Zyy (see  Appendix) and so cannot be rotated to explore two-dimensionality. They also lack vertical magnetic field (tipper) responses, which provide complementary information about resistivity strike and lateral variation (Chave and Jones, 2012). Nevertheless, the pronounced drop in ρa and high φ values that the westernmost MT sites exhibit are coherent and unique to that region. These patterns must correspond to a first-order conductive structure, and, when interpreted in tectonic context, they provide new insights into the Penokean-Superior deep crustal transition.

A two-dimensional (2-D) resistivity model was produced from the observations in Figure 2 using a stabilized fitting (inversion) algorithm developed in-house that produces smoothed images of the subsurface similar to seismic tomography (Fig. 4). In particular, the inversion emphasizes fitting of the TM mode data subset ρyx and φyx values, because these data are understood to be more robust (usually) to departures in the 2-D assumption (Wannamaker, 1999). Close attention was paid to replicating the responses from the western MT sites into the CSAMT profile, because the Penokean suture is presumed to lie in this area. Further details on model construction are described in the  Appendix.

The most striking feature of the inversion cross section of Figure 4 is the pronounced low-resistivity zone under the MT stations on the East-Central Minnesota batholith in the Snake River area, which exhibits a moderate dip to the southeast from ∼5 km depth near profile km 30 to ∼40 km depth near km 100 under the Midcontinent rift system domain. The conductor appears to pinch out at its west end and not extend fully across the McGrath gneiss dome terrane to reach the Mille Lacs Group lithologies. It is conceivable that a large volume of low resistivities could be hidden somewhere at sufficient depth below the CSAMT profile and not influence the response. Such a geometry was entertained in the early forward modeling of Wunderman (1988). However, given the elevated ρa values of most CSAMT stations, which imply via skin depth that resistivities are high to ∼5 km depth, the low values of φyx, the rising ρyx for T > 10 s at the northwesternmost MT site N03, and the tight data fit of the inversion model, some pinchout of the conductor under the McGrath gneisses seems likely.

The resistivity section implies that much of the Penokean conductive zone underlies the Yavapai-aged (post-Penokean) Isle and Warman granites of the East-Central Minnesota batholith, which rose buoyantly and mushroomed over this apparent suture. This is where the higher frequencies and denser lateral sampling of the recovered MT/CSAMT data offer improvement in resolution over that from the EarthScope MT TA data, having 70 km average spacing and ∼10 s as the shortest wave period. In particular, this large conductor at depth tends to show up centered under the McGrath gneiss dome or even the Malmo discontinuity in the regional TA models (Yang et al., 2015; Bedrosian, 2016). The apparent suturing conductor is distinct from the large, quasi-horizontal ABA conductor of Bedrosian (2016) to the north of our study, which lies in the upper crust, probably in the younger, more distant foreland zone of the Animikie Basin.

Moving southeast to the main axis of the Midcontinent rift system in the St. Croix State Park to Clam Falls area, from profile km 85 to km 125, the model conductor visually steps downward and is overlain by a compact resistor reaching 5000 Ω·m in the 5–20 km depth range. This resistor results from the low values of φyx in the 0.1–2 s range, shown in Figures 2 and 3C. Moreover, the deep crustal conductor appears to terminate below the southeast side of the Midcontinent rift system near profile km 125 as the high in φyx also diminishes (Fig. 4). Toward the far southeast end of the profile, the limited drop in model resistivity to the 500–700 Ω·m range in the deepest crust and uppermost mantle follows from the moderate high in φyx from 1 to 100 s. This structure is considered tenuous and, among other things, could be influenced by the tectonically related, crustal-scale Flambeau low-resistivity anomaly to the northeast of our profile (Sternberg and Clay, 1977; Yang et al., 2015; Bedrosian, 2016). However, early MT soundings over Wisconsin imply that such mild decreases in resistivity toward the deepest crust may be widespread under the Pembine-Wausau terrane and at odds with whole-crustal resistivities of 10,000 Ω·m or more thought to be characteristic of cold, solely crystalline shield lithologies (Dowling, 1970).

The computed response in ρyx and φyx from the inversion model appears also in Figures 2 and 3. These response variations are smoother than the observed data because statistical noise (scatter) in the calculations is absent, and the inversion balances closeness of fit against model smoothness. However, the primary features of the observed data are reproduced well, for example, the high φyx and drop in ρyx defining the strong, east-dipping crustal conductor under the Snake River area of Figure 4. The low φyx near T = 10 s at the westernmost MT site and the high ρyx of the southern CSAMT sites again imply that the rocks under the McGrath gneiss dome should largely be resistive. The compact midcrustal resistor under the central Midcontinent rift system reproduces the low φyx response there well. Very high ρyx values toward the far southeast end are reproduced by the high model resistivities there in the upper 20 km of the model.

Coincidentally, the agreement between the observed and computed responses in the TE mode ρxy and φxy in Figure 2 is also reasonably good, even though these data were not formally inverted out of concern for more severe three-dimensional (3-D) effects (see  Appendix). An exception is the compact high ρxy anomaly near profile km 130, just east of the bounds of the Midcontinent rift system, which persists to the longest periods. This behavior is characteristic of the finite strike length of a structure and cannot be reproduced in a purely 2-D response in ρxy and φxy (Wannamaker, 1999). Similarly, the low computed φxy at the longest periods (T > 100 s) for the four western MT soundings in the Snake River area, unlike the observed values, can be explained by the finite strike of the large conductive body. Third, the anomalous high in φxy is at least as great as that of φyx at the southeast end of the profile under Rice Lake and Flambeau Ridge, suggesting low resistivity exists somewhere beneath the profile. However, 3-D structural variation can locally amplify such anomalies to imply higher than actual contrasts from a 2-D analysis (e.g., Wannamaker et al., 2017), and so we do not wish to overemphasize those observations.

Despite data limitations, the MT resistivity model section provides new insights with regard to Paleoproterozoic terrane assembly and subsequent late Mesoproterozoic rift-related modification (Fig. 5). We maintain that the Penokean suture is marked by a massive, southeast-dipping conductor that extends below the post-Penokean East-Central Minnesota batholith and westernmost Keweenawan volcanic cover in east-central Minnesota and is most prominent in the depth range 5–20 km. As interpreted regionally, and in other settings globally, the highly conducting material of the suture volume is taken to be mainly metamorphic graphite with possible sulfide contributions (Boerner et al., 1996; Wannamaker, 2005; Yang et al., 2015; Bedrosian, 2016). The carbonaceous material generally is of biogenic origin, originally laid down in reduced rift basins or in sediment-starved compressional foreland basins (cf. Giles and Dickinson, 1995), subsequently overridden and possibly thickened through imbrication in convergence. Specifically, the suture would be represented by the interface between the highly conducting metasediments associated primarily with the Superior Province margin and the overlying Penokean arc rocks (Pembine-Wausau terrane) accreted from the southeast.

The most likely representative of the carbon-bearing lithology in our study appears to be early Proterozoic Mille Lacs and North Range schists, and possibly Little Falls schist, containing rifted Superior margin reduced sediments plus perhaps initial Animikie foreland basin material at the outset of Penokean convergence (Southwick et al., 1988; Ojakangas et al., 2001; Schulz and Cannon, 2007; Boerboom et al., 2011). The exposed schists now lie in a series of north-verging, low-angle thrust sheets (Fig. 1), which were likely created mainly during the middle Penokean orogeny, when the Marshfield terrane collided along the Eau Pleine shear zone, but which were thoroughly remetamorphosed by post-Penokean events (Holm et al., 2007; Schulz and Cannon, 2007; see also Fig. 1). They are interpreted to predate and unconformably underlie the main Penokean Animikie foreland basin (Boerboom et al., 2011). Graphite is abundant in these schists, and its increasing crystallinity southward with metamorphic grade (McSwiggen and Morey, 1989) should promote conductivity (Wannamaker and Doerner, 2002; Bedrosian, 2016). Most occurrences originally were carbonaceous mudstones, commonly holding abundant pyrite, which could contribute (with graphite) to the low resistivity. The 13C/12C isotope ratios of all graphitic samples studied strongly imply a biogenic origin (McSwiggen and Morey, 1989). Bedrosian (2016) also proposed that graphite of the Mille Lacs and North Range schists could be responsible for the conductor in this general area resolved in the EarthScope MT TA data.

Some samples classify as siliceous carbonites, reaching 44 wt% carbon, and appear to be chemical precipitates lacking a significant clastic component (McSwiggen and Morey, 1989). Fluid mobilization and reprecipitation of ordered graphite during metamorphism are widely recognized, producing vein-like textures amenable to long-range electrical connection (Luque et al., 1998; Wannamaker, 2000, 2005). Circulation of fluids from subduction closure or orogenic collapse is suggested to be abundant within the suture (Vallini et al., 2007). However, the bulk of the original carbonaceous material today lies deep in the suture, according to the resistivity model, so that correlation with surface units is tentative. The limited extent of low-resistivity material resolved by the CSAMT data across the Mille Lacs Group rocks exemplifies the heterogeneity possible in the rock properties, despite abundant graphite occurrences. Local airborne and ground EM surveys, plus outcrop studies, have revealed the sheet-like, anisotropic fabric of the graphite and sulfide conductors at finer scales here and near the Niagara fault zone (Sternberg and Clay, 1977; Young, 1977; McSwiggen, 1987). We suggest that the deep conductor may have such a high magnitude because of thrust imbrication, a greater quantity of original graphite or sulfides along this portion of the subducted paleomargin, or fluidized redistribution during Penokean and post-Penokean events.

The Penokean suture was not destroyed per se in our study area (cf. Schulz and Cannon, 2007), but post-Penokean emplacement from below the Archean-cored McGrath Dome gneisses and Yavapai-aged (1775–1750 Ma) East-Central Minnesota batholith (Figs. 1 and 5) disrupted and obscured it at the level of the present-day bedrock surface. Intrusive-cored gneiss doming replaces prior lithology with highly resistive crystalline rocks, as seen elsewhere along the orogen (Young and Rogers, 1985). The foreland rocks overlying the gneiss belt that were not downfolded and incorporated would have been uplifted and eroded away, separating the surviving Mille Lacs and North Range Groups from the conductive lithologies detected down the former trench. These events would complicate an interpretation of the Malmo discontinuity as the original Penokean suture due to deformation and the degree of later Proterozoic injected material (gneiss plus pluton). Our interpretation compares reasonably with that of Chandler et al. (2007, 2008), where they viewed the Malmo discontinuity as a structure diminished by subsequent events and suggested that one should perhaps seek the Penokean suture south of the McGrath gneiss dome in our study area. Boerboom et al. (2011) viewed the discontinuity today as a Yavapai-aged (post-Penokean) boundary along which doming and exhumation of gneiss (i.e., the McGrath) took place.

We consider the large crustal conductor resolved here to be correlative with the Flambeau anomaly and related conductors to the northeast of our study along the Niagara fault zone (Sternberg and Clay, 1977; Wunderman, 1988; Yang et al., 2015; Bedrosian, 2016). The 70 km spacing of the EarthScope MT TA sites again is an issue for precise resolution, but the conductors also appear to approximately correspond with organic carbon– or sulfide-bearing schists metamorphosed to moderately high grades in a similar setting to ours. The coarse TA site sampling leaves questionable the degree to which conductors extend southward down the Niagara fault zone, and regional erosion since early Proterozoic time may have caused southward shifts in the surface trace of the suture. Comparable MT studies in the Slave craton of northwest Canada have resolved conductors dipping to ∼150 km depth, interpreted to be ancient subduction material, through the full range of the graphite stability field (Jones et al., 2003).

The prominent conductor in Figure 4 dips southeastward, as expected for such a setting, but it becomes obscure below the superimposed Midcontinent rift system. The conductor appears to step downward at ∼85 km along the profile and is overlain by a compact resistor reaching 5000 Ω·m in the 5–20 km depth range from 90 to 120 km along profile. The resistor properties are compatible with intruded (synrift and near postrift) mafic rocks of the Midcontinent rift system event as imaged in detail in the seismic reflection and gravity lines across Lake Superior to the northeast and elsewhere (Hinze et al., 1997; Stein et al., 2015). Adams (1990) showed that possibly related rocks of the outcropping Duluth complex had resistivities often exceeding 10,000 Ω·m. Lateral limits are compatible with principal control by the Douglas and Lake Owens regional normal fault systems (Fig. 1). Asymmetry of the crustal rift structure is not resolved in our MT data, but we presume the controlling geometry is predominantly down on the east side, as our transect lies across the Brules second-order rift segment of the Superior rift zone (Dickas and Mudrey, 1997; see also Fig. 5). Upper-crustal rift strata contribute in only a modest way to the MT response, and there is no observable manifestation in the MT of structural inversion from later (e.g., Grenville) compression that formed the St. Croix horst (Hinze et al., 1997). However, an unknown proportion of the Pembine-Wausau arc rocks (Xtg) labeled as overlying the main conductor west of the rift could alternatively be rift volcanics.

Seismic P-wave receiver function (RF) imaging close to the northeast of our profile in the Superior Province Rifting EarthScope Experiments (SPREE) project (Zhang et al., 2016; their line SN) revealed the Moho deepened to >55 km under the Midcontinent rift system from an average regional value of ∼40 km (Keller, 2013; see also Fig. 5). In addition, a second imaged interface overlies this thickened crustal zone with a depth to top as shallow as 35 km, centered beneath the Midcontinent rift system. RF polarity signifies increased velocity from crustal background, and the zone is interpreted as mainly dense magmatic underplate, explaining much of the positive Bouguer gravity anomaly characterizing the majority of the Midcontinent rift system (Merino et al., 2013; Zhang et al., 2016). This structure, projected along strike of the Bouguer anomaly, does not show in the MT structure, as it should be highly resistive like the upper mantle immediately below. It also would be obscured by overlying, modestly lower resistivities of the suture, which are smoothed spatially in the inversion (Fig. 5). However, some of the gravity anomaly could be due to the middle-crustal resistive body noted earlier herein, as considered by Hinze et al. (1997). Speculatively, higher-resistivity upper mantle below the rift (Mdz of Fig. 5) could represent a possible Keweenawan melt depletion zone, perhaps corresponding to the high seismic Vs body of Shen et al. (2013) near there or, alternatively, underthrust Archean lithosphere.

Farther southeast toward the Penokean interior, the upper half of the crust is generally resistive and made up of high-grade meta-igneous rocks of the Pembine-Wausau arc and subsequent intrusives. At lower-crustal levels and into the upper mantle nominally, resistivities faintly fall below 1000 Ω·m, particularly below the Rice Lake area. This might represent trace amounts of solid phase conductors such as graphite introduced during precollisional subduction or redistributed during protracted metamorphism (e.g., Wannamaker, 2005), and possibly locally hydrated upper mantle at 100 km and beyond (Bedrosian, 2016). However, it also may include 3-D effects (sideswipe) from an offline structure such as the Flambeau anomaly or metasediments around the Eau Pleine shear zone (Marshfield suture) to the south, such that the true depth extent may be less than that in the image. Note that the color contours in Figures 4 and 5 are logarithmic, so the necessary volume of grain boundary conductors below the southeast model area is very small compared to the main suture zone. Nevertheless, the early MT soundings of Dowling (1970) spanning much of Wisconsin imply that such mild decreases in resistivity toward the deepest crust may be widespread. Clearly, additional 3-D data coverage and modeling are warranted.

The MT transect data improve resolution of the Penokean and Keweenawan structure in east-central Minnesota and western Wisconsin from that provided by the 70-km-spaced, long-period (T ≥ 10 s) EarthScope MT TA soundings. In both the Yang et al. (2015) and Bedrosian (2016) models, in the ∼10–30 km depth range, the main western Penokean “suture” conductor lies between four regional MT stations and roughly centered upon the Malmo discontinuity. Our results suggest that the main conductivity lies farther southeast under gneiss and plutonic rocks and, presumably, Penokean arc rocks covered by Keweenawan rift strata (Figs. 1 and 5). This is in accord with tectonic history at this locale involving doming, intrusion, and suture disruption. Our model also offers tighter geometric bounds for the crustal intrusive units of the Midcontinent rift system, which affect the MT data almost entirely at periods shorter than the 10 s limit (frequency >0.1 Hz) of the MT TA coverage. This volume likely was emplaced in the main magmatic stage at 1101–1094 Ma during crustal-scale graben formation (Miller and Nicholson, 2013).

Older geophysical data, even if not necessarily of modern-day standards, can yield fresh insights into Earth structure and processes. The MT/CSAMT results we interpreted here are somewhat reconnaissance in their coverage, but first-order features are strongly suggestive of convergent and rift structures. The recovered MT data in east-central Minnesota and western Wisconsin confirm an extension of Penokean orogen deep structure into eastern Minnesota from the comparable terranes in Wisconsin and upper Michigan. They improve upon the resolution of the regional EarthScope MT TA coverage by showing that the bulk of the deeply underthrust, former margin and early foreland organic-rich, sulfide-bearing sediments probably underlie the eastern portion of the uplifted gneiss belt and post-Penokean intrusives and the westernmost Keweenawan volcanic cover. To an extent, these carbonaceous volumes have been disconnected from their possible correlatives in the Mille Lacs and North Range schists by subsequent uplift. Additional MT data of similar frequency range to the current study would be needed to further resolve details of the transition from the Mille Lacs Group across the McGrath gneisses and East-Central Minnesota batholith, and of the suture geometry in other parts of the orogen such as across the Niagara fault zone. The suture zone conductive structure appears to be displaced below the later igneous Midcontinent rift system by formation of a compact resistive body in the 5–20 km depth range approximately bounded laterally by the regional Douglas and Lake Owens fault zones. The bounds of this resistor agree well with similar interpretations from reflection seismic imaging and gravity modeling elsewhere along the Superior zone of the Midcontinent rift system. No distinct Moho-level, fossil magmatic underplating zone beneath the rift is resolved in the MT data, but it would be difficult to distinguish such a deep resistor from generally high overall resistivities of the uppermost mantle today. The influence of prior (Penokean) lithospheric structure on the geometry of the Midcontinent rift system emplacement locally is obscure because the two systems converge at an acute angle through our study area.

Institutional support for Wunderman included Michigan Technological University, Phoenix Geophysics (Denver), and a Minnesota Geological Survey Grant-in-Aid to students. Magnetotelluric (MT) data collection field assistance was provided by Mark Kitchen, Tim Seiss, Alan Koivunen, Eric Hepp, and Josie O’Brien. Individuals aiding in the early data collection and analysis together with published abstracts are listed in Wunderman (1988). Recovery and inversion of the MT and controlled-source audio (CSA) MT data set were supported under U.S. National Science Foundation grant EAR11-48081 to Wannamaker, as a component of the Superior Province Rifting EarthScope Experiments (SPREE) project (Stein et al., 2011, 2016). Anders Jeppson, formerly of Electromagnetic Surveys, Inc., kindly loaned us an archival copy of Minnesota Department of Natural Resources Report 8714 in order to allow high-resolution scans of the plates for CSAMT data recovery. Suzan van der Lee is thanked for providing locations of the SPREE SN line passive seismic sites. Wannamaker wishes to acknowledge numerous insightful discussions with Professor Charles R. Bentley of the University of Wisconsin, who in particular pointed out the early MT work of Dowling (1970). Greg Nash and Virginie Maris produced several of the graphics.


In the MT method, the observed horizontal electric (E) and magnetic (H) field components in the frequency domain are related through the complex tensor impedance Z (Chave and Jones, 2012), expanded as:



A similar tensor relation can be defined for the observed vertical magnetic field as:


The xy and yx MT impedances typically are transformed through simple arithmetic to apparent resistivities (ρa) and impedance phases (φ) for display and analysis (Chave and Jones, 2012). For example, converting MT impedance to ρa over a uniform half-space returns the resistivity of that half-space, so ρa is intended to depict a smoothed version of resistivity structure toward depth as period increases (frequency decreases).

The 14 MT soundings were measured using instrumentation constructed by Phoenix Geophysics, Ltd., and owned by Michigan Technical University. As is standard for MT soundings in our wave period range (∼0.01 to 300 s), magnetic fields were recorded using induction coils buried several inches below surface for mechanical and thermal stability. Electric fields were obtained as voltage differences across bipoles spanning from 25 to 125 m, grounded at the end points with quiet Pb-PbCl2 chemical electrodes. Field time series were digitized to 16 bits accuracy and then processed to the Fourier (frequency) domain assuming an exp(+iωt) time dependence on a Hewlett Packard 9845-B mini-computer. Fourier components were computed via cascade decimation with quality sorting according to field multiple coherence (Stodt, 1983). Instrument calibrations were carried out before and after field deployment.

The original MT data exist as paper plotting records only for the off-diagonal impedance quantities (xy, yx). Neither the on-diagonal (xx, yy) impedance values nor tipper elements are available. Hence, good resolution plots of the xy and yx components of ρa and φ (ρxy, ρyx, φxy, φyx), including error intervals from Wunderman (1988, his Appendix A), were scanned at 300 dots per inch (dpi), and a graphical grid was overlain for manual digitization. We estimate the digitized ρa values and error intervals to be represented to ∼5%, with phases to 1°–2°. These are less than the data error floors we applied for inversion. Example soundings of the off-diagonal apparent resistivities and impedance phases are presented in Figure 3.

The CSAMT method substitutes EM fields at discrete frequencies from a man-made transmitter located outside a survey area for the natural fields of MT (Goldstein and Strangway, 1975; Wannamaker, 1997). This can be a logistical advantage for shallower, focused surveying. The original CSAMT recordings comprise a dense profile of 45 sites, ∼17.5 km long, running slightly west of north from the northeast margin of Mille Lacs Lake (profile A of EMS, Inc., 1985; see also Fig. 1). A crossed-bipole transmitter with legs of ∼2 km length was located ∼12 km to the west of center of the profile. The crossed bipoles broadcast source fields with E dominantly normal to the profile (Zxy), and separately with E dominantly along the profile (Zyx), which is a scalar recording setup. Only ρa values were recorded (no φ) at seven discrete frequencies over a useful range of 1000 to 10 Hz. Values of ρxy and ρxy for the original 45 sites at seven discrete frequencies from 1000 to 10 Hz were entered from paper records (EMS, Inc., 1985), converted to |Zxy| and |Zyx|, and averaged laterally using a triangle function in accord with Electromagnetic Array Profiling (EMAP) principles (Torres-Verdin and Bostick, 1992) to yield eight equivalent binned sites with ∼2 km separation (Fig. 1). These binned soundings well represent the major trends in the original 45 sites. Given that the CSAMT recording sensors were oriented along and across the profile, the CSAMT ρxy and ρyx values have a uniform x orientation of ∼N080E for all sites and frequencies.

The digitized CSAMT and MT data were combined as represented in Figure 2 for input to inversion imaging. When carrying out 2-D inversion, generic 3-D model studies suggest that useful resistivity cross sections can be recovered in the face of along-strike (3-D) variations in structure if the ρyx and φyx data defined with fixed axes having y approximately parallel to the profile trend are inverted as the equivalent 2-D transverse magnetic (TM) mode (electric current flow in the profile direction; Wannamaker, 1999; Ledo, 2006). Fortunately, the paper records of the MT soundings were so defined, with x = N060E throughout. This is somewhat, though not seriously, oblique to the profile trend and is nearly along strike of the rift and suspected suturing structures. The CSAMT data were defined with x = N080E, approximately, which deviates somewhat from the MT, but we do not consider this a large angle, and the CSAMT soundings are not strongly anisotropic (Fig. 2). Thus, it is feasible to merge these sites with the MT in a single profile.

Inversion of ρyx and φyx in Figure 2 as TM mode data was performed using the University of Utah finite-element algorithm described in Wannamaker et al. (2009, 2014). This stabilized (regularized) inversion code minimizes an objective function, which is a weighted sum of data misfit and model roughness defined according to lateral model slope in y and z (cf. Chave and Jones, 2012). An error floor of 7% for ρyx and 2 degrees for φyx was applied for the MT data, with a 20% floor applied to the CSAMT ρyx. The mesh was 178 elements columns wide by 60 element rows deep. The thinnest row at the surface was 25 m thick, and deeper rows grew in thickness geometrically to a total model depth of ∼400 km. Horizontally, there were two columns of element parameters per MT/CSAMT station, except where there were larger gaps between clusters of sites, in which case there were more parameter columns. The starting model was a uniform half-space of 1000 Ω·m. The inversion converged monotonically in eight iterations from an initial normalized root-mean-square misfit (nRMS) of 14.1 to nRMS = 1.94, which is considered a good fit. The model roughness matrix norm was 0.3 times the norm of the error-weighted Jacobian matrix, and decreasing this factor did not improve fit significantly but did generate needlessly rougher models (see, e.g., deGroot-Hedlin and Constable, 1990). Close attention was paid to fitting the responses from the western MT sites into the CSAMT profile, as the Penokean suture is presumed to lie in this area.

1GSA Data Repository Item 2018112, Text file of MT/CSAMT response functions analyzed in this paper, with terms defined in the header information, is available at http://www.geosociety.org/datarepository/2018, or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.