In the continental interior, isobases of glacio-isostatic rebound are usually defined by the elevations of a differentially uplifted shoreline of a former lake (Goldthwait, 1907, 1910; Leverett and Taylor, 1915; Hough, 1958; Walcott, 1972; Lewis and Anderson, 1989; Schaetzl et al., 2002). Sets of isobases selected for this study for basins within the Great Lakes watershed are illustrated in Figure 2A and 2B in which the lowest and highest elevations, trends of isobases, the name and uncalibrated radiocarbon age are shown for each lake. Although the Algonquin and Nipissing highstands were confluent in three of the basins about 10.6 ka BP (12.6 cal ka BP) and 5.0 ka BP (5.7 cal ka BP), respectively, most sets of isobases are confined to a single basin. These former, once-level lake surfaces are all warped upward in a north to northeasterly direction (Fig. 2A-B) as a result of differential glacio-isostatic recovery of Earth’s crust following deglaciation of the Laurentide Ice Sheet. The sources of isobase and related deglacial information for basins within and adjacent to the Great Lakes are listed in Table I.

Whereas isobases portray the present configuration of a rebounding surface, it is often desirable to describe the uplift of a given location through time, or to construct surfaces at intermediate times, especially when reference shorelines are widely spaced in time. Following Andrews (1970) and others, the exponential function is adopted to describe relative uplift (U_t) with time where time is expressed as age t (cal years BP) for a landscape previously loaded by an ice sheet. From Peltier (1994, 1998) we use:

where A and τ (tau) are parameters of the equation. A is a site-specific amplitude factor, and is evaluated as:

for a known relative uplift, age and relaxation time. Tau (τ) is the relaxation time or period in years for which decelerating uplift is reduced by 1/exp (1/2.7183 or 36.8%) in successive periods. It is evaluated by solving (Equation 1) at sites where U_t is known for at least two sets of isobases of different ages as shown below. As a first-order approximation, τ is assumed to be time invariant and similar in value throughout the Great Lakes region. Figure 3 shows the reasonably good fit of the exponential uplift curve to Great Lakes relative rebound data for the Nipissing, Algonquin and Iroquois shorelines between sites C and D for trial values of relaxation time between 3000 and 5000 years (see Fig. 2 for location and isobase values).

The 12 GIS paleogeographic reconstructions based on the empirical model of isostatic adjustment were used to measure land and water areas, and lake volumes for studies of the paleohydrology and meltwater flow through the Great Lakes system (Moore et al., 2000). Variations in the basin attributes in the 11.4 to 5.0 ka BP (13.3 to 5.7 cal ka BP) period were substantial (Fig. 8). Maximum variation in individual basins under overflow conditions ranges from +72% to ‑95% for lake area, and from +200% to ‑97% for lake volume, compared to the present Great Lakes.

Because of their connection via the Indian River lowland and Straits of Mackinac (Fig. 1), high water levels were always at common elevations in the Huron and Michigan basins, as at present, following retreat of ice from their northern regions (Eschman and Karrow, 1985; Hansel et al., 1985). Consequently, the history of lake-level variation above the Michigan basin sill defined in the Huron basin is applicable to the Michigan basin as shown in Figure 11. Only evidence of lowstands and other unique indicators of Michigan basin lake elevations are shown and discussed here.

Submerged in situ tree stumps in southern Lake Michigan (the Olson Forest of Chrzastowski et al., 1991) have been traditionally interpreted as having been drowned in the Nipissing transgression (Chrzastowski and Thompson, 1992, 1994; Colman et al., 1994b). However, in this analysis, the Olson tree stumps (sites 74a, 74b at 8.2-8.4 ka BP, 9.2-9.4 cal ka BP, in Figs. 1 and 11) appear tens of metres above the North Bay outlet, and could not be affected by the Nipissing transgression. A rise of water level to the final highstand of Lake Mattawa follows the Olson Forest after 8.2 ka BP (Fig. 11). This flooding event seems to have been the cause of forest drowning.

As in the Huron and Georgian Bay basins, an extreme decline in lake level following the last Mattawa highstand is suggested in the Michigan basin by the occurrence of shallow-water sediments and molluscan fauna in the deepwater environment of central Lake Michigan (51a, 51b, 51c) (Lewis et al., in press). A final recovery of lake level to an overflowing condition at the North Bay outlet probably occurred about 8.1-7.8 cal ka BP, allowing for tree growth in northern Michigan basin (54). From this time forward, relative lake level rose under control of the North Bay outlet as the Nipissing transgression.

During the period 8.95-8.3 cal ka BP (8.05-7.4 ka BP) lake levels (late Stanley in Huron basin, late Hough in Georgian Bay basin, and late Chippewa in Michigan basin) were below the sill of the North Bay outlet, the lowest possible overflow outlet at the time (Figs. 10B‑11). Thus the late Stanley, late Chippewa, and late Hough phases were hydrologically closed lakes at their lowest level, and as such may have resulted from the impact of a severe dry climate in which evaporative water losses exceeded water inflows by precipitation and runoff. The inference of dry climate is supported by the presence of thecamoebians that indicate a more saline lake environment, and by climate transfer function analysis of pollen assemblages which suggests less precipitation and warmer temperatures at 7.7 ka BP than today in the Georgian Bay area (Blasco, 2001). Other lowstands around 9.3 and 9.8 ka BP (10.5 and 11.2 cal ka BP) may have also been closed briefly.

Possible shorelines for the dry climate-induced closed lowstands at 7.8 ka BP (8.6 cal ka BP) are illustrated in a paleogeographic reconstruction of the Huron, Erie and Ontario basins on Figure 12. The northern Huron lowstand shore is tied to the Stanley unconformity beneath northwestern Lake Huron (Hough 1962; Lewis et al., in press). The Georgian Bay lowstand shore is tied to the low-level Flowerpot beach in the entrance to Georgian Bay (Blasco, 2001). Water bodies in these basins and those beneath southern Lake Huron were isolated closed lowstands, well offshore from the present lake boundaries. The Erie lowstand shore is tied to a submerged beach in eastern Lake Erie (Fig. 9E; Coakley and Lewis, 1985; Lewis et al., 2004) at a lake level that would have extended into the central Lake Erie area before that sub-basin became mostly infilled with sediment (Sly and Lewis, 1972). The Ontario basin is assumed to have been impacted by the same dry climate that affected the other basins. Under these conditions, its lowstand shoreline was inferred at a level that made its basin area-to-lake area ratio equal to that of the Huron basin (Bengsston and Malm, 1997).

With knowledge or estimates of the elastic and viscous properties of Earth’s lithosphere and mantle, respectively, geophysical models compute the crustal isostatic response following a known or inferred history of ice sheet loading and pattern of deglaciation. Clark et al. (1994) used this approach to evaluate isostatic movements during deglaciation of the Great Lakes basin. They successfully illustrated the increasing amplitude of postglacial isostatic rebound towards the north and northeast in the direction of past thicker ice and ice retreat in accordance with the evidence of deformed paleo-lake shorelines. They applied an innovative approach to constrain and calibrate their model by using geological knowledge of large-lake drainage transfers from northern to southern outlets as basins were differentially tilted during the Algonquin and Nipissing phases. Lake-level history was derived by tracking the upward movement of overflow sills through time. Lake-level indicators were not tracked separately from the sills with the result that possible episodes of hydrologic closure could not be detected.

Gravitationally self-consistent models of glacio-isostatic adjustment on a global scale have been progressively developed and improved over the past few decades (Peltier, 1998). These models are constrained and calibrated to relative sea-level histories at ocean-continent boundaries. Model estimates of isostatic adjustment are produced for continental interiors, and one of these, the ICE‑3G model (Tushingham and Peltier, 1991), has been compared favourably with short-term evidence of tilting of the Great Lakes basins derived from trend analysis of lake-level gauge records (Tushingham, 1992). However, over the longer term, the geophysical model results are not in complete agreement with past events in the watershed. The relative performances of two geophysical models and the empirical model were assessed by comparing the elevations of the northern and southern outlets from the upper Great Lakes basins during the Algonquin and Nipissing drainage transfers which are known to occur at about 10.5 ka BP and 5 ka BP, respectively, from independent geological evidence (Karrow et al., 1975; Karrow, 1980; Monaghan et al., 1986). Outflows were transferred from the Kirkfield outlet to the Port Huron and possibly Chicago outlets during the Algonquin transfer (Fig. 7A), and from the North Bay to the same southern outlets during the Nipissing transfer. Estimates of the elevations of the overflow sills and timing for these drainage transfers when northern outlets rose isostatically above southern outlets were provided for the ICE‑4G glacio-isostatic model (Peltier, 1995, 1996) courtesy of W.R. Peltier. Differences between Port Huron and the other outlets which should have been near zero at the times of the drainage transfers are shown in Figure 13. The empirical model performs best with northern and southern sills being within two metres of each other during transfers. The Clark et al. (1994) model is similar, although the best agreement for the Algonquin transfer was obtained at 10 ka BP, rather than the expected 10.5 ka BP, and the Nipissing-age transfer at Chicago differed by 3 m. Variances in sill elevations during transfers for the ICE‑4G model are larger, and ranged from 3.5 to 10 m.

Near and beyond the maximum margins of the ice sheets, the geophysical models predict crustal uplift as a forebulge of modest relief compared with the amplitude of depression beneath the centre of the ice load. This has been demonstrated in ice-marginal Atlantic and Arctic coastal regions (Barnhardt et al., 1995; Dyke, 1998). The forebulge effect has not been recognized from empirical evidence in the continental region adjoining or south of the Great Lakes basin. However, ice marginal areas in this region could have undergone uplift and subsidence associated with the growth, migration and decay of a glacio-isostatic forebulge (Colman et al., 1994a), and some evidence for subsidence exists, for example, anomalies in modern tilting of the Erie and southern Michigan basins (Mainville and Craymer, 2005).

Benefits of geophysical models are their ability to compute both vertical and horizontal earth movements associated with the isostatic process, and to relate these estimates to an absolute datum such as present sea level. Confidence in predictions of past crustal movements in the Great Lakes region by geophysical models, especially deformation over several millennia, will be increased markedly when these models are calibrated and constrained by the empirical observations of glacio-isostatic adjustments within the same region.

An extreme lowering of lake level in the Huron and Michigan basins after the Algonquin phases was predicted by Stanley (1936) when the ice sheet margin receded from the isostatically-depressed lowland drainage route past North Bay, Ontario, to the Mattawa and Ottawa river valleys. Sediment unconformities, discovered by Hough (1955, 1962), confirmed lake-level lowstands in each of the Michigan and Huron basins; these lowstands were named Chippewa and Stanley, respectively. Hough (1962) inferred a delay in isostatic adjustment so that, in his model, the North Bay sill remained at its depressed Algonquin level, and Lake Stanley could be interpreted as an open overflowing water body (Lewis et al., in press).

Although a regional advance (Marquette readvance) of the Laurentide Ice Sheet about 500 km wide reached the southern coast of the Superior basin at 10 ka BP (11.5 cal ka BP) (Lowell et al., 1999), and ice retreat was slowed during the Younger Dryas in eastern Ontario and western Québec (Simard et al., 2003), there is no evidence of ice advance or of delayed isostatic recovery in the Georgian Bay-North Bay outlet region, about 300-500 km east of the Superior basin (Dyke et al., 2003). As a result, it is highly unlikely that rebound was delayed significantly and, at the time of late Lake Stanley (Fig. 10), the recovered outlet was above water level, as was also interpreted by Lewis et al. (in press). Closed lake conditions for late Lake Stanley are also consistent with the preliminary paleoecological findings based on thecamoebian analysis of somewhat higher salinity in late Lake Hough, the equivalent closed water body in the Georgian Bay basin, as reported by Blasco (2001).

Major shifts in Great Lakes water levels have long been understood in terms of overflowing lakes. Lake elevations change as a result of shifting ice dams during glacial retreat or advance, outlet erosion, or by differential glacio-isostatic adjustment of the outlets (Eschman and Karrow, 1985; Hansel et al., 1985; Larsen, 1987; Barnett, 1992; Larson and Schaetzl, 2001). These mechanisms apply to the Algonquin highstand and the transgression to the Nipissing Great Lakes. The intervening Mattawa high phases (Lewis and Anderson, 1989) are also considered to be overflowing lakes, but with variable levels controlled by resistance at hydraulic constrictions downstream of North Bay to variable high-discharge flows from upstream Lake Agassiz or subglacial drainage.

The discovery of low lake levels below the lowest possible overflow outlet in their basins implies a controlling factor of excess evaporation (water loss) over precipitation and runoff (water supply) during a period of dry climate. These lowstands, the late Chippewa, late Stanley, and late Hough lakes in the Michigan, northern Huron, and Georgian Bay basins, respectively, are thought to reflect a phase of climatically-driven hydrologic closure that has not been recognized previously. The earlier pre‑9.5 ka BP (pre‑10.9 cal ka BP) lowstands, middle and early lakes Chippewa, Stanley and Hough, are inferred to be close in level to, or slightly below, their basin outlets, and thus were possibly also subject to enhanced evaporative losses of water. During these early, low lake phases, sustaining upstream inflows were not available, as Lake Agassiz discharge was diverted away from the Great Lakes by advances of ice in the Superior and Nipigon basins (Thorleifson and Kristjansson, 1993; Lewis et al., 1994). Dry air from the glacial atmospheric circulation over the nearby Laurentide Ice Sheet (Bryson and Wendlund, 1967; David, 1988; Anderson and Lewis, 2002; Wolfe et al., 2004) probably exerted continuous evaporative stress on nearby water surfaces, causing drawdown of these lakes during phases of reduced inflow.

By the time of the latest and longest-duration lowstands (late Chippewa, Stanley and Hough) the effects of glacial atmospheric circulation would have diminished greatly owing to the small size of the remaining Laurentide Ice Sheet, and to its more distant location relative to the Great Lakes basin. At this time (8 ka BP, 8.9 cal ka BP), meltwater drainage from the merged glacial lakes Agassiz-Ojibway began bypassing the upper Great Lakes basins directly into the Ottawa Valley (Veillette, 1994; Teller et al., 2002; Teller and Leverington, 2004). Thus, the Great Lakes watershed suddenly lost a significant source of water supply and was susceptible to the early Holocene dry climate (Edwards et al., 1996).

Atmospheric water supply today can be thought of as a function of the relative time spent over the Great Lakes basin by three major air masses over North America. These masses are the Arctic air from the north (dry and cold), the Pacific air from the west (dry and warm), and the Maritime Tropical air bringing moist, warm air north from the Gulf of Mexico (Bryson and Hare, 1974; Bradbury and Dean, 1993). Once glacial lake drainage began bypassing the Great Lakes, and as the Laurentide Ice Sheet downwasted over Hudson Bay, southward incursions of dry Arctic air, previously blocked by the high Laurentide Ice Sheet, were likely becoming increasingly frequent (Yu and Wright, 2001), and this dry air may have initiated or intensified draw down of lake levels to the late Chippewa, Stanley, and Hough lowstands. The closed lowstands may have been maintained by enhanced evaporation into increasingly strong flows of dry, warm Pacific air from the west. These flows are indicated by abundant evidence of vegetation shifts to drought-resistant plant taxa west of the Michigan basin (Baker et al., 1992; Wright et al., 2004). By 7 ka BP (7.8 cal ka BP), increasing incursions of the Maritime Tropical air mass were delivering sufficient precipitation, indicated by the appearance of mesic forest species in pollen diagrams of the Great Lakes region (Webb III et al., 1998) to convert the Michigan, Huron and Georgian Bay water bodies to open, overflowing lakes, as at present. These changes are consistent with paleovegetation maps which show a rapid northward migration through the Great Lakes region of the mixed-boreal forest biome boundary between 8.0 and 7.0 ka BP (8.9 and 7.8 cal ka BP) (Dyke et al., 2004).

The foregoing hypothesized climatic history is consistent with changes inferred for small Elk Lake in Minnesota, based on a comprehensive study of proxy climatic and limnological indicators (Bradbury and Dean, 1993). Similarly, this climatic history is supported by a coeval depletion in the ¹⁸O composition of precipitated carbonate in Deep Lake, Minnesota, attributed by Yu and Wright (2001) to the blocking of southern air masses by more frequent presence of dry Arctic air.

The phase of reduced and closed lakes at the onset of the present hydrologic regime offers a unique opportunity to evaluate the sensitivity of the Great Lakes system to high-amplitude climatic change. Understanding the sensitivity of the lakes to high-amplitude and long-duration change would be a distinct benefit in the light of the need to project and adapt to future changes under global warming which may drive the lakes below instrumentally-observed variability (Mortsch et al., 2000).

As the evaporation process favours concentration of water molecules containing the heavier ¹⁸O isotope, previous findings of high concentrations of the lighter ¹⁶O isotope, similar to that of glacial meltwater, in fossil valves of Huron basin benthic ostracodes at about 7600 ka BP (Rea et al., 1994a, 1994b; Dettman et al., 1995) appear to contradict the coeval presence of evaporatively-driven lowstands as postulated in this paper and in Lewis et al. (in press). Although surface lake water isotopic composition undoubtedly became concentrated in ¹⁸O during seasonal periods of rapid evaporation, it apparently did not greatly influence bottom water composition, similar to the isotopic stratification found for glacial Lake Agassiz (Buhay, 1998; Birks et al., in press). Alternatively, adjustments in the chronology of isotopic events in the Huron and Michigan basins suggest that the low ¹⁸O inflow occurred prior to the lowstands (Breckenridge and Johnson, 2005; Breckenridge, in press), and may have remained as bottom water during the evaporative phase. Additional research is needed to resolve the origin of bottom water in the Huron and Michigan basins.