Archean Crust

The Lake Superior region is largely underlain by Archean rocks of the Superior Province, except in the Penokean orogen in Michigan and Wisconsin and in areas where the crust was intruded by magma of the Midcontinent Rift System (fig. 1). The Superior Province is divided into subprovinces based on distinctive rock types and tectonic setting (Card, 1990; Percival and others, 2012). Located in the Lake Superior region are the Wawa-Abitibi, Minnesota River Valley, Quetico, and Wabigoon subprovinces (fig. 1).

The Wawa-Abitibi subprovince in Ontario, Canada, and northern Minnesota is a granite-greenstone terrane, comprising several distinct greenstone belts contained within a widespread region of metamorphosed granitic plutonic rocks (Williams and others, 1991; Jirsa and others, 2012). The granitic rocks range in composition from metamorphosed tonalite to granodiorite, are foliated to gneissic, and are commonly intruded by more potassic, granitoid rocks. The greenstone belts are composed of metamorphosed volcanic packages and associated metasedimentary rocks.

North of the Great Lakes tectonic zone in Michigan (fig. 1), Archean granitic rocks containing greenstone belts are correlated to the Wawa-Abitibi subprovince in Ontario (Sims and Day, 1993; Southwick, 1996; Bickford and others, 2007). Between the Niagara fault and the Great Lakes tectonic zone, exposures of Archean granitic gneiss are characteristic of the Minnesota River Valley subprovince (Sims and Day, 1993; Sims, 1996). The principal part of the Minnesota River Valley subprovince is located outside of the study area in southwestern Minnesota (Sims, 1996; Jirsa and others, 2012).

The Quetico subprovince in Ontario and northern Minnesota is predominantly metasedimentary. In Ontario, the major rock type is feldspathic metagraywacke, with lesser volumes of metamorphosed iron formation, conglomerate, and ultramafic graywacke and siltstone (Williams, 1991). In northern Minnesota, graywacke metamorphosed to biotite-plagioclase schist is accompanied by granitoid intrusions and migmatite (Jirsa and others, 2012).

The Wabigoon subprovince, divided in the study area into the Winnipeg River and Marmion terranes by Percival and others (2012) on the basis of age, consists mainly of tonalitic gneiss interspersed with greenstone belts. Rock types are similar to those of the Wawa-Abitibi subprovince.

Paleoproterozoic Rocks

The Niagara fault zone (labeled “NFZ” in fig. 1) represents the north-verging collision zone at the southern margin of the Superior Province that developed during the Paleoproterozoic Penokean orogeny (Sims and Schulz, 1996; Schulz and Cannon, 2007). The fault zone divides the Wisconsin magmatic terranes on the south from a domain of deformed Paleoproterozoic supracrustal rocks overlying Archean crust to the north (fig. 1, Schulz and Cannon, 2007; Drenth and others, 2021). Within the study area, the Wisconsin magmatic terranes consist mainly of volcanic and plutonic rocks of oceanic island arc affinity (Sims and Schulz, 1996; Schulz and Cannon, 2007). Rocks north of the Niagara fault zone record multiple episodes of sedimentary deposition and deformation that reflect an active continental margin tectonic regime (Schulz and Cannon, 2007; Craddock and others, 2013; Cannon and others, 2018, 2019; Drenth and others, 2021). The section includes graywackes, shales, argillites, quartzites, dolomite, iron formations, and lesser volumes of volcanic and intrusive rocks (Ojakangas and others, 2001b, Schulz and Cannon, 2007). Although the banded iron formations, for which the Lake Superior region is famous, cover a wide region (fig. 3), they are volumetrically small compared to the other Paleoproterozoic rocks.

Thick successions of unmetamorphosed to metamorphosed sandstones, graywackes, and shales form the upper part and greatest volume of the Paleoproterozoic sedimentary section (Morey, 1996; Ojakangas and others, 2001b; Schulz and Cannon, 2007; Craddock and others, 2013). These rocks include the Animikie Group (which is preserved southwest and northeast of the Duluth Complex in Minnesota and Ontario) and metasedimentary rocks in the Marquette Range Supergroup in Michigan and Wisconsin (fig. 3; Schulz and Cannon, 2007). These rocks were once part of the same Penokean basin system that extended over a large region that includes present-day Lake Superior (Schulz and Cannon, 2007; Craddock and others, 2013).

Several episodes of post-Penokean collisional events affected the region before the onset of magmatism related to the Midcontinent Rift System (Schulz and Cannon, 2007; Drenth and others, 2021). Both Archean and Paleoproterozoic rocks were metamorphosed and faulted.

Mesoproterozoic Sibley Group

The lower Mesoproterozoic Sibley Group is a roughly 1-km thick package of mainly siliciclastic sedimentary rocks exposed in the Nipigon region in the northern part of the study area (figs. 1, 3). Only about 300 meters (m) of the basal Sibley Group is preserved on land; a much thicker upper portion was discovered through drilling in Nipigon Bay (fig. 3; Rogala and others, 2007). The lower half of the section is dominated by fine-grained clastic rocks deposited in lacustrine and fluvial environments, whereas the upper half is dominated by sandstone deposited in eolian and fluvial environments (Rogala and others, 2007). Carbonate units are present near the base but compose only a small percentage of the section.

The Sibley Group has a small exposed footprint and is largely confined by faults or hidden by water. Its original extent is unknown, but evidence of a significant change in provenance and depositional environment in the middle of the stratigraphic section implies that these rocks have greater tectonic importance than previously realized (Rogala and others, 2007). Moreover, recent seismic interpretation indicating a sedimentary basin intruded by sills in northern Lake Superior (Grauch and others, 2020) opens the possibility that a greater thickness of the Sibley Group, possibly overlying Paleoproterozoic sedimentary rocks, is preserved under the lake.

Keweenawan Supergroup and Related Igneous Rocks

The Keweenawan Supergroup refers to layered volcanic and sedimentary rocks that formed as part of the 1.1 Ga Midcontinent Rift System and lie unconformably over much older Proterozoic and Archean rocks (fig. 4). Following the terminology adopted by Woodruff and others (2020), the history of the Midcontinent Rift System can be visualized as five transitional stages: initial widespread, plateau-stage eruptions over basal sandstones; voluminous basaltic magmatism of the rift stage; a transition between magmatism and clastic sedimentation of the late stage; thick accumulations of sediments of the post-rift stage of thermal subsidence; and a final transition to a compressional stage as sedimentation waned (Cannon and Hinze, 1992; Hinze and others, 1997; Ojakangas and others, 2001a; Woodruff and others, 2020). Ojakangas and others (2001a) consider the basal sandstones and quartzites as a record of the initiation of rifting, but they are relatively thin (roughly 100 m) across the region, so they are not considered further in this study.

Phanerozoic Rocks

Cambrian and Ordovician sedimentary rocks overlie the older rocks of interest for modeling of the Midcontinent Rift System in the eastern Upper Peninsula of Michigan (fig. 1 and fig. 3). About 200 m of these Paleozoic rocks were encountered in two boreholes near Munising, Michigan (Ojakangas and Dickas, 2002); the rocks consist of Ordovician sandy dolomite overlying friable Cambrian quartz arenite. No bedrock younger than Proterozoic has been found within the lake (Johnson, 1980). Examination of the sources of till in sediment cores suggests that bedrock underlying Lake Superior is mainly composed of Mesoproterozoic Jacobsville Sandstone and Bayfield Group (Dell, 1975; Wold and others, 1982).

Pleistocene and Holocene surficial deposits cover most of the Lake Superior region. These deposits are typically less than 60 m thick onshore, but are more than 100 m thick in many areas of Lake Superior and reach more than 500 m thickness adjacent to the northwestern shore (Wold and others, 1982; Soller, 1998; Soller and others, 2012). Onshore deposits consist mainly of sandy glacial till except on the Bayfield Peninsula, Wisconsin, and in the area of the Ashland syncline (fig. 3), where stratified coarse-grained deposits predominate (Soller, 1992). Lake Superior sediments are characterized by sandy or clayey tills overlain by as much as 30–60 m of glaciolacustrine and postglacial sediments, which vary from interbedded sand, silt, and clay to mostly clay (Lineback and others, 1979; Wold and others, 1982). The tills fill the deep troughs. The deepest troughs (more than 200 m) extend from the north side of Isle Royale all along the north shore of Minnesota (Johnson, 1980; Wold and others, 1982).

Empirical Velocity-Density Relations

Several studies have used large numbers of velocity and density measurements of rocks at various depths to solve for equations that give the best fit between the two variables. Advantages to this approach are that (1) the natural variability of mineralogy and other physical parameters are accommodated, (2) there is no need to separately consider variations of the physical properties with depth, and (3) the relations can be used in similar tectonic terranes where samples have not been collected. The disadvantage is that one empirical curve is commonly used to represent velocity versus density, which implies that velocity and density have a one-to-one correspondence that does not adequately represent the wide scatter of the data.

The most widely used velocity-density relations for sediments and sedimentary rocks are the Nafe-Drake curve (Talwani and others, 1959; Nafe and Drake, 1963; Ludwig and others, 1970) and Gardner’s relation (Gardner and others, 1974). The Nafe-Drake curve was originally presented as a polyline drawn to fit through a plot of a large number of widely scattered, measured velocity-density points from marine sediments and sedimentary rocks. The data were supplemented to include metamorphic and igneous rock types by Ludwig and others (1970). The data and the polyline fit to the data were originally in graphical form, later tabulated by others (for example, Barton, 1986), and more recently characterized by way of two different equations as follows. One equation is by Onizawa and others (2002) for velocities between 3 and 6 km/s. Their equation, recast with density as a function of velocity, is shown as equation 1 and is designated as the Nafe-Drake curve of Onizawa and others.

The other equation is by Brocher (2005) for velocities between 1.5 and 8.5 km/s. It is designated as the Nafe-Drake curve of Brocher and is shown as equation 2.

Gardner’s relation is based on measurements of rocks in sedimentary basins (Gardner and others, 1974). It is given in equation 3.

This equation is valid for velocities between 1.5 and 6.1 km/s. Castagna and others (1993) refined Gardner’s relation by promoting use of equations based on different lithologic types. The equations for shale and sandstone are applicable to the Lake Superior region. The equation for shale gives density as shown in equation 4.

The equation for sandstone gives density as shown in equation 5.

Both equations are valid for velocities between 1.5 and 5.0 km/s.

Velocity-density relations for oceanic basalts were investigated by Christensen and Salisbury (1975) using a large quantity of dredge and core samples from various sites within ocean basins. In particular, they tested the velocity-density relations from 77 tholeiitic basalt samples obtained from the Deep-Sea Drilling Project (DSDP), which represent the composition of the top 1 km of the igneous oceanic crust. They did not provide the range of velocities for these 77 samples. Equation 6 shows their linear density equation, derived by comparing velocities measured at 100 megapascals (MPa) to saturated bulk density for the DSDP samples.

Velocity-density relations for Icelandic basalts differ somewhat from those for the oceanic basalts, which were collected closer to the Mid-Atlantic Ridge. Icelandic basalts were examined by Christensen and Wilkens (1982) from more than a hundred samples of tholeiitic basalt flows and tholeiitic, basaltic to andesitic dikes from 1.9 km of core. They found that a regression equation relating compressional velocity measurements at 100 MPa to saturated bulk density gave equation 7.

This equation is valid for the velocity range of 3.6 to 6.7 km/s. Velocity-density data for these flows and dikes show considerable overlap, so this equation can apply to both basalt and diabase.

The Nafe-Drake curve is commonly used to represent crystalline crust (as well as its sedimentary cover) in conjunction with continent-scale earth models, including some for North America (Magistrale and others, 2000; Musacchio and others, 2004). Brocher (2005) judged the Nafe-Drake curve (eq. 2) as representative of velocity-density relations for continental crust by comparing data from northern California and considered the results to be widely applicable. Barton (1986) discussed the shortcomings of using the Nafe-Drake curve for developing gravity models in general. Maceira and Ammon (2009) used the Nafe-Drake curve of Onizawa and others (eq. 1) to represent sedimentary rocks and sediments of low density. Additionally, Maceira and Ammon (2009) used a form of Birch’s law to represent denser and deeper rocks for input to a joint inversion model of sedimentary basins underlain by Proterozoic and Archean rocks in central Asia.

Velocity-Density Relations from Models

Another approach to developing velocity-density equations is based on models that incorporate laboratory measurements of specific rock types (Christensen and Fountain, 1975; Christensen and Mooney, 1995). Velocity-density equations are developed by examining the data for specific rock types that are considered to represent a particular crustal depth or tectonic setting. Christensen and Mooney (1995) compiled a large set of velocity and density measurements for 29 rock types. These measurements were calibrated to different depths and temperatures; the authors then fit regression curves to the average values of various combinations of these rock types to represent continental crust in different settings. Their linear equation developed for crystalline crust that excludes volcanic and monomineralic rock types and assumes an average crustal heat flow at 10 km depth is the most widely used equation among contemporary geophysicists to represent the velocity and density relations within the crust (Brocher, 2005). It is designated as the global crustal line and is shown as equation 8.

The equation is valid for velocities between 5.5 and 7.5 km/s.

In developing the crustal model LITH05 for the Canadian shield regions, Perry and others (2002) adopted a crustal model of average petrology proposed by Christensen and Mooney (1995) but the model was not reported with an associated equation. Assuming a mixture of five different rock types that varies with depth (fig. 5), Christensen and Mooney (1995) used their laboratory results for each rock type to develop velocity and density curves plotted against depth. Plotting velocity versus density yields a nonlinear velocity-density relation (fig. 5B), which is estimated from a second-order polynomial regression. This regression, shown in equation 9, represents velocity-density relations for a crust of average petrology.

Velocity-density relations derived from rock-type assumptions have been used in conjunction with other approaches to develop continent-scale or global crustal models. Mooney and Kaban (2010) used the velocity-density equations of Christensen and Mooney (1995) to compute the gravity effects of crystalline crust from velocities compiled from seismic studies. Musacchio and others (2004) used a combination of measured surface densities, Nafe-Drake relations, and a linear equation developed for granulite facies rocks (Christensen and Fountain, 1975) to find the best-fit crustal model for seismic and gravity data in the Superior Province north of Lake Superior. Nitescu and others (2006) used both the Nafe-Drake curve and the linear equation of Christensen and Mooney (1995) for crystalline crust (the global crustal line of eq. 8) to investigate density contrasts predicted by velocity discontinuities found from seismic data in the Superior Province.

An additional approach to investigating density and velocity relations involves using petrophysical modeling to directly calculate the physical properties from presumed mineralogy and estimated pressure-temperature conditions (Tassara, 2006; Sowers and Boyd, 2019; Boyd, 2020). In this approach, mineral assemblages and their physical properties are calculated by rock type at specified pressure and temperature conditions. Sowers and Boyd (2019) used velocities and densities computed from mineral assemblages representing several common rock types at standard temperature and pressure (0 °C and 0.1 MPa) to compare the Nafe-Drake curve of Brocher (eq. 2). They found that their results were consistent with Brocher’s (2005) assessment of the scatter exhibited by field and laboratory measurements about the Nafe-Drake curve. The advantage of both model-based approaches is that physical properties can be predicted at depths where they cannot be measured. The disadvantage is the subjectivity of the models themselves.

Laboratory Studies of Seismic Velocity

Velocity is commonly determined in the laboratory by measuring the time it takes an input pulse to travel through a rock core in comparison to the travel time of a similar pulse through a standard medium, such as mercury (Birch, 1960). The temperature and confining pressure—or both—of fluid surrounding the rock core is varied to simulate conditions at depths that are appropriate for understanding seismic surveys of the Earth (Christensen and Mooney, 1995). The reader is referred to the original studies for detailed descriptions of the procedures and measurement uncertainties.

Velocity and density were measured from hand and core samples representing different rock types in the Lake Superior region by Halls (1969), Lippus (1988), and Allen (1994). Halls (1969) included measurements from a wide spectrum of sedimentary and igneous rock types as well as analysis of velocity-density and velocity-porosity relations and velocity variations with depth. The author considered the precision of the velocity measurements for all the rock types analyzed to be on the order of 1 percent for pressures over 100 MPa. Lippus (1988) measured velocities and densities for basalt samples from boreholes located in the northern Keweenaw Peninsula (fig. 3) to investigate variations in velocity and velocity-density relations with stratigraphic position. From hysteresis analysis, the author showed a typical reproducibility error on the order of 0.1 km/s for the lowest pressure (20 MPa) and errors less than 0.01 km/s for the highest pressures (200 and 220 MPa). Allen (1994) measured velocities and densities of core samples from a deep natural gas test well (Minnesota unique well number 207022), which was drilled through a thick section of Keweenawan basalt near the city of Osseo in southeastern Minnesota (Allen and Chandler, 1993). The original velocity and density values plotted by Allen (1994) were not available, so the data are reproduced from digitized figures from the author’s thesis. The samples were measured from the interiors of basalt flows at approximate depths of 300–1,220 m in the borehole. Allen considered the velocities (reported only for a confining pressure of 200 MPa) to represent basalt velocities at depth below Lake Superior (Allen and others, 1997). Velocities and densities for a few samples of mafic igneous rocks from the Lake Superior region were also included in early rock-property studies by Hughes and Maurette (1957) and Birch (1960). As for most of the early studies, details and location information about the rock samples are lacking.

All the laboratory studies that included samples from the Lake Superior region (except the one by Allen, 1994) systematically measured velocities at intervals of increasing confining pressures to simulate conditions at different depths. However, densities were only measured at atmospheric conditions. This practice does not properly sample the relations between the physical properties at greater pressures. The earlier laboratory studies also did not account for changes in sample volume caused by increasing pressure (Christensen and Mooney, 1995). These shortcomings represent a significant problem for measuring porous clastic rocks, which compact with increasing pressure, but the studies’ oversights are not as important for measurements of less porous crystalline rocks (Schön, 2015b). For this reason, velocity measurements of sedimentary rocks are discussed further here, whereas laboratory results for velocity and density of igneous rocks are considered together in the section “Velocity-Density Relations for the Lake Superior Region.”

Plots of velocity versus confining pressure allow for investigations of uncertainty and variability of the data. Such plots for Keweenawan sedimentary rocks from the Halls study (fig. 6) show moderately well-behaved correlations and consistent variations between samples for pressures greater than 100 MPa. The sharp increases and inconsistencies of the data observed for pressures less than 50 MPa are likely caused by experimental problems, which were described by Birch (1960) for the laboratory set-up that Halls (1969) used. Velocity curves for three samples are inconsistent with the other sample curves: sample BJ1 of Jacobsville Sandstone (fig. 6A) and mudstone samples CF22 and CF23 of Freda Sandstone (fig. 6B). An explanation for the problem with sample BJ1 for pressures less than 150 MPa is not obvious from Halls (1969, supplemental data therein), but the pattern suggests that the sample was not saturated before measurement, despite how it was reported. If so, the constant differences between measured values for BJ1 and BJ2 at pressures greater than 150 MPa suggest that the results for sample BJ1 may be valid at these pressures. Both samples BJ1 and BJ2 were collected from Caribou Island, Ontario (fig. 3). A likely explanation for the problems with samples CF22 and CF23 is that they were not saturated before measurement, unlike the other mudstone samples (Halls, 1969, his supplemental data). These samples are not considered in subsequent analysis because they are inconsistent with the velocities for the other mudstone samples even at the highest pressures measured. None of the samples of Copper Harbor Conglomerate were saturated before measurement, but the values show consistent, nearly level values at pressures greater than 100 MPa (fig. 6B).

Although the data from figure 6 show a wide range in values for all samples, there are general trends that can be distinguished between groups. For the units overlying the Oronto Group (fig. 6A), well indurated Jacobsville Sandstone samples have the highest velocities, moderately friable sandstone samples that do not have quartz overgrowths have the lowest velocities, and samples intermediate between these two types (except for samples BJ1 and BJ2) have velocities in between. For Oronto Group samples (fig. 6B), arenaceous and argillaceous samples show a significant difference in ranges of velocities measured, although the ranges for both groups are large. These relations are revisited in the “Velocity-Density Relations for the Lake Superior Region” section.

Sonic Log from Terra-Patrick 7–22 Borehole

Only a few deep-exploration boreholes (more than 1.5 km total depth) have been drilled in the region covered by figure 1: the Terra-Patrick 7–22 borehole near Ashland, Wisconsin (Dickas and Mudrey, 1999) and the St. Amour and Hickey Creek boreholes near Munising, Michigan (Ojakangas and Dickas, 2002; fig. 3). Sonic logs are available only for the Terra-Patrick borehole from the Wisconsin Geological and Natural History Survey (WGNHS, 2019), which was drilled through 1,425 m of Oronto Group in the Ashland Syncline located about 19 km west of Ashland, Wisconsin (fig. 3). Although numerous shallower boreholes (less than 1.5 km total depth) have been drilled in the region for mineral exploration and development, sonic logs were typically not acquired.

Compressional-wave velocities were measured for the Terra-Patrick 7–22 borehole using both compensated short-spaced and long-spaced sonic tools. Figure 7 shows the sonic log from the long-spaced tool. With transmitter-to-receiver spacings of 2.4, 3.1, and 3.7 m, the long-spaced tool allows for wider investigation of the borehole walls than can be achieved with the short-spaced tool, which uses spacings of 1.0 and 1.5 m (Lindblom, 1999). The sonic-log data are noisy in places, especially where hole rugosity was extreme, as indicated by the abrupt and wide variations in caliper spread (fig. 7).

The sonic log was recorded from 600 to 1,500 m depth through Freda Sandstone, Nonesuch Formation, and Copper Harbor Conglomerate, as picked by Daniels (1999). Daniels also identified transition zones between the Nonesuch Formation and the overlying and underlying formations. Although all three formations of the Oronto Group were drilled, the top 2.7 km of Freda Sandstone was missing from the section and the borehole penetrated only 260 m of Copper Harbor Conglomerate, which was estimated to be almost twice this thickness from predrilling seismic interpretations (Dickas and Mudrey, 1999).

A median filter was applied to the noisy log data (fig. 7) before computing statistics for each of the geologic units defined by Daniels (1999; table 1). The data are divided following the lithologic division of arenaceous versus argillaceous Oronto Group described in the “Geologic Setting” section. The Freda Sandstone has velocity values that remain fairly close to 3.60 km/s in the top 600 m of the borehole. Below this, velocity gradually increases, reaching approximately 4.50 km/s for the Nonesuch Formation below 1,000 m. Below depths of 1,200 m, the transitional interval and Copper Harbor Conglomerate show large velocity fluctuations, which are likely an expression of the heterogeneity of grain size in these units when compared with the overlying rocks.

Seismic Refraction Studies of Seismic Velocity

Seismic refraction studies examine the arrival times of seismic waves that are refracted at interfaces to determine velocities above and below the interfaces and the depth to the interface. First and second arrivals of refracted waves are recorded by receivers that are located along profiles and separated from shot locations at distances commensurate with depths that are to be investigated. Many seismic refraction surveys have been conducted in the Lake Superior region, reported by Ocola and Meyer (1973) and Perry and others (2002). Several surveys were conducted along long-range profiles extending for hundreds of kilometers, which are best for investigating properties of the crust and mantle. Others used shorter profiles (5–25 km) in comparison to mapped geology to investigate the velocities of rocks within the top 5–15 km of the crust (Thiel, 1956; Steinhart and Meyer, 1961; Mooney and others, 1970a; theses described in Ocola and Meyer, 1973). Some studies used other constraints on the refraction velocities to interpret subsurface geology, such as physical-property information (Halls and West, 1971; Luetgert and Meyer, 1982) or seismic-reflection profiles (Tréhu and others, 1991; Shay and Tréhu, 1993). These multiple studies, conducted over a large span of time using a wide range of technology and a variety of approaches, differ in the uncertainty of their results that is dependent upon details of the study. On the other hand, the largest errors for all the surveys are likely caused by the three-dimensional nature of the subsurface that is difficult to represent using linear refraction profiles.

Several short-range refraction studies conducted in areas where geology is well known provide velocity information about bedrock lying just under cover. In separate studies, Thiel (1956) and Steinhart and Meyer (1961) determined the velocities of Bayfield Group sandstones on the Bayfield Peninsula (fig. 8) from shallow explosions and from examination of short-range recordings from a long-range profile, respectively. Steinhart and Meyer (1961) also examined short-range recordings offshore of the northern tip of the Keweenaw Peninsula (fig. 3), where Copper Harbor Conglomerate dips gently into the lake (Cannon and Nicholson, 2001). Results from the short-range investigations by Steinhart and Meyer (1961) are described in both Halls (1969) and Ocola and Meyer (1973). Ocola and Meyer (1973) also characterized the velocities of geologic units from analysis of several published and unpublished short-range surveys from the Lake Superior region.

Mooney and others (1970a, 1970b) provide the most extensive study of short-range refraction profiles to characterize geologic units by velocity, but their geologic correlations require updates before their data can be included in modern velocity compilations. They analyzed 10 profiles in the Lake Superior region (fig. 8) and 23 more to the southwest, tracking the geophysical expression of the Midcontinent Rift. They relied on pre-1970 geologic maps for their analysis and did not have access to exploration cores drilled in the late 1950s within the Ashland syncline (fig. 8), which are now stored at the WGNHS (WGNHS, 2009).

The work contained in Mooney and others (1970a, 1970b) can be updated by comparing the refraction models to the now-available exploration cores and modern geologic maps. The locations of 10 refraction-model profiles in the Lake Superior region from their study and the shallowest sub-drift velocities derived from them are presented with respect to mapped geology on figure 8. Figure 9 compares the velocities of the entire depth range of four of these profile models to the lithologic picks for nearby cores and the Terra-Patrick 7–22 well. As described earlier, the Freda Sandstone in this well is generally argillaceous because of the considerable siltstone content. The Freda Sandstone in cores DO–5 and DO–14 is also mainly argillaceous because of the large proportion of fine material that is incorporated with the sandstone, mostly in the form of flaser bedding and mud drapes (E. Stewart, WGNHS, written commun., 2021). The Copper Harbor Conglomerate in well DO–3 is arenaceous. It is composed of medium- to very coarse-grained sandstone with carbonate cement (E. Stewart, WGNHS, written commun., 2021).

Velocities for profiles 16 and 30, previously correlated by Mooney and others (1970b) to the Freda Sandstone, properly match the depths to the top of the Freda Sandstone encountered in nearby cores (fig. 9), but interpretations for profiles 17 and 18 require some revision. The low velocity (3.23 km/s) determined for the top 550 m of bedrock on profile 17 was tentatively interpreted as Orienta Sandstone of the Bayfield Group by Mooney and others (1970b), based on comparisons to low velocities for Orienta Sandstone elsewhere. The location of profile 17 falls in an area of no outcrop inferred as Oronto Group—not Bayfield Group—by White (1966b) and Nicholson and others (2006). If the 3.23 km/s velocity layer corresponds to Freda Sandstone, it is significantly lower than velocities of 3.81–3.84 km/s for profiles 16 and 30 that are well correlated to Freda Sandstone in nearby wells, and also lower than averages of 3.59 and 4.03 km/s found for two intervals of Freda Sandstone in the Terra-Patrick 7–22 well (table 1). A velocity of 3.23 km/s for Freda Sandstone is plausible, however, given the velocities of 3.0–3.5 km/s for this unit at depths shallower than 250 m in the well (fig. 7). Thus, the model for profile 17 below 122 m of overburden could either represent Bayfield Group preserved over folded Freda Sandstone in the core of the Ashland syncline or Freda Sandstone that has lower velocity than expected in the top 550 m. The latter interpretation is preferred because preliminary resistivity images from an airborne electromagnetic survey about 7 km to the northeast, within a deeper part of the Ashland syncline, show a moderately resistive unit overlying a conductive unit at depths of 200–400 m. This relation is consistent with the interpretation of 200–400 m of overburden directly overlying Freda Sandstone (P. Bedrosian, USGS, written commun., 2021).

The high velocity (5.18 km/s) correlated to the Copper Harbor Conglomerate on profile 18 (fig. 9) was previously interpreted to indicate basalt by Mooney and others (1970b). The new correlation is substantiated by Copper Harbor Conglomerate in nearby well DO–3 and the profile’s location within the axis of the Ashland Syncline (fig. 8), where Oronto Group thickens above the basalt. This new velocity interpretation for the Copper Harbor Conglomerate supports the conclusions of Halls (1969), who called attention to overlapping, measured velocity ranges for Copper Harbor Conglomerate and unmetamorphosed basalt flows.

Velocity Data Compiled for Sedimentary Rocks

Velocities are compiled for units overlying the Oronto Group in table 2, divided by formation, if known. Laboratory measurements of velocity for all units overlying the Oronto Group vary widely, which is also demonstrated in the velocity versus pressure plots (fig. 6A). Laboratory measurements for the Bayfield Group are generally higher (more than 3.6 km/s) than those determined from refraction studies (less than 3.6 km/s). The difference does not appear to depend on which formation of the Bayfield Group was sampled (table 2), suggesting that velocities for the Bayfield Group are in the general range of 3–4 km/s but are highly variable. Velocities for the other units overlying the Oronto Group, which are only represented by laboratory measurements, are largely higher than those of the Bayfield Group and are also highly variable. They fall in the general range of 4.0–4.5 km/s, except for well indurated Jacobsville Sandstone, which has a higher range of 4.5–5.0 km/s.

Velocities are compiled for the Oronto Group in table 3, divided between argillaceous and arenaceous rock types. Velocities for the argillaceous and arenaceous groups are generally 3.3–4.5 km/s and 4.5–5.9 km/s, respectively, showing a significant separation in range as was evident on the velocity-pressure plots (fig. 6B). The velocity range of 3.3–4.5 km/s for the argillaceous rocks falls well within the total range determined for all the units overlying the Oronto Group (3–5 km/s). This overlap was noted earlier by Halls (1969).

Velocity measurements for the Paleozoic sandstones and dolomites in the Lake Superior region are not apparent from published data. Typical velocities of Paleozoic sandstones and dolomites range from 2.0 to 6.0 km/s and 2.5 to 6.5 km/s, respectively (Kearey and Brooks, 1991).

Velocities for unconsolidated deposits in the Lake Superior region and analogous regions nearby are summarized in table 4. The velocities vary widely and are mainly dependent on porosity and degree of saturation (Silver and Lineback, 1972; Kearey and Brooks, 1991). Shallow refraction data over unconsolidated deposits from Minnesota (Petersen, 2001) show that materials above the water table have the lowest velocities (on the order of 0.2–1.0 km/s) and those below have the highest (on the order of 1.5–2.3 km/s). At similar saturations, velocities for clay are generally higher than those for sand. Velocities of unconsolidated lake sediments commonly are 4–6 percent faster or slower than the velocity of lake water (Silver and Lineback, 1972; Johnson, 1980). The velocity of Lake Superior water is approximately 1.42 km/s (Halls and West, 1971). Glacial tills generally have somewhat higher velocities, on the order of 1.6–1.8 km/s (table 4). Wold and others (1982) estimated the velocity for unconsolidated sediments in the central part of Lake Superior at 5,000 feet per second, or ~1.5 km/s, which is likely a reasonable representation of the velocity of a combination of lacustrine sediments and glacial till.

Laboratory Studies of Density

Density is determined in the laboratory by dividing the measured mass of the sample by its volume. Mass and volume were determined in a variety of ways, as described in the original references. For porous rocks, such as clastic sedimentary samples, it is important to saturate the rock with water before measuring. This practice simulates the conditions of rocks located below the water table.

Campaigns to measure density from rock samples in the Lake Superior region date as far back as 1896 (Dutch and others, 1995). The most recent are published as databases of numerous density measurements for Minnesota (Chandler and Lively, 2021), the Iron Mountain-Menominee region in Michigan and Wisconsin (Anderson and others, 2020), and specific gravity measurements for Ontario (Rainsford and others, 2018). Specific gravity is the unitless ratio of the rock’s bulk density to that of water, whereas bulk density is an absolute quantity that takes the density of water into account. Assuming a density of 1.0 gram per cubic centimeter (g/cm³) for water, the specific gravity values can be equated to bulk density in units of g/cm³ and then converted to kg/m³ by multiplying by 1,000. Densities from these or preliminary versions of these more recent databases are summarized for specific regions in Minnesota by Chandler (2002) and Chandler and others (2007), in Michigan and Wisconsin by Drenth and others (2021), and in southern Ontario by Musacchio and others (2004).

Most of the samples measured for density were taken from outcrops (Thiel, 1956; Bacon, 1966; Klasner and others, 1979; Dutch and others, 1995; Rainsford and others, 2018; Anderson and others, 2020); some also include samples from core or mines (Birch, 1954; Halls, 1969; Lippus, 1988; Allen, 1994; Chandler and Lively, 2021). Thomas and Teskey (1994) included bedrock samples dredged from the lake bottom. With the exception of the modern compilations (Rainsford and others, 2018; Anderson and others, 2020; Chandler and Lively, 2021), the densities for particular formations are reported without details of the sampling method or sample location.

Density Data Compiled for Sedimentary Rocks

Densities are compiled for units overlying the Oronto Group in table 5, divided by formation and lithology (if known). The measurements for these units are fairly well clustered. The general ranges (rounded to the nearest 50 kg/m³) are 2,150–2,450 kg/m³ for the Bayfield Group and 2,300–2,550 kg/m³ for both the Fond du Lac and Hinckley Formations. The Jacobsville Sandstone samples show the widest variation, with densities of roughly 2,550 kg/m³ for the well-indurated samples and a range of 2,150–2,400 kg/m³ (rounded to the nearest 50 kg/m³) for the saturated samples. Bacon (1966) found an average value of 2,250 kg/m³ for density measurements compiled from multiple studies. However, the author considered this average too low to represent the Jacobsville Sandstone for gravity modeling and used a density of 2,400 kg/m³ instead. If valid, this supposition also brings into question the average value of 2,150 kg/m³ reported by Birch (1954). The low values from Birch’s study might be explained by the failure to saturate the samples before measurement.

Densities are summarized for the Oronto Group in table 6, divided by lithologic type following table 1 if known. The range of average densities for argillaceous rock types is approximately 2,350–2,550 kg/m³. The range for arenaceous rock types is approximately 2,600–2,700 kg/m³ if the study by Birch (1954) is discounted. This omission is justified because (1) the measured values in this study are much lower than the those of the other studies, and (2) the author acknowledged difficulties in obtaining representative samples for measuring the conglomerates.

Density measurements for the Paleozoic sandstones and dolomites in the Lake Superior region are not apparent from published data. Typical densities measured for saturated Paleozoic sandstones and dolomites have ranges of 2,350–2,550 kg/m³ and 2,280–2,900 kg/m³, respectively (Kearey and Brooks, 1991).

Densities for unconsolidated materials are difficult to measure from hand samples. Those reported by Chandler and Lively (2021) for Minnesota are inferred from seismic refraction velocities using Gardner’s relation (eq. 3) rather than measured. Densities reported for saturated sand and clay elsewhere are on the order of 1,500–2,500 kg/m³ (Kearey and Brooks, 1991). Densities obtained for unconsolidated materials in cores from Lake Erie, which may be similar to those of Lake Superior, are generally in the range 1,250–2,100 kg/m³ (Morgan, 1969).

Density Data Compiled for Felsic Igneous Rocks

Only limited density data are available for rhyolites and felsic intrusions related to the Midcontinent Rift System, despite their significance within the Duluth Complex and North Shore Volcanic Group (Green, 2002) and their abundance within volcanic centers in the Porcupine Mountains (Hubbard, 1975; Cannon and others, 1995) and on Michipicoten Island (Annells, 1974). Felsic rocks of unknown origin, which are commonly described as “felsite” in the older literature, compose a significant fraction of the pebbles contained within the Copper Harbor conglomerate (Daniels, 1982), suggesting that these rocks were once widespread across the area and have been eroded away. Although data are sparse, densities compiled in table 7 mainly fall in the range of 2,600–2,800 kg/m³; these values are consistent with densities determined for felsic rock types elsewhere (Kearey and Brooks, 1991).

Density Compiled for Rocks Older than the Midcontinent Rift System

A large number of density measurements for rocks older than the Midcontinent Rift System in the Lake Superior region, mainly on hand samples of igneous and metamorphic rocks, have been reported in the literature and collected in databases. These data are compiled in appendix 1, table 1.1 by rock type and region. Paleoproterozoic iron formations are excluded from the compilation, even though they commonly have high densities (more than 3,100 kg/m³; Klasner and Cannon, 1974; Dutch and others, 1995; Chandler and others, 2007; Drenth and others, 2021). The small volumes of the Paleoproterozoic iron formations relative to other Paleoproterozoic rocks makes them important only for geophysical modeling where they are present onshore. For geophysical models focused on the sections under the lake, their gravity effects would be negligible at the depths they would be expected to occur. Densities for the older metasedimentary and sedimentary rocks are discussed in the section “Velocity-Density Relations for the Lake Superior Region.” Densities for crystalline rocks from table 1.1 show that general variations mainly correlate with silica content. Felsic rocks are the least dense and mafic and ultramafic rocks are the densest, a common observation noted worldwide (Kearey and Brooks, 1991; Schön, 2015a).

Compiling densities by terrane rather than by rock type is a more useful exercise to understand density variations for use with gravity models of the heterogeneous crystalline crust. The modern databases for Ontario (Rainsford and others, 2018), Minnesota (Chandler and Lively, 2021), and a limited area near the Michigan-Wisconsin border (Anderson and others, 2020) contain sufficiently large sample populations to allow such statistical comparisons, shown in figure 10. The three histograms for rock densities for the Superior Province (fig. 10A, B, C) all display wide variability that is somewhat bimodal, with dominant peaks in the ranges 2,620–2,720 kg/m³ and 2,900–3,050 kg/m³ representing felsic and mafic rocks, respectively. The felsic range includes 2,700 kg/m³, which is commonly used for gravity modeling of the Archean crust (Bacon, 1966; Hutchinson and others, 1990; Hinze and others 1990, 1997; Mariano and Hinze, 1994; Thomas and Teskey, 1994; Allen and others, 1997; Chandler and Lively, 1998; and Drenth and others, 2021). The asymmetric, bimodal character of densities measured for the Superior Province in Ontario (fig. 10A) likely reflects the mixture of felsic and mafic rocks of the granite-greenstone terrane. The larger peaks in the felsic versus mafic ranges for the Superior Province in Minnesota, Wisconsin, and Michigan (fig. 10B) indicate a greater predominance of granitic rocks. The histogram for densities sampled from Paleoproterozoic crystalline rocks in Wisconsin and Michigan (fig. 10C), which are mostly from the Wisconsin magmatic terranes, also shows a separation in density ranges between felsic and mafic rocks, with the felsic range dominant. The bimodal and broad distribution of density values displayed by the histograms emphasize the importance of examining the dominant rock types of the terrane that is to be modeled rather than using one average density value for all crystalline crust.

Keweenawan Sedimentary Rocks

The only measurements of both velocity and density for the same samples of Keweenawan sedimentary rocks come from Halls (1969). Velocity data from this study for the Oronto Group and overlying units at confining pressure of 55 MPa are plotted against density in figure 11A. Velocities measured at 55 MPa (to simulate a depth of 1–2 km) from Halls’ study provide the best data (fig. 6) that most closely match the Earth’s surface conditions represented by the density measurements. Matching the measurement pressures as closely as possible is important for sedimentary rocks because both velocity and density vary greatly with depth.

In figure 11, the velocity-density data points from Halls (1969) are compared to the established velocity-density equations commonly used for sedimentary rocks: Gardner’s relation (eq. 3), the two versions of the Nafe-Drake curve (eqs. 1 and 2), and relations for shale and sandstone (eqs. 4 and 5). The shale and sandstone curves are purported to generally bracket the range of velocity-density relations expected for clastic rocks in sedimentary basins (Castagna and others, 1993). These rock types should be comparable to the Keweenawan sedimentary rocks. The scatter of the Halls data in relation to these curves points out the variability and possible problems with the data, issues with using a single curve to represent empirical relations that are variable, and difficulties in using limited data for determining representative velocities and densities and their relations.

The three well indurated Jacobsville Sandstone samples are anomalous compared to the other samples from the units overlying the Oronto Group; they plot closer to samples from the arenaceous Oronto Group (fig. 11A). These samples are all from one location near Sault Ste. Marie, Ontario (Halls, 1969), indicating that (1) the rocks in this location have anomalous lithologies compared to other areas, (2) Jacobsville Sandstone physical properties have much wider variation than the values from the other localities imply, or (3) the rocks should be correlated with the older Oronto Group. These possibilities suggest that velocity and density relations may differ among the units overlying the Oronto Group, but there is not enough data available to examine these differences.

Geophysical analyses from previous work can help address the problem of limited velocity-density data. Table 8 pairs velocities and densities that others have separately considered representative of particular geologic units or lithologies at some depth below Lake Superior to help understand velocity-density relations. For example, Thiel (1956) constructed density-depth curves from surface measurements and theoretical considerations of compaction versus lithology for Bayfield and Oronto Group units to best represent these densities in the subsurface of the lake. These values at depths of 1,500 m for Bayfield Group and 5,000 m for Oronto Group can be compared to velocities of the same lithologies available from Thiel’s refraction results from the same area (fig. 8) and measurements of similar argillaceous Oronto Group rocks from Halls (1969). Other inferred, representative velocity-density pairs were derived from one or more authors who provided the reasoning behind their choices of representative velocity or density for particular rock types. For example, White (1966a) used both detailed geophysical analysis and geologic considerations to derive representative lithologies and related subsurface densities for the Keweenawan sedimentary units. Halls (1969) developed theoretical curves of velocity versus depth for different lithologies found for these same rocks. Both Halls (1969) and White (1966a) reasoned that rocks of Freda Sandstone and Copper Harbor Conglomerate below the lake would likely be dominated by fine-grained sandstone. From this conclusion, an inferred, representative velocity-density pair can be formed from two separate sources. As another example, Allen (1994) provided detailed analysis on his development of velocity-density pairs for converting seismic-reflection profiles to depth (reported with less detail in Allen and others, 1997). Because this analysis included evaluation of velocities and densities measured by others, the results are naturally similar to several other authors.

The inferred, representative velocity-density pairs and ranges plot near the measured data points for their corresponding Keweenawan sedimentary units and within the two endmember curves for sandstone and shale (fig. 11). They provide additional evidence for comparing the observed data to the established velocity-density curves to see which most closely matches each rock suite. Both representations of the Nafe-Drake curve pass near the cluster of data points, inferred ranges, and inferred point for the units overlying the Oronto Group for velocities less than 4.5 km/s (fig. 11A), suggesting that the Nafe-Drake curve best represents this rock suite. The anomalous, well indurated Jacobsville Sandstone data points fall closer to the Nafe-Drake curve of Onizawa and others (2002) and Gardner’s relation at velocities greater than 4.5 km/s. Despite the sparse data at these velocities, the Nafe-Drake curve of Onizawa and others (2002) is chosen to represent this rock suite to maintain consistency between all the units overlying the Oronto Group. However, this choice remains tentative and could change if more data were available for the Jacobsville Sandstone. Gardner’s relation is chosen to represent the Oronto Group because this curve passes closer to the mudstone data compared to the other curves and provides a seamless way to represent the transition from argillaceous to arenaceous Oronto Group rocks (fig. 11B).

Igneous Rocks Related to the Midcontinent Rift System

Previous workers have investigated velocity-density relations sufficiently for basalts, diabase, and gabbro related to the Midcontinent Rift System in the Lake Superior region to present these properties together in this section. These data are compiled along with measurements for other rock types known or suspected to be related to the Midcontinent Rift System in table 9. Associated felsic intrusions may be present at depth, as evidenced by rhyolite flows and rhyolite breccias exposed at the surface (see the “Geologic Setting” section). Although only limited data are available for these rocks, they are included for completeness. Anorthosite may also be present at depth, because large xenolithic blocks of anorthosite are included in diabase related to the Beaver Bay Complex on the north shore of Minnesota (fig. 3; Morrison and others, 1983; Miller and Chandler, 1997). Although a paucity of density data and no velocity data are available for anorthosite in the Lake Superior region, other geophysical observations led Grauch and Heller (2021) to speculate that a considerable volume of these rocks might be present at depth in western Lake Superior. Therefore, velocity and density data for this rock type compiled from elsewhere are included here for comparison.

The laboratory studies of Halls (1969) and Lippus (1988) provide the best data to examine velocity-density relations for basalts of the Keweenawan Supergroup. Halls (1969) includes samples from a variety of basalts around the region (table 1.3), whereas Lippus (1988) thoroughly sampled a section of Portage Lake Volcanics near Eagle Harbor at the northern tip of the Keweenaw Peninsula (fig. 3). Both Halls (1969) and Lippus (1988) reported velocity measurements for individual samples at increasing confining pressures to a maximum of approximately 220 MPa. Density is plotted versus velocity at low (55–60 MPa) and high (193–200 MPa) confining pressures for both datasets (fig. 12) to demonstrate the scatter of the velocity data. As explained earlier, velocity measured at 55 MPa is the most reliable data that most closely matches the conditions of the density measurements. However, comparing density measurements at low pressures to velocity at higher pressures is not as problematic for crystalline rocks as these comparisons are for sedimentary rocks because the densities and velocities of crystalline rocks do not vary as much with depth (Schön, 2015a, b). Moreover, velocity is commonly measured and discussed for a confining pressure of 200 MPa in other rock property studies (Hughes and Maurette, 1957; Birch, 1960; Allen, 1994), which provides the impetus for plotting velocity at this pressure versus density in subsequent figures.

Velocity-density relations for all the basalts measured at 193–200 MPa show a similar relation to equation 11 (fig. 12B). The tighter clustering of the velocity-density points about the Halls line for the higher confining pressure suggests that the line is a suitable representation of relations as deep as 7–8 km below the subsurface, which is estimated considering depth versus overburden pressure (Halls, 1969; Lippus, 1988; Castagna and others, 1993). Thus, for simplicity, the Halls line (eq. 11) is shown in subsequent figures as a proxy for basalts of the Keweenawan Supergroup rather than plotting the individual data points.

Figure 12 shows that the Halls line representing basalts of the Keweenawan Supergroup is similar to the one for oceanic basalts sampled by the DSDP (Christensen and Salisbury, 1975; eq. 6) and has a slope that matches the one for Icelandic basalts (Christensen and Wilkens, 1982; eq. 7). Christensen and Wilkens (1982) attribute the difference between the velocity-density relations from DSDP and Iceland to greater iron content in the basalts from Iceland. Velocity versus density data compiled for basalt by Christensen and Mooney (1995; fig. 12B) plot on the line for the Icelandic basalts and do not capture the variability of basalt data that is demonstrated by Christensen and Wilkens (1982), although they do acknowledge that a wide range of velocities in basalt is expected as a function of glass content, porosity, and secondary alteration.

Velocities and densities for basalts of the Keweenawan Supergroup include higher velocities (more than 6.0 km/s) and densities (more than 2,900 kg/m³) than those reported by Christensen and Mooney (1995). The basalts from the borehole near Osseo (Allen and Chandler, 1993; Allen, 1994) provide an example of these high values (fig. 12B), with measured averages for velocity and density of 6.5 km/s and 2,950 kg/m³, respectively. These samples were obtained from the interiors of flows at about 1 km depth. Allen (1994) and Allen and others (1997) considered these samples to represent the conditions of basalts of the Keweenawan Supergroup in the subsurface because the high porosities of flow tops are reduced at depth.

Halls (1969) found that velocity and density generally varied in concert as a function of mineralogy and porosity. Amphibole-rich, diabasic, and fresh basalt samples have the highest densities and velocities (generally more than 2,900 kg/m³ and more than 6.5 km/s, respectively), whereas samples with high porosity have the lowest densities and velocities (generally less than 2,500 kg/m³ and less than 5.3 km/s, respectively).

Velocity and density data from Lippus (1988), plotted in figure 13, provide an interesting test of how velocity and density increase with depth within a thick section of Portage Lake Volcanics. Samples were obtained from three drillholes that collectively represent different stratigraphic levels within the 4.8-km thick basalt section near Eagle Harbor in the northern Keweenaw Peninsula (fig. 3). Both velocity and density increase down section, as seen from the progression of data points from the top, middle, and base of the section (fig. 13; drillholes EH2, C27, and D57, respectively). Lippus (1988) argued that the increase in both velocity and density at lower stratigraphic levels is explained by intensifying metamorphic grade, although the author did not investigate the mineralogy of the samples in detail.

The velocity-density data from Lippus (1988), using the velocity measured at 200 MPa, cluster well about a linear regression line (R² of 0.88) corresponding to equation 12.

In figure 14, the Halls line representing basalt of the Keweenawan Supergroup is compared to data for other igneous rocks related to the Midcontinent Rift System and to velocity-density data or model computations for similar rock types from other studies. The Christensen and Mooney (1995) data are plotted for multiple depths to demonstrate variability at deeper conditions than what can be achieved from the limited range of pressures used for measuring velocity from the Halls (1969) and Lippus (1988) studies. Their study also has the advantage that densities and velocities were measured at matching confining pressures.

Plots for diabase samples generally follow the Halls line for velocities greater than 6.5 km/s (fig. 14). The wide scatter of two of the Halls (1969) diabase samples compared to the tight clusters from the other studies suggests additional investigation is required. In contrast, velocity and density values for most of the gabbro samples plot significantly away from the Halls line, including the six samples of gabbro related to the Midcontinent Rift System from the Mellen Intrusive and Duluth Complexes (fig. 14). The similar densities for diabase and gabbro yet higher velocity for gabbro have been noted in previous rock physics studies (Birch, 1961; Christensen and Mooney, 1995). Christensen and Mooney (1995) attribute the difference to higher olivine content in gabbro. The similarity in densities for diabase and gabbro of the Midcontinent Rift System is corroborated by histograms of density measured on these two rock types from the Duluth Complex and related rocks that show comparable populations and median values both close to 2,920 kg/m³ (fig. 15). Although high-density outliers do exist (with averages greater than 3,000 kg/m³), they likely represent distinct gabbro intrusions (Chandler, 2002). Thus, the velocity-density plots indicate that a separate velocity-density equation is required to represent the bulk of gabbros in the Midcontinent Rift System. A line qualitatively honoring the gabbro data points that plot away from the Halls line is designated as the gabbro line and is shown as equation 13.

Only one sample of felsic rock related to the Midcontinent Rift System has been measured for both velocity and density (table 9). The single velocity-density point is plotted on figure 14 in combination with the range of density measurements for various felsic igneous rocks from table 7. The data for this sample fall in between those for granitic rocks compiled globally (Birch, 1960; Christensen and Mooney, 1995) and the density and velocity values for granite computed from mineral physics (Sowers and Boyd, 2019), which provide support for using the data from these other studies as a guide for what to expect for felsic igneous rocks related to the Midcontinent Rift System.

As explained at the beginning of this section, velocity-density data for anorthosite from the Lake Superior region are not plotted in figure 14 owing to lack of information. As seen from the plots from other geologic settings, anorthosite has unique properties compared to other igneous rock types, with relatively low densities and high velocities. Two samples of anorthosite xenoliths from the Duluth Complex yield measured densities of 2,640 and 2,680 kg/m³ (Chandler and Lively, 2021; Grauch and Heller, 2021). These densities are substantially lower than what is shown in figure 14, suggesting that further study of the Lake Superior anorthosites may be warranted.

Crystalline Crust

Because geophysical modeling of the Midcontinent Rift System must also deal with the crystalline country rock (fig. 2), developing velocity-density relations for rock types of the Archean and Paleoproterozoic crust is also important. The general rock types that immediately surround and possibly underlie large parts of Lake Superior, in order of regional extent, include Archean crystalline rocks of the Superior Province, Paleoproterozoic crystalline and metasedimentary rocks, and lower Mesoproterozoic sedimentary rocks (figs. 1, 3). The metasedimentary and sedimentary rocks are considered in the section “Older Metasedimentary and Sedimentary Rocks.”

Velocity-density relations for geophysical modeling of Archean and Proterozoic crystalline crust in the Lake Superior region can be examined by way of crustal models for North America, gravity and long-range refraction models for the region, and rock physics studies that measure relevant rock types. Values for these models are compiled for upper, middle, and lower crust and upper mantle in table 10 and plotted against various previously established equations for crystalline crust in figure 16, as well as points for relevant rock types from laboratory studies (Birch, 1961; Christensen and Mooney, 1995).

All the previously established velocity-density equations for crystalline crust pass through the rectangles representing the ranges of velocity and density values from the crustal models (fig. 16). However, the one that best matches both the upper and middle crust is equation 9 that was developed for the crust of average petrology from Christensen and Mooney (1995). This curve is also the only one that matches velocity-density points for granitic rocks, which most previous workers consider as representative of the upper Archean crust in the Lake Superior area (for example, Thomas and Teskey, 1994; Hinze and others, 1997; Drenth and others, 2021). Moreover, Perry and others (2002) considered the average petrology model as representative in predicting densities for a crustal model of the Canadian shield.

On the other hand, a simple crustal model that uses single densities for the upper, middle, and lower crust may be preferred where velocity information is unknown or poorly constrained. This approach has worked successfully in multiple gravity models for the Lake Superior region (references cited in table 10). From an inspection of table 10, figure 10, figure 16, and in consideration of the previous gravity models, a reasonable approach is to assign densities of 2,700 kg/m³ to the upper crust, 2,870 kg/m³ to the middle crust, 3,000 kg/m³ to the lower crust, and 3,310 kg/m³ to the upper mantle, with the depths between these layers determined subjectively or from seismic interpretation. Several passive seismic studies provide updated models of the crust and shallow mantle (Shen and others, 2013; Pollitz and Mooney, 2016; Zhang and others, 2016).

Older Metasedimentary and Sedimentary Rocks

Establishing velocity-density relations for metasedimentary rocks that predate the initiation of the Midcontinent Rift System is important to consider separately from the crystalline crust because (1) metasedimentary rock types are likely to have different physical properties than crystalline rocks, and (2) the foreland basins associated with the Penokean orogeny were regional-scale features that are preserved on either side of Lake Superior and presumably underlie the lake as well (fig. 1; Schulz and Cannon, 2007).

Investigating velocity-density relations for these older rocks is difficult because of the lack of velocity data and the heterogeneity of rock types. Figure 17 compares histograms of density data for rocks from the Lake Superior region (fig. 17A) to velocity-density data for similar rock types from elsewhere and the Nafe-Drake curve of Brocher (eq. 2; fig. 17B). A histogram of density measurements for sandstones and siltstones of the lower Mesoproterozoic Sibley Group are included for completeness (fig. 17A), although the volumetric significance of these rocks is unclear. The velocity-density plots for three rock types relevant to the Lake Superior region (metagraywacke, slate, and quartzite) on figure 17B are from the global compilation of Christensen and Mooney (1995). One inferred velocity-density point representing rocks of the Animikie Group comes from Halls and West (1971; fig. 17B). They estimated a typical density of 2,700 kg/m³ for Animikie Group metasedimentary rocks and reported velocities for two samples of slate and graywacke of 5.8 and 6.0 km/s at 30 MPa, respectively, which is represented by a value of 5.9 km/s on figure 17B.

The histograms of measured densities for Paleoproterozoic metasedimentary rocks show wide variation (fig. 17A). This heterogeneity is underscored by the different spread of values displayed by the histograms for Minnesota compared to the one for Michigan and Wisconsin and by the multiple peaks in the histogram plots. The multiple peaks can be partially explained by dissimilar densities for the different metasedimentary lithologies, as demonstrated by the plot of medians or means of subpopulations separated out by rock type (fig. 17C). A predominance of quartzite or metagraywacke can explain density peaks near 2,650 kg/m³ or 2,750 kg/m³, respectively, but the spread of densities for other rock types is large. Moreover, the densities of the Lake Superior rocks differ from those of comparable rock types in the global compilation of Christensen and Mooney (1995; fig. 17B), except for quartzite. The differences in density by rock type raise questions about how representative the rock types used by Christensen and Mooney (1995) are and confirm how the classification of rock types commonly fails to capture the heterogeneity of mineralogy, texture, and the degree of metamorphism that is important for determining physical properties

Faced with heterogeneous density data to compare to a paucity of velocity data, the Nafe-Drake curve is a reasonable choice for representing the Paleoproterozoic metasedimentary rocks. The marine setting of the rocks sampled for the Nafe-Drake curve matches the marine sedimentary environments common during Penokean cratonic margin basin formation (Schulz and Cannon, 2007). Use of the Nafe-Drake curve is buoyed by its proximity to velocity-density points for the inferred Animikie Group and metagraywacke and quartzite from the global compilation (fig. 17B).

Choosing one of the established velocity-density equations for use with the data-poor Sibley Group by comparing geologic settings breaks down because the Sibley Group includes a wide variety of depositional environments, including fluvial, lacustrine, sabkha, delta, and eolian (Rogala and others, 2007). In this case, the Nafe-Drake curve is a practical choice because it matches the relation used for the Paleoproterozoic rocks and thus provides a continuum between all the sedimentary and metasedimentary rocks older than the Midcontinent Rift System.

The wide range of densities for the older metasedimentary and sedimentary rocks of the Midcontinent Rift System makes identifying objective, quantifiable associations between density and velocity challenging. Despite this uncertainty, using the Nafe-Drake curve to represent these rocks helps visualize the range of physical properties that can be expected in relation to other rock types.

Summary of Velocity-Density Relations

Empirical velocity-density relations for the Lake Superior region can be simplified to specific equations and velocity-density criteria based on geologic unit or rock suite. These measures are for ease of use and should be used as guidelines for geophysical modeling. The equations or criteria for each of the groups are summarized in table 11 along with the expected velocity and density ranges. Equations could not be developed for felsic igneous rocks, which may be present in significant volumes at depth under the lake, or for surficial deposits, which cover wide areas. Instead, only target densities and velocities for these materials and for those of the upper mantle are listed. Moreover, two alternate guidelines are presented for crystalline crust: one that relies heavily on accurate velocity information and assumes a certain progression of rock types with depth, and another that views the crust as three discrete layers.

The graphical view of the velocity-density relations (fig. 18) underscores the dissimilarities of relations between the rock suites and geologic groups as well as the overlaps in velocity and density ranges. This view of the relations and the expected ranges for each group from table 11 can help guide geophysical modeling and geologic interpretation of seismic and gravity data. For example, if a velocity of 5.5 km/s is observed, the most likely candidates based on velocity range (table 11) and their corresponding densities would be basalt if 2,730 kg/m³, Oronto Group if 2,670 kg/m³, or older metasedimentary or sedimentary rocks if 2,620 kg/m³. The possibilities could presumably be tested by gravity modeling. The ranges of velocity and density also can be used to determine the most likely candidates for geologic interpretation. For example, if velocities are greater than 7.0 km/s, the most likely candidate is gabbro. If the velocities are less than 3.3 km/s or densities less than 2,300 kg/m³, the most likely candidate is one of the units overlying the Oronto Group. Velocities greater than 5.9 km/s are not likely related to Keweenawan sedimentary rocks.

Uncertainties in the Inferred Velocity-Density Relations

The guides presented herein are intended to be useful as representations of the physical properties of the different rock suites, but the representations are limited for several reasons. First, there is not a one-to-one correspondence between geologic units and density or velocity. For example, sedimentary units commonly vary in lithology, grain size, and degree of induration or compaction over a wide geographic area, even if they are of similar age. An example from the Lake Superior area comes from the St. Amour well near Munising, Michigan (fig. 3), wherein the stratigraphy below the Freda Sandstone could not be correlated with Oronto Group known elsewhere (Ojakangas and Dickas, 2002). Another example are the anomalous, well-indurated samples of Jacobsville Sandstone that were collected near Sault Ste. Marie that are discussed in the section on velocity-density relations for sedimentary rocks of the Keweenawan Supergroup. This problem is a fundamental limitation to geophysical modeling of geologic units and is not unique to this study.

Second, the lack of velocity information limits the ability to confidently determine velocity-density relations, especially for metasedimentary and sedimentary rocks. For example, the limited velocity data for the units overlying the Oronto Group leads to tentative conclusions about velocity-density relations and prevents investigation of differences between the formations. As another example, the expected velocity and density ranges for metasedimentary and sedimentary rocks older than the Midcontinent Rift System span almost the entire ranges of all rock types combined. These substantial ranges are in large part due to the wide variety of rock types and degree of metamorphism within the rock suite. Moreover, in the velocity range of 5.8–6.9 km/s, the velocity-density relation inferred to represent the older rocks is similar to that of the model of crystalline crust estimated from average petrology. There is enough uncertainty in both models that crust composed of crystalline versus metasedimentary rock cannot be confidently distinguished using these guidelines.

Finally, how well velocity and density measurements from rock samples represent the bulk physical properties of large volumes measured by geophysical methods is unclear. Samples collected only from outcrops may not provide an adequate representation of variability with depth. Moreover, velocity varies within one basalt flow (Christensen and Wilkens, 1982; Planke and others, 2000), which begs the question of whether (1) enough samples are present to represent the breadth of this variation, (2) if this variability is maintained at depth, and (3) how the variability translates into a bulk average at depth when geophysical methods cannot resolve the variations. Velocity also can vary due to microfractures or seismic anisotropy without an accompanying difference in density (Mavko and others, 2020). Variations in average velocity or density with depth also may be due to mixtures of rock types rather than variability within one type or may be due to rock types that are not anticipated from our current understandings.

Additional velocity measurements on rock samples and perhaps short-range seismic refraction experiments where geology is well known could help reduce the shortage of velocity data. Rock types having the least data are felsic igneous rocks and anorthosite related to the Midcontinent Rift System as well as older metasedimentary and sedimentary rocks. Although sufficient density data are available for the Paleoproterozoic rocks, additional velocity measurements need to be tied directly to measured density to make sense of velocity-density relations. An investigation into the physical properties of the Jacobsville Sandstone samples found near Sault Ste. Marie could help understand why they are anomalous, but probably would have to be accompanied by detailed study of lithologic variations of the formation. Finally, the unique properties and presence of anorthosite, both within layers of the Duluth Complex and as separate bodies brought to the surface as xenoliths, may warrant further investigation.

Implications for Seismic Interpretation

The table summarizing typical velocity ranges by rock type (table 11) can be used to address the differences in interpretation of basalt versus gabbro between earlier seismic refraction and reflection studies (presented in the “Introduction” section). Basalt, diabase, and gabbro have expected velocity ranges of 5.1–6.8, 5.9–7.0, and 6.7–7.4 km/s, respectively, which overlap. Using these ranges, velocities of 6.6–6.9 km/s consistently found for depths from approximately 10 km down to the Moho across Lake Superior from earlier seismic refraction studies (Berry and West, 1966; Smith and others, 1966; O’Brien, 1968; Ocola and Meyer, 1973; Luetgert and Meyer, 1982) could be interpreted as average velocities over large volumes of basalt, diabase, and gabbro. Basalt, possibly mixed with diabase, is more likely where these velocities correspond to strong, semihorizontal layering in the GLIMPCE seismic-reflection sections (Cannon and others, 1989; Green and others, 1989). However, velocities of 6.9–7.1 km/s associated with similar layered intervals at depths of 10–20 km for GLIMPCE Line A (fig. 3) are too high to represent basalt, as was presented by Tréhu and others (1991) and Shay and Tréhu (1993). The velocities are more consistent with an interpretation of diabase and gabbro. Metamorphism of basalt at these depths can only explain velocities as high as 6.9 km/s, assuming basalt in this region metamorphoses through amphibole replacement and velocity increases with greater amphibole content (Halls, 1969). If so, the maximum velocity for metamorphosed basalt would be represented by the velocity of amphibolite, which ranges from 6.7 to 6.9 km/s at depths of 10–20 km (Christensen and Mooney, 1995). Alternatively, basalt flows that are thick enough to develop a pegmatitic internal zone may resemble diabase, as Huber (1973, 1975) observed for several basalt flows on Isle Royale. If so, perhaps the coarse-grained layers would increase basalt velocities to 6.9–7.0 km/s without complete replacement by amphibole.

Thus, geologic interpretations of layered intervals with velocities greater than 6.9 km/s in seismic-reflection profiles may need to be reassessed. If the intervals represent mafic intrusions rather than basalt, then their layered character becomes significant and lacks current explanation. If they represent amphibolite in places, the metamorphic history is poorly understood, because amphibolite is not observed in outcrop.

Geologic interpretation of the seismic refraction study of Mooney and others (1970a, 1970b) could be reassessed as well. Overlaps in velocities between Freda Sandstone and Bayfield Group and between Copper Harbor Conglomerate and basalt led to misinterpretation of geologic units for at least one of their refraction models in the Lake Superior region (fig. 9). These ambiguities may have affected additional interpretations in their study area farther to the southwest.