Integrated geophysical analysis provides an alternate interpretation of the northern margin of the North American Midcontinent Rift System, Central Lake Superior

The Midcontinent Rift System (MRS) is a 1.1 Ga sequence of voluminous basaltic eruptions and multiple intrusions followed by widespread sedimentation that extends across the Midcontinent and northern Great Lakes region of North America. Previous workers have commonly used seismic-reflection data (Great Lakes International Multidisciplinary Program on Crustal Evolution [GLIMPCE] line A) to demonstrate that the northern rift margin in central Lake Superior developed as a normal growth fault that was structurally inverted to a reverse fault during a compressional event after rifting had ended. A prominent, curvilinear aeromagnetic anomaly that extends from Isle Royale, Michigan, to Superior Shoal in central Lake Superior, Ontario (the IR-SS anomaly), is commonly presented as a manifestation of this reverse fault. We have integrated multidisciplinary geophysical analyses (seismic-reflection, seismic-refraction, aeromagnetic, and gravity), physical-property information (density, magnetic susceptibility and remanence, and compressional-wave velocity), and geologic concepts to develop an alternate interpretation of the rift margin along GLIMPCE line A, where it intersects the IR-SS anomaly. Our new model indicates that a normal fault is the dominant structure at the northern rift margin along line A, contrary to the original rift-margin paradigm, which asserts that compressional structures are the dominant features preserved today. Integral to this alternate model is a newly interpreted, prerift sedimentary basin intruded by sills in northern Lake Superior. Our alternate model of the northern rift margin has implications for interpreting the style, scale, and timing of extension, rift-related intrusion, and compression during development of the MRS.

The North American Midcontinent Rift System (MRS) is a major magmatic continental rift that stretches approximately 2200 km from Kansas northeast to the Lake Superior region, where it curves southward and is inferred to extend into lower Michigan (Figure ¹; Ojakangas et al., 2001a). Because much of the MRS lies beneath Phanerozoic and surficial cover, it was first recognized by a series of prominent, mostly positive Bouguer gravity anomalies (Figure ¹) that were correlated to the 1.1 Ga basalts exposed in the Lake Superior region (Thiel, 1956; Hinze et al., 1997). The gravity expression of the MRS reflects the thick accumulation (tens of kilometers) of rift basalts overlain by an equally substantial clastic sedimentary section that developed over 60 million years (Figure ²) (Van Schmus and Hinze, 1985; Hinze et al., 1997; Ojakangas et al., 2001a).

The Lake Superior region has undergone extensive geoscientific investigation driven by mineral exploration and related geologic studies, academic studies of Proterozoic and Archean crust and plate-tectonic evolution, and, during a brief period in the late 1980s, oil and gas exploration. These investigations commonly used geophysical methods because of the poor rock exposure. Publicly available seismic and aeromagnetic data were collected as part of the Great Lakes International Multidisciplinary Program on Crustal Evolution (GLIMPCE) in the late 1980s (Green et al., 1989; Teskey et al., 1991). More recently, industry seismic-reflection data acquired over the lake in the 1980s have become publicly available as images of the processed seismic sections (McGinnis and Mudrey, 2003).

Paradigms about the nature and timing of the tectonic evolution of the MRS in the Lake Superior region have evolved over the decades of study. Modern workers commonly demonstrate their concepts through the interpretation of the GLIMPCE seismic-reflection lines, including line A in central Lake Superior (Figure ³), which crosses the northern margin of the rift near the central part of the lake (Cannon et al., 1989). This margin is marked by a prominent, narrow, aeromagnetic high-low pair that extends from Isle Royale to Superior Shoal (IR-SS aeromagnetic anomaly) (Figure ³). The IR-SS anomaly has been interpreted as an extension of a reverse fault inferred at Isle Royale, Michigan (Hinze et al., 1982; Teskey and Thomas, 1994). A reverse fault has also been interpreted from the seismic section and gravity data for line A, with near-vertical dip and approximately 3 km of throw (Thomas and Teskey, 1994). The commonly accepted paradigm asserts that the fault began as a normal growth fault that was later structurally inverted during postrift compression to a reverse fault (Cannon et al., 1989). This scenario is commonly applied to explain the overall tectonic style of basin margins within Lake Superior (championed in Ojakangas et al. [1997], for example).

In a review of previous reverse-fault models along line A, we discovered that the models were inconsistent and satisfied only some of the available geophysical data, as presented in the “Geophysical background” section. This observation led to questions about the validity of the reverse-fault interpretation for line A and an effort to develop a new, better integrated model. To approach the problem, we developed 2D potential-field models constrained by seismic information that test geologic concepts regarding reverse versus normal faults. Integral to these tests was a reinterpretation of the sedimentary section north of the fault and recognition of seismic patterns consistent with an extensive sill complex. The final integration of the geophysical data and geologic constraints provides an alternate concept to the commonly held interpretation of a reverse fault along line A. This has interesting geologic implications that require a reexamination of the northern margin of the MRS in central Lake Superior.

The 1.1 Ga MRS is characterized by voluminous, basaltic magmatism driven by the arrival of a mantle plume, followed by thick accumulations of sediments (Figure ²; Van Schmus and Hinze, 1985; Hinze et al., 1997; Ojakangas et al., 2001a). Despite extensive magmatism similar to that commonly associated with continental breakup, complete separation of the continental crust was never achieved. Instead, as magmatism waned, the rift transitioned to a sedimentary regime, with clastic sedimentation attributed to postrifting thermal subsidence. Faulting within a compressional regime heralded by the Grenville orogeny to the southeast (in current coordinates) created the geometry of the rift still preserved today (Cannon, 1994). Only minor tectonic activity, marine incursions, and glaciation have followed in the subsequent billion years (Huber, 1983).

Following the terminology adopted by Woodruff et al. (2019), magmatism developed in three main stages (plateau-, rift-, and late-rift-stages), recognized on the basis of physical volcanology, magnetic polarity (age), and geochemistry (Figure ²). The earliest plateau stage (approximately 1112–1105 Ma) was characterized by subaerial basalts erupted from fissures and low-angle shield volcanoes across a broad area, forming a volcanic plateau as much as 10 km thick. Paleomagnetic studies show that these rocks have strong, reversed-polarity remanent magnetization (for example, Books, 1972; Halls, 1972; Halls and Pesonen, 1982), commensurate with the ages of eruption during a reversal in the earth’s magnetic field (Figure ²). As demonstrated subsequently, this reversed-polarity magnetization is a distinctive attribute of plateau-stage basalts that produces large-amplitude negative aeromagnetic anomalies and can be used to generalize basalt packages in the geologic map and in geophysical models (Figure ³). Volcanic rocks in the Osler Group (“Osler volcanics”) are important representatives of plateau-stage basalts in the study area. Osler volcanics exposed on the land that stretches between the lake and Black and Nipigon Bays correspond to strong negative aeromagnetic anomalies (Figure ³).

Rift-stage magmatism (approximately 1102–1090 Ma) was characterized by voluminous, subaerial basalt flows that accumulated within rapidly subsiding basins in what is now Lake Superior. Basins exhibit synclinal shapes attributed to synmagmatic and postrift structural processes (White, 1966; Huber, 1983; Allen et al., 1997). Within the study area, synclinal axes generally follow broad aeromagnetic and gravity lows south of Isle Royale that curve to the southeast around the tip of the Keweenaw Peninsula and to the northeast next to Superior Shoal (Figure ³). Rift-stage volcanism began soon after the earth’s field returned from a reversed- to normal-polarity epoch (Figure ²), bestowing normal-polarity remanent magnetization to rocks of this and subsequent stages.

Mafic and ultramafic intrusions accompanied the plateau- and rift-stage volcanism and are exposed throughout the Lake Superior region (Figure ³). The best known is the expansive, mineral-rich, Duluth Complex in Minnesota (Figure ¹; Miller et al., 2002). In the study area, dikes from both magmatic stages are abundant along the northern shores of Lake Superior (Hollings et al., 2010). Sills that were emplaced along generally flat-lying, prerift sedimentary layers during plateau-stage volcanism are widely exposed north of Black and Nipigon Bays and south and southwest of the town of Thunder Bay (Figure ³; Hart and MacDonald, 2007; Hollings et al., 2010).

The late-rift stage of magmatism (approximately 1090–1083 Ma) was a time of waning volcanism and increasing sedimentation, and it is recorded mainly in the Copper Harbor Conglomerate of the lower part of the Oronto Group. Sporadic eruptions resulted in laterally discontinuous flows interspersed within clastic sediments, emplaced at ever-increasing intervals of time (Figure ²). The youngest dated volcanic rocks are approximately 1083 Ma (Fairchild et al., 2017).

After a widespread lacustrine incursion represented by the approximately 1081 Ma Nonesuch Shale, the locus, maturity, and provenance of sedimentation evolved (Ojakangas et al., 2001a) through to the final compressional stage (approximately 1060–1040 Ma). During the postrift stage, immature sediments of the Freda Sandstone of the upper part of the Oronto Group were derived from proximal MRS volcanic rocks and deposited mainly over the earlier volcanic basins. During the compressional stage, mature sandstones of the Bayfield Group and Jacobsville Sandstone were derived from a variety of rock types and deposited along the flanks of the volcanic basins. Although rift-related sedimentary rocks are as much as 9 km thick beneath the lake (Ojakangas et al., 2001a), only conglomerates and sandstones of the Copper Harbor Conglomerate (lowermost part of the Oronto Group) are exposed in our study area (Figure ³).

After the GLIMPCE seismic data became available in the 1980s, work to understand the structural architecture of the MRS drew from studies of the modern East African Rift to interpret several MRS rift segments partitioned by bounding normal faults (Cannon et al., 1989; Dickas and Mudrey, 1997; Hinze et al., 1997). These normal faults were interpreted as the locus of major reverse reactivation during the final compressional stage based on observations from outcrop and geophysics, especially along the western arm of the MRS that extends southwestward toward Iowa (Figure ¹; Chandler et al., 1989). The south-verging, low- to high-angle, reverse Keweenaw fault is considered to be a key example of such a normal-inverted-to-reverse fault within Lake Superior (Cannon et al., 1989). The Keweenaw fault, which carries rift-stage basalts over the younger Jacobsville Sandstone, has long been recognized on the Keweenaw Peninsula from drilling and outcrop (Lane, 1911; Butler and Burbank, 1929). It can be traced, generally parallel to the strike of the basalt flows, for the length of the Keweenaw Peninsula into Wisconsin where it dies out, a distance of approximately 230 km. An analogous, north-verging reverse fault is inferred (but not exposed) along the north shore of Isle Royale. The inference arises from several lines of evidence (Lane, 1911; Halls and West, 1971; Hinze et al., 1982; Huber, 1983), including (1) linear ridges of southeast-dipping, stacked basalt flows that correlate with and mirror the geology of the Keweenaw Peninsula, (2) a steep scarp along the island’s north shore that suggests faulting, (3) sharp, linear aeromagnetic gradients that follow the channel between Isle Royale and the mainland, and (4) geophysical interpretation of basalt underneath approximately 1 km of rift-related sedimentary section from seismic refraction studies in the Isle Royale channel.

Sedimentary rocks deposited in two separate tectonic events well before the onset of MRS volcanism also are important to the geologic setting of the study area: the Paleoproterozoic Animikie Group and the early Mesoproterozoic Sibley Group (Figure ³; Ojakangas et al., 2001b; Rogala et al., 2007). Both packages consist of partially eroded, siliciclastic, and chemical sedimentary rocks that are intruded or overlain by MRS igneous rocks. The Animikie Group comprises a relatively thin lower quartzite unit overlain by iron formation and a thick succession of black shale, argillite, and graywacke (the Gunflint Iron and Rove Formations in the study area) that were deposited in a broad foreland basin (Ojakangas et al., 2001b). Although the Animikie Group is less than 1 km thick in the study area, it thickens to the southwest; more than 1 km of section has been drilled in northeastern Minnesota without reaching the basement (Lucente and Morey, 1983; Ojakangas et al., 2001b). The Sibley Group is an approximately 1 km thick package of sedimentary rocks deposited in lacustrine, fluvial, and eolian environments (Rogala et al., 2007). Only approximately 300 m of the Sibley Group is preserved on land; the other portion was discovered through drilling in Nipigon Bay (Rogala et al., 2005; Figure ³). The Animikie and Sibley Groups were deposited on the peneplaned Archean Superior Province craton, which is composed of a series of generally northeast-trending belts of metasedimentary and metavolcanic rocks intruded by granites (Santaguida, 2001a, 2001b).

Geophysical interpretations have been critical to understanding the geology and tectonic evolution of the MRS in the Lake Superior region because so much of the rift is covered. MRS basalts commonly have high compressional-wave (P-wave) velocities and high densities, and they are very magnetic (Tables 1 and 2). These high values contrast well with the weak values of the same physical properties held by MRS sedimentary rocks, a situation allowing for strong geophysical expressions related to geology. Thus, much of the structure within the lake previously has been interpreted from geophysical data, sometimes presented as though it were known geology.

The commonly accepted Isle Royale fault is a geophysically inferred structure that is commonly mapped as one of the known, normal-inverted-to-reverse faults, which extends from Isle Royale, Michigan, to Superior Shoal, a shallowly submerged bedrock high approximately 80 km south of Terrace Bay, Ontario (Figure ³). The fault was first interpreted from vintage aeromagnetic data by connecting prominent, sharp anomalies between Isle Royale and Superior Shoal (Hinze et al., 1966; Wold and Ostenso, 1966). These anomalies appear as a prominent, fairly continuous, narrow, aeromagnetic high-low pair in the GLIMPCE aeromagnetic data (the IR-SS anomaly, on Figure ³), which extends from the northeast tip of Isle Royale to a strong, negative anomaly over Superior Shoal. These early workers recognized linear aeromagnetic patterns that followed the northern edge of the northeast-striking, upturned basalt layers along the rib of Isle Royale and extrapolated the interpreted origin as a reverse fault from Isle Royale to Superior Shoal. Thus, the Isle Royale fault was anticipated when seismic data for line A became available. It was intuitive to interpret a prominent, high-impedance reflection north of Superior Shoal, which we call the North Lake (NL) reflection, as the top of the MRS basalt section on the north side of the fault (Cannon et al., 1989; Figure ⁴). Prominent reflections higher in the section south of the IR-SS anomaly thus appeared to support the presence of the Isle Royale fault as a steeply dipping, reverse fault with approximately 3 km of throw in that area (Cannon et al., 1989; Green et al., 1989; Dickas and Mudrey, 1997). The reverse-fault interpretation also has been applied to other geophysical data along GLIMPCE line A, including seismic-refraction modeling (Shay and Tréhu, 1993) and gravity modeling (Thomas and Teskey, 1994; Figure ⁴). The interpretation of a prominent reversed fault along line A that is an inverted normal fault is still generally accepted (Stein et al. [2015], for example).

Irrespective of the presence of a reverse fault along the north shore of Isle Royale, the descriptions and models of a reverse fault extending all the way to Superior Shoal are not consistent in the literature. For example, previous mapping (based on the vintage aeromagnetic data) place the reverse fault north of Superior Shoal, whereas the prominent positive peak of the IR-SS aeromagnetic anomaly passes to the south (Figure ³). Teskey et al. (1991) point out how the arcuate segments of the IR-SS anomaly strongly suggest low-angle thrusting, yet steeply dipping, large displacement reverse faults are consistently interpreted for line A (Cannon et al., 1989; Green et al., 1989). A gravity model based on seismic-reflection interpretations (Thomas and Teskey, 1994) shows a steeply dipping reverse fault under the IR-SS anomaly that truncates a block of subhorizontal basalt layers (Figure ⁴). In contrast, Teskey and Thomas (1994) model the IR-SS aeromagnetic anomaly approximately 25 km west of line A as a narrow protrusion emanating upward from a basalt layer at approximately 3 km depth. Moreover, the gravity model of Thomas and Teskey (1994) does not account for magnetic bedrock at the lake bottom near Superior Shoal, which is related to the IR-SS anomaly (Figure ⁴). In this paper, we focus on developing a revised model of the source of the IR-SS anomaly along GLIMPCE line A that better integrates all known geophysical and geologic data.

A large amount of diverse geophysical and physical-property data has been acquired and compiled by numerous workers since the 1960s. Data useful for providing constraints for our study are described in this section and include geophysical data (seismic-reflection, seismic-refraction, and aeromagnetic and gravity) and physical-property information (density, magnetic susceptibility and remanence, and P-wave velocity). A heat-flow study of lake sediments (Hart et al., 1994) also provides qualitative constraints. Magnetotelluric data have been collected onshore, but they only poorly resolve features in the central lake (Bedrosian, 2016). Onshore and near-onshore areas have been extensively drilled for metallic mineral exploration since the 1840s, but most boreholes do not penetrate more than 1 km. Only a handful of deeper boreholes far outside the study area have geophysical logs. Thus, boreholes and drill core provide only limited information for interpreting geophysical data within the central lake.

Physical-property measurements and values used for geophysical modeling by previous workers for the central Lake Superior region are presented in Table 1 (densities and P-wave velocities) and Table 2 (magnetic susceptibilities and natural remanent magnetizations). Previous workers have recognized a general correspondence between densities and velocities of the volcanic and sedimentary rock suites in the Lake Superior region (Steinhart and Smith, 1966; Halls, 1969). As expected, basalts and mafic intrusions have high density and velocity values compared to sedimentary rocks. However, limited information for sedimentary rocks of the Rove Formation of the Animikie Group and the Copper Harbor Conglomerate, which is the basal section of the Oronto Group, indicates that these rocks have velocities in the range of 5.0–5.5 km/s, which overlaps with velocities of some MRS basalts (Table 1; Halls, 1969). This example emphasizes that using physical properties as proxies for geologic units must be applied with caution.

Extensive paleomagnetic studies of MRS basalts in the Lake Superior region have consistently shown that primary remanent magnetization is the dominant component contributing to the total magnetizations of the rocks (Table 2), with typical Köenigsberger ratios of remanent to induced components of 1–5 (Halls, 1972; Halls and Pesonen, 1982; Chandler, 2002). Remanent directions are commonly generalized for plateau-stage basalts as reversed with inclination of –60° and declination of 110°–120° and rift-stage basalts as normal with inclination of 40°–45° and declination 290° (Halls, 1972; Halls and Pesonen, 1982; Mariano and Hinze, 1994; Teskey and Thomas, 1994). Recent studies that investigate the migration of paleopoles with respect to plate movements show more nuanced variations in remanent magnetization orientations that depend on age (Swanson-Hysell et al. [2019], for example). Thus, we use remanent directions for our modeling that are generalized from samples measured in the northern Lake Superior area (Halls, 1974): a declination of 296.5° and inclination of 39.5° for normal-polarity basalts and a declination of 115° and inclination of –75.5° for reversed-polarity basalts. The earth’s field (and direction of the induced component) had a declination of –3° and inclination of 75° when the GLIMPCE aeromagnetic survey was acquired (see the “Aeromagnetic data” section).

Seismic-reflection sections considered in this study include the publicly available GLIMPCE lines A and G in northern Lake Superior (Figure ³; Green et al., 1989) and scanned, geolocated images of industry seismic-reflection sections for the western lake (McGinnis and Mudrey, 2003; LS-8 and LS-25 on Figure ³). The deep crustal, marine, multichannel GLIMPCE data were collected under contract to the U.S. Geological Survey and the Geological Survey of Canada in 1986 with a recording length of 20 s. Data collection used an approximately 128 L air-gun array, recorded with a 120-channel, towed streamer (Lee et al., 1988; Green et al., 1989). Conventional marine seismic and additional poststack processing were applied (Lee et al., 1988), and then the data were migrated (Milkereit et al., 1988). We obtained the digital data from a web site hosted by the Canadian government (Natural Resources Canada, 1987). The trace interval was 50 m, and the vertical sampling was 8 ms. Based on a range of velocities of 4.0–6.5 km/s and a frequency range of 20–30 Hz, we estimate the vertical resolution to be approximately 120 m, increasing with depth. We selected the version of the data to which an f-x decon filter (Canales, 1984; Gulunay, 1986) optimized for random noise attenuation and then a poststack time migration had been applied. Detailed descriptions of the acquisition, processing, migration, and filtering of the data are available from Agena et al. (1988), Lee et al. (1988), and Milkereit et al. (1988). The northern half of the migrated time section is shown in Figure ⁴ for 5 s of two-way traveltime (TWT).

Two-dimensional industry seismic-reflection lines with a recording length of 8 s were collected in 1985 by Grant-Norpac, Inc., for oil and gas exploration. Data were collected using a 20-gun tuned air-gun array and 120-channel, towed receivers and processed using a conventional marine processing sequence (McGinnis and Mudrey, 2003). We geolocated published images of the industry seismic lines by assigning the x-, y-coordinates to pixels along the length of each line. We then converted the images to SEG-Y files for input into geologic interpretation software using image conversion software IMAGE2SEGY (Farran, 2008).

Several seismic-refraction studies conducted in Lake Superior provide constraints on the velocity structure of the upper 10 km of the crust. Models derived from scattered sonobuoy profiles 15–30 km in length provide local constraints on modeling in the northern lake (Halls and West, 1971). However, only the model from profile 30 is close enough to project onto line A (Figure ³). Ray-tracing models of large-aperture data observed during GLIMPCE acquisition (Shay and Tréhu, 1993) provide comprehensive P-wave velocity models coincident with line A that can be used for modeling and for time-to-depth conversion of the seismic data. The northern half of the velocity model of Shay and Tréhu (1993) for line A, converted to TWT and superposed on the reflection section, is shown in Figure ⁵.

To do a time-to-depth conversion of line A, we generalized the velocity model of Shay and Tréhu (1993) by smoothing the transitions between boundaries (Figure ⁵). In particular, we disregarded details of the two-pronged, abrupt shallowing of velocities modeled near Superior Shoal (distance 70–78 km) to stabilize the depth conversion in that area. We converted the seismic section from time to depth by applying our smoothed velocity models to the legacy, migrated time sections. We did not have access to the original digital data for the industry lines, and data were missing for the GLIMPCE lines, so prestack depth migration was not possible.

Aeromagnetic data come from a recent compilation that merges GLIMPCE aeromagnetic data acquired over Lake Superior with magnetic grids from other previously published aeromagnetic surveys over onshore areas (Anderson and Grauch, 2018). Most of our study area is covered by the GLIMPCE aeromagnetic data, which was flown in 1987 with north–south flight lines spaced 1.9 km apart at a 305 m nominal altitude above the lake level or ground surface (Teskey et al., 1991). East–west tie lines were flown 20 km apart. A geomagnetic reference model of the earth’s magnetic field at the time that the survey was flown (DGRF) was removed from the total-field observations to obtain the residual total-field magnetic values, or the total-field magnetic anomaly (Anderson and Grauch, 2018). The total-field magnetic anomaly map, displayed using OASIS Montaj software, is shown in Figure ³.

We used flight-line and gridded data (500 m cell size) in this study. The gridded data were analytically continued to 300 m above the lake level offshore and ground surface onshore to best maintain the anomaly shapes commensurate with the resolution of the GLIMPCE survey then filtered to reduce flight-line noise (Anderson and Grauch, 2018). Two-dimensional modeling of the original data along flight lines instead of extracting data from the magnetic anomaly grid maintains the highest data resolution and reduces the possibilities of errors arising from gridding of the flight-line data, so we gave preference to modeling north–south profiles that are well situated to anomaly strike.

After several tests to separate short-wavelength from long-wavelength anomalies (representing shallow versus deep sources, respectively), we found that magnetic depth estimation was better suited to understanding the magnetic sources along line A. Based on the principle that wavelengths increase with the increasing source depth, magnetic depth estimation methods analyze gradients of the field data to estimate depths to the tops of magnetic sources, such as dikes and near-vertical contacts (Blakely, 1995). We present two methods for this study, which were applied in 2D and 3D from flight-line and gridded magnetic data, respectively. Two-dimensional methods applied to the high-resolution, flight-line data help us to understand the details of the gradients arising from primarily shallow sources, but they may be inaccurate where the lines are at an oblique angle to the trend of the anomaly gradient. Three-dimensional methods can handle anomaly gradients of variable trends, but they are limited by interference from multiple anomalies and by gridding artifacts arising from poor flight-line resolution. As is common with magnetic depth estimation, solutions often look scattered due to noise in the data or to magnetic contrasts that do not consistently follow geologic contacts because of the natural variability of magnetic rock properties.

We applied the 2D multisource Werner method (Hansen and Simmonds, 1993) as implemented by the program PDEPTH (Phillips, 1997, 2018) on flight-line data and the 3D extended Euler method (Reid et al., 1990) as implemented by Phillips (2002) on gridded data. Both of these methods analyze gradients within a sliding window given a choice of source type. The source type, i.e., structural index (SI), is set to 0 (contact) or 1 (sheet) for each run of the application. The window size is varied to test for shallow solutions (a small window size) versus deep solutions (a large window size). The most reliable results are obtained with a window size that captures the entire gradient without also containing nearby interfering anomalies. Each solution is intended to represent the locations, depths, and dips of the top edges of steeply dipping magnetic contrasts that match the assumptions of source type, such as dikes and dipping thin basalt layers for sheet solutions and faulted or steeply sloping basalt layers for contact solutions. Given a set window size, solutions returned for sheet geometries tend to be much deeper than those for contact geometries; the two types of solutions can be viewed as providing a maximum and minimum estimate, respectively, of the depth to the magnetic contrasts.

The 2D multisource Werner method was applied to the magnetic profile along GLIMPCE flight line 11160, which intersects seismic line A at the crossing of the IR-SS anomaly. The magnetic depth solutions were then projected onto line A for approximately 30 km north and south of the crossing of line 11160 and line A, a projected lateral distance of no more than 5 km. The method, which is a reformulation of the standard Werner deconvolution algorithm using the analytic signal (Hansen and Simmonds, 1993), looks for a specified number of source edges (e.g., the top and/or bottom of the source[s]) within each window as it moves over the profile. Parameters returned are the location of the source edge, dip of the contact or sheet, magnetization contrast, and number of solutions that cluster at one location. The more solutions in a cluster, the greater confidence that a source edge is located at the cluster location. Whereas dips are also estimated, these may be somewhat inaccurate in our study area because the analytic signal used by the method does not account for variable magnetization directions that could be due to postdepositional tectonic tilting (Nabighian, 1972).

The 3D extended Euler method was applied to the magnetic anomaly grid by testing window sizes of 4 km (11 × 11 points), 7.6 km (19 × 19 points), and 10 km (25 × 25 points) using contact and sheet geometries (SI = 0 or 1). Parameters returned are the source location averaged from all successful solutions for each window and the error estimated for the average (percentage of depth). The minimum source depth was constrained by the bathymetric surface (Anderson and Grauch, 2018).

Profile models were constructed using the program PDEPTH (Phillips, 2018) and GMsys2D, which is part of Geosoft’s OASIS software package. Both programs use the polygonal modeling approach of Talwani (1965) and allow for bodies to be truncated perpendicular to the profile (known as a 2.5D modeling). We treated all modeled bodies as 2D (extending infinitely in both directions perpendicular to the profile), except where noted. Physical properties are assigned to the bodies, and the magnetic and gravity effects of each body are computed and summed to give calculated magnetic and gravity curves, respectively. The corners of the polygonal bodies are moved by the user, and/or the physical properties are varied. The model effects are recalculated, and the process is repeated until a good fit between the observed and calculated curves is achieved. During this process, a constant value is automatically removed from the calculated curve to improve the comparison between the observed and calculated curves. The value of the constant can be attributed to the unknown effects of very deep sources and is not significant for the upper crustal focus of our modeling.

Three-dimensional forward modeling was accomplished using an unpublished program based on subroutines given by Blakely (1995). The program defines each grid interval of a topographic grid as a vertical prism extending infinitely downward from the topographic surface. Using uniform remanent and induced magnetizations, which can have separate orientations, the magnetic field of each prism is calculated and then summed together to give the computed magnetic field for the entire grid.

Simple Bouguer gravity station and gridded data were obtained from a recent compilation of publicly available databases, covering the same area as the aeromagnetic and bathymetric compilation (Anderson and Grauch, 2018). Data at individual station locations from several government data sources were combined, duplicates were removed, and the data were reprocessed to obtain Bouguer anomaly values using a reduction density of 2670 kg/m³. To maintain consistency between onshore and offshore stations, the Bouguer correction was calculated differently depending on whether the station was located on land, at the lake surface, or on the lake bottom and whether the lake-bottom stations were above or below mean sea level. Details can be found in Anderson and Grauch (2018). Note that Bouguer and free-air gravity values are very similar across our study area because the lake bottom varies no more than 200 m from sea level, the elevation at which the Bouguer correction is zero.

In our study area, gravity stations are spaced 5 km apart on average. This spacing is good for regional study, but it lacks sufficient resolution to resolve geologic features near the lake bottom that span less than 10 km. In the lake area, stations were mainly collected on the lake bottom. No correction was made to the Bouguer values equivalent to terrain correction on land. From inspection of gravity values across high-relief bathymetry in low-gradient areas of the Bouguer anomaly map of Lake Superior outside our study area (Anderson and Grauch, 2018), we estimate that the effects due to bathymetric relief are generally less than 4 mgal.

Bouguer anomaly data were gridded using a cell size of 1000 m and were low-pass filtered to produce a smooth surface (Anderson and Grauch, 2018). The final Bouguer anomaly map, displayed using OASIS Montaj software, is shown in Figure ³. Profile modeling of the gravity data followed the procedure described in the “Aeromagnetic data” section. Observed data for the profile models were extracted from the grid to better characterize the shapes of the gravity curve, but locations of the original data points are also included to emphasize the data limitations.

In this section, we present the interpretations and analyses that are required to integrate the seismic information with the gravity and magnetic data. Topics include an initial interpretation of seismic line A, analysis of the negative magnetic anomaly at Superior Shoal to ascertain geologic sources there, and simple magnetic models to test the viability of different geologic concepts.

The unmigrated time (Agena et al., 1988) and migrated time and depth-converted seismic-reflection sections for line A all show an approximately 10 km wide disrupted zone that extends vertically below the lateral extent of the IR-SS anomaly (Figures ⁴ and ⁶). The disrupted seismic zone separates packages of subhorizontal reflections on either side that do not appear to connect across it, confirming that the zone occupies the site of a major crustal discontinuity and cannot be discounted solely as a data acquisition problem. Within the top 2.5 km of the refraction velocity model of Shay and Tréhu (1993), the disrupted zone is modeled as a vertical, moderately high-velocity (5.0 km/s) column with two separate prongs that reach upward to the lake bottom (Figure ⁵). Velocities on either side of this column increase with depth and range from 2.8 to 4.0 km/s, a range typical of rift-related sedimentary rocks (Table 1). Identifying the rock types with velocities of approximately 5.0 km/s is more uncertain because velocities measured for several geologic units fall within this range (Table 1). However, the tips of the two prongs of the disrupted zone both correspond to clusters of magnetic depth estimates that indicate shallow magnetic bedrock (Figure ⁶), which is consistent with the presence of basalt at Superior Shoal less than 10 km away. Thus, it is likely that the rocks within the disrupted zone are at least partially composed of basalt, as either upturned, faulted basalt layers, intrusive feeders to basalt at the lake bottom, or a combination of these. Hints of reflections that turn upward into the disrupted zone along its south side could support either interpretation but more likely are migration artifacts.

Surfaces extrapolated from correlations that we have developed from the network of seismic reflection lines in western and central Lake Superior (Stewart et al., 2019) are plotted on top of the line A depth section in Figure ⁶. The surfaces, which only have been interpreted south of the disrupted zone, are generalized to represent, from top to bottom, packages of late-rift interbedded sediments and volcanic flows, rift-stage (normal-polarity) basalts, and plateau-stage (reversed-polarity) basalts, following the expected chronology of the rock types of the rift (Figure ²). We interpret the reflections below the plateau-stage basalts as possible pre-MRS sedimentary rocks plus sheet intrusions or as an earlier package of basalts of unknown affinity. Our interpretation on the south side of Figure ⁶ is similar to that of the Thomas and Teskey (1994) gravity model of Figure ⁴. The interpreted location of the top surface of the rift-stage basalt section is corroborated by several clusters of the 2D magnetic depth estimates from flight line 11160 projected onto the section (Figure ⁶). The several clusters of sheet solutions (the red circles with vertical tails) from 82 to 92 km distance are likely finding local structural disruptions or abrupt differences in magnetizations within the subhorizontal layers indicated by the seismic reflectivity.

Several features all occur at an approximate distance of 95–100 km on line A, south of the disrupted zone (Figure ⁶): an abrupt southward downturn of reflections interpreted as rift-stage basalt and plateau-stage basalt, a cluster of contact-solution magnetic depth estimates (the blue circles with tails) where reflections from our interpreted late-rift-stage section turn upward, and a broad peak in the gravity data. In addition, a northward-thickening, wedge-shaped package of reflections (labeled as “sedimentary wedge” on Figure ⁶) overlies the interpreted basalts in this area and appears to truncate the late-rift-stage section. The cluster of magnetic depth solutions may be locating upturned, interbedded flows within the late-rift-stage section that were tilted then truncated by erosion.

To interpret line A north of the disrupted seismic zone and the IR-SS anomaly, we pan to a view of the northern half of line A (Figure ⁷). The high-impedance NL reflection recognized from the time section (Figure ⁴) is prominent at depths of 2.5–4.0 km. This reflection is also prominent and extensive on crossing GLIMPCE line G (Figure ³; Cannon et al., 1989). Magnetic depth estimates projected onto the seismic section from the 3D Euler analysis do not support the NL reflection as the top of a thick basalt section. The solutions shown on Figure ⁷ assume a contact geometry, which generally represents the minimum depth to a magnetic contrast. They cluster 1–2 km below the NL reflection, suggesting that magnetic basalts lie well below the reflection. We do not consider that errors in our velocity model for this region, which would affect the depth to the NL reflection, would be large enough to account for this discrepancy. In addition, Shay and Tréhu (1993) point out that a velocity interface was not indicated in their refraction model at the NL reflection, which would be expected if it marked the top of a package of igneous layers. Instead, velocity interfaces were indicated above and below the NL reflection (Figure ⁵). These geophysical relations allow for an alternate, nonigneous interpretation of the NL reflection, such as a significant thickness of carbonate and/or chert layers overlain by clastic sediments.

Seismic patterns in line A differ above and below the NL reflection between approximately 25 and 70 km of horizontal distance. Above the NL reflection, the section is characterized by low-amplitude, subhorizontal reflections (except for possible multiples near the lake bottom). Refraction velocities within this interval range from 2.8 to 4.8 km/s (Figures ⁵ and ⁷), indicative of dominantly sedimentary rocks (Table 1). The absence of mafic rocks in this area is supported by a region of abnormally high heat flow localized over the northern part of line A (Hart et al., 1994). These authors argue that shallow sections of mafic rocks in Lake Superior correspond to abnormally low, not high, heat flow.

Below the NL reflection is an approximately 2–4 km thick zone characterized by low- and high-amplitude subhorizontal reflections. The locations of the magnetic depth estimates suggest that magnetic rocks lie within or below this zone (Figure ⁷). Below the zone of subhorizontal reflections are multiple, inclined and discontinuous reflections contained within an overall V-shape that appears to converge near a depth of 13 km at 40–45 km distance. Within the V, patterns show irregular, high-amplitude reflections that dip in opposite directions or in saucer shapes, cross-cut subhorizontal reflections, and interconnect in ramps or a staircase network (Figure ⁷). A similar succession of seismic patterns is evident for the northern part of line G (Cannon et al., 1989). These patterns are characteristic of those recognized as sill complexes intruded into sedimentary basins from modern seismic data in marine settings worldwide (Smallwood and Maresh, 2002; Planke et al., 2005; Polteau et al., 2008; Magee et al., 2016; Galland et al., 2018). The regional extent of the complex recognized on line A (approximately 20 km distance by approximately 10 km depth) is similar in appearance and scale to one documented in the Exmouth Plateau sedimentary basin of northwestern Australia (Rohrman, 2013). Thus, we interpret the strong reflections below the NL reflection as mafic sills and dikes cutting up into the sedimentary section of a pre-MRS sedimentary basin, possibly the prerift Sibley and/or Animikie basins (Figure ³ explanation) that are sporadically exposed on land or an unrecognized prerift basin. In contrast, previous workers have speculated that the dipping reflections on line A represent an early Proterozoic structural basin or a rift-related mafic igneous complex (Cannon et al., 1989; Green et al., 1989; Shay and Tréhu, 1993; Thomas and Teskey, 1994; Figure ⁴).

We interpret the eroded top of the pre-MRS sedimentary basin at the sharp reflection with a negative impedance contrast that aligns farther north with a velocity interface modeled from the refraction profile for sonobuoy 30 (Figure ⁷). The negative impedance contrast may reflect the remnants of basalt overlying the older sedimentary section (a decrease in velocity with depth), whereas the sonobuoy velocity interface may represent younger rift sediments directly overlying the older sedimentary section (an increase in velocity with depth). The rationale for this interpretation is further clarified by the modeling results in the “Integration” section.

Knowing the age and attitude of basalt flows at Superior Shoal is important for constraining geologic concepts during geophysical modeling. The shoal is an easterly trending, approximately 20 km long bathymetric bedrock high just below the water’s surface in central Lake Superior that is characterized by a strong negative aeromagnetic anomaly (Figures ³ and ⁸). It is composed of basalt on the south and sandstone on the north (Manson and Halls, 1991).

Previous studies of Superior Shoal give conflicting results on the age of the basalts, based on their magnetic polarities (Figure ²). Manson and Halls (1991) conclude from paleomagnetic measurements obtained from submersible sampling that the basalts have normal-polarity remanence and thus represent a rift-stage origin (Figure ²). Teskey et al. (1991) notice an inverse correlation of observed magnetic profiles with bathymetry, concluding instead that magnetizations of the Superior Shoal basalts are dominated by reverse-polarity remanence and thus represent an older, plateau-stage origin. Grauch and Schulz (2019) find that several normal-polarity samples from the Manson and Halls study have chemistries that support a plateau- and not a rift-stage origin, following Nicholson et al. (1997).

A 3D magnetic model using a typical reverse-polarity remanent direction assigned to the bathymetry (Figure ⁸) clearly shows the inverse correlation noted by Teskey et al. (1991), providing evidence that most of Superior Shoal is composed of the older, reverse-polarity (plateau-stage) basalts. Because the Manson and Halls (1991) paleomagnetic study appears to be valid (Grauch and Schulz, 2019), the normal-polarity results may evidence an excursion or short reversal in the earth’s magnetic field that has not been previously documented, or, less likely, an aberration in the regional evolution of magma chemistry observed throughout the rift (Nicholson et al., 1997). In any case, the magnetic effect of the normal-polarity rocks is negligible compared to the volume of reversely polarized plateau-stage rocks underneath (Figure ⁸). Thus, Superior Shoal provides an important geologic constraint that a significant volume of the older, plateau-stage basalts is present at very shallow depths at Superior Shoal. Moreover, Manson and Halls (1991) infer west–northwest strikes and steep northerly dips for the basalts adjacent to the IR-SS anomaly (Figure ⁸).

The complex geologic history in the Lake Superior region, the large number of unknowns at depth, and the overlap and variability of physical properties burden geophysical modeling with a high degree of ambiguity. Thus, before attempting a complicated geophysical model, we use simple magnetic models to explore the viability of different geologic concepts for the interpretation of the IR-SS anomaly and the disrupted seismic zone on line A (Figure ⁶). Ambiguity in the modeling is reduced by using simple polygonal bodies, fixing geometries south of the disrupted zone where we have confidence from the seismic interpretation, assigning consistent magnetizations to model bodies, using magnetic depth estimates as guides, and recognizing geologic constraints from Superior Shoal. Arriving at magnetic models that produce similar (but not necessarily exact) anomaly shapes to what is observed for the IR-SS anomaly provides an important criterion for determining the viability of the tested geologic concepts. Focusing on magnetic anomalies is appropriate because of (1) the sensitivity of the magnetic data to shallow sources, (2) the strong influence of remanent magnetizations of different polarities that indicate different age basalts, and (3) the independence of magnetic properties, whereas density and velocity are commonly correlated (Steinhart and Smith, 1966).

Figure ⁹ shows two geologic concepts that are tested for data along flight line 11110, which passes near line A but avoids the complicated anomalies crossed by the seismic line (Figures ³ and ⁸). The models represent (a) a reverse fault with upturned basalt layers that is thrust over a basalt layer that could be normal or reversed polarity and (b) a normal fault with only reversed-polarity basalt present in the footwall. In both scenarios, the separate layers labeled “R” are considered to be reverse-polarity, plateau-stage basalts that are older than the normal-polarity, rift-stage basalts represented by layers labeled “N,” in keeping with the age relations of rift basalts (Figure ²). Geologically, where the normal-polarity (N) layer is absent on top of the reversed-polarity (R) layer, we assume that normal-polarity basalt never accumulated there or was originally present but subsequently eroded away. Where the normal-polarity (N) layer is proposed without an underlying reverse-polarity (R) layer in the alternate for the scenario in Figure ⁹, we envision that the reversed-polarity basalts were eroded away before extrusion of the normal-polarity basalts.

Because the strong remanent magnetizations tilt in concert with any deformation of their associated basalt layers, it is important to account for dips of layers in the modeled magnetizations (Mariano and Hinze, 1994). Table 3 lists the magnetizations assigned to bodies in the magnetic models as well as the remanent directions computed after assuming specific tectonic tilts and a geologic strike of 93° from north. Each model was constructed so that reversed-polarity rocks reach near the surface, as dictated by the constraints at Superior Shoal. South of the fault in each case, the geometry of layers and their magnetizations are fixed, following the initial interpretation of line A (Figure ⁶). The northern reversed polarity body in the footwall of the normal fault model is forced to taper to zero thickness to the north, following the interpretation from seismic data that basalts are generally absent within the top 3 km of the section some distance north of Superior Shoal.

The results of the model tests (Figure ⁹) show that the anomaly produced by the normal fault model with reversed-polarity basalts in the footwall comes much closer to a fit of the observed anomaly than does the reverse fault with upturned basalt layers, no matter whether the footwall is normal or reversed polarity. The critical variable that makes the normal fault model work is the shallow wedge of strongly magnetized reversed-polarity basalts in the footwall that (1) is truncated and beveled on its south side and (2) extends outward toward the north. The wedge is also pleasing because it allows for a northerly tilt of the reversed-polarity basalts, consistent with the dip directions inferred by Manson and Halls (1991; Figure ⁸).

An alternate view of the scenario tested in Figure ⁹ is that it represents a south-verging reverse fault, where the fault that offsets the reversed-polarity layers (R-tilt15°N on the north and R on the south) is dipping to the north. Although this explanation is possible, it is challenging to explain geologically. A reverse fault indicates formation during the compressional stage (Figure ²) and as such, should verge to the north rather than to the south because of the north-directed stresses at this time (Cannon, 1994). A south-verging reverse fault is also incompatible with the north-verging reverse fault at Isle Royale, to which the IR-SS anomaly connects. The south-verging Keweenaw fault on the Keweenaw Peninsula is an exception because north-dipping basalt layers that formed during the rift stage provided a favorable orientation for the later fault displacement (Huber, 1983).

We now construct a more complicated gravity and magnetic model of flight line 11110 (Figure ¹⁰) to demonstrate the viability of the concept of a normal fault with reversed-polarity basalts on the north that resulted from the simple model tests of Figure ⁹. The densities and magnetic susceptibilities assigned to the bodies are listed in Figure ¹⁰; the remanent magnetizations are coded and explained in Table 3. Assigned densities follow the seismic refraction model of Shay and Tréhu (1993), using velocity and density relations observed by Steinhart and Smith (1966) for igneous rocks in the Lake Superior region and Gardner et al. (1974) for sedimentary basins. Geologic interpretations are derived from seismic interpretations (Stewart et al., 2019) and inspection of density ranges typical for rock types in the area (Table 1).

Aside from the depiction of a normal fault at Superior Shoal that marks a major difference in geometry of the rift basalts on either side, several other concepts from the initial analyses and interpretations are incorporated into the integrated magnetic and gravity model of Figure ¹⁰. The body labeled “sill complex” at less than −5000 m elevation on the north end of the profile is intended to represent the sill complex intruded into the lower part of a pre-MRS sedimentary basin that we postulated from the north half of line A in the “Initial interpretation of seismic line A” section. Bodies labeled “Older Sed” represent the upper part of the older sedimentary basin that is above the sill complex and devoid of significant igneous intrusion. “Carbonate/chert layers” represent our interpretation of the NL reflection as a sedimentary impedance contrast rather than the top of the basalt. Magma that formed the sill complex likely originated north and northwest of the Figure ¹⁰ model, where a possible magmatic source is indicated by broad, large highs in the gravity and magnetic data in northern Lake Superior (Figures ³ and ⁴). The extent of the sills to the south of the sill complex body is uncertain, as is the geologic identification of the body labeled “prerift crust.” This latter body could represent more of the older sedimentary basin intruded by a smaller volume of sills, Archean gneiss that is intruded by sills or has a seismically reflective fabric, or an unknown combination of these scenarios. Because all of these possibilities are somewhat magnetic (Table 2), we have assigned the prerift crust body a magnetic susceptibility of 0.01 SI.

The body labeled “pervasive intrusions?” below approximately −6000 m in elevation between the distances of 80 and 96 km is intended to represent multiple mafic intrusions that extend upward from the south, from below the extent of the model. This interpretation is supported by the broad gravity high that peaks at approximately distance 97 km and with seismic velocities reaching greater than 6.3 km/s below and to the south of our model (Figure ⁵; Shay and Tréhu, 1993), which may indicate the magmatic source. The “pervasive intrusions?” body, as well as the less-intruded layer labeled “intruded sed” above it, may be analogous to the older sedimentary basin intruded by a sill complex that we interpreted on the north; an interpretation that is supported by the layered seismic reflectivity in line A below approximately −4500 m in elevation between the distances 85 and 95 km. Alternatively, these two bodies may represent varying concentrations of sheet-like intrusions that followed metamorphic fabric in the Archean basement, or they may be different portions of a thick section of pre-plateau-stage basalts, as yet unrecognized in the evolution of the rift (Figure ²).

Bodies labeled “rift sed” at the top of the model across the whole distance are assumed to be effectively nonmagnetic clastic sedimentary rocks and are assigned densities in line with variations measured for rift sediments in the Lake Superior region (Table 1). A higher density within the expected range of densities (2520 kg/m³) is assigned to the wedge-shaped rift sed body on the south that thickens northward toward Superior Shoal. This body is intended to generally represent the sedimentary wedge interpreted from line A (Figure ⁶) as well as overlying sediments next to the fault at the IR-SS anomaly. The dipping layer labeled “Sed + volcanic” on Figure ¹⁰ represents the interbedded volcanic and sedimentary section interpreted from line A in conjunction with magnetic depth estimates, discussed in the “Initial interpretation of line A” section. The inference of magmatic material in this package prompted an assignment of a moderate magnetic susceptibility.

The apparent mismatch of approximately 3 mgal between the calculated gravity curve and the observed curve extracted from the Bouguer gravity grid between the 75 and 80 km distances can be explained by incomplete data coverage over Superior Shoal (Figure ¹⁰). Shallow basalt with moderate to high density at this location, surrounded by lower density rocks to the north and south, is dictated by the geology and velocities in the refraction model of Shay and Tréhu (1993) (Figure ⁵).

The final step toward an integrated interpretation is to apply the proof-of-concept gravity and magnetic model of Figure ¹⁰ to a geologic interpretation of seismic line A, shown in Figure ¹¹. An overall picture emerges of a normal fault in which basalts only from the older, plateau stage are preserved on the footwall. We drew the basalt wedge so that it thins out against the sharp reflection with a negative impedance contrast that aligns with a velocity interface modeled from the refraction profile from sonobuoy 30 (Halls and West, 1971) farther north (Figure ⁷). This is consistent with our interpretation in the “Initial interpretation of seismic line A” section that this reflection represents the top of the preexisting sedimentary basin. The basalt overlying the sedimentary section explains the negative impedance contrast. The shallowest strong reflection within this basin is likely a carbonate and/or chert layer or layers within the prerift sedimentary sequence. Many of the reflections below this layer we interpret as sheet intrusions related to the sill complex discussed for the northern part of line A. Assuming a simple picture of parallel normal faults, we interpret the carbonate/chert layer as extending southward into the seismic disrupted zone, where data problems, intrusions, or a combination of these may explain why the layer is not imaged in the seismic data.

As pointed out earlier, the seismic disrupted zone and IR-SS anomaly undoubtedly occupy the site of a major crustal discontinuity. In the model of flight line 11110 (Figure ¹⁰), we represent the discontinuity with a single normal fault. However, aeromagnetic patterns suggest that the discontinuity is more complicated where line A crosses the IR-SS anomaly (Figure ⁸). At line A, the IR-SS anomaly separates into two high-low pairs that overlap in an en echelon fashion. Because of the 3D characteristics of the en echelon patterns, we do not attempt to model the aeromagnetic data over seismic line A. Instead, we qualitatively infer that the compound anomalies represent multiple displacements within a zone of normal faulting (Figure ¹¹). Sheet intrusions also may be present within the seismic disrupted zone, as suggested by some of the upturned reflections between distances 75 and 82 km and the high refraction velocities in this zone (Figure ⁵).

Our final integrated interpretation is supported or permitted by all of the geophysical data and geologic constraints available to us and thus presents a persuasive argument that our alternate model for the northern margin of the MRS is viable. However, some uncertainty and ambiguities remain. These uncertainties and several geologic implications presented by the new ideas are discussed in this section.

Despite integration of multiple geophysical data sets and geologic constraints, the final gravity and magnetic model of Figure ¹⁰ is still ambiguous. Given that the thickness and physical properties of the modeled rock layers are not independent variables, both can be varied and still obtain a good model fit. However, the variations must still meet the overall constraints laid out in the previous sections. For example, densities on the south side of the model must have overall greater value than those on the north side to fit the gravity curve and agree with overall differences in observed P-wave velocities (Figure ⁵). Thickness of interpreted rift-stage and plateau-stage basalt layers on the south side of the section could be varied, as long as they agree with interpreted reflections observed from seismic reflection data.

Although the simple magnetic model tests and evidence for shallow basalt at Superior Shoal suggest a general configuration of geometries as shown in Figure ⁹, detailed variations of this configuration are still viable. For example, the geometry of the body representing strong, reversed-polarity rocks near Superior Shoal (marked R-tilt15°N under Superior Shoal on Figures ⁹ and ¹⁰) could be modeled differently depending on how one interprets its geologic origin. Although the overall northerly dip of this body and its shallow upward extent are well constrained, empirical tests indicate that the model is highly sensitive to slight variations of thickness or smoothness of the lower contact on the south side of this body. In addition, the presence of normal-polarity rocks measured by Manson and Halls (1991) near the top of Superior Shoal adds a degree of uncertainty. We have assumed that the magnetic effect of these rocks is negligible for purposes of modeling, but they may have more influence than anticipated.

Deep structures below the model extent, such as the lower crust and Moho, probably produce a smooth, regional gravity gradient from high in the north to low in the south over the model area, based on seismic refraction models for line A (Tréhu et al., 1991; Shay and Tréhu, 1993). However, these authors used the gravity data to constrain their interpretations on the north side of line A, which is heavily influenced by a large gravity high off axis to the seismic line (the gravity high that peaks approximately 20 km south of “NB” on Figure ³). As discussed by Thomas and Teskey (1994), gravity effects of the lower crust and Moho can be fit by reasonable variations in densities without significant changes to the geometry of modeled bodies in the upper crust. They suggest that the effects of the deeper section are difficult to untangle from those of the natural variability of densities in the rift section. Thus, we do not attempt to remove a regional field before modeling in the present study. In a way, the dense and magnetic bodies near the bottom of the model of Figure ¹⁰ (sill complex on the north and pervasive intrusions? on the south) may be viewed as the primary drivers of the regional magnetic and gravity fields across the model. Because these bodies are each intended to collectively represent multiple intrusions of uncertain number and extent, their shapes and properties are subject to alternate interpretations.

Having only legacy data available for modeling, no borehole information, and inadequate data coverage present the largest uncertainties in the results. Potential problems with side-swipe and misrepresentations in interpreting 2D seismic data may be large, given the complicated patterns indicated by the potential-field data. The wide spacing of gravity stations and flight-line magnetic data prevent accurate representation of anomalies related to shallow bathymetric features. These problems can only be solved with additional data collection.

Although a fault that places the older, plateau-stage (reversed-polarity) basalts at higher elevations on the north compared to the south accounts for the IR-SS anomaly at Superior Shoal, this representation does not fit the geology at Isle Royale, which begs the question of how to explain the extension of the IR-SS anomaly from Superior Shoal to Isle Royale. Industry seismic lines LS-8 and LS-25 suggest that the younger, normal-polarity rocks ramp up toward the IR-SS anomaly (McGinnis and Mudrey, 2003), although neither line crosses the anomaly (Figure ³). These relations are similar to Isle Royale, suggesting that, if the normal fault extends westward from Superior Shoal, it might be overridden by the thrusting of rift-stage basalts or have variable amounts of throw. In addition, preliminary magnetic modeling at different places along the anomaly suggests that intrusions as well as normal faults are required to properly fit the anomaly shapes. Perhaps a series of normal faults along the extent of the IR-SS anomaly that involve reversed-polarity plateau-stage basalts are locally overprinted by variable degrees of reverse faulting and/or intrusions that are contemporaneous and/or later than normal faulting.

Geophysical evidence supporting the absence of mafic volcanic rocks within the top 2–3 km of crust north of Superior Shoal from refraction models, magnetic depth estimates, and heat-flow data appears problematic considering geologic evidence from the northern shore. Reversed-polarity basalts (Osler volcanics) exposed near Black and Nipigon Bays (Figure ³) are as much as 3 km thick where they dip into the lake (McIlwaine and Wallace, 1976; Hollings et al., 2007; Swanson-Hysell et al., 2014). Several flows with normal-polarity magnetization lie at the top of the reversed-polarity Osler section and also dip into the lake (Halls, 1974). They total approximately 200 m in thickness and are correlated with the much younger rift-stage basalts using paleomagnetic arguments (Swanson-Hysell et al., 2019; Figure ²). Geochemically, they are unlike the normal-polarity rocks at Superior Shoal (Grauch and Schulz, 2019). Thus, thick sections of reversed- and normal-polarity Osler basalts likely underlie much of northernmost Lake Superior, as is first postulated by Halls (1972). Their absence along the northern part of line A and the presence of proximal gravity and magnetic lows (Figure ³) suggest a structural high composed of prerift crust, a hypothesis also first presented by Halls (1972). If so, the body labeled prerift crust in Figure ¹⁰ is likely composed mainly of Archean basement, upon which the postulated chart/carbonate layers were directly deposited well before rifting began.

Why older, plateau-stage basalts are locally preserved at Superior Shoal is unclear. One explanation is that the basalts represent the top of a plateau-stage shield volcano that was tectonically tilted to the north after eruption had ceased. This origin might explain the wedge shape that is preferred by the magnetic modeling and the somewhat limited, elongate bathymetric form of Superior Shoal that has withstood a billion years of erosion and glaciation.

Integral to our alternate view of line A is the interpretation of a preexisting sedimentary basin intruded by a sill complex on the northern part of line A (Figure ⁷). Based on our interpretation of plateau-stage lavas occurring high up in the seismic section, this sedimentary basin would have to be older than the MRS. Two candidates are represented by exposures on the Ontario shores north of the lake (Figure ³): the Paleoproterozoic Animikie Basin and the early Mesoproterozoic Sibley Basin. Stratigraphic sequences for both basins contain significant intervals of carbonate and chert (Ojakangas et al., 2001b; Rogala et al., 2007), which could explain the NL reflection (Figure ⁷). They are both intruded by numerous sills of somewhat differing chemistry and ages during plateau-stage magmatism (Hollings et al., 2010). The sills are aerially extensive and generally follow sedimentary layering but can also show field evidence of saucer shapes and dikes that may feed the sills (Hart and MacDonald, 2007; Hollings et al., 2010).

Although the proximity, stratigraphy, and presence of sills in Animikie and Sibley Group rocks suggest that these basins may be related to the postulated sedimentary basin that we interpret below northern Lake Superior, it is uncertain whether sedimentary deposits of either basin were likely to have extended or be preserved to the southeast under the lake. Sedimentary rocks from both basins have been thinned or never were deposited on the Archean basement of the northern lake shore (Santaguida, 2001a, 2001b). However, the Michigamme Formation in northern Michigan and the Virginia/Thomson Formations in Minnesota are equivalents to the Rove Formation of the Animikie Group in Ontario, indicating a pre-MRS spatially extensive basin (Ojakangas et al., 2001b). The Sibley Group has a much smaller footprint because it is largely confined by faults or hidden by the lake. Its original extent is unknown, but evidence of a significant change in provenance and depositional environment in the middle of the stratigraphic section tempts speculation that the Sibley Basin may have greater tectonic importance than was previously realized (Rogala et al., 2007). Alternatively, the deeper parts of the sill complex may intrude crystalline rock instead of a sedimentary basin, as was observed for some crustal-scale sills in the North Atlantic and Fennoscandia (White et al., 2008; Buntin et al., 2019).

The MRS is a 1.1 Ga sequence of voluminous basaltic eruptions and multiple intrusions followed by extensive sedimentation that extends across the Midcontinent and northern Great Lakes region of North America. Hidden by Lake Superior, tens of kilometers of basalts overlain by kilometers of sediments are mostly known from extensive geophysical studies tied to geologic mapping onshore. Previous workers have commonly used GLIMPCE seismic-reflection line A to demonstrate that the northern margin of the rift developed as a normal growth fault that was structurally inverted to a reverse fault during the compressional Grenville orogeny after rifting had ended. Key to this previous interpretation is the assumption that a widespread seismic horizon with high acoustic impedance contrast (NL reflection) on the north side of the fault is the top of the basalt section on the footwall of the reverse fault. A prominent, curvilinear aeromagnetic anomaly (the IR-SS anomaly) is commonly presented in the literature as a manifestation of this reverse fault. It extends northeastward from the northeastern tip of Isle Royale, Michigan, to Superior Shoal, a bathymetric high in the central lake. However, previous models of a reverse fault where the IR-SS anomaly intersects line A are inconsistent and do not integrate all of the available data.

To develop a more integrated interpretation of the fault along line A and the IR-SS anomaly, we first analyzed and investigated constraints provided by seismic-reflection, seismic-refraction, magnetic, gravity, and physical-property data and by geologic concepts. In particular, the MRS basalts present a unique opportunity for aeromagnetic interpretation because of their strong remanent magnetizations that chronicle a change in the earth’s magnetic field from reverse to normal polarity that generally corresponds to the transition from plateau-stage to rift-stage magmatism. From this initial work, we developed constraints for constructing 2D potential-field models and tested concepts of reversed versus normal faults. We then constructed a more complete potential-field model and applied the results to interpret the seismic section for line A.

The final interpretation largely satisfies all data sets and leads to a profoundly different concept of the rift margin crossed by line A. Instead of a steeply dipping reverse fault with a thick section of late rift sediments on the north, we propose that the margin is principally a normal fault, with a prerift sedimentary basin on the north. Several outcomes related to this alternate model are that (1) the older, reverse-polarity (plateau-stage) basalts are present at shallow levels north of the fault in a north-dipping wedge, (2) the NL reflection represents a sedimentary interval within the prerift basin rather than the top of rift basalts, and (3) the IR-SS anomaly has multiple geologic sources and is not a simple northeastward extension of a reverse fault at Isle Royale. Moreover, we recognize seismic patterns consistent with an extensive sill complex that intrudes the prerift basin. Two candidates for the prerift basin under the northern lake are represented by sedimentary sections and related sills exposed on the northern shore: the Paleoproterozoic Animikie and early Mesoproterozoic Sibley basins.

Our new model indicates that a normal fault is the dominant structure at the northern rift margin along line A, as opposed to the original rift-margin paradigm, which asserts that compressional structures are the dominant features preserved today. Apparently, the northern rift margin here has not experienced as much compression as previously thought. Although the normal fault may be locally modified by intrusion and/or compression along strike away from line A, this alternate model has implications for interpreting the style, scale, and timing of extension, rift-related intrusion, and compression during development of the MRS.

We greatly appreciate the constructive comments and reviews from I. Filina, I. Guerra, S. Haines, M. Hoggett, J. Phillips, and an anonymous reviewer, all of which have improved the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.

All data associated with this research is publicly available and access is cited within the text of this manuscript.

V. J. S. (“Tien”) Grauch received a B.A. (1975) in geology from Carleton College and a Ph.D. (1986) in geophysics from the Colorado School of Mines. She has been employed by the U.S. Geological Survey in Denver, Colorado since 1977, where she is currently a senior research geophysicist. Her current research interests include integration of geologic and geophysical data to understand 3D geology and interpretation of high-resolution aeromagnetic surveys. She is a member of SEG, ASEG, AGU, and a Fellow of the Geological Society of America (GSA).

Eric Anderson received a B.A. (1997) in geology from Augustana College, an M.Sc. (2003) in geology from the University of Colorado, Boulder, and a Ph.D. (2013) in geology from Colorado School of Mines. He has been employed by the U.S. Geological Survey in Denver since 2001, where he currently is a research geophysicist. His research interests include the integration of geology and geophysical data to understand the 3D framework of ore deposit formation. He is a member of SEG, the Society of Economic Geologists, the Society for Geology Applied to Mineral Deposits, and the Geological Society of America.

Samuel Heller received a B.S. (2001) in geology from Brigham Young University and an M.S. (2005) in geology from California State University, East Bay. He was employed by Lawrence Livermore National Laboratory (2001–2007) as a scientific associate and WesternGeco (Schlumberger) (2007–2015) as a geophysicist. He is now employed by the U.S. Geological Survey (2015–present) as a geophysicist where he processes land and marine reflection seismic projects in time and depth. He has a strong background in tomography, advanced migration, earth model building, anisotropy, AVO, noise attenuation, geometry, and statics.

Esther K. Stewart received a B.A. (2005) in geology from Mount Holyoke College and an M.Sc. (2008) in geology from Idaho State University. She worked for ExxonMobil in Houston as a petroleum geologist from 2008 to 2012 before joining the Wisconsin Geological and Natural History Survey, UW-Madison Division of Extension in 2013. In her current role as a Precambrian geologist, she focuses on mapping and interpreting Wisconsin’s Precambrian and Paleozoic geology. She is a member of Epsilon Sigma Phi and the Geological Society of America.

Laurel Woodruff received a B.S. (1973) in geology from the University of Michigan, an M.S. (1977) in geology from Michigan Technological University, and a Ph.D. (1986) in geology from the University of Chicago. She has been employed by the U.S. Geological Survey for 37 years and currently is a research geologist in St. Paul, Minnesota. Her current research interests include the tectonics, metallogeny, and geochemistry of the MRS.