Landscape evolution under the southern Laurentide Ice Sheet

Description

We compiled and refined bedrock surface data for northwestern Missouri, Kansas, Montana, and North Dakota. We digitized and edge-matched existing bedrock surface maps for northeastern Montana, southeastern Saskatchewan, and northwestern North Dakota (97). We georeferenced and digitized previously published bedrock contours for the rest of North Dakota (98). Similarly, existing bedrock surface maps of northeastern Kansas (99) and northwestern Missouri(100) were digitized, georeferenced, and attributed with elevation values. In northeastern Missouri, we interpolated a bedrock surface from water-well logs with bedrock contact designations (101).
For Manitoba, the existing bedrock surface map (102) contains sparse data along the Pembina Escarpment (G. Keller, personal communication, 17 December, 2015), therefore, we digitized and incorporated previously mapped bedrock contours (103). Additionally, eastern Manitoba is mapped as bedrock outcrop (104) and we used ground-surface elevations (Shuttle Radar Topography Mission data) for the bedrock surface in this region.
Bedrock surface resolution We interpolated the final bedrock surface interpolation to 250-m raster cell resolution. We chose this resolution because a 1:500,000 map scale corresponds this raster cell resolution (105) and the source data were predominantly mapped at the 1:500,000 scale or better. We compared the final DEM to a 250-m resample ground surface DEM derived from 1-arc second Shuttle Radar Topography Mission (SRTM) data (106), and cells with bedrock elevations above the ground surface were set to equal the SRTM elevation.
Technical validation of the bedrock surface We conducted a spatial statistical analysis to ensure that the final ANUDEM interpolation results, derived primarily from digital contour input data, compare well with point data throughout the study area. Geologic boring records with bedrock elevations were readily available from 8 states and provinces (Indiana, northern Kentucky, Manitoba, Michigan's Lower Peninsula, Minnesota, northeastern Missouri, Ohio, and southern Ontario) to conduct a statistical analysis of the final 250-m bedrock-surface topography DEM. The assessment points were paired with elevation values from the final interpolated DEM at the same locations. An ordinary least squares regression analysis indicated regression coefficients (R) between 0.95 and 0.99 for each of the analysis regions (Fig. S2). We also calculated root mean square error (RMSE) between ground-truth point data and the associated raster elevation cell in each of the state/provinces (Fig. S2) and obtained a mean RMSE of 12.5 m and a range from 8 to 17 m of RMSE values for individual states/provinces. This is close to the range of error in the SRTM specification (16 m), and error between the gridded data set and point measurements can also be produced by sub-grid topographic variability in bedrock surface elevation.
Residual glacial isostatic adjustment (GIA) The GIA correction produced an additional 3 to 96 meters of uplift, which we applied to the modern bedrock surface. The greatest amount of residual uplift lies in the northeastern portion of the study region, which hosted thicker ice that also retreated most recently. Because of the radial pattern of GIA around the former Laurentide Ice Sheet and increasing residual GIA magnitudes farther within its footprint, the Paleo-Missouri River and Paleo-Saskatchewan River long profiles (see below) were most strongly impacted by this residual GIA.
Erosional isostatic adjustment was not simulated because (1) insufficient data exist beyond the study area to constrain the spatial distribution of erosional unloading and (2) magnitudes of expected erosional isostatic deformation (~100 m) are similar to uncertainties among different tuned solid-Earth models for dynamic topography in Canada (56). Although we constrain mass redistribution within the study area, erosional isostatic response is impacted significantly by mass redistribution within one flexural wavelength (~500-1000 km) of the study region margins. Furthermore, local isostatic adjustment links solid-Earth rheology with both surface and subsurface loads, and as such, significant research on the deformation of mid-continental North America must be completed in order to solve this joint problem. We acknowledge that the isostatic effects of Quaternary erosion, as well as long-wavelength land-surface deformation due to dynamic topography, represent an important problem for future research. Therefore, we created an extensive data package from this work to support future research into surface deformation. Our geologically constrained river courses and paleotopography, alongside changes in surface loading from our reconstructed erosion and deposition depths, should help to limit and inform plausible geophysical solutions.
Strath terrace elevations and pre-glacial (end-Pliocene) surface reconstruction Multiple strath terrace levels are present in the modern Missouri River basin (107) and in the buried bedrock topography of the Southern Paleozoic Basins and Arches province ( Fig. 1) as indicated by "deep stage" bedrock incision (i.e., incision to the deepest point along a buriedvalley cross-sectional profile) events (65,108). The occurrence of multiple terrace levels that record former floodplain elevations is problematic when establishing consistent paleo channel elevation indicators along a valley. We overcome this issue by generally choosing prominent valley-bounding strath terraces (Fig. S6) and ensuring that the elevations along a channel coincide with a consistent gradient as shown by the sequence of profiles interpreted in Figures  S5 and S6. Furthermore, we integrate data from previous studies on Late Neogene geology when available such as the Flaxville and Souris River stratigraphic units(109) present along the section of the paleo Missouri River shown in Fig. S5. The robustness of this method is further supported by pre-glacial river longitudinal profiles (Fig. 5).

Reconstructing pre-glacial channel and interfluve patterns Paleo Saskatchewan and Missouri Rivers
The west-to-east orientation of valleys that routed eastward drainage across the Pliocene Canadian Plains (110,111) is reconstructed using low-gradient channels delineated from strath terraces along major valleys in the Rocky Mountain foothills that trend toward eastern modernday channels feeding Hudson Bay (Figs. 2, and S5). The northernmost valley in the Paleo Saskatchewan drainage network follows the trend of the modern Beaver River course and has headwaters consistent with those of the modern Peace River (Fig. S10). Others have documented the historic path of the Upper Peace River toward Hudson Bay (112). To the south, the current Upper Missouri River and the ancestral South Saskatchewan River(113) also flowed toward Hudson Bay. These reconstructions are consistent with paleoflow indicators in the Cenozoic gravels of Western Canada and the northern United States (109). Pre-glacial elevations in northern North Dakota were generally consistent with the modern Missouri Escarpment (36). Similarly, geomorphic analyses of terraces in the U.S. Black Hills(107) as well as Late Miocene fluvial benches of the northern U.S. and western Canada(109) provide guidance for inferring preglacial floodplain elevations and river longitudinal profiles in the region.

Paleo St. Lawrence River
Early work in the Great Lakes led to a proposed pre-glacial drainage that followed the axes of the lakes except where the pre-glacial "Laurentian River" crosscut southern Ontario between the Georgian Bay of Lake Huron and Toronto, along the northwestern shore of Lake Ontario (114). This interpretation is largely supported by later studies investigating the lithologic (differential erodibility of carbonate and shale rock units) and structural (bedrock jointing and rock unit dip) controls on Great Lakes Region bedrock valley orientation (114,115) and regional bedrocksurface morphology (60,83). More recently work by Carson et al. (1) provided evidence for the pre-glacial "Wyalusing River" with headwaters near the modern confluence of the Mississippi and Wisconsin Rivers in the mid-continent and which flowed northeastward towards the Lake Michigan basin. The Sawtooth Mountains to the north and west of Lake Superior are interpreted as the northern pre-glacial divide of the Paleo St. Lawrence River drainage as are bedrock highs north of Lakes Michigan and Huron. Pre-glacial tributaries in the southern portion of the St. Lawrence drainage, were sourced from topographic highs along a cuesta formed by resistant Paleozoic strata between Lake Michigan and Lake Erie (Fig. 1). Here, we used high bedrock surface elevations to interpolate a divide between this drainage and the ancestral Wabash River to the south.

Paleo Wabash (Teays) River
Adjacent to the Great Lakes, the current southwesterly path of the Ohio River toward the Mississippi River is widely known to be a result of Quaternary glacial diversion (59,65,(116)(117)(118)(119)(120), but its pre-glacial course is debated. Some researchers identify the Teays-Mahomet buried valley system, spanning from West Virginia in the east to Illinois in the west, as the pre-glacial drainage (120). Others suggest that Pliocene rivers east of Indiana flowed to the north into the Lake Erie basin (59), formerly occupied by the Erigan River (39,121). However, most researchers acknowledge the northwest trending entrenched Teays Valley in Ohio as the pre-glacial main channel that drained much of Ohio and northern Kentucky(108).
Preserved sediments and landforms provide evidence that the Old Kentucky and Licking Rivers in northern Kentucky (Fig. S10) once were the headwaters of a northwestwardly flowing system that was abandoned during the Quaternary diversion/capture of the modern Ohio River (117,119). Valley morphology and abandoned terraces capped with gravel deposits indicate that the valley containing the present-day South Fork of the Licking River once held the pre-glacial trunk stream (119) (122) indicate that early-to mid-Quaternary deposits exist along the current drainage divides of high-order watersheds suggesting that pre-glacial headwater channel elevations were closer to current ridgeline elevations in the region. This is consistent with slackwater sediments on terraces associated with the northward flowing ancient Pittsburgh River in western Pennsylvania and northeastern Ohio (123). These high-level terraces exist at ~330 m and were dated as early Quaternary (prior to 0.774 Myr BP; millions of years before present) based on their reversed magnetic polarity (123,124). Headwater channel elevations for the shortened ancestral Ohio River are estimated based on strath terraces near the pre-glacial drainage divide for the ancestral Ohio River basin near Madison, Indiana, where an undissected upland exists (120). The elevations are consistent with cosmogenic dating work in the Ohio River basin that suggest Quaternary incision rates of ~30 m/Myr, mostly during three primary pulses of river downcutting (125,126).
In the modern Wabash and Ohio River basins (Figs. S4 and S10), bedrock strath terrace elevations visible in the bedrock topography (Fig. 1) do not clearly support the previously suggested east-west course of the pre-glacial Teays that drained into the modern Mississippi channel. The interpretation presented herein, based on 180-215-m strath elevations in the central Paleo Wabash valley and the longstanding tectonic activity in the Wabash Valley Seismic Zone (127), is that the upper Teays was a tributary to the Paleo Wabash River and not to the modern Mississippi channel.

Paleo Platte River
In the western portion of the current Mississippi River basin, the present-day Lower Missouri River developed as an ice-marginal feature (128) and the paleo Platte River was the major Pliocene drainage to carry water from Nebraska, northern Kansas, and northern Missouri(5). Calculated post-depositional tilt for the Miocene-Pliocene Ogallala Group based on the preserved channel-scour geometry and grain size (62) indicates that the eastern portion of the Paleo Platte, which lies in our study area, did not experience uplift. Therefore, SRTM elevations of Cenozoic sediments exposed and mapped at the ground surface in central Nebraska are used to represent the trend of pre-glacial topography in the headwaters of the Paleo Platte. The bedrock surface data presented (Fig. 1) indicates that the Paleo Platte River drained a large portion of the mid-continent west of the Indiana and Illinois state border (Fig. S3).
Geologic mapping indicates that the precursor to the central Mississippi River course is the ancestral Illinois River(129), a southward flowing tributary of the Paleo Platte River (Fig. 2). North of this, bedrock elevation data indicate that the northern portion of the modern upper Mississippi River system drained towards the Paleo St. Lawrence River (1). Near the headwaters of the ancestral Illinois River, regionally correlated geomorphic surfaces consistent with the pre-0.774 ka Bridgeport strath terrace (58,124,130) along the pre-glacial Wyalusing River(1) provide minimum elevation constraints on the pre-glacial Gulf of Mexico / St. Lawrence River drainage divide (Fig. S9). The reconstructed position of this pre-glacial continental drainage divide as well as its northwards relocation during the Quaternary are consistent with Pliocene drainage-basin reconstructions based on the sedimentary record in the Mississippi embayment, which that indicate a major Pliocene-Quaternary reorganization of the Missouri River(5).      Colorbar is glacial sediment deposition in meters.

Fig. S8. Bedrock erosion with contour lines showing distance (km) from the Keewatin and Quebec-Labrador ice domes and connecting ice saddle (yellow area to the north of study area) derived from sediment dispersal patterns(21).
Colorbar is glacial sediment erosion in meters. The scale of the colorbar is limited to 500 m but there are erosion extremes that reach 1000 km in the Lake Superior basin (location shown in Fig.  1).

Fig. S9 Mississippi River valley bedrock topography showing diversion from end Pliocene drainage.
The dashed box on the inset map (same inset map from Fig. 1 of the main text) shows the location of this region relative to the full analysis extent. The Wyalusing River was dammed during an early Quaternary LIS advance(6) and flow was diverted through Military Ridge(1) to create the incised modern Mississippi River channel.